【文章內(nèi)容簡(jiǎn)介】 素的前一個(gè)元素 (同一行 ) //迭代初值 數(shù)值計(jì)算課程設(shè)計(jì) 17 h = bb aa。 y[0] = h*(f(aa) + f(bb))/。 m = 1。 n = 1。 ep = epsilon + 。 //迭代計(jì)算 while ((ep = epsilon) amp。amp。 (m MAXREPT)) { //復(fù)化積 分公式求 T2n(Romberg 計(jì)算表中的第一列 ),n初始為 1,以后倍增 p = 。 for (int i=0。 in。 i++)//求 Hn { x = aa + (i+)*h。 p = p + f(x)。 } p = (y[0] + h*p)/。//求 T2n = 1/2(Tn+Hn),用 p指示 //求第 m行元素 ,根據(jù) Romberg 計(jì)算表本行的前一 個(gè)元素 (p 指示 ), //和上一行左上角元素 (y[k1]指示 )求得 . s = 。 for (int k=1。 k=m。 k++) { s = *s。 q = (s*p y[k1])/(s )。 y[k1] = p。 p = q。 龍貝格求積分算法 18 } p = fabs(q y[m1])。 m = m + 1。 y[m1] = q。 n = n + n。 h = h/。 } return (q)。 } int main() { double a,b。 coutRomberg 積分 ,請(qǐng)輸入積分范圍 a,b:endl。 cinab。 cout積分結(jié)果 :Romberg(a, b)endl。 return 0。 } 數(shù)值計(jì)算課程設(shè)計(jì) 19 三次樣條插值算法 三次樣條插值算法說 明表 表 51 三次樣條插值算法說明表 策略描述 包含 0m 和 Nm 的方程 (i)三次緊壓樣條，確定 0()Sx? ， ()nSx? （如果導(dǎo)數(shù)已知，這是“最佳選擇”） 10 0 003 ( ( )) 2mm d S xh ?? ? ? 1113 ( ( ) ) 2NN N NN mm S x dh ??? ?? ? ? (ii)natural 三次樣條（一條“松弛曲線”） 0 0m? ， 0Nm? (iii)外掛 ()Sx?? 到端點(diǎn) 0 2 101 1()h m mmm h??? 1 1 21 2()N N NNN Nh m mmm h? ? ?? ???? (iv) ()Sx?? 是靠近端點(diǎn)的常量 01mm? ， 1NNmm?? (v)在每個(gè)端點(diǎn)處指定 ()Sx?? 00()m S x??? ， ()NNm S x??? 、 三次樣條插值算法（壓緊樣條）程序調(diào)試 我們將所編寫的程序三次樣條插值算法（壓緊樣條）程序進(jìn)行調(diào)試 圖 51三次樣條插值算法（壓緊樣條）程序 輸入界面、運(yùn)行結(jié)果 圖 52三次樣條插值算法程序 運(yùn)行結(jié)果（ a） 圖 52三次樣條插值算法程序運(yùn)行結(jié)果 （ b）三次樣條插值算法 20 運(yùn)行結(jié)果分析： ? ?? ? ? ? ? ?? ?? ?????????????????????????????????????32)2()2()3(21)1()2()1(10)0(333333xxxxxxxxxxxxxxxxf （ 51） 作圖程序（ Matlab）： x1=0::1。 y1=*(x1 1).^3 + *(x1 0).^3 *(x1 1) + *(x1 0)。 x2=1::2。 y2=*(x2 2).^3 *(x2 1).^3 *(x2 2) + *(x2 1)。 x3=2::3。 y3=*(x3 3).^3 + *(x3 2).^3 *(x3 3) + *(x3 2)。 X=[0 1 2 3]。 Y=[0 2 ]。 plot(x1,y1,x2,y2,x3,y3,X,Y,39。*39。) gtext(39。S139。) gtext(39。S239。) gtext(39。S339。) 圖形為： 圖 53 三次樣條插值算法 Matlab作圖分析 數(shù)值計(jì)算課程設(shè)計(jì) 21 、 三次樣條插值算法（壓緊樣條）代碼 includeiostream includeiomanip using namespace std。 const int max = 50。 float x[max], y[max], h[max]。 float c[max], a[max], fxym[max]。 float f(int x1, int x2, int x3) { float a = (y[x3] y[x2]) / (x[x3] x[x2])。 float b = (y[x2] y[x1]) / (x[x2] x[x1])。 return (a b)/(x[x3] x[x1])。 } //求差分 void cal_m(int n) { //用追趕法求解出彎矩向量 M?? float B[max]。 B[0] = c[0] / 2。 for(int i = 1。 i n。 i++) B[i] = c[i] / (2 a[i]*B[i1])。 fxym[0] = fxym[0] / 2。 for(i = 1。 i = n。 i++) fxym[i] = (fxym[i] a[i]*fxym[i1]) / (2 a[i]*B[i1])。 for(i = n1。 i = 0。 i) fxym[i] = fxym[i] B[i]*fxym[i+1]。 三次樣條插值算法 22 } void printout(int n)。 int main() { int n,i。 char ch。 do{ cout輸入 x的最大下標(biāo) :。 cinn。 for(i = 0。 i = n。 i++) { coutPlease put in Xi39。:39。 cinx[i]。 coutPlease put in Yi39。:39。 ciny[i]。 } for(i = 0。 i n。 i++) //求步長(zhǎng) h[i] = x[i+1] x[i]。 coutPlease 輸入邊界條件 \n 1 :已知兩端的一階導(dǎo)數(shù) \n 2 :兩端的二階導(dǎo)數(shù)已知 \n 默認(rèn) :自然邊界條件 \n。 int t。 float f0, f1。 cint。 switch(t) { case 1:cout輸入 Y0\39。 Yn\39。\n。 cinf0f1。 數(shù)值計(jì)算課程設(shè)計(jì) 23 c[0] = 1。 a[n] = 1。 fxym[0] = 6*((y[1] y[0]) / (x[1] x[0]) f0) / h[0]。 fxym[n] = 6*(f1 (y[n] y[n1]) / (x[n] x[n1])) / h[n1]。 break。 case 2:cout輸入 Y0\ Yn\\n。 cinf0f1。 c[0] = a[n] = 0。 fxym[0] = 2*f0。 fxym[n] = 2*f1。 break。 default:cout不可用 \n。//待定 }。 for(i = 1。 i n。 i++) fxym[i] = 6 * f(i1, i, i+1)。 for(i = 1。 i n。 i++) { a[i] = h[i1] / (h[i] + h[i1])。 c[i] = 1 a[i]。 } a[n] = h[n1] / (h[n1] + h[n])。 cal_m(n)。 cout\n輸出三次樣條插值函數(shù)： \n。 printout(n)。 coutDo you have another try ? y/n :。 cinch。 三次樣條插值算法 24 }while(ch == 39。y39。 || ch == 39。Y39。)。 return 0。} void printout(int n) {coutsetprecision(6)。 for(int i = 0。 i n。 i++) {co