《北航数值分析A大作业3.doc》由会员分享，可在线阅读，更多相关《北航数值分析A大作业3.doc（17页珍藏版）》请在金锄头文库上搜索。
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+x2+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;
