数值分析第三次作业.docxVIP

数值分析第三次作业一算法设计方案欲获得z=f(x,y)的值，需要根据已知方程组和数表，把每一对待求的x和y值代入方程组，得到X=(t,u,v,w)的非线性方程组，用牛顿法求解所得的非线性方程组，得到解向量X=(t,u,v,w)；对数表z、t、u进行分片二次插值，把解向量X=(t,u,v,w)中的t,u代入插值公式，即可求得z值。由（1）和（2）即可以得到数表x、y、f(x,y)，此数表由C程序输出后导入Excel中打印出来。对已求得的数表进行曲线拟合，方法采用教材中的曲线拟合方法，实质上就是矩阵的计算求解过程。拟合基函数取(x)=， (k=0,1…K，l=0,1…K)，K的初始值设为1，每循环一次K增加1，同时判断sigma值是否达到1e-7，如果达到精度要求，输出K和sigma的值以及此时的系数矩阵Crs，此时就得到了拟合曲面方程。校验过程：重新给定一组和的值，需要重复（1）和（2）过程，得到f(,)，再由步骤（4）所确定的函数求得拟合值P(,)，把两组值做一个对比，可以看出拟合结果的好坏。此数表仍然把C程序得到的数据输出导入Excel中打印出来，这样便于观察。在具体编写程序时没有采用过多的子程序，牛顿法求解非线性方程组、列主元素Gauss消去法、曲线拟合等程序均含在主程序中，只是把分片二次插值作为一个子程序进行循环调用。同时步骤（5）中的校验过程均是采用前面的程序计算，只是x和y的数据稍有不同，因此为了避免重复没有把具体的计算写入主程序中，而只是给出了最后的Excel数表以便观察拟合结果的好坏。程序虽然有点繁琐，但是易懂！！！二 C语言程序# includestdio.h# includemath.h# define ee 1e-12# define M 100# define N 4main (){ double f(double t,double u);int i,j,r,s,tt,k,k1,k2,row,kk,p; double x[11],y[21],X[N],dertx[N],F1[N][N],F[N]; double max,l1,m,mdx,mx,z,sigma,a1[11][11],b1[11],x2[11]; double Z[11][21],B[11][11],G[21][11],a[11][11],b[11],A[11][21],D[11][21]; double C[11][11],x1[11],fai[11],psai[11],H[11][11],ZZ[11][21];for(i=0;i11;i++)x[i]=0.08*i;for(i=0;i21;i++)y[i]=0.5+0.05*i; /*牛顿法解非线性方程组*/for(i=0;iN;i++) X[i]=1; for(r=0;r11;r++)for(s=0;s21;s++) {for(kk=0;kkM;kk++) {F1[0][0]=-0.5*sin(X[0]);F1[0][1]=1;F1[0][2]=1;F1[0][3]=1; F1[1][0]=1;F1[1][1]=0.5*cos(X[1]);F1[1][2]=1;F1[1][3]=1; F1[2][0]=0.5;F1[2][1]=1;F1[2][2]=-sin(X[2]);F1[2][3]=1; F1[3][0]=1;F1[3][1]=0.5;F1[3][2]=1;F1[3][3]=cos(X[3]); F[0]=-(0.5*cos(X[0])+X[1]+X[2]+X[3]-x[r]-2.67);F[1]=-(X[0]+0.5*sin(X[1])+X[2]+X[3]-y[s]-1.07);F[2]=-(0.5*X[0]+X[1]+cos(X[2])+X[3]-x[r]-3.74);F[3]=-(X[0]+0.5*X[1]+X[2]+sin(X[3])-y[s]-0.79);for(k=0;kN-1;k++){ max=fabs(F1[k][k]);row=k;for(i=k;iN;i++)if(fabs(F1[i][k])max){ max=fabs(F1[i][k]);row=i;}if(row!=k){ for(i=k;iN;i++){ l1=F1[k][i];F1[k][i]=F1[row][i];F1[row][i]=l1; } l1=F[k]; F[k]=F[row]; F[row]=l1; } for(i=k+1;iN;i++){ m=F1[i][k]/F1[k][k];for(j=k+1;jN;j++) F1[i][j