北航数值分析A大作业

1、一、算法设计方案1、解非线性方程组将各拟合节点（xi，yj）分别带入非线性方程组，求出与相对应的数组teij，ueij，求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。2、二元二次分偏插值对数表z(t,u)进行分片二次代数插值，求得对应（tij，uij）处的值，即为 的值。根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断(t ij, u ij ) 所在，的区域，再作此区域的插值，计算 z ij，相应的Lagrange形式的插值多项式为：其中 (k=m-1, m, m+1) (r=n-1, n, n+1)3、曲面拟合从

2、k=1开始逐渐增大k的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当时结束计算。拟合基函数r(x)s(y)选择为r(x)=xr，s(y)=ys。拟合系数矩阵c通过连续两次解线性方程组求得。， 其中，4、观察比较计算的值并输出结果，以观察逼近的效果。其中。二、全部源程序/ hean.cpp : 定义控制台应用程序的入口点。/#include stdafx.h#include #include #include void Set_non_JacobiA(double* A,double* x);/求题中非线性方程组对应自变量向量x的雅克比/矩阵void Set_non_B(double

3、* B,double* x,double a,double b);/求非线性方程组Newton迭代法的右/端式:-F(x)void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C);/矩阵相乘ABCvoid Transpose(double *A,int m,int n,double* AT);/转置void Doolittle(double *A,int n,int *M);void Solve_LU(double* A,int n,double* b,double* x);/解方程LUxbvoid Solve

4、_lin(double* A,int n,double* B,double* x,int m);/解线性方程组AxBdouble Vector_FanShu(double *V,int n);void Solve_non_Equation(double a,double b,double* x);/求解非线性方程组int Near_Index(double* v,double a,int n);/查找n维向量V中与常数a最接近的元素的下标double ChaZhi(double a,double b);/求数表z(t,u)在点（x，y）处的分片二次代数差值void NiHe(double*U,

5、double*x,double*y,int m,int n,int p,int q,double*C);/对mn数表U(x,y)进行二元多项式拟合double Pxy(double *C,int p,int q,double x,double y);/求多项式拟合函数P(x,y)在点（x，y）处的函数值void main()double non4; /非线性方程组变量向量（t,u,v,w）double f1121;/f(x,y)在拟合节点处的值double x11,y21;/拟合节点double *C;/拟合系数矩阵double det=1e-7;/拟合精度要求double p,d;int i

6、,j,k,t;/设置拟合节点for(i=0;i11;i+)xi=0.08*i;for(j=0;j21;j+)yj=0.5+0.05*j;/求拟合节点处的f(x,y)的值for(i=0;i11;i+)for(j=0;j21;j+)Solve_non_Equation(xi,yj,non);fij=ChaZhi(non0,non1);printf(n 数值分析第三次大作业n);/打印数表f(xi,yi)printf(n f(xi,yi) n);for(t=0;t21/3;t+)printf(-n);printf(| f(x,y) |);for(i=3*t;i3*t+3;i+)printf( y=%

7、4.2f |,yi);printf(n);for(j=0;j11;j+)printf(| x=%4.2f |,xj);for(i=3*t;i3*t+3;i+)printf( %+.12e |,fji);printf(n);printf(-n);printf(n);k=0;printf( 不同k对应的精度 );printf(n-);doC=(double*) calloc(k+1)*(k+1),sizeof(double);NiHe(*f,x,y,11,21,k+1,k+1,C);/求在当前k值拟合的节点处误差d=0;for(i=0;i11;i+)for(j=0;j=det)free(C);el

8、seprintf(n-);printf(n 故可知k=k=%d%.0e，满足题设要求,det);break;while(+k11);/打印拟合系数矩阵cprintf(nnn 当k=%d时拟合函数p(x,y)的系数矩阵c：n,k);printf(-n);for(t=0;t(k+1)/3;t+)for(i=0;i(k+1);i+)for(j=3*t;j3*t+3;j+)printf( %+.12e ,Ci*(k+1)+j);printf(n);printf(-n);/打印(x*,y*)处的拟合逼近情况printf(nn观察以下各点的逼近情况n);printf( | (x*,y*) | f(x*,y

9、*) | p(x*,y*) | |f-p| |n);printf(-n);for(i=1;i=8;i+)for(j=1;j=5;j+)Solve_non_Equation(0.1*i,0.5+0.2*j,non);d=ChaZhi(non0,non1);p=Pxy(C,k+1,k+1,0.1*i,0.5+0.2*j);printf( | (%3.1f,%3.1f)| %+.12e | %+.12e | %f |n,0.1*i,0.5+0.2*j,d,p,abs(d-p);printf(-nnn);/求题中非线性方程组对应自变量向量x的雅克比矩阵void Set_non_JacobiA(doub

10、le* A,double* x)/A为*4阶系数矩阵，x为自变量(t,u,v,w)A0=-0.5*sin(x0);A1=1;A2=1;A3=1;A4=1;A5=0.5*cos(x1);A6=1;A7=1;A8=0.5;A9=1;A10=-sin(x2);A11=1;A12=1;A13=0.5;A14=1;A15=cos(x3);/求非线性方程组Newton迭代法的右端式:-F(x)void Set_non_B(double* B,double* x,double a,double b)/a非线性方程的参数，对应题中的x/b非线性方程的参数，对应题中的yB0=-(0.5*cos(x0)+x1+x

11、2+x3-a-2.67);B1=-(x0+0.5*sin(x1)+x2+x3-b-1.07);B2=-(0.5*x0+x1+cos(x2)+x3-a-3.74);B3=-(x0+0.5*x1+x2+sin(x3)-b-0.79);/矩阵相乘ABC，其中A为ms阶，B为sn阶。void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C)int i,j,k;for(i=0;im;i+)for(j=0;jn;j+)Ci*n+j=0;for(k=0;ks;k+)Ci*n+j+=Ai*s+k*Bk*n+j;/将mn阶矩阵A的

12、转置为ATvoid Transpose(double *A,int m,int n,double* AT)int i,j;for(i=0;im;i+)for(j=0;jn;j+)ATj*m+i=Ai*n+j;/对n阶矩阵A进行选主元的Doolittle分解void Doolittle(double *A,int n,int *M)/将分解得到的L、U仍然存储在原来的矩阵A中int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double);for(k=0;kn;k+) for(i=k;in;i+)si=Ai*n+

13、k;for(t=0;tk;t+)si-= Ai*n+t * At*n+k;Maxs=abs(sk); Mk=k;for(i=k+1;in;i+)if(Maxsabs(si) Maxs=abs(si);Mk=i;if(Mk!=k)for(t=0;tn;t+)temp=Ak*n+t;Ak*n+t=AMk*n+t;AMk*n+t=temp;temp=sk;sk=sMk;sMk=temp;if(Maxs(1e-14) sk=1e-14;printf(%.16e方阵奇异n,Maxs);Ak*n+k=sk;for(j=k+1;(jn)&(kn-1);j+)for(t=0;tk;t+)Ak*n+j-=Ak*

14、n+t*At*n+j;Aj*n+k=sj/Ak*n+k;/解方程LUxbvoid Solve_LU(double* A,int n,double* b,double* x)int i,t;for(i=0;in;i+)xi=bi;for(t=0;t-1;i-)for(t=i+1;tn;t+)xi-=Ai*n+t*xt;xi/=Ai*n+i;/解线性方程组AxBvoid Solve_lin(double* A,int n,double* B,double* x,int m)/B为nm矩阵int* M,i,j;M=(int*) calloc(n,sizeof(int);double *BT,*xT,

15、temp;BT=(double*) calloc(n*m,sizeof(double);xT=(double*) calloc(n*m,sizeof(double);Transpose(B,n,m,BT); Doolittle(A,n,M);/将A三角分解for(i=0;im;i+)for(j=0;jn-1;j+)temp=BTi*n+j;BTi*n+j=BTi*n+Mj;BTi*n+Mj=temp;/将B转置，使得同一方程组对应的系数连续存储Solve_LU(A,n,BT+i*n,xT+i*n);Transpose(xT,m,n,x);/求n维向量V的无穷范数double Vector_Fa

16、nShu(double *V,int n)int i;double max=0;for(i=0;in;i+)if(maxabs(Vi)max=abs(Vi);return max;/求解非线性方程组void Solve_non_Equation(double a,double b,double* x)double A44,F4,detx4;double det=1e-12;/求解精度要求int i,k=0;int M=5000;/最大迭代次数/设定迭代初始值x0=0.5;x1=1;x2=1;x3=1;doSet_non_B(F,x,a,b);Set_non_JacobiA(*A,x);Solv

17、e_lin(*A,4,F,detx,1);if(Vector_FanShu(detx,4)/Vector_FanShu(x,4)det)return;for(i=0;i4;i+)xi+=detxi;k+;while(kM);printf(Newton法在该初值不收敛n);/查找n维向量V中与常数a最接近的元素的下标int Near_Index(double* v,double a,int n)int i,Index;double min;min=abs(v0-a);Index=0;for(i=1;iabs(vi-a)min=abs(vi-a);Index=i;return Index;/求数表

18、z(t,u)在点（x，y）处的分片二次代数差值double ChaZhi(double a,double b)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,-0.5;double x6=0,0.2,0.4,0.6,0.8,1;double y6=0,0.

19、4,0.8,1.2,1.6,2;double L3,L_3,Z;int i,j,k,r,t;i=Near_Index(x,a,6);j=Near_Index(y,b,6);/插值区域边界处插值节点的选取if(i=0)i=1;if(i=5)i=4;if(j=0)j=1;if(j=5)j=4;for(k=i-1;k=i+1;k+)Lk-i+1=1;for(t=i-1;t=i+1;t+)if(t!=k)Lk-i+1*=(a-xt)/(xk-xt);for(r=j-1;r=j+1;r+)L_r-j+1=1;for(t=j-1;t=j+1;t+)if(t!=r)L_r-j+1*=(b-yt)/(yr-y

20、t);Z=0;for(k=i-1;k=i+1;k+)for(r=j-1;r=j+1;r+)Z+=(Lk-i+1*L_r-j+1*zkr);return Z;/对mn数表U(x,y)进行二元多项式拟合void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C)/U为拟合数据点处的函数值int i,j;double *B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*CT;B=(double*) calloc(m*p,sizeof(double);BT=(double*) callo

21、c(p*m,sizeof(double);BTB=(double*) calloc(p*p,sizeof(double);G=(double*) calloc(n*q,sizeof(double);GT=(double*) calloc(q*n,sizeof(double);GTG=(double*) calloc(q*q,sizeof(double);BTUG=(double*) calloc(p*q,sizeof(double);temp=(double*) calloc(p*n,sizeof(double);for(i=0;im;i+)for(j=0;jp;j+)Bi*p+j=pow(xi,j);for(i=0;in;i+)for(j=0;jq;j+)Gi*q+j=pow(yi,j);Tr