数值分析作业最后一题文档格式.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
升级会员 