数值分析作业最后一题.docx
《数值分析作业最后一题.docx》由会员分享,可在线阅读,更多相关《数值分析作业最后一题.docx(33页珍藏版)》请在冰点文库上搜索。
数值分析作业最后一题
数值分析作业3
一、算法设计方案:
1.使用牛顿迭代法(简单迭代法不收敛),对原题中给出的
,
,(
)的11*21组
分别求出原题中方程组的一组解,于是得到一组和
对应的
。
2.对于已求出的
,使用分片二次代数插值法对原题中关于
的数表进行插值得到
。
于是产生了z=f(x,y)的11*21个数值解。
3.从k=1开始逐渐增大k的值,并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合,得到每次的
。
当
时结束计算,输出拟合结果。
4.计算
的值并输出结果,以观察
逼近
的效果。
其中
。
二、源代码:
/******************************************
*计算实习大作业第三题源程序
*作者:
理学院系统与控制科学系35093120吴豪
*本程序在Dev-C++4.9.9.0版本下调试通过
*最后修改:
2008.06.11
******************************************/
#defineN6//矩阵的阶;
#defineM4//方程组未知数个数;
#defineEPSL1.0e-12//迭代的精度水平;
#defineMAXDIGIT1.0e-219
#defineSIGSUM1.0e-7
#defineL500//迭代最大次数;
#defineK10//最小二乘法拟合时的次数上限;
#defineX_NUM11
#defineY_NUM21
#defineOUTPUTMODE1//输出格式:
0--输出至屏幕,1--输出至文件;
#include
#include
doublemat_u[N]={0,0.4,0.8,1.2,1.6,2.0},
mat_t[N]={0,0.2,0.4,0.6,0.8,1.0},
mat_ztu[N][N]={
{-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},
},
mat_val_z[X_NUM][Y_NUM]={0},
init_tuvw[4]={1,2,1,2},
mat_crs[K][K]={0};
FILE*p;
intmain()//主程序入口
{
inti,j,x,y,kmin,tmp=0;
doubletmpval;
intgetval_z(double,double,int,int);
intleasquare();
voidresult_out(int);
if(OUTPUTMODE)
{
p=fopen("c:
\\Result.txt","w+");
fprintf(p,"计算结果:
\n");
fclose(p);
}
for(i=0;ifor(j=0;jif(!
tmp)printf("迭代求解z=f(x,y)完毕。
\n");
elseprintf("迭代超过最大次数,计算结果可能不准确!
\n");
if(kmin=leasquare())
{
printf("近似表达式已求出!
\n");
for(i=0;i<8;i++)
for(j=0;j<5;j++)tmp+=getval_z(0.1*(i+1),0.5+0.2*(j+1),i,j);
if(!
tmp)
{
printf("迭代求解z=f(x*,y*)完毕。
\n");
for(x=0;x<8;x++)
{
for(y=0;y<5;y++)
{
tmpval=0;
for(i=0;i{
for(j=0;j{
tmpval=tmpval+mat_crs[i][j]*pow(0.1*(x+1),i)*pow(0.5+0.2*(y+1),j);
}
}
if(OUTPUTMODE)
{
p=fopen("c:
\\Result.txt","at+");
fprintf(p,"x*[%d]=%f,y*[%d]=%f:
f(x*[%d],y*[%d])=%.12le,p(x*[%d],y*[%d])=%.12le\n",x+1,0.1*(x+1),y+1,0.5+0.2*(y+1),x+1,y+1,mat_val_z[x][y],x+1,y+1,tmpval);
fclose(p);
}
else
printf("x*[%d]=%f,y*[%d]=%f:
f(x*[%d],y*[%d])=%.12le,p(x*[%d],y*[%d])=%.12le\n",x+1,0.1*(x+1),y+1,0.5+0.2*(y+1),x+1,y+1,mat_val_z[x][y],x+1,y+1,tmpval);
}
}
}
printf("结果见C:
\\Result.txt!
");
getchar();
}
else
{
printf("近似表达式未求出,请增大K值后重试!
\n");
getchar();
}
}
intgetval_z(doublex,doubley,inti,intj)//使用牛顿迭代法迭代求解方程组f(x,y,t,u,v,w)。
{
doublevec_tuvw[M],vec_delta[M],fdao[M][M],f[M],mo_k,mo_delta_k,shang,t,u;
intn=0,k,flag_max,flag_stat;
voidsolve(doublefdao[M][M],doublevec_delta[M],doublef[M]);
doubleinterpolation(double,double);
flag_stat=1;
for(k=0;kdo
{
f[0]=-1.0*(0.5*cos(vec_tuvw[0])+vec_tuvw[1]+vec_tuvw[2]+vec_tuvw[3]-x-2.67);
f[1]=-1.0*(vec_tuvw[0]+0.5*sin(vec_tuvw[1])+vec_tuvw[2]+vec_tuvw[3]-y-1.07);
f[2]=-1.0*(0.5*vec_tuvw[0]+vec_tuvw[1]+cos(vec_tuvw[2])+vec_tuvw[3]-x-3.74);
f[3]=-1.0*(vec_tuvw[0]+0.5*vec_tuvw[1]+vec_tuvw[2]+sin(vec_tuvw[3])-y-0.79);
fdao[0][0]=-0.5*sin(vec_tuvw[0]);
fdao[0][1]=1;
fdao[0][2]=1;
fdao[0][3]=1;
fdao[1][0]=1;
fdao[1][1]=0.5*cos(vec_tuvw[1]);
fdao[1][2]=1;
fdao[1][3]=1;
fdao[2][0]=0.5;
fdao[2][1]=1;
fdao[2][2]=-1*sin(vec_tuvw[2]);
fdao[2][3]=1;
fdao[3][0]=1;
fdao[3][1]=0.5;
fdao[3][2]=1;
fdao[3][3]=cos(vec_tuvw[3]);
solve(fdao,vec_delta,f);
flag_max=0;
for(k=1;kfabs(vec_delta[flag_max]))flag_max=k;
mo_delta_k=vec_delta[flag_max];
flag_max=0;
for(k=1;kfabs(vec_tuvw[flag_max]))flag_max=k;
mo_k=vec_tuvw[flag_max];
shang=fabs(mo_delta_k/mo_k);
if(shang{
t=vec_tuvw[0]+vec_delta[0],u=vec_tuvw[1]+vec_delta[1];
flag_stat=0;
break;
}
else
for(k=0;k}
while(n++<=L);
if(!
flag_stat)
{
if(OUTPUTMODE)
{
p=fopen("c:
\\Result.txt","at+");
fprintf(p,"x%d=%f,y%d=%f:
f(x%d,y%d)=%.12le\n",i,0.08*i,j,0.5+0.05*j,i,j,mat_val_z[i][j]=interpolation(u,t));
fclose(p);
}
else
printf("x%d=%f,y%d=%f:
f(x%d,y%d)=%.12le\n",i,0.08*i,j,0.5+0.05*j,i,j,mat_val_z[i][j]=interpolation(u,t));
}
returnflag_stat;
}
voidsolve(doublefdao[M][M],doublevec_delta[M],doublef[M])//牛顿消元法解线性方程组;
{
inti,j,k,ik;
doubletmp,mik;
for(k=0;k{
ik=k;
for(i=k;ifor(j=k;j{
tmp=fdao[k][j];fdao[k][j]=fdao[ik][j];fdao[ik][j]=tmp;
}
tmp=f[k];f[k]=f[ik];f[ik]=tmp;
for(i=k+1;i{
mik=fdao[i][k]/fdao[k][k];
for(j=k;jf[i]-=mik*f[k];
}
}
vec_delta[M-1]=f[M-1]/fdao[M-1][M-1];
for(k=M-2;k>=0;k--)
{
tmp=0;
for(j=k+1;jvec_delta[k]=(f[k]-tmp)/fdao[k][k];
}
}
doubleinterpolation(doubleu,doublet)//分片代数二次插值;
{
intr,i,j,k;
doublecoe_u[3],coe_t[3],z;
z=0;
r=int(fabs((t/0.2)+0.5));
k=int(fabs((u/0.4)+0.5));
if(r=0)r=1;
if(r=5)r=4;
if(k=0)k=1;
if(k=5)k=4;
coe_u[0]=(u-mat_u[k])*(u-mat_u[k+1])/(mat_u[k-1]-mat_u[k])/(mat_u[k-1]-mat_u[k+1]);
coe_u[1]=(u-mat_u[k-1])*(u-mat_u[k+1])/(mat_u[k]-mat_u[k-1])/(mat_u[k]-mat_u[k+1]);
coe_u[2]=(u-mat_u[k-1])*(u-mat_u[k])/(mat_u[k+1]-mat_u[k-1])/(mat_u[k+1]-mat_u[k]);
coe_t[0]=(t-mat_t[r])*(t-mat_t[r+1])/(mat_t[r-1]-mat_t[r])/(mat_t[r-1]-mat_t[r+1]);
coe_t[1]=(t-mat_t[r-1])*(t-mat_t[r+1])/(mat_t[r]-mat_t[r-1])/(mat_t[r]-mat_t[r+1]);
coe_t[2]=(t-mat_t[r-1])*(t-mat_t[r])/(mat_t[r+1]-mat_t[r-1])/(mat_t[r+1]-mat_t[r]);
for(i=r-1;i<=r+1;i++)
for(j=k-1;j<=k+1;j++)z+=coe_t[i-r+1]*coe_u[j-k+1]*mat_ztu[i][j];
returnz;
}
intleasquare()//最小二乘法曲面拟合;
{
intx,y,i,j,k,k_max,k_now;
doublevec_x[X_NUM],vec_y[Y_NUM],mat_btb[K][K]={0},mat_gtg[K][K]={0},mat_btbbt[K][X_NUM],mat_btbbtu[K][Y_NUM],tmp,sigma;
voidinv(doublematrix[K][K],int);
k_now=1;
for(i=0;ifor(j=0;jdo
{
for(i=0;i{
for(j=0;j{
tmp=0;
for(k=0;kmat_btb[i][j]=tmp;
}
}
inv(mat_btb,k_now);
for(i=0;i{
for(j=0;j{
tmp=0;
for(k=0;kmat_gtg[i][j]=tmp;
}
}
inv(mat_gtg,k_now);
for(i=0;i{
for(j=0;j{
tmp=0;
for(k=0;kmat_btbbt[i][j]=tmp;
}
}
for(i=0;i{
for(j=0;j{
tmp=0;
for(k=0;kmat_btbbtu[i][j]=tmp;
}
}
for(i=0;i{
for(j=0;j{
tmp=0;
for(k=0;kmat_btb[i][j]=tmp;
}
}
for(i=0;i{
for(j=0;j{
tmp=0;
for(k=0;kmat_crs[i][j]=tmp;
}
}
sigma=0;
for(x=0;x{
for(y=0;y{
tmp=0;
for(i=0;i{
for(j=0;j{
tmp+=mat_crs[i][j]*pow(0.08*x,i)*pow(0.5+0.05*y,j);
}
}
sigma+=(tmp-mat_val_z[x][y])*(tmp-mat_val_z[x][y]);
}
}
if(OUTPUTMODE)
{
p=fopen("c:
\\Result.txt","at+");
fprintf(p,"%d%.12le\n",k_now,sigma);
fclose(p);
}
else
printf("%d%.12le\n",k_now,sigma);
if(sigma{
if(OUTPUTMODE)
{
p=fopen("c:
\\Result.txt","at+");
fprintf(p,"Crs[%d][%d]=\n{\n",k_now,k_now);
for(i=0;i{
fprintf(p,"{");
for(j=0;jfprintf(p,"%.12le",mat_crs[i][k_now-1]);
fprintf(p,"},\n");
}
fprintf(p,"}\n");
fclose(p);
}
else
{
printf("Crs[%d][%d]=\n{\n",k_now,k_now);
for(i=0;i{
printf("{");
for(j=0;jprintf("%.12le",mat_crs[i][k_now-1]);
printf("},\n");
}
printf("}\n\n");
}
returnk_now;
}
}
while(k_now++return0;
}
voidinv(doublematrix[K][K],intrank)//矩阵求逆:
使用LU分解法。
{
inti,j,k,t,n;
doublemat_tmp[K][K]={0},matrix_bak[K][K],vec_tmp[K],vec_x[K]={0},vec_dx[K]={0},vec_r[K]={0},tmp,max_x,max_dx;
voidpreprocess(double[K][K],double[K],int);
voidLUDecomposition(double[K][K],int);
voidLUSolve(double[K][K],double[K],double[K],int);
n=0;
for(i=0;i{
for(j=0;jvec_tmp[i]=0;
}
LUDecomposition(mat_tmp,rank);
for(i=0;i{
if(i==0)vec_tmp[i]=1;
else
{
vec_tmp[i-1]=0;
vec_tmp[i]=1;
}
LUSolve(mat_tmp,vec_x,vec_tmp,rank);
while(n++<=3)
{
for(j=0;j{
tmp=0;
for(t=0;t