北航数值分析A大作业3_第1页
北航数值分析A大作业3_第2页
北航数值分析A大作业3_第3页
北航数值分析A大作业3_第4页
北航数值分析A大作业3_第5页
已阅读5页,还剩20页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

一、算法设计方案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

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论