《北航数值分析第三次大作业.docx》由会员分享，可在线阅读，更多相关《北航数值分析第三次大作业.docx（36页珍藏版）》请在三一办公上搜索。
1、北航数值分析第三次大作业数值分析第三次大作业 一、算法的设计方案： 、总体方案设计： 解非线性方程组。将给定的(xi,yi)当作已知量代入题目给定的非线性方程组，求得与(xi,yi)相对应的数组tij,uij。 分片二次代数插值。通过分片二次代数插值运算，得到与数组t1121,u1121对应的数组z1121，得到二元函数z=f(xi,yi)。 曲面拟合。利用xi,yj,z1121建立二维函数表，再根据精度的要求选择适当k值，并得到曲面拟合的系数矩阵Crs。 观察和p(xi,yi)的逼近效果。观察逼近效果只需要重复上面和的过程，得到与新的插值节点(xi,yi)对应的f(xi,yi)，再与对应的p
2、(xi,yi)比较即可，这里求解p(xi,yi)可以直接使用中的Crs和k。 具体算法设计： 解非线性方程组 牛顿法解方程组F(x)=0的解x，可采用如下算法： 1）在x附近选取x*(0)*D，给定精度水平e0和最大迭代次数M。 2）对于k=0,1,LM执行 计算F(x(k)和F(x(k)。 求解关于Dx(k)的线性方程组 F(x(k)Dx(k)=-F(xk( ) 若Dx(k)x(k)*(k)e，则取xx，并停止计算；否则转。 计算x(k+1)=x(k)+Dx(k)。 若kxn-1-h2，插值节点对应取i=1或i=n-1， ，插值节点对应取j=1或i=m-1。 t2或yyn-1-t2hhx-x
3、x+,2in-2ii22 若 tty-yy+,2jm-2jj22则选择(xk,yr)(k=i-1,i,i+1;r=j-1,j,j+1)为插值节点。 2）计算 i+1lk(x)=t=i-1tkj+1x-xtxk-xty-ytyr-ytk(=i-1i,i+,1)l%r(y)=r(=j-1j,j+,1)t=j-1tr插值多项式的公式为： i+1j+1 p(x,y)=k=i-1r=j-1lk(x)l%r(y)f(xk,yr) 注：本步进行插值运算的是(t,u)，利用(xi,yj)与(t,u)的对应关系就可以得到z与(xi,yj)的对应关系。 曲面拟合 根据插值得到的数表xi,yj,f(xi,yj)进行
4、曲面拟合的过程： 1） 根据拟合节点和基底函数写出矩阵B和G： (x0)00(x)1 B=M(x)0n(x0)(x1)M(xn)11LLOL1k(y0)0(x0)k0(x1)(y)1 G=MMk(y)0(xn)m(y0)(y1)M1LLO1(ym)1Lk(y0)k(y1) Mk(ym)2） 计算 C=(BTB)-1BTUG(GTG)-1。 在这里，为了简化计算和编程、避免矩阵求逆，记： A=(BB)BU，DT=G(GTG)-1 对上面两式进行变形，得到如下两个线性方程组： (BB)A=BU，(GG)D=G TTTTT-1T通过解上述两个线性方程组，则有：C=AD kkT3） 对于每一个(xi,
5、yj)， p(xi,yj)=4） 拟合需要达到的精度条件为： nm*r=0s=0Crs(xi)(yj)。 rs s=pi=0j=0(xi,yj)-uij102-。 7 其中uij对应着插值得到的数表xi,yj,f(xi,yj)中f(xi,yj)的值。 5） 让k逐步增加，每一次重复执行以上几步，直到 nm*s= pi=0j=0(xi,yj)-uij102-7 成立。此时的k值就是要求解最小的k。 二、源程序： #include #include #include #include #include #include #define Epsilon1 1e-12 /*解线性方程组时近似解向量的精
6、度*/ #define M 200 /*解线性方程组时的最大迭代次数*/ #define N 10 /*求解迭代次数时假设的k的最大值，用于定义包含k的存储空间*/ void Newton; /*牛顿法求解非线性方程组子程序*/ void fpeccz; /*分片二次代数插值子程序*/ void qmnh; /*曲面拟合子程序*/ void duibi; /*对比𝑓和p逼近效果的子程序*/ double x11,y21,t1121,u1121;/*定义全局变量*/ double z1121,C1010; double kz; void Newton(double x11,dou
7、ble y21)/*牛顿法求解非线性方程组子程序*/ double X4,dx4,F4,dF44,temp,m,fx,fX; int i,j,k,l,p,ik,n; for(i=0;i=10;i+) for(j=0;j=20;j+) X0=1; /*选取迭代初始向量，四个分别代表t，u，v，w*/ X1=1; X2=1; X3=1; n=0; loop1: F0=0.5*cos(X0)+X1+X2+X3-xi-2.67; F1=X0+0.5*sin(X1)+X2+X3-yj-1.07; F2=0.5*X0+X1+cos(X2)+X3-xi-3.74; F3=X0+0.5*X1+X2+sin(X
8、3)-yj-0.79; /*求解F*/ dF00=-0.5*sin(X0); /*求解F*/ dF01=1; dF02=1; dF03=1; dF10=1; dF11=0.5*cos(X1); dF12=1; dF13=1; dF20=0.5; dF21=1; dF22=-sin(X2); dF23=1; dF30=1; dF31=0.5; dF32=1; dF33=cos(X3); /*高斯选主元消去法求解x*/ for(k=0;k3;k+) ik=k; for(l=k;l=3;l+) if(dFikkdFlk) ik=l; /*选主元*/ temp=0; temp=Fik; Fik=Fk;
9、 Fk=temp; for(l=k;l=3;l+) temp=0; temp=dFikl; dFikl=dFkl; dFkl=temp; for(l=k+1;l=3;l+) m=dFlk/dFkk; Fl=Fl-m*Fk; for(p=k+1;p=0;k-) temp=0; for(l=k+1;l=3;l+) temp=temp+dFkl*dxl/dFkk; dxk=-Fk/dFkk-temp; temp=0; for(l=0;l=3;l+) /*求解矩阵范数，用无穷范数*/ if(tempfabs(dxl) temp=fabs(dxl); fx=temp; temp=0; for(l=0;l
10、=3;l+) if(tempfabs(Xl) temp=fabs(Xl); fX=temp; if(fabs(fx/fX)Epsilon1) /*判断Dx(k) tij=X0; uij=X1; goto loop4; else for(l=0;l=3;l+) Xl=Xl+dxl; n=n+1; goto loop3; loop3:if(nM) /*判断是否超出规定迭代次数*/ goto loop1; else printf(迭代不成功n); goto loop4; loop4:continue; void fpeccz(double t1121,double u1121)/*分片二次代数插值子
11、程序*/ int s1121,r1121; int i,j,i1,j1,m; double z066=-0.5,-0.34,0.14,0.94,2.06,3.5, x(k)e是否成立*/ -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; double t06=0,0.2,0.4,0.6,0.8,1.0; double
12、u06=0,0.4,0.8,1.2,1.6,2.0; double temp1,temp2; for(i=0;i=10;i+) /*选取插值节点*/ for(j=0;j=20;j+) if(tij0.7) sij=4; else for(m=2;m0.2*m-0.1)&(tij=0.2*m+0.1) sij=m; for(i=0;i=10;i+) for(j=0;j=20;j+) if(uij1.4) rij=4; else for(m=2;m0.4*m-0.2)&(uij=0.4*m+0.2) rij=m; for(i=0;i=10;i+) /*插值运算*/ for(j=0;j=20;j+)
13、 zij=0; for(i1=sij-1;i1=sij+1;i1+) for(j1=rij-1;j1=rij+1;j1+) temp1=1.0; for(m=sij-1;m=sij+1;m+) if(m!=i1) temp1*=(tij-t0m)/(t0i1-t0m); temp2=1.0; for(m=rij-1;m1e-7) for(i=0;i=10;i+) for(j=0;j=k;j+) Bij=pow(xi,j); Btji=Bij; for(i=0;i=k;i+) for(j=0;j=k;j+) temp=0; for(l=0;l=10;l+) temp+=Btil*Blj; BtB
14、ij=temp; for(i=0;i=k;i+) for(j=0;j=20;j+) temp=0; for(l=0;l=10;l+) temp+=Btil*zlj; BtUij=temp; for(l=0;l=20;l+) for(i=0;i=k;i+) for(j=0;j=k;j+) BtB1ij=BtBij; for(m=0;m=k-1;m+) ik=m; for(i=m;i=k-1;i+) if(fabs(BtB1im)fabs(BtB1i+1m) ik=i+1; else ; if(ik!=m) for(i=m;i=k;i+) temp=BtB1mi; BtB1mi=BtB1iki;
15、BtB1iki=temp; temp=BtUml; BtUml=BtUikl; BtUikl=temp; for(i=m+1;i=k;i+) q=BtB1im/BtB1mm; for(j=m;j=0;m-) temp=0; for(i=m+1;i=k;i+) temp+=Ail*BtB1mi; Aml=(BtUml-temp)/BtB1mm; for(i=0;i=20;i+) for(j=0;j=k;j+) Gij=pow(yi,j); Gtji=Gij; for(i=0;i=k;i+) for(j=0;j=k;j+) temp=0; for(m=0;m=20;m+) temp+=Gtim*G
16、mj; GtGij=temp; for(l=0;l=20;l+) for(i=0;i=k;i+) for(j=0;j=k;j+) GtG1ij=GtGij; for(m=0;m=k-1;m+) ik=m; for(i=m;i=k-1;i+) if(fabs(GtG1im)fabs(GtG1i+1m) ik=i+1; else ; if(ik!=m) for(i=m;i=k;i+) temp=GtG1mi; GtG1mi=GtG1iki; GtG1iki=temp; temp=Gtml; Gtml=Gtikl; Gtikl=temp; for(i=m+1;i=k;i+) q=GtG1im/GtG
17、1mm; for(j=m;j=0;m-) temp=0; for(i=m+1;i=k;i+) temp+=Dil*GtG1mi; Dml=(Gtml-temp)/GtG1mm; 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=0;/*归零，开始计算sigma值*/ 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(xi,i1)*po
18、w(yj,j1); sigma+=pow(pij-zij,2); printf(k=%d sigma=%.12en,k,sigma); k=k+1; k-; printf(达到精度要求时的k和sigma值：n); printf(k=%d sigma=%.12en,k,sigma); kz=k; /*定义为全局变量，便于duibi调用*/ void duibi /*对比𝑓和p逼近效果的子程序*/ int i,j,i1,j1; double p85; for(i=0;i=7;i+) for(j=0;j=4;j+) xi=0.1*(i+1); yj=0.5+0.2*(j+1); /*
19、重新输入节点*/ Newton(x,y); fpeccz(t,u); /*求解𝑓*/ for(i=0;i=7;i+) /*求解p*/ for(j=0;j=4;j+) pij=0; for(i1=0;i1=kz;i1+) for(j1=0;j1=kz;j1+) pij+=Ci1j1*pow(xi,i1)*pow(yj,j1); printf(x%d=%.6f y%d=%.6fn,i+1,xi,j+1,yj); printf(f(x%d,y%d)=%.12en,i+1,j+1,zij); printf(p(x%d,y%d)=%.12en,i+1,j+1,pij); printf(=
20、%.12en,zij-pij); /*数表xi,yi,𝑓 (xi,yi),p*/ void main int i,j; for(i=0;i=10;i+) for(j=0;j=20;j+) xi=0.08*i; yj=0.5+0.05*j; /*输入节点*/ Newton(x,y); fpeccz(t,u); * for(i=0;i=10;i+) for(j=0;j=20;j+) printf(x%d=%.6f n,i,xi,j,yj,i,j,zij); /*数表：xi,yi, 𝑓 (xi,yi)*/ y%d=%.6f z%d%d=%.12e qmnh; for(
21、i=0;i=kz;i+) for(j=0;j=kz;j+) printf(C%d%d=%.12e n,i,j,Cij); /*数表：Crs*/ duibi; 三、运行结果输出 1、数表数表：xi,yi,𝑓(xi,yi) x0=0.000000 y0=0.500000 z00=4.465040184807e-001 x0=0.000000 y1=0.550000 z01=3.246832629277e-001 x0=0.000000 y2=0.600000 z02=2.101596866827e-001 x0=0.000000 y3=0.650000 z03=1.03043608
22、3160e-001 x0=0.000000 y4=0.700000 z04=3.401895562675e-003 x0=0.000000 y5=0.750000 z05=-8.873581363800e-002 x0=0.000000 y6=0.800000 z06=-1.733716327497e-001 x0=0.000000 y7=0.850000 z07=-2.505346114666e-001 x0=0.000000 y8=0.900000 z08=-3.202765063876e-001 x0=0.000000 y9=0.950000 z09=-3.826680661097e-0
23、01 x0=0.000000 y10=1.000000 z010=-4.377957667384e-001 x0=0.000000 y11=1.050000 z011=-4.857589414438e-001 x0=0.000000 y12=1.100000 z012=-5.266672548835e-001 x0=0.000000 y13=1.150000 z013=-5.606384797965e-001 x0=0.000000 y14=1.200000 z014=-5.877965387677e-001 x0=0.000000 y15=1.250000 z015=-6.082697790
24、899e-001 x0=0.000000 y16=1.300000 z016=-6.221894528764e-001 x0=0.000000 y17=1.350000 z017=-6.296883781856e-001 x0=0.000000 y18=1.400000 z018=-6.308997600028e-001 x0=0.000000 y19=1.450000 z019=-6.259561525454e-001 x0=0.000000 y20=1.500000 z020=-6.149885466094e-001 x1=0.080000 y0=0.500000 z10=6.380152
25、265113e-001 x1=0.080000 y1=0.550000 z11=5.066117551467e-001 x1=0.080000 y2=0.600000 z12=3.821763692774e-001 x1=0.080000 y3=0.650000 z13=2.648634911537e-001 x1=0.080000 y4=0.700000 z14=1.547802002848e-001 x1=0.080000 y5=0.750000 z15=5.199268349094e-002 x1=0.080000 y6=0.800000 z16=-4.346804020490e-002
26、 x1=0.080000 y7=0.850000 z17=-1.316010567885e-001 x1=0.080000 y8=0.900000 z18=-2.124310883088e-001 x1=0.080000 y9=0.950000 z19=-2.860045510580e-001 x1=0.080000 y10=1.000000 z110=-3.523860789794e-001 x1=0.080000 y11=1.050000 z111=-4.116554565222e-001 x1=0.080000 y12=1.100000 z112=-4.639049115188e-001
27、 x1=0.080000 y13=1.150000 z113=-5.092367247005e-001 x1=0.080000 y14=1.200000 z114=-5.477611179623e-001 x1=0.080000 y15=1.250000 z115=-5.795943883391e-001 x1=0.080000 y16=1.300000 z116=-6.048572588895e-001 x1=0.080000 y17=1.350000 z117=-6.236734213318e-001 x1=0.080000 y18=1.400000 z118=-6.361682484133e-001 x1=0.080000 y19=1.450000 z119=-6.424676566901e-001 x1=0.080000 y20=1.500000 z120=-6.426971026996e-001 x2=0.160000 y0=0.500000 z20=8.400813957666e-001 x2=0.160000 y1=0.550000 z21=6.997641656732e-001 x2=0.160000 y2=0.600000 z22=5.66061442351
