数值分析大作业三

1、数值分析B计算实习第三题算法的设计方案1.使用牛顿迭代法(简单迭代法不收敛，对原题中给出的，(的11*21组分别求出原题中方程组的一组解，于是得到一组和对应的。2.对于已求出的，使用分片二次代数插值法对原题中关于的数表进行插值得到。于是产生了z=f(x,y的11*21个数值解。3.从k=1开始逐渐增大k的值，并使用最小二乘法曲面拟合法对z=f(x,y进行拟合，得到每次的。当时结束计算，输出拟合结果。4.计算的值并输出结果，以观察逼近的效果。其中。源程序代码#include #include #define N 4double e=1e-12;double gauss(double aNN+1,

2、double xN /gauss求线性方程组int i,j,k;double p=0;for(k=0;k if(akk=0printf("a%d%d=0n",k,k;return 1;for(i=k+1;i p=aik/akk;for(j=k;j aij=aij-p*akj;if(aN-1N-1=0printf("a%d%d=0n",N-1,N-1;return 1;xN-1=aN-1N/aN-1N-1;for(i=N-2;i>=0;i-p=0;for(j=N-1;j>i;j-p=p+aij*xj;xi=(aiN-p/aii;return 0

3、;void A_ni(double A1010,int n /求逆矩阵 double a1010,b1020,c1010,t;int i,j,m;for(i=0;i for(j=0;j bij=Aij;for(i=0;i for(j=n;j<2*n;j+bij=0;for(i=0;i bin+i=1;for(m=0;m t=bmm; i=m; while(bmm=0 bmm=bi+1m;i+;if(i>mbim=t; for(j=0;j t=bmj;bmj=bij;bij=t;for(j=m+1;j<2*n;j+ t=bmj;bmj=bij;bij=t;for(i=m+1;i

4、 for(j=2*n-1;j>=m;j-bij-=bim*bmj/bmm; for(j=2*n-1;j>=m;j-bmj/=bmm; m=n-1;while(m>0for(i=0;i for(j=2*n-1;j>=m;j- bij-=bim*bmj;m-; for(i=0;i for(j=0;j cij=bin+j;for(i=0;i for(j=0;j Aij=cij; void Newton(double x,double y,double tu4 /Newton迭代法解非线性方程组int i,j;double A45,result4=1,2,1,2,result1

5、4;double num1,num2;while(1for(i=0;i<4;i+result1i=resulti;A00=-0.5*sin(result0;A01=1;A02=1;A03=1;A10=1;A11=0.5*cos(result1;A12=1;A13=1;A20=0.5;A21=1;A22=-1*sin(result2;A23=1;A30=1;A31=0.5;A32=1;A33=cos(result3;A04=-1*(0.5*cos(result0+result1+result2+result3-x-2.67;A14=-1*(result0+0.5*sin(result1+r

6、esult2+result3-y-1.07;A24=-1*(0.5*result0+result1+cos(result2+result3-x-3.74;A34=-1*(result0+0.5*result1+result2+sin(result3-y-0.79;gauss(A,result;num1=num2=0;for(i=0;i<4;i+num1+=resulti*resulti;num2+=result1i*result1i;num1=sqrt(num1;num2=sqrt(num2;if(num1/num2<=etu0=result10;tu1=result11;tu2=

7、result12;tu3=result13;break;for(i=0;i<4;i+resulti+=result1i;double Z(double tt,double uu /二元二次插值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,

8、0.46, -0.26, -0.66, -0.74, -0.5;double l13,l23,t6=0,0.2,0.4,0.6,0.8,1.0,u6=0,0.4,0.8,1.2,1.6,2.0;int i,j;double result=0;if(tt<=0.3i=1;else if(0.3 i=2;else if(0.5 i=3;elsei=4;if(uu<=0.6j=1;else if(0.6 j=2;else if(1.0 j=3;elsej=4;l10=(tt-ti*(tt-ti+1/(ti-1-ti*(ti-1-ti+1;l11=(tt-ti-1*(tt-ti+1/(ti

9、-ti-1*(ti-ti+1;l12=(tt-ti-1*(tt-ti/(ti+1-ti-1*(ti+1-ti;l20=(uu-uj*(uu-uj+1/(uj-1-uj*(uj-1-uj+1;l21=(uu-uj-1*(uu-uj+1/(uj-uj-1*(uj-uj+1;l22=(uu-uj-1*(uu-uj/(uj+1-uj-1*(uj+1-uj;result=l10*l20*zi-1j-1+l10*l21*zi-1j+l10*l22*zi-1j+1+l11*l20*zij-1+l11*l21*zij+l11*l22*zij+1+l12*l20*zi+1j-1+l12*l21*zi+1j+l1

10、2*l22*zi+1j+1;return result;double y_x(double y,int x /幂函数int i;double result=1;for(i=1;i<=x;i+result=result*y;return result;double Pk(double x,double y,double C1010,int k /曲面拟合函数int i,j,m;double num=0,result=0;double B1110,G2110,BB1010,GG1010,U1121,UU2121,tu4;double X11=0,0.08,0.16,0.24,0.32,0.4

11、0,0.48,0.56,0.64,0.72,0.80;double Y21=0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.05,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5; for(i=0;i<11;i+Bi0=1;for(j=1;j<=k;j+Bij=y_x(Xi,j; for(i=0;i<21;i+Gi0=1;for(j=1;j<=k;j+Gij=y_x(Yi,j;for(i=0;i<11;i+for(j=0;j<21;j+Newton(Xi,Yj,tu;Uij=

12、Z(tu0,tu1;for(i=0;i<=k;i+ for(j=0;j<=k;j+ for(m=0;m<11;m+num=num+Bmi*Bmj;BBij=num;num=0;for(i=0;i<=k;i+ for(j=0;j<=k;j+num=0;for(m=0;m<21;m+num=num+Gmi*Gmj;GGij=num;A_ni(BB,k+1;A_ni(GG,k+1;num=0;for(i=0;i<=k;i+for(j=0;j<21;j+for(m=0;m<11;m+num=num+Bmi*Umj;UUij=num;num=0;fo

13、r(i=0;i<=k;i+for(j=0;j<=k;j+for(m=0;m<21;m+num=num+UUim*Gmj;Cij=num;num=0;for(i=0;i<=k;i+for(j=0;j<=k;j+for(m=0;m<=k;m+num=num+BBim*Cmj;UUij=num;num=0;for(i=0;i<=k;i+for(j=0;j<=k;j+for(m=0;m<=k;m+num=num+UUim*GGmj;Cij=num;num=0;for(i=0;i<=k;i+for(j=0;j<=k;j+result=re

14、sult+Cij*y_x(x,i*y_x(y,j;return result; void main(double tu4,U1121,B1110,BB1010,C1010;double X11=0,0.08,0.16,0.24,0.32,0.40,0.48,0.56,0.64,0.72,0.80;double Y21=0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.05,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5;double x8=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8;doub

15、le y5=0.7,0.9,1.1,1.3,1.5;double num=0;int i,j,m,k=1;for(i=0;i<11;i+for(j=0;j<21;j+Newton(Xi,Yj,tu;Uij=Z(tu0,tu1; printf("x%d=%f,y%d=%f,f(x%d,y%d=%12.11en",i,Xi,j,Yj,i,j,Uij;while(1for(i=0;i<11;i+for(j=0;j<21;j+Newton(Xi,Yj,tu;Uij=Z(tu0,tu1;for(i=0;i<11;i+for(j=0;j<21;j+num=num+(Uij-Pk(Xi,Yj,C,k*(Uij-Pk(Xi,Yj,C,k;printf("k=%d,num=%12.11en",k,num;if(num<=1e-7break;k+;num=0;Pk(Xi,Yj,C,k;printf("数表Crs：n"for(i=0;i<=k;i+for(j=0;j<=k;j+printf("%