数值分析报告第三大题

《数值分析》计算实习报告第三大题学号：DY1305姓名：指导老师：一、题目要求1、试用数值方法求出f(x，y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个近似表达式,0(,)krsrsrspxycxy要求p(x,y)一最小的k值达到以下的精度10202700((,)(,))10ijijijfxypxy其中xi=0.08i，yj=0.5+0.05j。2、计算f(xi*,yj*),p(xi*,yj*)(i=1,2,…,8;j=1,2,…,5)的值，以观察p(x,y)逼近f(x,y)的效果，其中xi*=0.1i,yj*=0.5+0.2j。二、算法设计方案将给出的(,)iixy代入给定的非线性方程组中，采用是牛顿法解非线性方程组求出与(,)iixy对应的数组te[i][j]，ue[i][j]。再对求出的数组te[i][j]，ue[i][j]，通过分片二次代数插值运算得到与数组te[11][21]，ue[11][21]对应的数组ze[11][21]，从而得到二元函数(,)iizfxy，得到的x[i]，y[j]，ze[i][j]输出到一个文件table1.txt中。利用x[11]，y[21]，ze[11][21]建立二维函数表，进行曲面插值计算。逐步提高k值，计算其精度，看其是否满足要求，以此来确定循环结束的时刻，并得到曲面拟合的系数矩阵C[r][s]。循环中每一次的k值、误差平方和σ的值及系数C[r][s]输出到一个文件xishu.txt中。利用新给的点列(,)iixy重复上面的过程，得到与新的插值节点(,)iixy对应的(,)iifxy，再与对应的(,)iipxy比较，即可观测(,)iipxy逼近(,)iifxy的效果。逼近效果的结果输出到一个文件table2.txt中。这里求解(,)iipxy可以直接使用3中的C[r][s]和k。三、源程序#includestdio.h#includemath.hdoublete[11][21]={0};doubleue[11][21]={0};doubleze[11][21]={0};/doublex[11]={0};doubley[21]={0};doubleinv[10][10]={0};doubleC[10][10]={0};doublers=0;doublefanshu(double*p);//求向量的无穷范数voidcompare();//比较f(x*,y*)与p(x*,y*)voidinverse(doubleX[10][10],intN);//求逆矩阵voidnewton(doublex[11],doubley[11]);//牛顿法求解非线性方程组voideryuanchazhi(doublete[11][21],doubleue[11][21]);//分片二次代数插值，求z=f(x,y)voidqumianchazhi();//曲面插值，求p(x,y)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();}doublefanshu(double*p){inti=0;doublemax=fabs(p[1]);for(i=1;i=4;i++){if(fabs(p[i])max)max=fabs(p[i]);}return(max);}voidcompare(){doublepn[8+1][5+1]={0};inti,j;doubletemp=0;intii,jj;for(i=1;i=8;i++)x[i]=0.1*i;for(j=1;j=5;j++)y[j]=0.5+0.2*j;for(i=1;i=8;i++)for(j=1;j=5;j++){temp=0;for(ii=0;ii=rs;ii++)for(jj=0;jj=rs;jj++)temp+=C[ii+1][jj+1]*pow(x[i],ii)*pow(y[j],jj);pn[i][j]=temp;}newton(x,y);eryuanchazhi(te,ue);FILE*fp;fp=fopen(table2.txt,a);for(i=1;i=8;i++)for(j=1;j=5;j++){fprintf(fp,x*(%d)=%f\ty*(%d)=%f\t\n,i,x[i],j,y[j]);fprintf(fp,f(x*,y*)=%13.11e\np(x*,y*)=%13.11e\t\n,ze[i][j],pn[i][j]);fprintf(fp,σ=%13.11e\n\n,ze[i][j]-pn[i][j]);}}voidinverse(doubleX[10][10],intN){doublematrix[10][10];doubleI[10][10]={0};inti,j;intk,ik,cnt,count;doublem;doubletemp;for(i=1;i=N;i++)for(j=1;j=N;j++)matrix[i][j]=X[i][j];for(i=1;i=N;i++)I[i][i]=1;//单位阵for(k=1;k=N-1;k++)//列主元的高斯消元法求逆矩阵{ik=k;for(cnt=k;cnt=N;cnt++){if(fabs(matrix[cnt][k])fabs(matrix[ik][k])){ik=cnt;}}for(cnt=1;cnt=N;cnt++){temp=matrix[k][cnt];matrix[k][cnt]=matrix[ik][cnt];matrix[ik][cnt]=temp;temp=I[k][cnt];I[k][cnt]=I[ik][cnt];I[ik][cnt]=temp;}for(cnt=k+1;cnt=N;cnt++){m=matrix[cnt][k]/matrix[k][k];for(count=1;count=N;count++){matrix[cnt][count]=matrix[cnt][count]-m*matrix[k][count];I[cnt][count]=I[cnt][count]-m*I[k][count];}}}for(i=1;i=N;i++){inv[N][i]=I[N][i]/matrix[N][N];for(k=N-1;k=1;k--){temp=0;for(cnt=k+1;cnt=N;cnt++)temp+=matrix[k][cnt]*inv[cnt][i];inv[k][i]=(I[k][i]-temp)/matrix[k][k];}}}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(j