数值分析第三次大作业一、 算法的设计方案：（一）、总体方案设计：（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)