北航数值分析实习第三题.docx

1、北航数值分析实习第三题数值分析计算实习报告第三题所在班级A1 班学生姓名学生学号2015年11月1算法设计方案1.1 求对应于，的1.1.1 通过求解非线性方程组获得x，y与t，u的关系题目中给出的方程组建立了t，u和x，y的联系，而数表建立了t，u与z的联系。为了求得x，y与z的联系，我们先求解方程组，获得t，u和x，y的联系。将和代入非线性方程组中，用Newton法解出和。1.1.2 通过分片二次代数插值求得t，u和z的关系，进而得到x，y和z的关系上一步，已经获得对应于，的，现在只要再获得对应于，的z，就间接求得了对应于，的。根据题目给出的数表，进行分片二次代数插值，之后带入，获得了，它

2、对应于，。1.2 曲面插值上一步建立了数表,，利用该数表进行曲面插值，本质上就是求得系数矩阵，就可获得。k从1开始尝试，每个k都进行曲面插值，满足精度要求时停止计算。1.3 观察逼近的效果这里新给了点，和将其重复上面步骤，得到与新的插值节点对应的。再将带入上一步求得的即可得到。比较时可以比较相对误差和绝对误差。1.4 补充计算逆矩阵计算系数矩阵时，需要求解一些矩阵的逆矩阵。这里采用列主元Gauss消去法来求解，即将原矩阵初等变换为单位阵，同样的变换可以把单位阵变为原矩阵的逆矩阵。2C+程序#include stdio.h#include math.hdouble te1121=0;double

3、 ue1121=0; double ze1121=0; double x11=0; double y21=0; double inv1010=0; double C1010=0; double rs=0; double fanshu(double *p)/求向量的无穷范数 int i=0; double max=fabs(p1); for(i=1;imax) max=fabs(pi); return(max);void inverse(double X1010,int N)/采用列主元Gauss消去法求逆矩阵 double matrix1010; double I1010=0; int i,j

4、; int k,ik,cnt,count; double m; double temp; for(i=1;i=N;i+) for(j=1;j=N;j+) matrixij=Xij; for(i=1;i=N;i+) Iii=1;/单位阵 for (k=1;k=N-1;k+) ik=k; for(cnt=k;cntfabs(matrixikk) ik=cnt; for (cnt=1;cnt=N;cnt+)/交换 temp=matrixkcnt; matrixkcnt=matrixikcnt; matrixikcnt=temp; temp=Ikcnt; Ikcnt=Iikcnt; Iikcnt=te

5、mp; for (cnt=k+1;cnt=N;cnt+)/消去了 m=matrixcntk/matrixkk; for (count=1;count=N;count+) matrixcntcount=matrixcntcount-m*matrixkcount; Icntcount=Icntcount-m*Ikcount; for(i=1;i=1;k-) temp=0; for (cnt=k+1;cnt=N;cnt+) temp+=matrixkcnt*invcnti; invki=(Iki-temp)/matrixkk; void newton(double x11,double y11)/牛

6、顿法求解非线性方程组 double t=0; double u=0; double w=0; double v=0; double dF55=0;/F的导数 double F5=0; double d5=0; int i,j,k; int ik; int cnt,count; double temp=0; double m=0; double jie4=0; for(i=0;i11;i+) for(j=0;j21;j+) t=1; u=1; w=1; v=1; do dF11=-0.5*sin(t); dF12=1; dF13=1; dF14=1; dF21=1; dF22=0.5*cos(u

7、); dF23=1; dF24=1; dF31=0.5; dF32=1; dF33=-sin(v); dF34=1; dF41=1; dF42=0.5; dF43=1; dF44=cos(w); F1=-(0.5*cos(t)+u+v+w-xi-2.67);/这里实际是-F F2=-(t+0.5*sin(u)+v+w-yj-1.07); F3=-(0.5*t+u+cos(v)+w-xi-3.74); F4=-(t+0.5*u+v+sin(w)-yj-0.79); for (k=1;k=3;k+)/列主元Gauss消去法求解dF*x=-F ik=k; for(cnt=k;cntfabs(dFik

8、k) ik=cnt; temp=Fk; Fk=Fik; Fik=temp; for (cnt=k;cnt=4;cnt+) temp=dFkcnt; dFkcnt=dFikcnt; dFikcnt=temp; for (cnt=k+1;cnt=4;cnt+) m=dFcntk/dFkk; for (count=k+1;count=1;k-) temp=0; for (cnt=k+1;cnt1e-12);/引用求范数的函数 void eryuanchazhi(double te1121,double ue1121)/分片二次代数插值求z=f(x,y) int i,j,ii,jj,m; int r1

9、121=0; int k1121=0; double temp=0; double u6=0,0.4,0.8,1.2,1.6,2; double t6=0,0.2,0.4,0.6,0.8,1.0; double z66=-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,

10、-0.5; for(i=0;i=10;i+) for(j=0;j=20;j+) if(teij=0.7) rij=4; else for(m=2;m0.2*m-0.1)&(teij0.2*m+0.1) rij=m; if(ueij=1.4) kij=4; else for(m=2;m0.4*m-0.2)&(ueij0.4*m+0.2) kij=m; zeij=0; for(ii=kij-1;ii=kij+1;ii+) for(jj=rij-1;jj=rij+1;jj+)/这两个循环是求和层次的 temp=1; for(m=kij-1;m=kij+1;m+)/这两个循环是求l的求积层次的 if(

11、m!=ii) temp*=(teij-tm)/(tii-tm); for(m=rij-1;m=rij+1;m+) if(m!=jj) temp*=(ueij-um)/(ujj-um); zeij+=temp*ziijj; FILE *fp=fopen(shubiao-xyf.txt,a); for(i=0;i=10;i+) for(j=0;j=20;j+) fprintf(fp,x%d=%fty%d=%ftf(x,y)=%13.11en,i,xi,j,yj,zeij); fclose(fp);void qumianchazhi()/曲面插值，求p(x,y) int i,j; int ii,jj

12、; int k,t; double temp=0; double temp1=0; double sigma=0; double B11+110=0; double G21+110=0; double BB1010=0;/B转置*B double GG1010=0;/G转置*G double BBN1010=0;/(B转置*B)的逆 double GGN1010=0;/(G转置*G)的逆 double BBB1012=0;/(B转置*B)的逆*B的转置 double BBBU1022=0;/(B转置*B)的逆*B的转置*U double BBBUG1010=0;/(B转置*B)的逆*B的转置*

13、U*G double p=0; FILE *fp; fp=fopen(shubiao-ksic.txt,a); k=0+1;/这里先加上1是为了矩阵运算方便 do /首先把系数矩阵C求出来 for(i=1;i=11;i+) for(j=1;j=k;j+) Bij=pow(xi-1,j-1); for(i=1;i=21;i+) for(j=1;j=k;j+) Gij=pow(yi-1,j-1); for(i=1;i=k;i+)/求B转置*B for(j=1;j=k;j+) temp=0; for(t=1;t=11;t+) temp+=Bti*Btj; BBij=temp; for(i=1;i=k

14、;i+)/求G转置*G for(j=1;j=k;j+) temp=0; for(t=1;t=21;t+) temp+=Gti*Gtj; GGij=temp; inverse(BB,k); for(i=1;i=k;i+)/（B转置*B）的逆 for(j=1;j=k;j+) BBNij=invij; inverse(GG,k); for(i=1;i=k;i+)/（G转置*G）的逆 for(j=1;j=k;j+) GGNij=invij; for(i=1;i=k;i+)/（B转置*B）的逆*B的转置 for(j=1;j=11;j+) temp=0; for(t=1;t=k;t+) temp+=BBN

15、it*Bjt; BBBij=temp; for(i=1;i=k;i+)/(B转置*B)的逆*B的转置*U for(j=1;j=21;j+) temp=0; for(t=1;t=11;t+) temp+=BBBit*zet-1j-1; BBBUij=temp; for(i=1;i=k;i+)/(B转置*B)的逆*B的转置*U*G for(j=1;j=k;j+) temp=0; for(t=1;t=21;t+) temp+=BBBUit*Gtj; BBBUGij=temp; for(i=1;i=k;i+)/得到系数矩阵C for(j=1;j=k;j+) temp=0; for(t=1;t=k;t+

16、) temp+=BBBUGit*GGNtj; Cij=temp; /下面判断现在的k对应的C是否满足精度要求 for(i=0;i11;i+) for(j=0;j21;j+) temp=0; for(ii=0;iik;ii+) for(jj=0;jj=1e-7); k-;/减去最后一次判断循环条件前加的1 rs=k;/存储k值，用于比较函数中使用 fprintf(fp,k=%dn,k-1);/减去一开始为了矩阵运算加的1 for(i=1;i=k;i+) for(j=1;j=k;j+) fprintf(fp,C%d%d=%13.11en,i-1,j-1,Cij);/这里为了匹配矩阵下标和书上所示的

17、系数下标所以减1void compare()/比较f(x*,y*)与p(x*,y*) double pn8+15+1=0;/下标都加1，方便以后运算 int i,j; double temp=0; int ii,jj; for(i=1;i=8;i+) xi=0.1*i; for(j=1;j=5;j+) yj=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+=Cii+1jj+1*pow(xi,ii)*pow(yj,jj); pnij=temp; newt

18、on(x,y); eryuanchazhi(te,ue); FILE *fp; fp=fopen(shubiao-compare.txt,a); for(i=1;i=8;i+) for(j=1;j=5;j+) fprintf(fp,x*%d=%ft y*%d=%ft n,i,xi,j,yj); fprintf(fp,z*%d%d=%13.11et p*%d%d=%13.11etn,i,j,zeij,i,j,pnij); fprintf(fp,绝对误差e=f*(x%d,y%d)-p*(x%d,y%d)=%13.11en,i,j,i,j,zeij-pnij); fprintf(fp,相对误差er=

19、e/f*(x%d,y%d)=%13.11en,i,j,(zeij-pnij)/zeij); void main() int i,j; for(i=0;i11;i+) for(j=0;j21;j+) xi=0.08*i; yj=0.5+0.05*j; newton(x,y); eryuanchazhi(te,ue); qumianchazhi(); compare();3运行结果3.1 数表,x0=0.000000 y0=0.500000 f(x,y)=4.46504018480e-001x0=0.000000 y1=0.550000 f(x,y)=3.24683262927e-001x0=0.

20、000000 y2=0.600000 f(x,y)=2.10159686683e-001x0=0.000000 y3=0.650000 f(x,y)=1.03043608316e-001x0=0.000000 y4=0.700000 f(x,y)=3.40189556266e-003x0=0.000000 y5=0.750000 f(x,y)=-8.87358136380e-002x0=0.000000 y6=0.800000 f(x,y)=-1.73371632750e-001x0=0.000000 y7=0.850000 f(x,y)=-2.50534611467e-001x0=0.000000 y8=0.900000 f(x,y)=-3.20276506388e-001x0=0.000000 y9=0.950000 f(x,y)=-3.82668066110e-001x0=0.000000 y10=1.000000 f(x,y)=-4.37795