数值分析B第三题插值拟合Word格式.docx
- 文档编号:5816054
- 上传时间:2023-05-05
- 格式:DOCX
- 页数:37
- 大小:24.78KB
数值分析B第三题插值拟合Word格式.docx
《数值分析B第三题插值拟合Word格式.docx》由会员分享,可在线阅读,更多相关《数值分析B第三题插值拟合Word格式.docx(37页珍藏版)》请在冰点文库上搜索。
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)
f=t+0.5*sin(u)+v+w-y-1.07;
doubleF3(doublet,doubleu,doublev,doublew,doublex,doubley)
f=0.5*t+u+cos(v)+w-x-3.74;
doubleF4(doublet,doubleu,doublev,doublew,doublex,doubley)
f=t+0.5*u+v+sin(w)-y-0.79;
//牛顿法求解非线性方程组
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);
i<
i++)
x[i]+=dx[i][0];
a=Fanshu(y);
doubletdx[4];
tdx[i]=dx[i][0];
b=Fanshu(tdx);
k++;
if(k>
M)//k>
M则停止计算,M为最大迭代次数
break;
}while((b/a)>
=e);
//e为求解精度
res[i]=x[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;
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;
row1;
for(j=0;
j<
col2;
j++)
{
R[i][j]=0;
for(k=0;
k<
col1;
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];
n;
LU[i]=newdouble[n];
M=newint[n];
s=newdouble[n];
x=newdouble[n];
y=newdouble[n];
j<
j++)
LU[i][j]=arrin[i][j];
for(k=0;
k<
k++)//LU分解
{
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;
n-1;
if(fabs(s_temp)<
fabs(s[t+1]))
s_temp=s[t+1];
s_max=t+1;
}
M[k]=s_max;
doublelu_temp;
//第三步:
判断,如有必要,交换主元行元素
if(M[k]!
=k)
lu_temp=LU[k][t];
LU[k][t]=LU[M[k]][t];
LU[M[k]][t]=lu_temp;
for(t=k;
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);
temp=0;
for(t=0;
temp+=LU[k][t]*LU[t][j];
LU[k][j]=LU[k][j]-temp;
LU[j][k]=s[j]/LU[k][k];
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-1);
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--)
for(t=(j+1);
temp+=LU[j][t]*arrout[t][i];
arrout[j][i]=(y[j]-temp)/LU[j][j];
delete[]b;
for(i=0;
delete[]LU[i];
delete[]LU;
delete[]M;
delete[]s;
delete[]y;
delete[]x;
//存储已知的z,t,u的二维数表
voidstore(double**z,doubleu[6],doublet[6])
6;
u[i]=i*(0.4);
t[i]=i*(0.2);
z=newdouble*[6];
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;
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!
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;
=j+1;
=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;
i=(int)((T-t[0])/tao);
//通过差除以步长确定t的插值范围
if((t[0]+i*tao+tao/2)>
T)
i++;
intk,r;
k=0,r=0;
doublef=0;
for(k=j-1;
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确定
//曲面拟和,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];
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];
=10;
i++)//初始化B矩阵
=k;
B[i][j]=pow((0.08*i),j);
i++)//初始化BT矩阵
BT[i][j]=B[j][i];
=20;
i++)//初始化G矩阵
G[i][j]=pow((0.5+0.05*j),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;
p=0;
for(r=0;
for(s=0;
s<
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);
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()
doublex,y,ret;
double**U;
U=newdouble*[11];
11;
U[i]=newdouble[21];
doublee=1e-7;
doubleres[4];
doubleu[6],t[6];
double**z;
z=newdouble*[6];
z[i]=newdouble[6];
intM=100,m=11,n=21;
double**c;
c=newdouble*[n];
for(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 第三 题插值 拟合