数值分析B第三题插值拟合.docx
- 文档编号:3247189
- 上传时间:2023-05-05
- 格式:DOCX
- 页数:37
- 大小:24.78KB
数值分析B第三题插值拟合.docx
《数值分析B第三题插值拟合.docx》由会员分享,可在线阅读,更多相关《数值分析B第三题插值拟合.docx(37页珍藏版)》请在冰点文库上搜索。
数值分析B第三题插值拟合
北航2009级研究生
《数值分析B》计算实习题目
(第三题:
插值拟合)
设计文档与源程序
2009年12月12日
打印内容
1算法的设计方案
1.1运行与开发平台
1.2算法设计
2全部源代码
3输出结果:
3.1数表:
f(x,y):
3.2k值及精度σ
3.3系数Crs值
3.4数表(x,y),f(x,y),p(x,y)
1算法的设计方案
1.1运行与开发平台
操作系统:
Windows7;
开发平台:
VC++6.0;
工程类型:
Win32ConsoleApplication;
工程名:
ChazhiNihe;
1.2算法设计
设计思想如下:
1、牛顿法求解非线性方程组
从拟和节点出发,将每一对(x,y)带入所给的非线性方程中,则可以对应解出一组u和t。
解非线性方程组的具体方法是牛顿法,其中牛顿法中的每一次迭代过程中都要求雅可比矩阵,利用求雅可比矩阵逆的方法计算关于Δx的线性方程。
2、二元分片二次代数插值求拟和点的函数值
解方程组所求得的u和t,通过二元插值求得z的值,具体方法是根据已知的关于z,t,u的二维数表采用分片二次代数插值,选用拉格朗日插值基函数。
求得二元函数z=f(x,y)。
3、曲面拟和
利用已知的拟和节点、2中求得的函数值和已知的基函数进行曲面拟和。
利用拟和节点和函数值建立U阵,对于每一个K值,求出相应的B阵与G阵。
利用C=(BTB)-1BTUG(GTG)T,求出拟合系数矩阵C,则拟合函数为P(x,y)=∑∑Crsxrys,其中(0≤r≤k;0≤s≤k)。
4、将每一次拟和求出的p值和函数值z进行比较,用最小二乘法确定达到精度要求的最小的k值。
2、全部源代码
#include"stdafx.h"
#include
#include
voidMatrixReverse(double**arrin,double**arrout,intn);
voidMultiAB(double**A,double**B,introw1,intcol1,introw2,intcol2,double**R);
//非线性方程组
doubleF1(doublet,doubleu,doublev,doublew,doublex,doubley)
{
doublef;
f=0.5*cos(t)+u+v+w-x-2.67;
returnf;
}
doubleF2(doublet,doubleu,doublev,doublew,doublex,doubley)
{
doublef;
f=t+0.5*sin(u)+v+w-y-1.07;
returnf;
}
doubleF3(doublet,doubleu,doublev,doublew,doublex,doubley)
{
doublef;
f=0.5*t+u+cos(v)+w-x-3.74;
returnf;
}
doubleF4(doublet,doubleu,doublev,doublew,doublex,doubley)
{
doublef;
f=t+0.5*u+v+sin(w)-y-0.79;
returnf;
}
//牛顿法求解非线性方程组
voidNewton(doublee,intM,doublex1,doubley1,doubleres[4])
{
doubleFanshu(doublea[4]);
intk=0,i;
doublet,u,v,w,a,b;
doublex[4],y[4];
double**Array_J,**ni,**dx,**A;
Array_J=newdouble*[4];
ni=newdouble*[4];
dx=newdouble*[4];
A=newdouble*[4];
for(i=0;i<4;i++)
{
Array_J[i]=newdouble[4];
ni[i]=newdouble[4];
dx[i]=newdouble[1];
A[i]=newdouble[1];
}
x[0]=t=1;
x[1]=u=1;
x[2]=v=1;
x[3]=w=1;
do{
y[0]=x[0];
y[1]=x[1];
y[2]=x[2];
y[3]=x[3];
t=x[0];
u=x[1];
v=x[2];
w=x[3];
Array_J[0][0]=-0.5*sin(t);//Jacobi矩阵
Array_J[0][1]=1;
Array_J[0][2]=1;
Array_J[0][3]=1;
Array_J[1][0]=1;
Array_J[1][1]=0.5*cos(u);
Array_J[1][2]=1;
Array_J[1][3]=1;
Array_J[2][0]=0.5;
Array_J[2][1]=1;
Array_J[2][2]=-sin(v);
Array_J[2][3]=1;
Array_J[3][0]=1;
Array_J[3][1]=0.5;
Array_J[3][2]=1;
Array_J[3][3]=cos(w);
A[0][0]=(-1)*F1(t,u,v,w,x1,y1);
A[1][0]=(-1)*F2(t,u,v,w,x1,y1);
A[2][0]=(-1)*F3(t,u,v,w,x1,y1);
A[3][0]=(-1)*F4(t,u,v,w,x1,y1);
intn=4,m=1;
MatrixReverse((double**)Array_J,(double**)ni,n);
MultiAB((double**)ni,(double**)A,n,n,n,m,dx);
for(i=0;i<4;i++)
{
x[i]+=dx[i][0];
}
a=Fanshu(y);
doubletdx[4];
for(i=0;i<4;i++)
tdx[i]=dx[i][0];
b=Fanshu(tdx);
k++;
if(k>M)//k>M则停止计算,M为最大迭代次数
break;
}while((b/a)>=e);//e为求解精度
for(i=0;i<4;i++)
{
res[i]=x[i];//解向量
}
for(i=0;i<4;i++)
{
delete[]Array_J[i];
delete[]ni[i];
delete[]dx[i];
delete[]A[i];
}
delete[]Array_J;
delete[]ni;
delete[]dx;
delete[]A;
}
//求范数
doubleFanshu(doublea[4])
{
inti;
doubletemp_sum=0;
for(i=0;i<4;i++)
{
temp_sum+=a[i]*a[i];
}
temp_sum=sqrt(temp_sum);
returntemp_sum;
}
//计算矩阵A*B
voidMultiAB(double**A,double**B,introw1,intcol1,introw2,intcol2,double**R)
{
inti,j,k;
for(i=0;i { for(j=0;j { R[i][j]=0; for(k=0;k { R[i][j]+=A[i][k]*B[k][j]; } } } } //求逆矩阵 voidMatrixReverse(double**arrin,double**arrout,intn) { inti,j,k,t; int*M; doubletemp; double*s; double**LU; double*x,*y; LU=newdouble*[n]; for(i=0;i { LU[i]=newdouble[n]; } M=newint[n]; s=newdouble[n]; x=newdouble[n]; y=newdouble[n]; for(i=0;i { for(j=0;j LU[i][j]=arrin[i][j]; } for(k=0;k { for(i=k;i 计算中间量 { temp=0; for(t=0;t<=(k-1);t++) temp+=LU[i][t]*LU[t][k]; s[i]=LU[i][k]-temp; } ints_max=k;//第二步: 选主元,确定主元行号 doubles_temp=s[k]; for(t=k;t { if(fabs(s_temp) { s_temp=s[t+1]; s_max=t+1; } } M[k]=s_max; doublelu_temp;//第三步: 判断,如有必要,交换主元行元素 if(M[k]! =k) { for(t=0;t<=(k-1);t++) { lu_temp=LU[k][t]; LU[k][t]=LU[M[k]][t]; LU[M[k]][t]=lu_temp; } for(t=k;t { lu_temp=LU[k][t]; LU[k][t]=LU[M[k]][t]; LU[M[k]][t]=lu_temp; } s[M[k]]=s[k]; s[k]=s_temp; } LU[k][k]=s[k];//第四步: LU分解,计算l和u if(k<(n-1)) { for(j=(k+1);j { temp=0; for(t=0;t<=(k-1);t++) temp+=LU[k][t]*LU[t][j]; LU[k][j]=LU[k][j]-temp; LU[j][k]=s[j]/LU[k][k]; } } } for(i=0;i { double*b; b=newdouble[n]; b[i]=1.0; for(k=0;k { temp=b[k]; b[k]=b[M[k]]; b[M[k]]=temp; } y[0]=b[0]; for(j=1;j { temp=0; for(t=0;t<=(j-1);t++) { temp+=LU[j][t]*y[t]; } y[j]=b[j]-temp; } arrout[(n-1)][i]=y[n-1]/LU[n-1][n-1]; for(j=(n-2);j>=0;j--) { temp=0; for(t=(j+1);t temp+=LU[j][t]*arrout[t][i]; arrout[j][i]=(y[j]-temp)/LU[j][j]; } delete[]b; } for(i=0;i { delete[]LU[i]; } delete[]LU; delete[]M; delete[]s; delete[]y; delete[]x; } //存储已知的z,t,u的二维数表 voidstore(double**z,doubleu[6],doublet[6]) { inti; for(i=0;i<6;i++) { u[i]=i*(0.4); t[i]=i*(0.2); } z=newdouble*[6]; for(i=0;i<6;i++) { z[i]=newdouble[6]; } z[0][0]=-0.5; z[1][0]=-0.42; z[2][0]=-0.18; z[3][0]=0.22; z[4][0]=0.78; z[5][0]=1.5; z[0][1]=-0.34; z[1][1]=-0.5; z[2][1]=-0.5; z[3][1]=-0.34; z[4][1]=-0.02; z[5][1]=0.46; z[0][2]=0.14; z[1][2]=-0.26; z[2][2]=-0.5; z[3][2]=-0.58; z[4][2]=-0.5; z[5][2]=-0.26; z[0][3]=0.94; z[1][3]=0.3; z[2][3]=-0.18; z[3][3]=-0.5; z[4][3]=-0.66; z[5][3]=-0.66; z[0][4]=2.06; z[1][4]=1.18; z[2][4]=0.46; z[3][4]=-0.1; z[4][4]=-0.5; z[5][4]=-0.74; z[0][5]=3.5; z[1][5]=2.38; z[2][5]=1.42; z[3][5]=0.62; z[4][5]=-0.02; z[5][5]=-0.5; for(i=0;i<6;i++) { delete[]z[i]; } delete[]z; } //计算插值多项式 doubleL1(inti,intk,doubleU,doubleu[6],doublet[6]) { intt1=0; doublel1; inttempi=i; t1=tempi; l1=1; for(t1=tempi-1;t1<=(tempi+1);t1++) { if(t1! =k) { l1=l1*((U-u[t1])/(u[k]-u[t1])); } } returnl1; } doubleL2(intj,intr,doubleT,doubleu[6],doublet[6]) { intt1; doublel2; l2=1; for(t1=j-1;t1<=j+1;t1++) { if(t1! =r) { l2=l2*(T-t[t1])/(t[r]-t[t1]); } } returnl2; } //二元分片二次代数插值 doublechazhi(doubleU,doubleT,double**z,doubleu[6],doublet[6]) { doubleh=0.4; doubletao=0.2; inti,j; if(U<=(u[1]+h/2)) { j=1; } elseif(U>(u[4]-h/2)) { j=4; } else { j=(int)((U-u[0])/h);//通过差除以步长确定u的插值范围 if((u[0]+j*h+h/2)>U) { } else { j++; } } if(T<=(t[1]+tao/2)) { i=1; } elseif(T>(t[4]-tao/2)) { i=4; } else { i=(int)((T-t[0])/tao);//通过差除以步长确定t的插值范围 if((t[0]+i*tao+tao/2)>T) { } else { i++; } } intk,r; k=0,r=0; doublef=0; for(k=j-1;k<=j+1;k++) { for(r=i-1;r<=i+1;r++) { f+=z[r][k]*L1(j,k,U,u,t)*L2(i,r,T,u,t);//此处z[r][k]由R,K确定 } } returnf; } //曲面拟和,c为曲面拟和系数矩阵,U为函数值 doubleNihe(double**c,double**U,intUrow,intUcol,intn,intk) { inti,j,r,s; doubleret=0; double**B,**BT,**G,**GT,**temp,**R,**R1,**Rni,**R1ni; B=newdouble*[n]; BT=newdouble*[n]; G=newdouble*[n]; GT=newdouble*[n]; R=newdouble*[n]; R1=newdouble*[n]; Rni=newdouble*[n]; R1ni=newdouble*[n]; temp=newdouble*[n]; for(i=0;i { B[i]=newdouble[n]; BT[i]=newdouble[n]; G[i]=newdouble[n]; GT[i]=newdouble[n]; R[i]=newdouble[n]; R1[i]=newdouble[n]; Rni[i]=newdouble[n]; R1ni[i]=newdouble[n]; temp[i]=newdouble[n]; } for(i=0;i<=10;i++)//初始化B矩阵 { for(j=0;j<=k;j++) B[i][j]=pow((0.08*i),j); } for(i=0;i<=k;i++)//初始化BT矩阵 { for(j=0;j<=10;j++) BT[i][j]=B[j][i]; } for(i=0;i<=20;i++)//初始化G矩阵 { for(j=0;j<=k;j++) G[i][j]=pow((0.5+0.05*j),j); } for(i=0;i<=k;i++)//初始化BT矩阵 { for(j=0;j<=20;j++) GT[i][j]=G[j][i]; } intBTrow=k+1,BTcol=11,Brow=11,Bcol=k+1; MultiAB((double**)BT,(double**)B,BTrow,BTcol,Brow,Bcol,(double**)R); intGTrow=k+1,GTcol=21,Grow=21,Gcol=k+1; MultiAB((double**)GT,(double**)G,GTrow,GTcol,Grow,k+1,(double**)R1); MatrixReverse(R,Rni,k+1); MatrixReverse(R1,R1ni,k+1); MultiAB(G,R1ni,Grow,Gcol,k+1,k+1,R1); MultiAB(BT,U,BTrow,BTcol,Urow,Ucol,R); MultiAB(Rni,R,k+1,k+1,k+1,Ucol,temp); MultiAB(temp,R1,k+1,Ucol,Grow,k+1,c);//计算曲面拟合系数矩阵 doublep; doublep1=0; for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { p=0; for(r=0;r<=k;r++) { for(s=0;s<=k;s++) p+=c[r][s]*pow((0.08*i),r)*pow((0.5+0.05*j),s); } p1=U[i][j]-p; ret+=pow(p1,2); } } for(i=0;i { delete[]B[i]; delete[]BT[i]; delete[]G[i]; delete[]GT[i]; delete[]temp[i]; delete[]R[i]; delete[]R1[i]; delete[]Rni[i]; delete[]R1ni[i]; } delete[]B; delete[]BT; delete[]G; delete[]GT; delete[]temp; delete[]R; delete[]R1; delete[]Rni; delete[]R1ni; returnret; } voidmain() { inti,j,k; doublex,y,ret; double**U; U=newdouble*[11]; for(i=0;i<11;i++) { U[i]=newdouble[21]; } doublee=1e-7; doubleres[4]; doubleu[6],t[6]; double**z; z=newdouble*[6]; for(i=0;i<6;i++) { z[i]=newdouble[6]; } intM=100,m=11,n=21; double**c; c=newdouble*[n]; for(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 第三 题插值 拟合