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

数值分析第三次大作业一、算法的设计方案1、求解非线性方程组将题目中给出的(,)iixy当作已知量代入题目给定的非线性方程组，求出与(,)iixy相对应的数组te[i][j]，ue[i][j]，此处采用的是牛顿法解非线性方程组，其算法如书上91页所示。2、分片二次代数插值对所求出的数组te[i][j]，ue[i][j]，通过分片二次代数插值运算，得到与数组te[11][21]，ue[11][21]对应的数组ze[11][21]，从而得到二元函数z=(,)iifxy，此处采用如书上101页例2中所示的分片二次代数插值。3、曲面插值利用x[11]，y[21]，ze[11][21]建立二维函数表，进行曲面插值计算，逐步提高k值，计算其精度，看其是否满足要求，以此来确定循环结束的时刻，并得到曲面拟合的系数矩阵C[r][s]，此处的算法如书142页所示，只需将所需矩阵给出，然后按公式进行计算即可。4、比较观察和),(jiyxp逼近(,)iifxy的效果。观察逼近效果只需要利用新给的点列(,)iixy重复上面（1）和（2）的过程，得到与新的插值节点(,)iixy对应的(,)iifxy，再与对应的(,)iipxy比较即可，这里求解(,)iipxy可以直接使用（3）中的C[r][s]和k。5、几点说明分片二次插值的结果x[i]，y[j]，ze[i][j]输出到一个文件shubiao.txt中，方便结果的复制与粘贴。曲面插值的结果输出到一个文件xishu.txt中，包括循环中每一次的k值以及误差平方和sigma的值，还有最后满足误差要求时曲面插值的系数C[r][s]。观察逼近效果的结果输出到一个文件shubiao1.txt中，方便结果的复制与粘贴。二、源程序如下：#includestdio.h#includemath.hvoidnewton(doublex[11],doubley[11]);//牛顿法求解非线性方程组voideryuanchazhi(doublete[11][21],doubleue[11][21]);//分片二次代数插值，求z=f(x,y)voidqumianchazhi();//曲面插值，求p(x,y)voidcompare();//比较f(x*,y*)与p(x*,y*)voidinverse(doubleX[10][10],intN);//求逆矩阵doublefanshu(double*p);//求向量的无穷范数doublete[11][21]={0};//存储非线线方程组的解tdoubleue[11][21]={0};//存储非线线方程组的解udoubleze[11][21]={0};//存储分片二次插值的解zdoublex[11]={0};//xdoubley[21]={0};//ydoubleinv[10][10]={0};//存储某矩阵对应的逆矩阵doubleC[10][10]={0};//存储曲面插值的系数矩阵doublers=0;//存储k值，用于比较函数中使用voidmain(){inti,j;for(i=0;i11;i++)for(j=0;j21;j++){x[i]=0.08*i;y[j]=0.5+0.05*j;}newton(x,y);eryuanchazhi(te,ue);qumianchazhi();compare();}voidnewton(doublex[11],doubley[11]){doublet=0;doubleu=0;doublew=0;doublev=0;doubledF[5][5]={0};doubleF[5]={0};doubled[5]={0};inti,j,k;intik;intcnt,count;doubletemp=0;doublem=0;doublejie[4]={0};for(i=0;i11;i++)for(j=0;j21;j++){t=1;u=1;w=1;v=1;//初始化//===============求解非线性方程组=============//do{dF[1][1]=-0.5*sin(t);dF[1][2]=1;dF[1][3]=1;dF[1][4]=1;dF[2][1]=1;dF[2][2]=0.5*cos(u);dF[2][3]=1;dF[2][4]=1;dF[3][1]=0.5;dF[3][2]=1;dF[3][3]=-sin(v);dF[3][4]=1;dF[4][1]=1;dF[4][2]=0.5;dF[4][3]=1;dF[4][4]=cos(w);F[1]=-(0.5*cos(t)+u+v+w-x[i]-2.67);F[2]=-(t+0.5*sin(u)+v+w-y[j]-1.07);F[3]=-(0.5*t+u+cos(v)+w-x[i]-3.74);F[4]=-(t+0.5*u+v+sin(w)-y[j]-0.79);//对Newton法的系数矩阵进行赋值,注意这里是-F()//=========================列主元高斯消元法求解方程组====================//for(k=1;k=3;k++){ik=k;for(cnt=k;cnt=4;cnt++){if(fabs(dF[cnt][k])fabs(dF[ik][k])){ik=cnt;}}temp=F[k];F[k]=F[ik];F[ik]=temp;for(cnt=k;cnt=4;cnt++){temp=dF[k][cnt];dF[k][cnt]=dF[ik][cnt];dF[ik][cnt]=temp;}for(cnt=k+1;cnt=4;cnt++){m=dF[cnt][k]/dF[k][k];for(count=k+1;count=4;count++)dF[cnt][count]=dF[cnt][count]-m*dF[k][count];F[cnt]=F[cnt]-m*F[k];}}d[4]=F[4]/dF[4][4];for(k=4-1;k=1;k--){temp=0;for(cnt=k+1;cnt=4;cnt++)temp+=dF[k][cnt]*d[cnt];d[k]=(F[k]-temp)/dF[k][k];}t=t+d[1];u=u+d[2];v=v+d[3];w=w+d[4];te[i][j]=t;ue[i][j]=u;jie[1]=t;jie[2]=u;jie[3]=v;jie[4]=w;}while(fanshu(d)/fanshu(jie)1e-12);printf(te[%d][%d]=%13.11e\t,i,j,te[i][j]);printf(ue[%d][%d]=%13.11e\n,i,j,ue[i][j]);}}voideryuanchazhi(doublete[11][21],doubleue[11][21]){inti,j,ii,jj,m;intr[11][21]={0};intk[11][21]={0};doubletemp=0;doubleu[6]={0,0.4,0.8,1.2,1.6,2};doublet[6]={0,0.2,0.4,0.6,0.8,1.0};doublez[6][6]={{-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}};for(i=0;i=10;i++)for(j=0;j=20;j++){if(te[i][j]=0.3)r[i][j]=1;elseif(te[i][j]=0.7)r[i][j]=4;else{for(m=2;m=3;m++)if((te[i][j]0.2*m-0.1)&&(te[i][j]0.2*m+0.1))r[i][j]=m;}//t的插值点if(ue[i][j]=0.6)k[i][j]=1;elseif(ue[i][j]=1.4)k[i][j]=4;else{for(m=2;m=3;m++)if((ue[i][j]0.4*m-0.2)&&(ue[i][j]0.4*m+0.2))k[i][j]=m;}//u的插值点ze[i][j]=0;for(ii=k[i][j]-1;ii=k[i][j]+1;ii++)for(jj=r[i][j]-1;jj=r[i][j]+1;jj++){temp=1;for(m=k[i][j]-1;m=k[i][j]+1;m++)if(m!=ii)temp*=(te[i][j]-t[m])/(t[ii]-t[m]);for(m=r[i][j]-1;m=r[i][j]+1;m++)if(m!=jj)temp*=(ue[i][j]-u[m])/(u[jj]-u[m]);ze[i][j]+=temp*z[ii][jj];}}FILE*fp=fopen(shubiao.txt,a);for(i=0;i=10;i++)for(j=0;j=20;j++)fprintf(fp,x[%d]=%f\ty[%d]=%f\tf(x,y)=%13.11e\n,i,x[i],j,y[j],ze[i][j]);fclose(fp);}voidqumianchazhi(){inti,j;intii,jj;intk,t;doubletemp=0;doubletemp1=0;doublesigma=0;doubleB[11+1][10]={0};doubleG[21+1][10]={0};doubleBB[10][10]={0};doubleGG[10][10]={0};doubleBBN[10][10]={0};doubleGGN[10][10]={0};doubleBBB[10][12]={0};doubleBBBU[10][22]={0};doubleBBBUG[10][10]={0};doublep=0;FILE*fp;fp=fopen(xishu.txt,a);k=0+1;//为了以后矩阵运算方便，这里先加上1do{for(i=1;i=11;i++)for(j=1;j=k;j++){B[i][j]=pow(x[i-1],j-1);}//B阵赋值for(i=1;i=21;i++)for(j=1;j=k;j++){G[i][j]=pow(y[i-1],j-1);}//G阵赋值for(i=1;i=k;i++)for(j=1;j=k;j++){temp=0;for(t=1;t=11;t++)temp+=B[t][i]*B[t][j];BB[i][j]=temp;}//Bt*Bfor(i=1;i=k;i++)for(j=1;j=k;j++){temp=0;for(t=1;t=21;t++)temp+=G[t][i]*G[t][j];GG[i][j]=temp;}//Gt*Ginverse(BB,k);for(i=1;i=k;i++)for(j=1;j=k;j++)BBN[i][j]=inv[i][j];//Bt*B的逆inverse(GG,k);for(i=1;i=k;i++)for(j=1;j=k;j++)GGN[i][j]=inv[i][j];//Gt*G的逆for(i=1;i=k;i++)for(j=1;j=11;j++){temp=0;for(t=1;t=k;t++)temp+=BBN[i][t]*B[j][t];BBB[i][j]=temp;}//Bt*B的逆*Btfor(i=1;i=k;i++)for(j=1;j=21;j++){temp=0;for(t=1;t=11;t++)temp+=BBB[i][t]*