北航数值分析第三次大作业.doc

数值分析第三次大作业一、 算法的设计方案：（一）、总体方案设计：（1）解非线性方程组。将给定的当作已知量代入题目给定的非线性方程组，求得与相对应的数组tij,uij。（2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t1121,u1121对应的数组z1121，得到二元函数z=。（3）曲面拟合。利用xi,yj,z1121建立二维函数表，再根据精度的要求选择适当k值，并得到曲面拟合的系数矩阵Crs。（4）观察和的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点对应的，再与对应的比较即可，这里求解可以直接使用（3）中的Crs和k。（二）具体算法设计：（1）解非线性方程组牛顿法解方程组的解，可采用如下算法：1）在附近选取，给定精度水平和最大迭代次数M。2）对于执行 计算和。 求解关于的线性方程组 若，则取，并停止计算；否则转。 计算。 若，则继续，否则，输出M次迭代不成功的信息，并停止计算。（2）分片双二次插值给定已知数表以及需要插值的节点，进行分片二次插值的算法：设已知数表中的点为： ，需要插值的节点为。1) 根据选择插值节点：若或，插值节点对应取或，若或，插值节点对应取或。 若则选择为插值节点。2）计算 插值多项式的公式为： 注：本步进行插值运算的是，利用与的对应关系就可以得到与的对应关系。（3）曲面拟合根据插值得到的数表进行曲面拟合的过程：1） 根据拟合节点和基底函数写出矩阵B和G： 2） 计算 。在这里，为了简化计算和编程、避免矩阵求逆，记：，对上面两式进行变形，得到如下两个线性方程组：，通过解上述两个线性方程组，则有：3） 对于每一个， 。4） 拟合需要达到的精度条件为： 。 其中对应着插值得到的数表中的值。5） 让k逐步增加，每一次重复执行以上几步，直到 成立。此时的k值就是要求解最小的k。二、 源程序：#include#include#include #include #include#include #define Epsilon1 1e-12 /*解线性方程组时近似解向量的精度*/#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,double 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(X3)-yj-0.79; /*求解F（x）*/ dF00=-0.5*sin(X0); /*求解F（x）*/ 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;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=3;l+) if(tempfabs(Xl) temp=fabs(Xl); fX=temp; if(fabs(fx/fX)Epsilon1) /*判断是否成立*/ 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)/*分片二次代数插值子程序*/ int s1121,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; double t06=0,0.2,0.4,0.6,0.8,1.0; double 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+) 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; BtBij=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; 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*Gmj; 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/GtG1mm; 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)*pow(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); /*重新输入节点*/ Newton(x,y); fpeccz(t,u); /*求解（x*，y*）*/ for(i=0;i=7;i+) /*求解p（x*，y*）*/ 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(=%.12en,zij-pij); /*数表xi*,yi*, (xi*,yi*),p（xi*,yi*）*/ 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 y%d=%.6f z%d%d=%.12e n,i,xi,j,yj,i,j,zij); /*数表：xi,yi, (xi,yi)*/ qmnh(); for(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-001x0=0.000000 y1=0.550000 z01=3.246832629277e-001x0=0.000000 y2=0.600000 z02=2.101596866827e-001x0=0.000000 y3=0.650000 z03=1.030436083160e-001x0=0.000000 y4=0.700000 z04=3.401895562675e-003x0=0.000000 y5=0.750000 z05=-8.873581363800e-002x0=0.000000 y6=0.800000 z06=-1.733716327497e-001x0=0.000000 y7=0.850000 z07=-2.505346114666e-001x0=0.000000 y8=0.900000 z08=-3.202765063876e-001x0=0.000000 y9=0.950000 z09=-3.826680661097e-001x0=0.000000 y10=1.000000 z010=-4.377957667384e-001x0=0.000000 y11=1.050000 z011=-4.857589414438e-001x0=0.000000 y12=1.100000 z012=-5.266672548835e-001x0=0.000000 y13=1.150000 z013=-5.606384797965e-001x0=0.000000 y14=1.200000 z014=-5.877965387677e-001x0=0.000000 y15=1.250000 z015=-6.082697790899e-001x0=0.000000 y16=1.300000 z016=-6.221894528764e-001x0=0.000000 y17=1.350000 z017=-6.296883781856e-001x0=0.000000 y18=1.400000 z018=-6.308997600028e-001x0=0.000000 y19=1.450000 z019=-6.259561525454e-001x0=0.000000 y20=1.500000 z020=-6.149885466094e-001x1=0.080000 y0=0.500000 z10=6.380152265113e-001x1=0.080000 y1=0.550000 z11=5.066117551467e-001x1=0.080000 y2=0.600000 z12=3.821763692774e-001x1=0.080000 y3=0.650000 z13=2.648634911537e-001x1=0.080000 y4=0.700000 z14=1.547802002848e-001x1=0.080000 y5=0.750000 z15=5.199268349094e-002x1=0.080000 y6=0.800000 z16=-4.346804020490e-002x1=0.080000 y7=0.850000 z17=-1.316010567885e-001x1=0.080000 y8=0.900000 z18=-2.124310883088e-001x1=0.080000 y9=0.950000 z19=-2.860045510580e-001x1=0.080000 y10=1.000000 z110=-3.523860789794e-001x1=0.080000 y11=1.050000 z111=-4.116554565222e-001x1=0.080000 y12=1.100000 z112=-4.639049115188e-001x1=0.080000 y13=1.150000 z113=-5.092367247005e-001x1=0.080000 y14=1.200000 z114=-5.477611179623e-001x1=0.080000 y15=1.250000 z115=-5.795943883391e-001x1=0.080000 y16=1.300000 z116=-6.048572588895e-001x1=0.080000 y17=1.350000 z117=-6.236734213318e-001x1=0.080000 y18=1.400000 z118=-6.361682484133e-001x1=0.080000 y19=1.450000 z119=-6.424676566901e-001x1=0.080000 y20=1.500000 z120=-6.426971026996e-001x2=0.160000 y0=0.500000 z20=8.400813957666e-001x2=0.160000 y1=0.550000 z21=6.997641656732e-001x2=0.160000 y2=0.600000 z22=5.660614423517e-001x2=0.160000 y3=0.650000 z23=4.391716081176e-001x2=0.160000 y4=0.700000 z24=3.192421380408e-001x2=0.160000 y5=0.750000 z25=2.063761923874e-001x2=0.160000 y6=0.800000 z26=1.006385238914e-001x2=0.160000 y7=0.850000 z27=2.060740067837e-003x2=0.160000 y8=0.900000 z28=-8.935402476698e-002x2=0.160000 y9=0.950000 z29=-1.736269688648e-001x2=0.160000 y10=1.000000 z210=-2.507999561599e-001x2=0.160000 y11=1.050000 z211=-3.209322694446e-001x2=0.160000 y12=1.100000 z212=-3.840977350046e-001x2=0.160000 y13=1.150000 z213=-4.403821754175e-001x2=0.160000 y14=1.200000 z214=-4.898811523126e-001x2=0.160000 y15=1.250000 z215=-5.326979655338e-001x2=0.160000 y16=1.300000 z216=-5.689418792921e-001x2=0.160000 y17=1.350000 z217=-5.987265495151e-001x2=0.160000 y18=1.400000 z218=-6.221686297503e-001x2=0.160000 y19=1.450000 z219=-6.393865356972e-001x2=0.160000 y20=1.500000 z220=-6.504993507878e-001x3=0.240000 y0=0.500000 z30=1.051515091803e+000x3=0.240000 y1=0.550000 z31=9.029274308310e-001x3=0.240000 y2=0.600000 z32=7.605802668596e-001x3=0.240000 y3=0.650000 z33=6.247151981456e-001x3=0.240000 y4=0.700000 z34=4.955197560009e-001x3=0.240000 y5=0.750000 z35=3.731340427746e-001x3=0.240000 y6=0.800000 z36=2.576567488723e-001x3=0.240000 y7=0.850000 z37=1.491505594102e-001x3=0.240000 y8=0.900000 z38=4.764698677337e-002x3=0.240000 y9=0.950000 z39=-4.684932320146e-002x3=0.240000 y10=1.000000 z310=-1.343567603849e-001x3=0.240000 y11=1.050000 z311=-2.149133449274e-001x3=0.240000 y12=1.100000 z312=-2.885737006348e-001x3=0.240000 y13=1.150000 z313=-3.554063647857e-001x3=0.240000 y14=1.200000 z314=-4.154913964886e-001x3=0.240000 y15=1.250000 z315=-4.689182499695e-001x3=0.240000 y16=1.300000 z316=-5.157838831247e-001x3=0.240000 y17=1.350000 z317=-5.561910752001e-001x3=0.240000 y18=1.400000 z318=-5.902469305629e-001x3=0.240000 y19=1.450000 z319=-6.180615482412e-001x3=0.240000 y20=1.500000 z320=-6.397468392579e-001x4=0.320000 y0=0.500000 z40=1.271246751483e+000x4=0.320000 y1=0.550000 z41=1.115002018147e+000x4=0.320000 y2=0.600000 z42=9.646077272157e-001x4=0.320000 y3=0.650000 z43=8.203473694751e-001x4=0.320000 y4=0.700000 z44=6.824476781795e-001x4=0.320000 y5=0.750000 z45=5.510852085975e-001x4=0.320000 y6=0.800000 z46=4.263923859018e-001x4=0.320000 y7=0.850000 z47=3.084629956332e-001x4=0.320000 y8=0.900000 z48=1.973571296919e-001x4=0.320000 y9=0.950000 z49=9.310562085941e-002x4=0.320000 y10=1.000000 z410=-4.285992234033e-003x4=0.320000 y11=1.050000 z411=-9.483392529689e-002x4=0.320000 y12=1.100000 z412=-1.785729903640e-001x4=0.320000 y13=1.150000 z413=-2.555537790546e-001x4=0.320000 y14=1.200000 z414=-3.258401501575e-001x4=0.320000 y15=1.250000 z415=-3.895069883634e-001x4=0.320000 y16=1.300000 z416=-4.466382045995e-001x4=0.320000 y17=1.350000 z417=-4.973249517677e-001x4=0.320000 y18=1.400000 z418=-5.416640326994e-001x4=0.320000 y19=1.450000 z419=-5.797564797951e-001x4=0.320000 y20=1.500000 z420=-6.117062881476e-001x5=0.400000 y0=0.500000 z50=1.498321052482e+000x5=0.400000 y1=0.550000 z51=1.334998632066e+000x5=0.400000 y2=0.600000 z52=1.177125123739e+000x5=0.400000 y3=0.650000 z53=1.025024055020e+000x5=0.400000 y4=0.700000 z54=8.789600231744e-001x5=0.400000 y5=0.750000 z55=7.391451087035e-001x5=0.400000 y6=0.800000 z56=6.057448714871e-001x5=0.