数值分析作业最后一题文档格式.docx
- 文档编号:7824204
- 上传时间:2023-05-09
- 格式:DOCX
- 页数:33
- 大小:47KB
数值分析作业最后一题文档格式.docx
《数值分析作业最后一题文档格式.docx》由会员分享,可在线阅读,更多相关《数值分析作业最后一题文档格式.docx(33页珍藏版)》请在冰点文库上搜索。
2008.06.11
******************************************/
#defineN6//矩阵的阶;
#defineM4//方程组未知数个数;
#defineEPSL1.0e-12//迭代的精度水平;
#defineMAXDIGIT1.0e-219
#defineSIGSUM1.0e-7
#defineL500//迭代最大次数;
#defineK10//最小二乘法拟合时的次数上限;
#defineX_NUM11
#defineY_NUM21
#defineOUTPUTMODE1//输出格式:
0--输出至屏幕,1--输出至文件;
#include<
stdio.h>
math.h>
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;
i<
X_NUM;
i++)
for(j=0;
j<
Y_NUM;
j++)tmp+=getval_z(0.08*i,0.5+0.05*j,i,j);
if(!
tmp)printf("
迭代求解z=f(x,y)完毕。
elseprintf("
迭代超过最大次数,计算结果可能不准确!
if(kmin=leasquare())
printf("
近似表达式已求出!
8;
5;
j++)tmp+=getval_z(0.1*(i+1),0.5+0.2*(j+1),i,j);
tmp)
迭代求解z=f(x*,y*)完毕。
for(x=0;
x<
x++)
for(y=0;
y<
y++)
tmpval=0;
kmin;
j++)
tmpval=tmpval+mat_crs[i][j]*pow(0.1*(x+1),i)*pow(0.5+0.2*(y+1),j);
}
at+"
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);
else
结果见C:
\\Result.txt!
"
getchar();
近似表达式未求出,请增大K值后重试!
}
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;
k<
M;
k++)vec_tuvw[k]=init_tuvw[k];
do
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;
k++)if(fabs(vec_delta[k])>
fabs(vec_delta[flag_max]))flag_max=k;
mo_delta_k=vec_delta[flag_max];
k++)if(fabs(vec_tuvw[k])>
fabs(vec_tuvw[flag_max]))flag_max=k;
mo_k=vec_tuvw[flag_max];
shang=fabs(mo_delta_k/mo_k);
if(shang<
EPSL)
t=vec_tuvw[0]+vec_delta[0],u=vec_tuvw[1]+vec_delta[1];
flag_stat=0;
break;
k++)vec_tuvw[k]+=vec_delta[k];
while(n++<
=L);
flag_stat)
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;
M-1;
k++)
ik=k;
for(i=k;
i++)if(fabs(fdao[ik][k])<
fabs(fdao[i][k]))ik=i;
for(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;
mik=fdao[i][k]/fdao[k][k];
j++)fdao[i][j]=fdao[i][j]-mik*fdao[k][j];
f[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;
j++)tmp+=fdao[k][j]*vec_delta[j];
vec_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;
=r+1;
for(j=k-1;
=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;
i++)vec_x[i]=0.08*i;
j++)vec_y[j]=0.5+0.05*j;
k_now;
{
k++)tmp+=pow(vec_x[k],i)*pow(vec_x[k],j);
mat_btb[i][j]=tmp;
inv(mat_btb,k_now);
k++)tmp+=pow(vec_y[k],i)*pow(vec_y[k],j);
mat_gtg[i][j]=tmp;
inv(mat_gtg,k_now);
k++)tmp+=mat_btb[i][k]*pow(vec_x[j],k);
mat_btbbt[i][j]=tmp;
k++)tmp+=mat_btbbt[i][k]*mat_val_z[k][j];
mat_btbbtu[i][j]=tmp;
k++)tmp+=mat_btbbtu[i][k]*pow(vec_y[k],j);
k++)tmp+=mat_btb[i][k]*mat_gtg[k][j];
mat_crs[i][j]=tmp;
sigma=0;
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]);
%d%.12le\n"
k_now,sigma);
if(sigma<
SIGSUM)
fprintf(p,"
Crs[%d][%d]=\n{\n"
k_now,k_now);
{"
k_now-1;
j++)fprintf(p,"
%.12le,"
mat_crs[i][j]);
%.12le"
mat_crs[i][k_now-1]);
},\n"
}\n"
j++)printf("
}\n\n"
returnk_now;
while(k_now++<
K);
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;
rank;
j++)matrix_bak[i][j]=mat_tmp[i][j]=matrix[i][j];
vec_tmp[i]=0;
LUDecomposition(mat_tmp,rank);
if(i==0)vec_tmp[i]=1;
vec_tmp[i-1]=0;
vec_tmp[i]=1;
LUSolve(mat_tmp,vec_x,vec_tmp,rank);
while(n++<
=3)
for(t=0;
t<
t++)tmp+=matrix_bak[j][t]*vec_x[t
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 作业 最后