数值分析第三次大作业.docx
《数值分析第三次大作业.docx》由会员分享，可在线阅读，更多相关《数值分析第三次大作业.docx（34页珍藏版）》请在冰点文库上搜索。
数值分析第三次大作业
数值分析第三次大作业
姓名：
陈怡然
学号：
38091112
2010/6/16
数值分析第三次大作业
一、设计方案：
运用newton迭代法，求解已给方程的根
对应的
。
并根据题中所给的u,t,z的二维数表对z进行插值得到
，并使用最小二乘法曲面拟合法对z=f（x,y）进行拟合，在所给的
条件下停止运算。
将
带入z=f（x,y）和p（x,y）求出所求的值，并比较。
二、源程序：
#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）//牛顿迭代法求方程的解。
{
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\n",k_now,sigma）;
fclose（p）;
}
else
printf（"%d%.12le\n\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）//dolittle分解
{
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;tvec_r[j]=vec_tmp[j]-tmp;
}
LUSolve（mat_tmp,vec_dx,vec_r,rank）;
max_x=vec_x[0],max_dx=vec_dx[0];
for（j=0;j{
if（vec_x[j]>
升级会员 