北航数值分析A大作业3

1一、算法设计方案1、解非线性方程组将各拟合节点（xi，yj）分别带入非线性方程组，求出与(,)iixy相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。2、二元二次分偏插值对数表z(t,u)进行分片二次代数插值，求得对应（tij，uij）处的值，即为),(jiyxf的值。根据给定的数表，可将整个插值区域分成16个小的区域，故先判断tij,uij所在，的区域，再作此区域的插值，计算zij，相应的Lagrange形式的插值多项式为：112211(,)()()(,)mnkrkrkmrnptultluftu其中11()mwkwmkwwkttlttt(k=m-1,m,m+1)11()nwrwnrwwryyluyy(r=n-1,n,n+1)23、曲面拟合从k=1开始逐渐增大k的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当710时结束计算。拟合基函数φr(x)ψs(y)选择为φr(x)=xr，ψs(y)=ys。拟合系数矩阵c通过连续两次解线性方程组求得。[]rscC，11()()TTTCBBBUGGG其中0011101011[()]1kkrikxxxxxxxB，0011101011[()]1kksjkyyyyGyyy[(,)]ijfxyU4、观察比较计算)5,,2,1,8,,2,1)(,(),,(****jiyxpyxfjiji的值并输出结果，以观察),(yxp逼近),(yxf的效果。其中jyixji2.05.0,1.0**。二、全部源程序//hean.cpp:定义控制台应用程序的入口点。//#includestdafx.h#includestdio.h3#includestdlib.h#includemath.hvoidSet_non_JacobiA(double*A,double*x);//求题中非线性方程组对应自变量向量x的雅克比//矩阵voidSet_non_B(double*B,double*x,doublea,doubleb);//求非线性方程组Newton迭代法的右//端式:-F(x)voidArray_Mult_Array(double*A,double*B,intm,ints,intn,double*C);//矩阵相乘AB＝CvoidTranspose(double*A,intm,intn,double*AT);//转置voidDoolittle(double*A,intn,int*M);voidSolve_LU(double*A,intn,double*b,double*x);//解方程LUx＝bvoidSolve_lin(double*A,intn,double*B,double*x,intm);//解线性方程组Ax＝BdoubleVector_FanShu(double*V,intn);voidSolve_non_Equation(doublea,doubleb,double*x);//求解非线性方程组intNear_Index(double*v,doublea,intn);//查找n维向量V中与常数a最接近的元素的下标doubleChaZhi(doublea,doubleb);//求数表z(t,u)在点（x，y）处的分片二次代数差值voidNiHe(double*U,double*x,double*y,intm,intn,intp,intq,double*C);4//对m×n数表U(x,y)进行二元多项式拟合doublePxy(double*C,intp,intq,doublex,doubley);//求多项式拟合函数P(x,y)在点（x，y）处的函数值voidmain(){doublenon[4];//非线性方程组变量向量（t,u,v,w）doublef[11][21];//f(x,y)在拟合节点处的值doublex[11],y[21];//拟合节点double*C;//拟合系数矩阵doubledet=1e-7;//拟合精度要求doublep,d;inti,j,k,t;//设置拟合节点for(i=0;i11;i++)x[i]=0.08*i;for(j=0;j21;j++)y[j]=0.5+0.05*j;//求拟合节点处的f(x,y)的值for(i=0;i11;i++){for(j=0;j21;j++){Solve_non_Equation(x[i],y[j],non);f[i][j]=ChaZhi(non[0],non[1]);}5}printf(\n数值分析第三次大作业\n);//打印数表f(xi,yi)printf(\nf(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=%4.2f|,y[i]);}printf(\n);for(j=0;j11;j++){printf(|x=%4.2f|,x[j]);for(i=3*t;i3*t+3;i++)printf(%+.12e|,f[j][i]);printf(\n);}printf(-------------------------------------------------------------------------------------\n);printf(\n);6}k=0;printf(不同k对应的精度);printf(\n----------------------------------------------);do{C=(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;j21;j++){p=Pxy(C,k+1,k+1,x[i],y[j]);d+=((f[i][j]-p)*(f[i][j]-p));}}printf(\nk=%d对应的σ=%.12e,k,d);if(d=det)free(C);else{printf(\n----------------------------------------------);printf(\n故可知k=k=%d%.0e，满足题设7要求,det);break;}}while(++k11);//打印拟合系数矩阵cprintf(\n\n\n当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,C[i*(k+1)+j]);printf(\n);}printf(----------------------------------------------------------------------------\n);}//打印(x*,y*)处的拟合逼近情况printf(\n\n观察以下各点的逼近情况\n);printf(|(x*,y*)|f(x*,y*)|p(x*,y*)||f-p||\n);printf(-----------------------------------------------------------------------\n);8for(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(non[0],non[1]);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(-----------------------------------------------------------------------\n\n\n);}//求题中非线性方程组对应自变量向量x的雅克比矩阵voidSet_non_JacobiA(double*A,double*x)//A为*4阶系数矩阵，x为自变量(t,u,v,w){A[0]=-0.5*sin(x[0]);A[1]=1;A[2]=1;A[3]=1;A[4]=1;A[5]=0.5*cos(x[1]);A[6]=1;A[7]=1;A[8]=0.5;A[9]=1;A[10]=-sin(x[2]);A[11]=1;A[12]=1;A[13]=0.5;A[14]=1;A[15]=cos(x[3]);}9//求非线性方程组Newton迭代法的右端式:-F(x)voidSet_non_B(double*B,double*x,doublea,doubleb)//a非线性方程的参数，对应题中的x//b非线性方程的参数，对应题中的y{B[0]=-(0.5*cos(x[0])+x[1]+x[2]+x[3]-a-2.67);B[1]=-(x[0]+0.5*sin(x[1])+x[2]+x[3]-b-1.07);B[2]=-(0.5*x[0]+x[1]+cos(x[2])+x[3]-a-3.74);B[3]=-(x[0]+0.5*x[1]+x[2]+sin(x[3])-b-0.79);}//矩阵相乘AB＝C，其中A为m×s阶，B为s×n阶。voidArray_Mult_Array(double*A,double*B,intm,ints,intn,double*C){inti,j,k;for(i=0;im;i++)for(j=0;jn;j++){C[i*n+j]=0;for(k=0;ks;k++)C[i*n+j]+=A[i*s+k]*B[k*n+j];}}//将m×n阶矩阵A的转置为ATvoidTranspose(double*A,intm,intn,double*AT){10inti,j;for(i=0;im;i++)for(j=0;jn;j++)AT[j*m+i]=A[i*n+j];}//对n阶矩阵A进行选主元的Doolittle分解voidDoolittle(double*A,intn,int*M)//将分解得到的L、U仍然存储在原来的矩阵A中{inti,j,k,t;double*s;doubleMaxs,temp;s=(double*)calloc(n,sizeof(double));for(k=0;kn;k++){for(i=k;in;i++){s[i]=A[i*n+k];for(t=0;tk;t++)s[i]-=A[i*n+t]*A[t*n+k];}Maxs=abs(s[k]);M[k]=k;for(i=k+1;in;i++){if(Maxsabs(s[i])){11Maxs=abs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;tn;t++){temp=A[k*n+t];A[k*n+t]=A[M[k]*n+t];A[M[k]*n+t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=