數(shù)值分析實(shí)習(xí)三

北京航空航天大學(xué)實(shí)驗(yàn)報(bào)告課程名稱：數(shù)值分析A第三次實(shí)驗(yàn)實(shí)驗(yàn)名稱：數(shù)值分析A計(jì)算實(shí)習(xí)題目三院〔系〕：自動(dòng)化科學(xué)與電氣工程學(xué)院專業(yè)：控制科學(xué)與工程姓名：保然學(xué)號(hào)：SY1203302實(shí)驗(yàn)時(shí)間：2012年11月26日算法的設(shè)計(jì)方案利用已給出的關(guān)于x,y,t,u,v,w的非線性方程組，以及已經(jīng)給出的一組(x,y)，可以解非線性方程組，計(jì)算出與(x,y)相對(duì)應(yīng)的(t,u)的值；在此根底上，再利用關(guān)于t,u,z的數(shù)表，對(duì)上一步中求出的(t,u)值進(jìn)行分片二次插值，即可得到與(t,u)對(duì)應(yīng)的z的值〔因?yàn)椴逯岛瘮?shù)在插值點(diǎn)上是準(zhǔn)確值〕，從而得到了與(x,y)對(duì)應(yīng)的z=f(x,y)值〔準(zhǔn)確值〕，從而構(gòu)造出關(guān)于x,y,z的數(shù)表；接著利用x,y,z的數(shù)表進(jìn)行曲面擬合，得到滿足最小二乘要求的擬合函數(shù)p(x,y)，即求出了p(x,y)的所有系數(shù)；最后用一組給定的(x*,y*)，求解關(guān)于x,y,t,u,v,w的非線性方程組，得到與(x*,y*)對(duì)應(yīng)的(t*,u*)，利用數(shù)表t,u,z進(jìn)行分片二次插值，即可得到與(x*,y*)對(duì)應(yīng)的z*=f(x*,y*)的值，同時(shí)利用上一步中解得的系數(shù)Crs求解出p*=(x*,y*)，比擬z*=f(x*,y*)和p*=(x*,y*)的大小以檢驗(yàn)擬合結(jié)果。注意：a)求解非線性方程組時(shí)使用Newton法，給定精度水平，最大迭代次數(shù)M=20，在每次迭代前構(gòu)造和，當(dāng)時(shí)停止計(jì)算。b)分片二次插值時(shí)注意函數(shù)值點(diǎn)所落在的區(qū)域，同時(shí)注意行值和列值得區(qū)別，t對(duì)應(yīng)的數(shù)值確定行值，u對(duì)應(yīng)的數(shù)值確定列值。c)曲面擬合時(shí)，關(guān)于x的基函數(shù)為，關(guān)于y的基函數(shù)為其中k從0開(kāi)始取值，每做完一次曲面擬合后，計(jì)算，當(dāng)時(shí)，使k=k+1，再次進(jìn)行曲面擬合，直到時(shí)停止，記錄此時(shí)的k值以及p(x,y)中的系數(shù)Crs。源程序源代碼如下：#include<stdio.h>#include<math.h>#include<iostream>#include<iomanip>usingnamespacestd;voidgauss_elimination(doubleA[][4],doubleb[])//高斯消去法求解線性方程組的解向量，返回一個(gè)向量{ intk,i,j,ik; doublec,d,mik; doubleB[4][4]; for(i=0;i<=3;i++) { for(j=0;j<=3;j++) { B[i][j]=A[i][j];//將系數(shù)矩陣存儲(chǔ)用于計(jì)算，而不破壞其原始結(jié)構(gòu) } } for(k=0;k<=2;k++) { /*選主元*/ c=fabs(B[k][k]); for(i=k,ik=k;i<=3;i++) { if(c<fabs(B[i][k])) { c=fabs(B[i][k]); ik=i; } } for(j=k;j<=3;j++) { d=B[k][j]; B[k][j]=B[ik][j]; B[ik][j]=d; } d=b[k];b[k]=b[ik]; b[ik]=d; /*選主元結(jié)束*/ for(i=k+1;i<=3;i++) { mik=B[i][k]/B[k][k]; for(j=k;j<=3;j++) { B[i][j]=B[i][j]-mik*B[k][j]; } b[i]=b[i]-mik*b[k]; } } b[3]=b[3]/B[3][3]; for(k=2;k>=0;k--) { for(j=k+1,d=0;j<=3;j++) { d+=B[k][j]*b[j]; } b[k]=(b[k]-d)/B[k][k]; } return;}doublemax(doublea[])//求一個(gè)向量的無(wú)窮范數(shù)，并返回?zé)o窮范數(shù)的值{ inti; doubleb; for(i=0,b=fabs(a[0]);i<=3;i++) { if(b<fabs(a[i]))b=fabs(a[i]); } return(b);}voidnewton(doubleB[],doubleX[])//Newton法求解非線性方程組的解向量，返回一個(gè)向量{ inti,k,M=20;//設(shè)置迭代上限M doublex[4]={1,0.5,1,1},dx[4],a,b; doubleF[4],dF[4][4]={{1,1,1,1},{1,1,1,1},{0.5,1,1,1},{1,0.5,1,1}}; for(k=0;k<=M;k++) { F[0]=0.5*cos(x[0])+x[1]+x[2]+x[3];F[1]=x[0]+0.5*sin(x[1])+x[2]+x[3]; F[2]=0.5*x[0]+x[1]+cos(x[2])+x[3]; F[3]=x[0]+0.5*x[1]+x[2]+sin(x[3]); dF[0][0]=-0.5*sin(x[0]); dF[1][1]=0.5*cos(x[1]); dF[2][2]=-sin(x[2]); dF[3][3]=cos(x[3]); for(i=0;i<=3;i++) { dx[i]=B[i]-F[i]; } gauss_elimination(dF,dx); a=max(dx); b=max(x); if(a/b<=1e-12)break; else { for(i=0;i<=3;i++) { x[i]=x[i]+dx[i]; } } } for(i=0;i<=3;i++) { X[i]=x[i]; } return;}doublelk(doublex[],doublexm,inti,intk)//構(gòu)造Lagrange二次插值基函數(shù){ doublelk; if(k==i-1) lk=(xm-x[i])*(xm-x[i+1])/(x[k]-x[i])/(x[k]-x[i+1]); elseif(k==i)lk=(xm-x[i-1])*(xm-x[i+1])/(x[k]-x[i-1])/(x[k]-x[i+1]); else lk=(xm-x[i-1])*(xm-x[i])/(x[k]-x[i-1])/(x[k]-x[i]); return(lk);}voidercichazhi_1(doublex[][21],doubley[][21],doublez[][21])//分片二次插值求解二次函數(shù)的函數(shù)值{ inti,j,m,n; doublea,b,c,d,e,f; doubleu[6]={0,0.4,0.8,1.2,1.6,2},t[6]={0,0.2,0.4,0.6,0.8,1}; doubleg[6][6]={{-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}}; for(m=0;m<=10;m++) { for(n=0;n<=20;n++) { i=int(fabs(x[m][n]/0.4)+0.5);//確定數(shù)據(jù)點(diǎn)落在哪一列 if(i==5)i=4; if(i==0)i=1;a=lk(u,x[m][n],i,i-1); b=lk(u,x[m][n],i,i); c=lk(u,x[m][n],i,i+1); j=int(fabs(y[m][n]/0.2)+0.5);//確定數(shù)據(jù)點(diǎn)落在哪一行 if(j==5)j=4; if(j==0)j=1;d=lk(t,y[m][n],j,j-1); e=lk(t,y[m][n],j,j); f=lk(t,y[m][n],j,j+1);z[m][n]=a*d*g[j-1][i-1]+a*e*g[j][i-1]+a*f*g[j+1][i-1]+b*d*g[j-1][i]+b*e*g[j][i]+b*f*g[j+1][i]+c*d*g[j-1][i+1]+c*e*g[j][i+1]+c*f*g[j+1][i+1]; } } return;}voidgauss_elimination_2(doubleA[][6],doubleb[],intn)//高斯消去法求系數(shù)矩陣階數(shù)可變的方程組的解向量{intk,i,j,ik; doublec,d,mik; doubleB[6][6]={}; for(i=0;i<=n-1;i++) { for(j=0;j<=n-1;j++) { B[i][j]=A[i][j];//將系數(shù)矩陣存儲(chǔ)用于計(jì)算，而不破壞其原始結(jié)構(gòu) } } for(k=0;k<=n-2;k++) { /*選主元*/ c=fabs(B[k][k]); for(i=k,ik=k;i<=n-1;i++) { if(c<fabs(B[i][k])) { c=fabs(B[i][k]); ik=i; } } for(j=k;j<=n-1;j++) { d=B[k][j]; B[k][j]=B[ik][j]; B[ik][j]=d; } d=b[k];b[k]=b[ik]; b[ik]=d; /*選主元結(jié)束*/ for(i=k+1;i<=n-1;i++) { mik=B[i][k]/B[k][k]; for(j=k;j<=n-1;j++) { B[i][j]=B[i][j]-mik*B[k][j]; } b[i]=b[i]-mik*b[k]; } } b[n-1]=b[n-1]/B[n-1][n-1]; for(k=n-2;k>=0;k--) { for(j=k+1,d=0;j<=n-1;j++) { d+=B[k][j]*b[j]; } b[k]=(b[k]-d)/B[k][k]; } return;}voidintial(doubleB[][6],doubleG[][6],intn)//初始化曲面擬合法中的矩陣B和G{ inti,j; for(i=0;i<=10;i++) for(j=0;j<=n;j++) B[i][j]=pow(0.08*i,j); for(i=0;i<=20;i++) for(j=0;j<=n;j++) G[i][j]=pow(0.5+0.05*i,j); return;}voidqumiannihe(doubleB[][6],doubleG[][6],doubleU[][21],doubleC[][6],intn)//曲面擬合求解擬合系數(shù){ inti,j,k; doubleBTB[6][6]={},GTG[6][6]={},BTU[6][21]={},GT[6][21]={},A[6][21]={},D[21][6]={},DT[6][21]={}; doubleBT[6][11]={}; doublea[6]={},d[6]={};for(i=0;i<=n;i++) for(j=0;j<=10;j++) BT[i][j]=B[j][i];//初始化BT for(i=0;i<=n;i++) for(j=0;j<=20;j++) GT[i][j]=G[j][i];//初始化GT for(i=0;i<=n;i++) for(j=0;j<=n;j++) for(k=0;k<=10;k++) BTB[i][j]+=BT[i][k]*B[k][j];//求解BT*Bfor(i=0;i<=n;i++) for(j=0;j<=n;j++) for(k=0;k<=20;k++) GTG[i][j]+=GT[i][k]*G[k][j];//求解GT*G for(i=0;i<=n;i++) for(j=0;j<=20;j++) for(k=0;k<=10;k++) BTU[i][j]+=BT[i][k]*U[k][j];//求解BT*U for(i=0;i<=20;i++) {for(j=0;j<=n;j++) { a[j]=BTU[j][i]; d[j]=GT[j][i]; } gauss_elimination_2(BTB,a,n+1); gauss_elimination_2(GTG,d,n+1); for(j=0;j<=n;j++) { A[j][i]=a[j]; DT[j][i]=d[j]; } } for(i=0;i<=20;i++) for(j=0;j<=n;j++) D[i][j]=DT[j][i];for(i=0;i<=n;i++) for(j=0;j<=n;j++) for(k=0;k<=20;k++) C[i][j]+=A[i][k]*D[k][j];//求解系數(shù)矩陣C return;}voidercichazhi_2(doublex[][5],doubley[][5],doublez[][5])//分片二次插值求解二次函數(shù)的函數(shù)值{ inti,j,m,n; doublea,b,c,d,e,f; doubleu[6]={0,0.4,0.8,1.2,1.6,2},t[6]={0,0.2,0.4,0.6,0.8,1}; doubleg[6][6]={{-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}}; for(m=0;m<=7;m++) { for(n=0;n<=4;n++) { i=int(fabs(x[m][n]/0.4)+0.5);//確定數(shù)據(jù)點(diǎn)落在哪一列 if(i==5)i=4; if(i==0)i=1;a=lk(u,x[m][n],i,i-1); b=lk(u,x[m][n],i,i); c=lk(u,x[m][n],i,i+1); j=int(fabs(y[m][n]/0.2)+0.5);//確定數(shù)據(jù)點(diǎn)落在哪一行 if(j==5)j=4; if(j==0)j=1;d=lk(t,y[m][n],j,j-1); e=lk(t,y[m][n],j,j); f=lk(t,y[m][n],j,j+1);z[m][n]=a*d*g[j-1][i-1]+a*e*g[j][i-1]+a*f*g[j+1][i-1]+b*d*g[j-1][i]+b*e*g[j][i]+b*f*g[j+1][i]+c*d*g[j-1][i+1]+c*e*g[j][i+1]+c*f*g[j+1][i+1]; } } return;}voidmain(){ inti,j,n=5,k=0,r,s; doublex[11],y[21],X[4],z[11][21],u[11][21],t[11][21]; doubleB[11][6]={},G[21][6]={},C[6][6]={}; doubleBM[4],p,eps=0,derta[11][21]; doublex1[8],y1[5],u1[8][5],t1[8][5],f[8][5],p1[8][5]={}; for(i=0;i<=10;i++) { x[i]=0.08*i;//初始化xi } for(j=0;j<=20;j++) { y[j]=0.5+0.05*j;//初始化yi } /*Newton法求解非線性方程組的解*/ for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { BM[0]=2.67+x[i]; BM[1]=1.07+y[j]; BM[2]=3.74+x[i]; BM[3]=0.79+y[j]; newton(BM,X); u[i][j]=X[1];//提取出u t[i][j]=X[0];//提取出t } } /*Newton法求解非線性方程組的解結(jié)束*/ ercichazhi_1(u,t,z);//用與x,y相對(duì)應(yīng)的t,u值進(jìn)行分片二次插值，得到相應(yīng)的函數(shù)值z(mì)=f(x,y) /*輸出數(shù)表x,y,z*/ for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { cout<<"x"<<i<<"="<<x[i]<<""<<"y"<<j<<"="<<y[j]<<""; cout<<"f(x"<<i<<",y"<<j<<")="<<setprecision(12)<<z[i][j]<<endl; } } cout<<endl; /*輸出數(shù)表x,y,z結(jié)束*/for(n=0;n<=5;n++){ for(i=0;i<=5;i++) { for(j=0;j<=5;j++) { C[i][j]=0; } } intial(B,G,n);qumiannihe(B,G,z,C,n);//用曲面擬合法求解擬合曲面，得到滿足最小二乘要求的擬合函數(shù)p〔x,y〕for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { for(r=0,p=0;r<=n;r++) { for(s=0;s<=n;s++) { p+=C[r][s]*pow(x[i],r)*pow(y[j],s); } } derta[i][j]=(z[i][j]-p)*(z[i][j]-p); } } for(i=0,eps=0;i<=10;i++) for(j=0;j<=20;j++) eps+=derta[i][j]; cout<<"當(dāng)擬合函數(shù)次數(shù)k為"<<n<<"時(shí):";//輸出迭代次數(shù)k cout<<"函數(shù)精度σ為:"<<eps<<setprecision(12)<<endl;//輸出精度 if(eps<=1e-7)break;//判斷是否到達(dá)精度要求 }cout<<endl; /*輸出到達(dá)精度要求時(shí)的擬合系數(shù)Crs*/ cout<<"當(dāng)函數(shù)精度σ到達(dá)"<<eps<<setprecision(12)<<"時(shí)p(x,y)中的系數(shù)Crs為:"<<endl; for(i=0;i<=n;i++) { for(j=0;j<=n;j++) { cout<<"C("<<i<<","<<j<<")="<<setprecision(12)<<C[i][j]<<""; } cout<<endl; } cout<<endl; /*輸出到達(dá)精度要求時(shí)的擬合系數(shù)Crs結(jié)束*/for