北航数值分析A大作业3

一、算法设计方案1、解非线性方程组将各拟合节点（xi，yj）分别带入非线性方程组，求出与相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。2、二元二次分偏插值对数表z(t,u)进行分片二次代数插值，求得对应（tij，uij）处的值，即为的值。根据给定的数表可将整个插值区域分成16个小的区域，故先判断tij,uij所在，的区域，再作此区域的插值，计算zij，相应的Lagrange形式的插值多项式为：其中(k=m-1,m,m+1)(r=n-1,n,n+1)3、曲面拟合从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<stdio.h>#include<stdlib.h>#include<math.h>voidSet_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);//对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;i<11;i++) x[i]=0.08*i; for(j=0;j<21;j++) y[j]=0.5+0.05*j; //求拟合节点处的f(x,y)的值 for(i=0;i<11;i++) { for(j=0;j<21;j++) { Solve_non_Equation(x[i],y[j],non); f[i][j]=ChaZhi(non[0],non[1]); } } printf("\n数值分析第三次大作业\n"); //打印数表f(xi,yi) printf("\nf(xi,yi)\n"); for(t=0;t<21/3;t++) { printf("-------------------------------------------------------------------------------------\n"); printf("|f(x,y)|"); for(i=3*t;i<3*t+3;i++) { printf("y=%4.2f|",y[i]); } printf("\n"); for(j=0;j<11;j++) { printf("|x=%4.2f|",x[j]); for(i=3*t;i<3*t+3;i++) printf("%+.12e|",f[j][i]); printf("\n"); } printf("-------------------------------------------------------------------------------------\n"); printf("\n"); } 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;i<11;i++) { for(j=0;j<21;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，满足题设要求",det); break; } }while(++k<11); //打印拟合系数矩阵c printf("\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;j<3*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"); 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(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]);}//求非线性方程组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;i<m;i++) for(j=0;j<n;j++) { C[i*n+j]=0; for(k=0;k<s;k++) C[i*n+j]+=A[i*s+k]*B[k*n+j]; }}//将m×n阶矩阵A的转置为ATvoidTranspose(double*A,intm,intn,double*AT){ inti,j; for(i=0;i<m;i++) for(j=0;j<n;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;k<n;k++) { for(i=k;i<n;i++) { s[i]=A[i*n+k]; for(t=0;t<k;t++) s[i]-=A[i*n+t]*A[t*n+k]; } Maxs=abs(s[k]);M[k]=k; for(i=k+1;i<n;i++) { if(Maxs<abs(s[i])) { Maxs=abs(s[i]); M[k]=i; } } if(M[k]!=k) { for(t=0;t<n;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]]=temp; } if(Maxs<(1e-14)) { s[k]=1e-14; printf("%.16e方阵奇异\n",Maxs); } A[k*n+k]=s[k]; for(j=k+1;(j<n)&&(k<n-1);j++) { for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j]; A[j*n+k]=s[j]/A[k*n+k]; } }}//解方程LUx＝bvoidSolve_LU(double*A,intn,double*b,double*x){ inti,t; for(i=0;i<n;i++) { x[i]=b[i]; for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t]; } for(i=n-1;i>-1;i--) { for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t]; x[i]/=A[i*n+i]; }}//解线性方程组Ax＝BvoidSolve_lin(double*A,intn,double*B,double*x,intm)//B为n×m矩阵{ int*M,i,j; M=(int*)calloc(n,sizeof(int)); double*BT,*xT,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;i<m;i++) { for(j=0;j<n-1;j++) { temp=BT[i*n+j]; BT[i*n+j]=BT[i*n+M[j]]; BT[i*n+M[j]]=temp; }//将B转置，使得同一方程组对应的系数连续存储 Solve_LU(A,n,BT+i*n,xT+i*n); } Transpose(xT,m,n,x);}//求n维向量V的无穷范数doubleVector_FanShu(double*V,intn){ inti; doublemax=0; for(i=0;i<n;i++) if(max<abs(V[i])) max=abs(V[i]); returnmax;}//求解非线性方程组voidSolve_non_Equation(doublea,doubleb,double*x){ doubleA[4][4],F[4],detx[4]; doubledet=1e-12; //求解精度要求 inti,k=0; intM=5000; //最大迭代次数 //设定迭代初始值 x[0]=0.5; x[1]=1; x[2]=1; x[3]=1; do { Set_non_B(F,x,a,b); Set_non_JacobiA(*A,x); Solve_lin(*A,4,F,detx,1); if(Vector_FanShu(detx,4)/Vector_FanShu(x,4)<det) return; for(i=0;i<4;i++) x[i]+=detx[i]; k++; }while(k<M); printf("Newton法在该初值不收敛\n");}//查找n维向量V中与常数a最接近的元素的下标intNear_Index(double*v,doublea,intn){ inti,Index; doublemin; min=abs(v[0]-a); Index=0; for(i=1;i<n;i++) { if(min>abs(v[i]-a)) { min=abs(v[i]-a); Index=i; } } returnIndex;}//求数表z(t,u)在点（x，y）处的分片二次代数差值doubleChaZhi(doublea,doubleb){ doublez[6][6]={-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}; doublex[6]={ 0, 0.2, 0.4, 0.6, 0.8, 1 }; doubley[6]={ 0, 0.4, 0.8, 1.2, 1.6, 2 }; doubleL[3],L_[3],Z; inti,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++) { L[k-i+1]=1; for(t=i-1;t<=i+1;t++) { if(t!=k) L[k-i+1]*=((a-x[t])/(x[k]-x[t])); } } 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-y[t])/(y[r]-y[t])); } } Z=0; for(k=i-1;k<=i+1;k++) { for(r=j-1;r<=j+1;r++) { Z+=(L[k-i+1]*L_[r-j+1]*z[k][r]); } } returnZ;}//对m×n数表U(x,y)进行二元多项式拟合voidNiHe(double*U,double*x,double*y,intm,intn,intp,intq,double*C)//U为拟合数据点处的函数值{ inti,j; double*B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*CT; B=(double*)calloc(m*p,sizeof(double)); BT=(double*)calloc(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*)callo