《北航数值分析大作业.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业.docx（35页珍藏版）》请在三一办公上搜索。
1、北航数值分析大作业一、题目： 关于x, y, t, u, v, w的下列方程组 0.5cost+u+v+w-x=2.67t+0.5sinu+v+w-y=1.07 0.5t+u+cosv+w-x=3.74t+0.5u+v+sinw-y=0.79以及关于z, t, u的下列二维数表 z u t 0 0.2 0.4 0.6 0.8 1.0 0 -0.5 -0.42 -0.18 0.22 0.78 1.5 0.4 -0.34 -0.5 -0.5 -0.34 -0.02 0.46 0.8 0.14 -0.26 -0.5 -0.58 -0.5 -0.26 1.2 0.94 0.3 -0.18 -0.5 -
2、0.66 -0.66 1.6 2.06 1.18 0.46 -0.1 -0.5 -0.74 2.0 3.5 2.38 1.42 0.62 -0.02 -0.5 确定了一个二元函数z = f(x, y)。 1、试用数值方法求出f(x, y)在区域 D=(x,y)|0x0.8,0.5y1.5上的一个近似表达式 p(x,y)=要求p(x,y)一最小的k值达到以下的精度 r,s=0ckrsxrys s=(f(xi,yj)-p(xi,yj)210-7 i=0j=01020其中，xi=0.08i,yj=0.5+0.05j。 *2、计算f(xi,yj),p(xi,yj)的值，以观察p(x,y)逼近f(x,y
3、)的效果，其中，xi*=0.1i, y*j=0.5+0.2j。 说明：1、用迭代方法求解非线性方程组时，要求近似解向量x(k)满足 |x(k)-x(k-1)|/|x(k)|10-12 2、作二元插值时，要使用分片二次代数插值。 3、要由程序自动确定最小的k值。 4、打印以下内容： 算法的设计方案。 全部源程序。 数表：xi,yj,f(xi,yj) 。 选择过程的k,s值。 达到精度要求时的k,s值以及p(x,y)中的系数crs。 *数表：xi*,y*。 j,f(xi,yj),p(xi,yj)5、采用f型输出xi,yj,xi*,y*j的准确值，其余实型数采用e型输出并且至少显示12位有效数字。
4、二、算法设计方案 1.使用牛顿迭代法，对原题中给出的xi=0.08i，yj=0.5+0.05j，,0j20)的11*21组xi,yj分别求出原题中方程组的一组解，于是得(0i10到一组和xi,yi对应的ui,tj。 2.对于已求出的ui,tj，使用分片二次代数插值法对原题中关于z,t,u的数表进行插值得到zij。于是产生了z=f(x,y)的11*21个数值解。 3.从k=1开始逐渐增大k的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的k,s。当sM时仍未达到迭代精度，则迭代计算失0.5*cos(t) + u + v + w - x - 2.67t + 0.5*sin(u)
5、 + v + w - y - 1.07， F(x)=0.5*t + u + cos(v) + w - x - 3.74t + 0.5*u + v + sin(w) - y - 0.79111-0.5*sin(t)10.5*cos(u)11， F(x)=0.51-sin(v)110.51cos(w)分片双二次插值： 根据题目给出的表格， ti=ih (i=0,1,L,5)ui=jt (j=0,1,L,5)对于给定的(t,u),如果(t,u)满足 t=0.4） 中取i=1或i=4; 如果uu1+或uu4+,2222则在式中取u=1或u=4。 在区域D=(x,y)|0x0.8,0.5y1.5上，将(
6、xi,yj) (xi=0.08i,i=0,1,10;yj=0.5+0.2j,j=0,1,20)代入到非线性方程组中，用Newton法解出(ti,j,ui,j)，再由分片双二次插值得zi,j=h(ti,j,ui,j)，则有zi,j=f(xi,yj)，即求出了z=f(x,y)。 2、求p(x,y)： 乘积型最小二乘拟合曲面： 求系数矩阵C： C=(BTB)-1BTUG(GTG)-1 其中， 1x01x1B=MM1x101y01y1G=MM1y20x02x12Mx102y02y12Mx202Lx0kLx1kOMLx10kLy0kLy1k OMLy20kz0,0z1,0U=Mz10,0z0,1z1,1
7、Mz10,1z0,20z1,20,z=f(x,y) (i=0,1,L,10;j=0,1,L,20) ijOMi,jLz10,20LL计算中涉及到对矩阵求逆，接着在后面将会具体说明列主元的高斯消去法求矩阵的逆的方法。 确定最小的k值，拟合曲面： p(x,y)=设s1020r,s=0rscxrsyk=(f(xi,yj)-p(xi,yj)2，给定精度水平e2=10-7和最大迭代次数N，i=0j=0则确定最小k值的迭代格式为： k(k)rsp(x,y)=cxy (i=0,1,L,10;j=0,1,L,20)ijrsijr,s=0(k)1020(k)2 s=(f(xi,yj)-p(xi,yj)i=0j=
8、0k=0,1,2,L迭代终止条件为s(k)e2，若kN时仍未达到迭代精度，则迭代计算失败。 待确定满足精度条件的最小k值后，就可以进行曲面拟合计算了。 3、关于列主元的高斯消去法求矩阵的逆： 设非奇异矩阵ACnn，且B=A ，则 -1AB=I， 对B和I列分块，有 A(b1,L,bj,L,bn)=(e1,L,ej,L,en) 即 Abj=ej, j=1,2,L,n 其中， Tej=(0,L,0,1,0,L,0) 14243j应用列主元的高斯消去法线性解方程组 则B=(b1,L,bj,L,bn)即为A的逆。 注：若A不可逆，则此算法失效。 四、源程序 #include Abj=ej, j=1,2
9、,L,n， / zhhfshzhfx3.cpp : 定义控制台应用程序的入口点。 / #include #include const int M= 500;/迭代最大次数 const double E=1.0e-12;/确定牛顿迭代精度水平 const double E1=1.0e-7;/确定拟合精度水平 const int kmax=9;/k的最大值 const double matrix66= -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,
10、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 const double mat_t6=0, 0.2, 0.4, 0.6, 0.8, 1.0; const double mat_u6=0, 0.4, 0.8, 1.2, 1.6, 2.0; double U1121;/对f(x,y)的值进行存储 double Ckmax+1kmax+1;/拟合系数矩阵 /两数绝对值取大/ double ma
11、x(double x,double y) /高斯消元法解线性方程组/ void linear_solution(double f4,double ff44, double (&delta)4) int i, j, k, ik; double tmp, mik; for(k=0; k3; k+) delta3=f3/ff33; ik=k; for(i=k;i4;i+) tmp=ffkj; ffkj=ffikj; ffikj=tmp; tmp=fk; fk=fik; fik=tmp; for(i=k+1;i4;i+) mik=ffik/ffkk; for(j=k;j4;j+) ffij=ffij-
12、mik*ffkj; if(fabs(ffikk)fabs(y)?fabs(x):fabs(y); for(j=k;j=0;k-) tmp=0; tmp+=ffkj*deltaj; for(j=k+1;j4;j+) deltak=(fk-tmp)/ffkk; /Newton法求解非线性方程组/ void equation_solution(double x, double y,double &u,double &t) double v=1.0,w=1.0; u=1.0; t=1.0; double delta4,f4,ff44; int k=0; for(k=0;k=M;k+) f0=-1.0*
13、(0.5*cos(t)+u+v+w-x-2.67); f1=-1.0*(t+0.5*sin(u)+v+w-y-1.07); f2=-1.0*(0.5*t+u+cos(v)+w-x-3.74); f3=-1.0*(t+0.5*u+v+sin(w)-y-0.79); ff00=-0.5*sin(t); ff01=1.0; ff02=1.0; ff03=1.0; ff10=1.0; ff11=0.5*cos(u); ff12=1.0; ff13=1.0; ff20=0.5; ff21=1.0; ff22=-sin(v); ff23=1.0; ff30=1.0; ff31=0.5; ff32=1.0;
14、 ff33=cos(w); linear_solution(f,ff,delta); if(max(delta3,max(delta2,max(delta0,delta1)/max(w,max(v,max(u,t)=E) t+=delta0; break; else u+=delta1; v+=delta2; w+=delta3; if(k=M) printf(非线性方程组迭代不成功!n); /分片双二次插值/ double interpolation(double a,double b) int i,j,k,r,t1,t2; double z=0.0; i=int(fabs(a/0.2)+0
15、.5); j=int(fabs(b/0.4)+0.5); if(i=0) i=1; if(i=5) i=4; if(j=0) j=1; if(j=5) j=4; /矩阵求逆/ void inverse(double (&matrix)kmax+1kmax+1,int k) double matrixTkmax+1kmax+1=0.0; double bkmax+1kmax+1=0.0,tmp,mik; int i,j,t,p,it; for(i=0;i=kmax;i+) for(j=0;j=kmax;j+) if(i=j) bij=1.0; else bij=0.0; for(k=i-1;k=
16、i+1; k+) for(r=j-1;r=j+1;r+) double sum=1.0; sum*=matrixkr; for(t1=i-1;t1=i+1;t1+) if(t1!=k) sum*=(a-mat_tt1)/(mat_tk-mat_tt1); for(t2=j-1;t2=j+1;t2+) if(t2!=r) sum*=(b-mat_ut2)/(mat_ur-mat_ut2); z+=sum; return z; for(t=0; tk; t+) for(i=0;i=k;i+) for(j=0;j=k;j+) matrixij=matrixTij; matrixTkp=bkp/mat
17、rixkk; for(tmp=0.0,j=i+1;j=0;i-) for(p=0;p=k;p+) it=t; for(i=t+1;i=k;i+) for(p=0;p=k;p+) tmp=btp; btp=bitp; bitp=tmp; for(i=t+1;i=k;i+) mik=matrixit/matrixtt; for(j=t+1;j=k;j+) matrixij-=mik*matrixtj; bip-=mik*btp; tmp=matrixitj; matrixitj=matrixtj; matrixtj=tmp; if(fabs(matrixitt)fabs(matrixit) it=
18、i; for(j=t;j=k;j+) for(p=0;p=k;p+) /针对每个k值进行曲面拟合/ double surface_fitting(int k, double x, double y) double B11kmax+1=0.0,BTkmax+111=0.0; double G21kmax+1=0.0,GTkmax+121=0.0; double BTBkmax+1kmax+1=0.0,BTB1kmax+1kmax+1=0.0,BTB2kmax+1kmax+1=0.0,GTGkmax+1kmax+1=0.0; for(i=0;i=k;i+) for(j=0;j=k;j+) for(
19、sum=0.0,t=0;t=10;t+) sum+=BTit*Btj; BTBij=sum; for(i=0;i=10;i+) for(j=0;j=k;j+) Bij=pow(0.08*i,j); double matrix1kmax+111=0.0,matrix321kmax+1=0.0; double matrix2kmax+121=0.0; double p,sum; int i,j,t,r,s; for(i=0;i=20;i+) for(j=0;j=k;j+) Gij=pow(0.5+0.05*i,j); for(i=0;i=10;i+) for(j=0;j=k;j+) BTji=Bi
20、j; for(i=0;i=20;i+) for(j=0;j=k;j+) GTji=Gij; inverse(BTB,k); for(i=0;i=k;i+) for(j=0;j=k;j+) for(sum=0.0,t=0;t=20;t+) sum+=GTit*Gtj; GTGij=sum; inverse(GTG,k); for(i=0;i=k;i+) for(j=0;j=10;j+) for(j=0;j=k;j+) for(j=0;j=20;j+) for(j=0;j=k;j+) for(sum=0.0,t=0;t=20;t+) sum+=matrix2it*matrix3tj; Cij=su
21、m; for(sum=0.0,t=0;t=10;t+) sum+=matrix1it*Utj; matrix2ij=sum; for(sum=0.0,t=0;t=k;t+) sum+=Git*GTGtj; matrix3ij=sum; for(sum=0.0,t=0;t=k;t+) sum+=BTBit*BTtj; matrix1ij=sum; for(i=0;i=20;i+) for(i=0;i=k;i+) for(i=0;i=k;i+) p=0.0; for(r=0;r=k;r+) for(s=0;s=k;s+) p+=Crs*pow(x,r)*pow(y,s); return p; /主
22、函数/ int _tmain(int argc, _TCHAR* argv) double x,y,t,u,sum,sigma; int i,j,k,kmin; printf(数表：xi,yi,f(xi,yi)(i=0,1,2.10;j=0,1,2.,20):n); for(i=0;i=10;i+) for(j=0;j=20;j+) x=0.08*i; y=0.5+0.05*j; equation_solution(x,y,u,t); printf(x%d=%f, y%d=%f, f(x%d, y%d)=%.12len,i,x,j,y,i,j,Uij); Uij=interpolation(t
23、,u); printf(选择过程的k,值分别为：n); for(k=0;k=kmax;k+) printf(数表：x*i,y*j,f(x*i,y*j),p(x*i,y*j)(i=1,2,.,8;j=1,2,.,5):n); for(i=1;i=8;i+) for(j=1;j=5;j+) equation_solution(0.1*i,0.5+0.2*j,u,t); printf(f(x*%d,y*%d)=%.12le,i,j,interpolation(t,u); printf(p(x*%d,y*%d)=%.12len,i,j,surface_fitting(kmin,0.1*i,0.5+0.
24、2*j); sum=0.0; for(i=0;i=10;i+) for(j=0;j=20;j+) sum+=pow(surface_fitting(k,0.08*i,0.5+0.05*j)-Uij,2); printf(%d %.12len,k,sum); if(sum=E1) else if(k=kmax) printf(无法达到曲面拟合精度要求!); kmin=k; sigma=sum; printf(达到精度要求的k,值分别为：nk=%d,=%.12len,kmin,sigma); printf(p(x,y)中的系数Crs(r=0,1,.,k;s=0,1,.,k):n); for(i=0
25、;i=k;i+) for(j=0;j=k;j+) printf(Crs%d%d=%.12len,i,j,Cij); break; printf(x*%d=%f,y*%d=%fn,i,0.1*i,j,0.5+0.2*j); return 0; 五、程序输出 - 1. 数表：xi,yi,f(xi,yi)(i=0,1,2.10;j=0,1,2.,20): x0=0.000000, y0=0.500000, f(x0, y0)=4.465040184807e-001 x0=0.000000, y1=0.550000, f(x0, y1)=3.246832629277e-001 x0=0.000000,
26、 y2=0.600000, f(x0, y2)=2.101596866827e-001 x0=0.000000, y3=0.650000, f(x0, y3)=1.030436083160e-001 x0=0.000000, y4=0.700000, f(x0, y4)=3.401895562675e-003 x0=0.000000, y5=0.750000, f(x0, y5)=-8.873581363800e-002 x0=0.000000, y6=0.800000, f(x0, y6)=-1.733716327497e-001 x0=0.000000, y7=0.850000, f(x0
27、, y7)=-2.505346114666e-001 x0=0.000000, y8=0.900000, f(x0, y8)=-3.202765063876e-001 x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x0=0.000000, x1=0.080000, x1=0.080000, x1=0.080000, x1=0.080000, x1=0.080
28、000, x1=0.080000, x1=0.080000, x1=0.080000, x1=0.080000, x1=0.080000, x1=0.080000, x1=0.080000, y9=0.950000, y10=1.000000, y11=1.050000, y12=1.100000, y13=1.150000, y14=1.200000, y15=1.250000, y16=1.300000, y17=1.350000, y18=1.400000, y19=1.450000, y20=1.500000, y0=0.500000, y1=0.550000, y2=0.600000
29、, y3=0.650000, y4=0.700000, y5=0.750000, y6=0.800000, y7=0.850000, y8=0.900000, y9=0.950000, y10=1.000000, y11=1.050000, f(x0, y9)=-3.826680661097e-001 f(x0, y10)=-4.377957667384e-001 f(x0, y11)=-4.857589414438e-001 f(x0, y12)=-5.266672548835e-001 f(x0, y13)=-5.606384797965e-001 f(x0, y14)=-5.877965
30、387677e-001 f(x0, y15)=-6.082697790899e-001 f(x0, y16)=-6.221894528764e-001 f(x0, y17)=-6.296883781856e-001 f(x0, y18)=-6.308997600028e-001 f(x0, y19)=-6.259561525454e-001 f(x0, y20)=-6.149885466094e-001 f(x1, y0)=6.380152265113e-001 f(x1, y1)=5.066117551467e-001 f(x1, y2)=3.821763692774e-001 f(x1,
31、y3)=2.648634911537e-001 f(x1, y4)=1.547802002848e-001 f(x1, y5)=5.199268349094e-002 f(x1, y6)=-4.346804020490e-002 f(x1, y7)=-1.316010567885e-001 f(x1, y8)=-2.124310883088e-001 f(x1, y9)=-2.860045510580e-001 f(x1, y10)=-3.523860789794e-001 f(x1, y11)=-4.116554565222e-001 x1=0.080000, y12=1.100000, f
32、(x1, y12)=-4.639049115188e-001 x1=0.080000, y13=1.150000, f(x1, y13)=-5.092367247005e-001 x1=0.080000, y14=1.200000, f(x1, y14)=-5.477611179623e-001 x1=0.080000, y15=1.250000, f(x1, y15)=-5.795943883391e-001 x1=0.080000, y16=1.300000, f(x1, y16)=-6.048572588895e-001 x1=0.080000, y17=1.350000, f(x1,
33、y17)=-6.236734213318e-001 x1=0.080000, x1=0.080000, x1=0.080000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x2=0.160000, x
34、2=0.160000, x2=0.160000, x2=0.160000, y18=1.400000, y19=1.450000, y20=1.500000, y0=0.500000, y1=0.550000, y2=0.600000, y3=0.650000, y4=0.700000, y5=0.750000, y6=0.800000, y7=0.850000, y8=0.900000, y9=0.950000, y10=1.000000, y11=1.050000, y12=1.100000, y13=1.150000, y14=1.200000, y15=1.250000, y16=1.
35、300000, y17=1.350000, y18=1.400000, y19=1.450000, y20=1.500000, f(x1, y18)=-6.361682484133e-001 f(x1, y19)=-6.424676566901e-001 f(x1, y20)=-6.426971026996e-001 f(x2, y0)=8.400813957666e-001 f(x2, y1)=6.997641656732e-001 f(x2, y2)=5.660614423517e-001 f(x2, y3)=4.391716081176e-001 f(x2, y4)=3.19242138
36、0408e-001 f(x2, y5)=2.063761923874e-001 f(x2, y6)=1.006385238914e-001 f(x2, y7)=2.060740067837e-003 f(x2, y8)=-8.935402476698e-002 f(x2, y9)=-1.736269688648e-001 f(x2, y10)=-2.507999561599e-001 f(x2, y11)=-3.209322694446e-001 f(x2, y12)=-3.840977350046e-001 f(x2, y13)=-4.403821754175e-001 f(x2, y14)=-4.898811523126e-001 f(x2, y15)=-5.326979655338e-001 f(x2, y16)=-5.689418792921e-001 f(x2, y1
