計算實習(xí)報告3

1、 1、算法設(shè)計：整體思路這一次的作業(yè)和前面兩次相比，感覺要難很多，最明顯的就是書上沒有完整的算法實現(xiàn)流程。再加上插值和逼近這一段在課上也沒聽太懂什么意思，所以剛開始思考這個程序時感覺很吃力。不過經(jīng)過后來好好看書和向同學(xué)請教，漸漸對插值和分片元二次擬合有點懂了，也大概知道了解決本次問題分別要完成的子問題。要完成本次任務(wù)，主體上有三部分問題要解決：（1）解非線性方程組。將d內(nèi)的)當(dāng)作已知量代入題目給定的非線性方程組，得到與相對應(yīng)的tij,uij。（2）分片雙二次插值。采用分片雙二次插值，得到與t1121,u1121對應(yīng)的z1121，也即建立與的對應(yīng)關(guān)系，得到二元函數(shù)。（3）曲面擬合。利用xi,yj

2、,z1121建立二維函數(shù)表，再根據(jù)精度的要求選擇適當(dāng)k值，并得到曲面擬合的系數(shù)矩陣ckk。（4）觀察逼近效果。觀察逼近效果只需要重復(fù)上面（1）和（2）的過程，得到與新的插值節(jié)點對應(yīng)的，再與對應(yīng)的比較即可。程序示意圖如下所示：各部分問題的算法實現(xiàn)（1）解非線性方程組非線性方程組常用數(shù)值解法有簡單迭代法和牛頓法。牛頓法收斂快，一般都能達(dá)到平方收斂，因此這里選擇newton法解非線性方程組。牛頓法解方程組的解，可采用如下算法：1）在附近選取，給定精度水平和最大迭代次數(shù)m。2）對于執(zhí)行 計算和。 求解關(guān)于的線性方程組 若，則取，并停止計算；否則轉(zhuǎn)。 計算。 若，則繼續(xù)，否則，輸出m次迭代不成功的信息，

3、并停止計算。注：第步中的求取，用到了解線性方程組，使用的是列主元素gauss消去法。 另外，本次作業(yè)中解非線性方程組實際求取的是解向量和，方程組中的、 是已知值，為當(dāng)前節(jié)點的值。（2）分片雙二次插值給定已知數(shù)表以及需要插值的節(jié)點，進(jìn)行分片二次插值的算法：設(shè)已知數(shù)表中的點為： ，需要插值的節(jié)點為。1) 根據(jù)選擇插值節(jié)點：若或，插值節(jié)點對應(yīng)取或，若或，插值節(jié)點對應(yīng)取或。 若則選擇為插值節(jié)點。2）計算 插值多項式的公式為： 注：在（1）中通過解非線性方程組得到的解向量和與插值節(jié)點有著一一對應(yīng)關(guān)系，而本題題目給出的函數(shù)表為關(guān)于的函數(shù)關(guān)系。因此，為得到關(guān)于的函數(shù)關(guān)系，本題中應(yīng)先根據(jù)給定數(shù)表對進(jìn)行插值，再

4、利用與的一一對應(yīng)關(guān)系得到與的對應(yīng)關(guān)系。（3）曲面擬合根據(jù)插值得到的數(shù)表進(jìn)行曲面擬合的過程：1） 根據(jù)擬合節(jié)點和基底函數(shù)寫出矩陣b和g： 2） 計算 。在這里，為了簡化計算和編程、避免矩陣求逆，記：，對上面兩式進(jìn)行變形，得到如下兩個線性方程組：，通過解上述兩個線性方程組，則有：3） 對于每一個， 。4） 擬合需要達(dá)到的精度條件為： 。 其中對應(yīng)著插值得到的數(shù)表中的值。5） 讓k逐步增加，每一次重復(fù)執(zhí)行以上幾步，直到 成立。此時的k值就是所需最小的k。注：曲面擬合過程中用到的數(shù)表為前面插值得到的數(shù)表，即。在2)中，求矩陣a和d的過程用的是列主元素gauss消去法。（4）觀察逼近效果。觀察逼近效果主

5、要要做的是，通過新的插值節(jié)點 建立新的插值數(shù)表，同時求出對應(yīng)的，比較即可。2、程序源代碼：#include stdio.h#include conio.h#include math.h#define epsilon1 1e-12 /解線性方程組時近似解向量的精度#define maxm 200 /解線性方程組時的最大迭代次數(shù)#define k 10 /預(yù)估的能達(dá)到精度要求的k值，k=k#define givenx 11 /已知要求的插值節(jié)點xi個數(shù)#define giveny 21 /已知要求的插值節(jié)點yi個數(shù)#define testx 8 /觀察逼近效果時，插值節(jié)點xi個數(shù)#define t

6、esty 5 /觀察逼近效果時，插值節(jié)點yi個數(shù)/*/* 函數(shù)功能描述： */*/* norm()：求向量x的無窮范數(shù) */* newton_nonlinear()：牛頓法解非線性方程組的主干部分 */* gausselemitation_select()：線性方程組的列主元素gauss消去法 */* interplate()：分片二元二次插值 */* surfacefit()：曲面擬合，函數(shù)返回滿足精度要求的最小k值 */* calculate_a()：根據(jù)擬合曲面的系數(shù)矩陣c=a(dt)，求矩陣a */* calculate_d()：根據(jù)擬合曲面的系數(shù)矩陣c=a(dt)，求矩陣d */*

7、test_observe()：觀察逼近效果子函數(shù) */*/double norm(double x4);int newton_nonlinear(double xstar4,double xi,double yi);void gausselemitation_select(double df44,double f4,double delta_x4);void interplate(double t1121,double u1121,double z1121,int numi,int numj);int surfacefit(double z1121,double x11,double y21,

8、double ckk);void calculate_a(double ak21,double z1121,double x11,int k);void calculate_d(double dk21,double z1121,double y11,int k);void test_observe(int k,double ckk,double z1121);void main()/*/* 變量描述： */*/* x11、y21：插值節(jié)點(xi,yi) */* t1121、u1121：與節(jié)點(xi,yi)對應(yīng)的二元組(tij,uij) */* gausselemitation_select()

9、：線性方程組的列主元素gauss消去法 */* z1121：存放插值得到的數(shù)表 */* ckk：存放曲面擬合得到的系數(shù)矩陣，其中k是預(yù)估的最小k值 */*/ double t1121,u1121,x11,y21,z1121,ckk; double xstar4; int i,j,k,flag;/*/* 解非線性方程組 */*/ for(i=0;i=givenx-1;i+) /通過解非線性方程組，建立(x,y)與(t,u)的對應(yīng)關(guān)系 for(j=0;j=giveny-1;j+) xi=0.08*i; yj=0.5+0.05*j; xstar0=1; xstar1=1; xstar2=1; xst

10、ar3=1; /選取迭代初始值x(0) flag=newton_nonlinear(xstar,xi,yj); /牛頓法解非線性方程組 if(flag=-1) goto end; /如果非線性方程組求解失敗，中止程序 tij=xstar0; uij=xstar1; /*/* 插值并輸出插值得到的數(shù)表 */*/ interplate(t,u,z,11,21); /分片二次插值 printf(1、插值得到的數(shù)表為：n); for(i=0;i=10;i+) /輸出插值得到的數(shù)表：xi,yi,f(xi,yi) for(j=0;j=20;j+) printf(x%d= %.6f ,i,xi); prin

11、tf(y%d= %.6f ,j,yj); printf(z%d%d= %.12en,i,j,zij); printf(n);/*/* 曲面擬合并輸出擬合得到的矩陣c */*/ k=surfacefit(z,x,y,c); /曲面擬合，返回滿足精度要求的最小k值 printf(系數(shù)矩陣c為：n); for(i=0;i=k;i+) /輸出曲面擬合得到的系數(shù)矩陣c for(j=0;j=k;j+) printf(c%d%d=%.12e ,i,j,cij); if(j!=0&j%2=1) printf(n); /兩個一行輸出 printf(n); /*/* 觀察插值的逼近效果 */*/ printf(3

12、、觀察插值的逼近效果：n); test_observe(k,c,z); /觀察逼近效果end: ; /作為求解非線性方程組迭代次超限的出口 getch();/* 牛頓法解非線性方程組 */int newton_nonlinear(double xstar4,double xi,double yj) int k,l; double delta_x4,f4,df44; k=0; while(1) /* 計算f(x(k)的值 */ f0=0.5*cos(xstar0)+xstar1+xstar2+xstar3-xi-2.67; f1=xstar0+0.5*sin(xstar1)+xstar2+xst

13、ar3-yj-1.07; f2=0.5*xstar0+xstar1+cos(xstar2)+xstar3-xi-3.74; f3=xstar0+0.5*xstar1+xstar2+sin(xstar3)-yj-0.79; /* 計算f(x(k)的值，f(x(k)為一個4x4矩陣*/ df00=-0.5*sin(xstar0); df01=1; df02=1; df03=1; df10=1; df11=0.5*cos(xstar1); df12=1; df13=1; df20=0.5; df21=1; df22=-sin(xstar2); df23=1; df30=1; df31=0.5; df

14、32=1; df33=cos(xstar3); gausselemitation_select(df,f,delta_x); /列主元素gauss消去法解線性方程組 if(norm(delta_x)/norm(xstar)=epsilon1) /達(dá)到精度要求，跳出 break; for(l=0;l=maxm) /迭代次數(shù)超限，返回-1作為后面程序是否執(zhí)行的標(biāo)志 printf(迭代次數(shù)k=%d超出限制！n,k); return(-1); /* 列主元素gauss消去法解線性方程組 */void gausselemitation_select(double df44,double f4,doubl

15、e delta_x4) int i,j,k,line_flag; double m_ik,temp; for(k=0;k=2;k+) /*選行號*/ line_flag=k; for(i=k;i=2;i+) if(fabs(dfik)fabs(dfi+1k) line_flag=i+1; else ; if(line_flag!=k) /第k個主元不在第k行時，交換兩行元素 for(i=k;i=3;i+) temp=dfki; dfki=dfline_flagi; dfline_flagi=temp; temp=fk; fk=fline_flag; fline_flag=temp; /*消元*

16、/ for(i=k+1;i=3;i+) m_ik=dfik/dfkk; for(j=k;j=0;k-) temp=0; for(j=k+1;j=3;j+) temp+=delta_xj*dfkj; delta_xk=(-fk-temp)/dfkk; /* 求向量x的無窮范數(shù) */double norm(double x4) int i; double temp=0; for(i=0;i=3;i+) if(tempfabs(xi) temp=fabs(xi); return(temp);/* 分片二元二次插值 */void interplate(double t1121,double u1121

17、,double z1121,int numi,int numj) int k1121,r1121; int i,j,i1,j1,m; double z066=-0.5,-0.34,0.14,0.94,2.06,3.5, -0.42,-0.5,-0.26,0.3,1.18,2.38, -0.18,-0.5,-0.5,-0.18,0.46,1.42, 0.22,-0.34,-0.58,-0.5,-0.1,0.62, 0.78,-0.02,-0.5,-0.66,-0.5,-0.02, 1.5,0.46,-0.26,-0.66,-0.74,-0.5; /已知數(shù)表 double t06=0,0.2,0.

18、4,0.6,0.8,1.0; double u06=0,0.4,0.8,1.2,1.6,2.0; /原始已知節(jié)點 double temp1,temp2; /* 選擇插值節(jié)點，xi、yi分開選擇 */ for(i=0;i=numi-1;i+) /選擇插值節(jié)點xi for(j=0;j=numj-1;j+) if(tij0.7) kij=4; else for(m=2;m0.2*m-0.1)&(tij=0.2*m+0.1) kij=m; /end if /end for /end xi for(i=0;i=numi-1;i+) /選擇插值節(jié)點yi for(j=0;j=numj-1;j+) if(ui

19、j1.4) rij=4; else for(m=2;m0.4*m-0.2)&(uij=0.4*m+0.2) rij=m; /end if /end for /end yi /* 插值 */ for(i=0;i=numi-1;i+) for(j=0;j=numj-1;j+) zij=0; /初始化節(jié)點zij for(i1=kij-1;i1=kij+1;i1+) for(j1=rij-1;j1=rij+1;j1+) temp1=1.0; for(m=kij-1;m=kij+1;m+) /計算 if(m!=i1) temp1*=(tij-t0m)/(t0i1-t0m); temp2=1.0; for

20、(m=rij-1;m1e-7) calculate_a(a,z,x,k); /求矩陣a calculate_d(d,z,y,k); /求矩陣d /*計算矩陣c，c=a(dt)*/ for(i=0;i=k;i+) for(j=0;j=k;j+) temp=0; for(m=0;m=20;m+) temp+=aim*djm; cij=temp; /*計算sigma，選擇滿足精度要求的最小k*/ sigma=0; for(i=0;i=10;i+) for(j=0;j=20;j+) pij=0; for(i1=0;i1=k;i1+) for(j1=0;j1=k;j1+) pij+=ci1j1*pow(

21、xi,i1)*pow(yj,j1); /end p(xi,yj) sigma+=pow(pij-zij,2); /end sigma printf(k=%d =%.12en,k,sigma); /輸出當(dāng)前k值和 k+; k-; printf(曲面擬合達(dá)到精度要求，此時：n); printf(k=%d =%.12en,k,sigma); return k;/*/* 函數(shù)名：calculate_a */*/* 功 能：根據(jù)公式 (bt)b)a=(bt)u，通過解線性方程組得到矩陣a */*/void calculate_a(double ak21,double z1121,double x11,i

22、nt k) int i,j,m,l,line_flag; double b11k,btk11,btbkk,btuk21; /b和bt的階數(shù)分別為11xk和kx11 double btb_1kk,temp; /a的階數(shù)為kx21 double m_ik; /* 計算矩陣b */ for(i=0;i=10;i+) for(j=0;j=k;j+) bij=pow(xi,j); btji=bij; /end b /* 計算矩陣bt和b相乘，并將結(jié)果存在矩陣btb中 */ for(i=0;i=k;i+) for(j=0;j=k;j+) temp=0; for(m=0;m=10;m+) temp+=bti

23、m*bmj; btbij=temp; /end btb /* 計算矩陣bt和u相乘，并將結(jié)果存在矩陣btu中 */ for(i=0;i=k;i+) for(j=0;j=20;j+) temp=0; for(m=0;m=10;m+) temp+=btim*zmj; btuij=temp; /end btu /* 列主元素gauss消元法解矩陣a */ for(l=0;l=20;l+) /*復(fù)制矩陣btb*/ for(i=0;i=k;i+) for(j=0;j=k;j+) btb_1ij=btbij; /*gauss消元法*/ for(m=0;m=k-1;m+) /*選行號*/ line_flag

24、=m; for(i=m;i=k-1;i+) if(fabs(btb_1imbtb_1i+1m) line_flag=i+1; else ; if(line_flag!=m) /第k個主元不在第k行時，交換兩行元素 for(i=m;i=k;i+) temp=btb_1mi; btb_1mi=btb_1line_flagi; btb_1line_flagi=temp; temp=btuml; btuml=btuline_flagl; btuline_flagl=temp; /*消元*/ for(i=m+1;i=k;i+) m_ik=btb_1im/btb_1mm; for(j=m;j=0;m-)

25、temp=0; for(i=m+1;i=k;i+) temp+=ail*btb_1mi; aml=(btuml-temp)/btb_1mm; /end 回代 /end 求矩陣a/*/* 函數(shù)名：calculate_d */*/* 功 能：根據(jù)公式 (gt)g)d=gt，通過解線性方程組得到矩陣d */*/void calculate_d(double dk21,double z1121,double y11,int k) int i,j,m,l,line_flag; double g21k,gtk21,gtgkk; /g和gt的階數(shù)分別為11xk和kx11 double gtg_1kk,tem

26、p; /g的階數(shù)為kx21 double m_ik; /* 計算矩陣g */ for(i=0;i=20;i+) for(j=0;j=k;j+) gij=pow(yi,j); gtji=gij; /end g /* 計算矩陣gt和g相乘，并將結(jié)果存在矩陣gtg中 */ for(i=0;i=k;i+) for(j=0;j=k;j+) temp=0; for(m=0;m=20;m+) temp+=gtim*gmj; gtgij=temp; /end gtg /* 列主元素gauss消元法解矩陣d */ for(l=0;l=20;l+) /*復(fù)制矩陣gtg*/ for(i=0;i=k;i+) for(

27、j=0;j=k;j+) gtg_1ij=gtgij; /*gauss消元法*/ for(m=0;m=k-1;m+) /*選行號*/ line_flag=m; for(i=m;i=k-1;i+) if(fabs(gtg_1imgtg_1i+1m) line_flag=i+1; else ; if(line_flag!=m) /第k個主元不在第k行時，交換兩行元素 for(i=m;i=k;i+) temp=gtg_1mi; gtg_1mi=gtg_1line_flagi; gtg_1line_flagi=temp; temp=gtml; gtml=gtline_flagl; gtline_flag

28、l=temp; /*消元*/ for(i=m+1;i=k;i+) m_ik=gtg_1im/gtg_1mm; for(j=m;j=0;m-) temp=0; for(i=m+1;i=k;i+) temp+=dil*gtg_1mi; dml=(gtml-temp)/gtg_1mm; /end 回代 /end 求矩陣d/* 觀察p(x,y)逼近f(x,y)效果 */void test_observe(int k,double ckk) double testx11,testy21,t1121,u1121; double xstar4,p1121,z1121; int i,j,i1,j1; for(

29、i=0;i=testx-1;i+) /通過解非線性方程組，建立(x,y)與(t,u)的對應(yīng)關(guān)系 for(j=0;j=testy-1;j+) testxi=0.1*(i+1); testyj=0.5+0.2*(j+1); xstar0=1; xstar1=1; xstar2=1; xstar3=1; /選取迭代初始值x(0) newton_nonlinear(xstar,testxi,testyj); tij=xstar0; uij=xstar1; interplate(t,u,z,8,5); /分片二次插值，求與(xi,yi)對應(yīng)的f(xi,yi)的值 for(i=0;i=testx-1;i+

30、) /求與(xi,yi)對應(yīng)的p(xi,yi)的值，并與f(xi,yi)比較輸出 for(j=0;j=testy-1;j+) pij=0; for(i1=0;i1=k;i1+) /求 p(xi,yj) 的值 for(j1=0;j1=k;j1+) pij+=ci1j1*pow(testxi,i1)*pow(testyj,j1); printf(x%d=%.6f y%d=%.6fn,i+1,testxi,j+1,testyj); printf(ttf(x%d,y%d)=%.12en,i+1,j+1,zij); printf(tttp(x%d,y%d)=%.12en,i+1,j+1,pij); pr

31、intf(tt=%.12en,zij-pij); 3、程序運行結(jié)果：1、插值得到的數(shù)表為：x0=0.000000 y0=0.500000 z00=4.465040184796e-001x0=0.000000 y1=0.550000 z01=3.246832629271e-001x0=0.000000 y2=0.600000 z02=2.101596866821e-001x0=0.000000 y3=0.650000 z03=1.030436083154e-001x0=0.000000 y4=0.700000 z04=3.401895561932e-003x0=0.000000 y5=0.750

32、000 z05=-8.873581363889e-002x0=0.000000 y6=0.800000 z06=-1.733716327506e-001x0=0.000000 y7=0.850000 z07=-2.505346114672e-001x0=0.000000 y8=0.900000 z08=-3.202765063879e-001x0=0.000000 y9=0.950000 z09=-3.826680661098e-001x0=0.000000 y10=1.000000 z010=-4.377957667384e-001x0=0.000000 y11=1.050000 z011=

33、-4.857589414438e-001x0=0.000000 y12=1.100000 z012=-5.266672548843e-001x0=0.000000 y13=1.150000 z013=-5.606384797965e-001x0=0.000000 y14=1.200000 z014=-5.877965387680e-001x0=0.000000 y15=1.250000 z015=-6.082697790900e-001x0=0.000000 y16=1.300000 z016=-6.221894528764e-001x0=0.000000 y17=1.350000 z017=