北航数值分析A大作业.docx
- 文档编号:3067634
- 上传时间:2022-11-17
- 格式:DOCX
- 页数:10
- 大小:45.03KB
北航数值分析A大作业.docx
《北航数值分析A大作业.docx》由会员分享,可在线阅读,更多相关《北航数值分析A大作业.docx(10页珍藏版)》请在冰豆网上搜索。
北航数值分析A大作业
一、算法设计方案
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、观察比较
计算
的值并输出结果,以观察
逼近
的效果。
其中
。
二、全部源程序
2e|",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);
2e",k,d);
if(d>=det)
free(C);
else
{
printf("\n----------------------------------------------");
printf("\n故可知k=k=%d<%.0e,满足题设要求",det);
break;
}
}while(++k<11);
2e",C[i*(k+1)+j]);
printf("\n");
}
printf("----------------------------------------------------------------------------\n");
}
2e|%+.12e|%f|\n",*i,+*j,d,p,abs(d-p));
}
printf("-----------------------------------------------------------------------\n\n\n");
}
6e方阵奇异\n",Maxs);
}
A[k*n+k]=s[k];
for(j=k+1;(j { for(t=0;t A[k*n+j]-=A[k*n+t]*A[t*n+j]; A[j*n+k]=s[j]/A[k*n+k]; } } } //解方程LUx=b voidSolve_LU(double*A,intn,double*b,double*x) { inti,t; for(i=0;i { x[i]=b[i]; for(t=0;t } for(i=n-1;i>-1;i--) { for(t=i+1;t x[i]-=A[i*n+t]*x[t]; x[i]/=A[i*n+i]; } } //解线性方程组Ax=B voidSolve_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 { for(j=0;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 if(max 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]=;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) return; for(i=0;i<4;i++) x[i]+=detx[i]; k++; }while(k 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 { 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]={,,,,,, ,,,,, ,,,,, ,,,,, ,,,,, ,,,,}; doublex[6]={0,,,,,1}; doubley[6]={0,,,,,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*)calloc(p*q,sizeof(double)); temp=(double*)calloc(p*n,sizeof(double)); for(i=0;i { for(j=0;j } for(i=0;i { for(j=0;j } Transpose(B,m,p,BT); Transpose(G,n,q,GT); Array_Mult_Array(BT,B,p,m,p,BTB); Array_Mult_Array(GT,G,q,n,q,GTG); Array_Mult_Array(BT,U,p,m,n,temp); Array_Mult_Array(temp,G,p,n,q,BTUG); free(B); free(BT); free(G); free(GT); free(temp); temp=(double*)calloc(p*q,sizeof(double)); Solve_lin(BTB,p,BTUG,temp,q)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业