北航数值分析实习第三题

数值分析计算实习报告第三题所在班级A1班学生姓名学生学号2015年11月北航数值分析计算实习报告第1页1算法设计方案1.1求对应于ix，jy的(,)ijijzfxy1.1.1通过求解非线性方程组获得x，y与t，u的关系题目中给出的方程组建立了t，u和x，y的联系，而数表建立了t，u与z的联系。为了求得x，y与z的联系，我们先求解方程组，获得t，u和x，y的联系。将0.08(0,1,,10)ixii和0.50.05(0,1,,20)jyjj代入非线性方程组中，用Newton法解出ijt和iju。1.1.2通过分片二次代数插值求得t，u和z的关系，进而得到x，y和z的关系上一步，已经获得对应于ix，jy的ijt，iju，现在只要再获得对应于ijt，iju的z，就间接求得了对应于ix，jy的ijz。根据题目给出的数表，进行分片二次代数插值，之后带入ijt，iju，获得了(,)ijijzfxy，它对应于ix，jy。1.2曲面插值上一步建立了数表(,,(,))ijijxyfxy,(0,1,,10;0,1,,20)ij，利用该数表进行曲面插值，本质上就是求得系数矩阵rsc，就可获得(,)iipxy。k从1开始尝试，每个k都进行曲面插值，满足精度要求时停止计算。1.3观察),(jiyxp逼近(,)iifxy的效果这里新给了点(*,*)iixy，*0.1(1,2,,8)ixii和*0.50.2(1,2,,5)jyjj将其重复上面步骤，得到与新的插值节点(*,*)iixy对应的(*,*)iifxy。再将(*,*)iixy带入上一步求得的(,)iipxy即可得到。比较时可以比较相对误差和绝对误差。1.4补充—计算逆矩阵计算系数矩阵rsc时，需要求解一些矩阵的逆矩阵。这里采用列主元Gauss消去法来求解，即将原矩阵初等变换为单位阵，同样的变换可以把单位阵变为原矩阵的逆矩阵。北航数值分析计算实习报告第2页2C++程序#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)//求向量的无穷范数{inti=0;doublemax=fabs(p[1]);for(i=1;i=4;i++){if(fabs(p[i])max)max=fabs(p[i]);}return(max);}voidinverse(doubleX[10][10],intN)//采用列主元Gauss消去法求逆矩阵{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++)北航数值分析计算实习报告第3页{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++)//把上对角矩阵化为单位阵的相同步骤就可把经过变形的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];}}}北航数值分析计算实习报告第4页voidnewton(doublex[11],doubley[11])//牛顿法求解非线性方程组{doublet=0;doubleu=0;doublew=0;doublev=0;doubledF[5][5]={0};//F的导数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);//这里实际是-FF[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);北航数值分析计算实习报告第5页for(k=1;k=3;k++)//列主元Gauss消去法求解dF*△x=-F{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;北航数值分析计算实习报告第6页jie[2]=u;jie[3]=v;jie[4]=w;}while(fanshu(d)/fanshu(jie)1e-12);//引用求范数的函数}}voideryuanchazhi(doublete[11][21],doubleue[11][21])//分片二次代数插值求z=f(x,y){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)//选择t的插值点，考察是边界还是中间的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;}if(ue[i][j]=0.6)//选择u的插值点，考察是边界还是中间的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;北航数值分析计算实习报告第7页}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++)//这两个循环是求l的求积层次的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-xyf.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()//曲面插值，求p(x,y){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};//B转置*BdoubleGG[10][10]={0};//G转置*GdoubleBBN[10][10]={0};//(B转置*B)的逆doubleGGN[10][10]={0};//(G转置*G)的逆doubleBBB[10][12]={0};//[(B转置*B)的逆]*B的转置doubleBBBU[10][22]={0};//{[(B转置*B)的逆]*B的转置}*UdoubleBBBUG[10][10]={0};//{[(B转置*B)的逆]*B的转置}*U*Gdoublep=0;北航数值分析计算实