计算实习报告3.doc

1、算法设计：整体思路这一次的作业和前面两次相比，感觉要难很多，最明显的就是书上没有完整的算法实现流程。再加上插值和逼近这一段在课上也没听太懂什么意思，所以刚开始思考这个程序时感觉很吃力。不过经过后来好好看书和向同学请教，渐渐对插值和分片元二次拟合有点懂了，也大概知道了解决本次问题分别要完成的子问题。要完成本次任务，主体上有三部分问题要解决：（1）解非线性方程组。将d内的)当作已知量代入题目给定的非线性方程组，得到与相对应的tij,uij。（2）分片双二次插值。采用分片双二次插值，得到与t1121,u1121对应的z1121，也即建立与的对应关系，得到二元函数。（3）曲面拟合。利用xi,yj,z1121建立二维函数表，再根据精度的要求选择适当k值，并得到曲面拟合的系数矩阵ckk。（4）观察逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点对应的，再与对应的比较即可。程序示意图如下所示：各部分问题的算法实现（1）解非线性方程组非线性方程组常用数值解法有简单迭代法和牛顿法。牛顿法收敛快，一般都能达到平方收敛，因此这里选择newton法解非线性方程组。牛顿法解方程组的解，可采用如下算法：1）在附近选取，给定精度水平和最大迭代次数m。2）对于执行 计算和。 求解关于的线性方程组 若，则取，并停止计算；否则转。 计算。 若，则继续，否则，输出m次迭代不成功的信息，并停止计算。注：第步中的求取，用到了解线性方程组，使用的是列主元素gauss消去法。 另外，本次作业中解非线性方程组实际求取的是解向量和，方程组中的、 是已知值，为当前节点的值。（2）分片双二次插值给定已知数表以及需要插值的节点，进行分片二次插值的算法：设已知数表中的点为： ，需要插值的节点为。1) 根据选择插值节点：若或，插值节点对应取或，若或，插值节点对应取或。 若则选择为插值节点。2）计算 插值多项式的公式为： 注：在（1）中通过解非线性方程组得到的解向量和与插值节点有着一一对应关系，而本题题目给出的函数表为关于的函数关系。因此，为得到关于的函数关系，本题中应先根据给定数表对进行插值，再利用与的一一对应关系得到与的对应关系。（3）曲面拟合根据插值得到的数表进行曲面拟合的过程：1） 根据拟合节点和基底函数写出矩阵b和g： 2） 计算 。在这里，为了简化计算和编程、避免矩阵求逆，记：，对上面两式进行变形，得到如下两个线性方程组：，通过解上述两个线性方程组，则有：3） 对于每一个， 。4） 拟合需要达到的精度条件为： 。 其中对应着插值得到的数表中的值。5） 让k逐步增加，每一次重复执行以上几步，直到 成立。此时的k值就是所需最小的k。注：曲面拟合过程中用到的数表为前面插值得到的数表，即。在2)中，求矩阵a和d的过程用的是列主元素gauss消去法。（4）观察逼近效果。观察逼近效果主要要做的是，通过新的插值节点 建立新的插值数表，同时求出对应的，比较即可。2、程序源代码：#include stdio.h#include conio.h#include math.h#define epsilon1 1e-12 /解线性方程组时近似解向量的精度#define maxm 200 /解线性方程组时的最大迭代次数#define k 10 /预估的能达到精度要求的k值，k=k#define givenx 11 /已知要求的插值节点xi个数#define giveny 21 /已知要求的插值节点yi个数#define testx 8 /观察逼近效果时，插值节点xi个数#define testy 5 /观察逼近效果时，插值节点yi个数/*/* 函数功能描述： */*/* norm()：求向量x的无穷范数 */* newton_nonlinear()：牛顿法解非线性方程组的主干部分 */* gausselemitation_select()：线性方程组的列主元素gauss消去法 */* interplate()：分片二元二次插值 */* surfacefit()：曲面拟合，函数返回满足精度要求的最小k值 */* calculate_a()：根据拟合曲面的系数矩阵c=a(dt)，求矩阵a */* calculate_d()：根据拟合曲面的系数矩阵c=a(dt)，求矩阵d */* test_observe()：观察逼近效果子函数 */*/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,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：插值节点(xi,yi) */* t1121、u1121：与节点(xi,yi)对应的二元组(tij,uij) */* gausselemitation_select()：线性方程组的列主元素gauss消去法 */* z1121：存放插值得到的数表 */* ckk：存放曲面拟合得到的系数矩阵，其中k是预估的最小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)的对应关系 for(j=0;j=giveny-1;j+) xi=0.08*i; yj=0.5+0.05*j; xstar0=1; xstar1=1; xstar2=1; xstar3=1; /选取迭代初始值x(0) flag=newton_nonlinear(xstar,xi,yj); /牛顿法解非线性方程组 if(flag=-1) goto end; /如果非线性方程组求解失败，中止程序 tij=xstar0; uij=xstar1; /*/* 插值并输出插值得到的数表 */*/ interplate(t,u,z,11,21); /分片二次插值 printf(1、插值得到的数表为：n); for(i=0;i=10;i+) /输出插值得到的数表：xi,yi,f(xi,yi) for(j=0;j=20;j+) printf(x%d= %.6f ,i,xi); printf(y%d= %.6f ,j,yj); printf(z%d%d= %.12en,i,j,zij); printf(n);/*/* 曲面拟合并输出拟合得到的矩阵c */*/ k=surfacefit(z,x,y,c); /曲面拟合，返回满足精度要求的最小k值 printf(系数矩阵c为：n); for(i=0;i=k;i+) /输出曲面拟合得到的系数矩阵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、观察插值的逼近效果：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+xstar3-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; df32=1; df33=cos(xstar3); gausselemitation_select(df,f,delta_x); /列主元素gauss消去法解线性方程组 if(norm(delta_x)/norm(xstar)=epsilon1) /达到精度要求，跳出 break; for(l=0;l=maxm) /迭代次数超限，返回-1作为后面程序是否执行的标志 printf(迭代次数k=%d超出限制！n,k); return(-1); /* 列主元素gauss消去法解线性方程组 */void gausselemitation_select(double df44,double f4,double 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; /*消元*/ 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的无穷范数 */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,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; /已知数表 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; /* 选择插值节点，xi、yi分开选择 */ for(i=0;i=numi-1;i+) /选择插值节点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+) /选择插值节点yi for(j=0;j=numj-1;j+) if(uij1.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; /初始化节点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(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(xi,i1)*pow(yj,j1); /end p(xi,yj) sigma+=pow(pij-zij,2); /end sigma printf(k=%d =%.12en,k,sigma); /输出当前k值和 k+; k-; printf(曲面拟合达到精度要求，此时：n); printf(k=%d =%.12en,k,sigma); return k;/*/* 函数名：calculate_a */*/* 功 能：根据公式 (bt)b)a=(bt)u，通过解线性方程组得到矩阵a */*/void calculate_a(double ak21,double z1121,double x11,int k) int i,j,m,l,line_flag; double b11k,btk11,btbkk,btuk21; /b和bt的阶数分别为11xk和kx11 double btb_1kk,temp; /a的阶数为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相乘，并将结果存在矩阵btb中 */ for(i=0;i=k;i+) for(j=0;j=k;j+) temp=0; for(m=0;m=10;m+) temp+=btim*bmj; btbij=temp; /end btb /* 计算矩阵bt和u相乘，并将结果存在矩阵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+) /*复制矩阵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=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-) temp=0; for(i=m+1;i=k;i+) temp+=ail*btb_1mi; aml=(btuml-temp)/btb_1mm; /end 回代 /end 求矩阵a/*/* 函数名：calculate_d */*/* 功 能：根据公式 (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的阶数分别为11xk和kx11 double gtg_1kk,temp; /g的阶数为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相乘，并将结果存在矩阵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+) /*复制矩阵gtg*/ for(i=0;i=k;i+) for(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_flagl=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(i=0;i=testx-1;i+) /通过解非线性方程组，建立(x,y)与(t,u)的对应关系 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)对应的f(xi,yi)的值 for(i=0;i=testx-1;i+) /求与(xi,yi)对应的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); printf(tt=%.12en,zij-pij); 3、程序运行结果：1、插值得到的数表为：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.750000 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=-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=