数值分析B第二次大作业.docx

数值分析B第二次大作业.docx

《数值分析B第二次大作业.docx》由会员分享，可在线阅读，更多相关《数值分析B第二次大作业.docx（31页珍藏版）》请在冰豆网上搜索。

数值分析B第二次大作业

一、算法设计方案

1．解非线性方程

使用valuez（）进行牛顿迭代，对

，

，（

）的11行21列数组

分别求出原题中方程组的一组解，得到一组与

一一对应的

。

2．分片法插值

应用Sharding（）函数（包括inv（）,kchoose（）等函数）对于已求出的

，使用分片二次代数插值法对原题中关于

的数表进行插值得到与每个

对应的

。

于是产生了z=f（x,y）的11*21个数值解。

3．最小二乘法曲面拟合

应用preprocess（），LUbreakdown（），LUsolving（）函数从k=1开始逐渐增大k的值，并使用最小二乘法曲面拟合法对z=f（x,y）进行拟合，得到每次的

。

当

时结束计算，输出拟合结果。

4．计算误差

调用函数fvalue（），shuchu（）计算

的值并输出结果，对

逼近

的效果进行观察。

其中

。

5．主函数

Main（）主要进行子函数的调用，放于程序尾部。

包括kchoose（）、fvalue（）、shuchu（）用于计算

，以及近似值和数表。

二、全部源程序如下：

#include

#include

#definess1.0e-12//精度水平ss；

#definesss1.0e-219

#defineL500//迭代最大次数；

#defineK10//确定最小二乘法拟合时的最多次数；

doubleaaa_u[6]={0,0.4,0.8,1.2,1.6,2.0},//输入矩阵数据表在z（u，t）

aaa_t[6]={0,0.2,0.4,0.6,0.8,1.0},

aaa_ztu[6][6]={

{-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},

},

aaa_val_z[11][21]={0},

init_tuvw[4]={1,2,1,2},

aaa_crs[K][K]={0};

intvaluez（doublex,doubley,inti,intj）//用牛顿迭代法计算x（i），y（i）对应的函数值z

{

doublevtuvw[4],det[4],ff[4][4],f[4],mo_k,mo_delta_k,shang,t,u;

intn=0,k,flag_max,flag_stat;

voidsolve（doubleff[4][4],doubledet[4],doublef[4]）;

doubleSharding（double,double）;

flag_stat=1;

for（k=0;k

do

{

f[0]=-1.0*（0.5*cos（vtuvw[0]）+vtuvw[1]+vtuvw[2]+vtuvw[3]-x-2.67）;//输入非线性方程组

f[1]=-1.0*（vtuvw[0]+0.5*sin（vtuvw[1]）+vtuvw[2]+vtuvw[3]-y-1.07）;

f[2]=-1.0*（0.5*vtuvw[0]+vtuvw[1]+cos（vtuvw[2]）+vtuvw[3]-x-3.74）;

f[3]=-1.0*（vtuvw[0]+0.5*vtuvw[1]+vtuvw[2]+sin（vtuvw[3]）-y-0.79）;

ff[0][0]=-0.5*sin（vtuvw[0]）;ff[0][1]=1;ff[0][2]=1;ff[0][3]=1;//对非线性方程组计算关于t,u,v,w的偏导数

ff[1][0]=1;ff[1][1]=0.5*cos（vtuvw[1]）;ff[1][2]=1;ff[1][3]=1;

ff[2][0]=0.5;ff[2][1]=1;ff[2][2]=-1*sin（vtuvw[2]）;ff[2][3]=1;

ff[3][0]=1;ff[3][1]=0.5;ff[3][2]=1;ff[3][3]=cos（vtuvw[3]）;

solve（ff,det,f）;

flag_max=0;

for（k=1;kfabs（det[flag_max]））flag_max=k;

mo_delta_k=det[flag_max];

flag_max=0;

for（k=1;kfabs（vtuvw[flag_max]））flag_max=k;

mo_k=vtuvw[flag_max];

shang=fabs（mo_delta_k/mo_k）;

if（shang

{

t=vtuvw[0]+det[0],u=vtuvw[1]+det[1];

flag_stat=0;

break;

}

else

for（k=0;k<4;k++）vtuvw[k]+=det[k];

}

while（n++<=L）;

if（!

flag_stat）

{

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,aaa_val_z[i][j]=Sharding（u,t））;

}

returnflag_stat;

}

voidsolve（doubleff[4][4],doubledet[4],doublef[4]）//牛顿消元法解线性方程组

{

inti,j,k,ik;

doublehe,mik;

for（k=0;k<4-1;k++）

{

ik=k;

for（i=k;i<4;i++）

if（fabs（ff[ik][k]）

for（j=k;j<4;j++）

{

he=ff[k][j];ff[k][j]=ff[ik][j];ff[ik][j]=he;

}

he=f[k];f[k]=f[ik];f[ik]=he;

for（i=k+1;i<4;i++）

{

mik=ff[i][k]/ff[k][k];

for（j=k;j<4;j++）

ff[i][j]=ff[i][j]-mik*ff[k][j];

f[i]-=mik*f[k];

}

}

det[4-1]=f[4-1]/ff[4-1][4-1];

for（k=4-2;k>=0;k--）

{

he=0;

for（j=k+1;j

det[k]=（f[k]-he）/ff[k][k];

}

}

doubleSharding（doubleu,doublet）//分片法进行代数二次插值

{

intr,i,j,k;

doubleduiyzh_u[3],duiyzh_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;

duiyzh_u[0]=（u-aaa_u[k]）*（u-aaa_u[k+1]）/（aaa_u[k-1]-aaa_u[k]）/（aaa_u[k-1]-aaa_u[k+1]）;

duiyzh_u[1]=（u-aaa_u[k-1]）*（u-aaa_u[k+1]）/（aaa_u[k]-aaa_u[k-1]）/（aaa_u[k]-aaa_u[k+1]）;

duiyzh_u[2]=（u-aaa_u[k-1]）*（u-aaa_u[k]）/（aaa_u[k+1]-aaa_u[k-1]）/（aaa_u[k+1]-aaa_u[k]）;

duiyzh_t[0]=（t-aaa_t[r]）*（t-aaa_t[r+1]）/（aaa_t[r-1]-aaa_t[r]）/（aaa_t[r-1]-aaa_t[r+1]）;

duiyzh_t[1]=（t-aaa_t[r-1]）*（t-aaa_t[r+1]）/（aaa_t[r]-aaa_t[r-1]）/（aaa_t[r]-aaa_t[r+1]）;

duiyzh_t[2]=（t-aaa_t[r-1]）*（t-aaa_t[r]）/（aaa_t[r+1]-aaa_t[r-1]）/（aaa_t[r+1]-aaa_t[r]）;

for（i=r-1;i<=r+1;i++）

for（j=k-1;j<=k+1;j++）z+=duiyzh_t[i-r+1]*duiyzh_u[j-k+1]*aaa_ztu[i][j];

returnz;

}

intkchoose（）//计算k的值，以及相应的误差

{

intx,y,i,j,k,kvalue;

doublevec_x[11],vec_y[21],aaa_btb[K][K]={0},aaa_gtg[K][K]={0},aaa_btbbt[K][11],aaa_btbbtu[K][21],he,wucha;

voidinv（doubleaaarix[K][K],int）;

kvalue=1;//赋初值

printf（"选择过程中的k与截断误差：

\n"）;

for（i=0;i<11;i++）vec_x[i]=0.08*i;

for（j=0;j<21;j++）vec_y[j]=0.5+0.05*j;

do

{

for（i=0;i

{

for（j=0;j

{

he=0;

for（k=0;k<11;k++）he+=pow（vec_x[k],i）*pow（vec_x[k],j）;

aaa_btb[i][j]=he;

}

}

inv（aaa_btb,kvalue）;

for（i=0;i

{

for（j=0;j

{

he=0;

for（k=0;k<21;k++）he+=pow（vec_y[k],i）*pow（vec_y[k],j）;

aaa_gtg[i][j]=he;

}

}

inv（aaa_gtg,kvalue）;

for（i=0;i

{

for（j=0;j<11;j++）

{

he=0;

for（k=0;k

aaa_btbbt[i][j]=he;

}

}

for（i=0;i

{

for（j=0;j<21;j++）

{

he=0;

for（k=0;k<11;k++）he+=aaa_btbbt[i][k]*aaa_val_z[k][j];

aaa_btbbtu[i][j]=he;

}

}

for（i=0;i

{

for（j=0;j

{

he=0;

for（k=0;k<21;k++）he+=aaa_btbbtu[i][k]*pow（vec_y[k],j）;

aaa_btb[i][j]=he;

}

}

for（i=0;i

{

for（j=0;j

{

he=0;

for（k=0;k

aaa_crs[i][j]=he;

}

}

wucha=0;

for（x=0;x<11;x++）

{

for（y=0;y<21;y++）

{

he=0;

for（i=0;i

{

for（j=0;j

{

he+=aaa_crs[i][j]*pow（0.08*x,i）*pow（0.5+0.05*y,j）;

}

}

wucha+=（he-aaa_val_z[x][y]）*（he-aaa_val_z[x][y]）;

}

}

printf（"%d%.12le\n\n",kvalue,wucha）;

if（wucha<1.0e-7）

{

printf（"Crs的系数矩阵：

\n"）;

printf（"Crs[%d][%d]=\n{\n",kvalue,kvalue）;

for（i=0;i

{

printf（"{"）;

for（j=0;j

printf（"%.12le",aaa_crs[i][kvalue-1]）;

printf（"},\n"）;

}

printf（"}\n\n"）;

returnkvalue;

}

}

while（kvalue++

return0;

}

voidinv（doubleaaarix[K][K],intrank）

{

inti,j,t,n;

doubleaaa_he[K][K]={0},aaarix_bak[K][K],vec_he[K],vec_x[K]={0},vec_dx[K]={0},vec_r[K]={0},he,max_x,max_dx;

voidpreprocess（double[K][K],double[K],int）;//对矩阵进行处理，使其条件数变小。

voidLUbreakdown（double[K][K],int）;//LU分解法

voidLUsolving（double[K][K],double[K],double[K],int）;//LU分解法求矩阵的逆

n=0;

for（i=0;i

{

for（j=0;j

vec_he[i]=0;

}

LUbreakdown（aaa_he,rank）;

for（i=0;i

{

if（i==0）vec_he[i]=1;

else

{

vec_he[i-1]=0;

vec_he[i]=1;

}

LUsolving（aaa_he,vec_x,vec_he,rank）;

while（n++<=3）

{

for（j=0;j

{

he=0;

for（t=0;t

vec_r[j]=vec_he[j]-he;

}

LUsolving（aaa_he,vec_dx,vec_r,rank）;

max_x=vec_x[0],max_dx=vec_dx[0];

for（j=0;j

{

if（vec_x[j]>max_x）max_x=vec_x[j];

if（vec_dx[j]>max_dx）max_dx=vec_dx[j];

}

for（j=0;j

}

for（j=0;j

}

}

voidpreprocess（doubleaaarix[K][K],doublevec_he[K],intrank）/对矩阵进行处理，使其条件数变小。

{

inti,k;

doublehe;

for（i=0;i

{

he=aaarix[i][0];

for（k=0;khe）he=aaarix[i][k];

for（k=0;k

vec_he[i]/=he;

}

}

voidLUbreakdown（doubleaaarix[K][K],intrank）//LU分解法

{

inti,j,k,t;

doublehe;

for（k=0;k

{

for（i=k;i

{

he=0;

for（t=0;t

he+=aaarix[i][t]*aaarix[t][k];

aaarix[i][k]-=he;

}

if（k

for（j=k+1;j

he=0;

for（t=0;t

he+=aaarix[k][t]*aaarix[t][j];

aaarix[k][j]=（aaarix[k][j]-he）/aaarix[k][k];

}

}

}

}

voidLUsolving（doubleaaarix[K][K],doublevec_x[K],doublevec_he[K],intrank）//LU分解法求矩阵的逆

{

inti,t;

doublehe;

vec_x[0]=vec_he[0]/aaarix[0][0];

for（i=1;i

{

he=0;

for（t=0;t

vec_x[i]=（vec_he[i]-he）/aaarix[i][i];

}

for（i=rank-2;i>=0;i--）

{

he=0;

for（t=i+1;t

vec_x[i]-=he;

}

}

voidfvalue（inthe）{

inti,j;

printf（"近似表达式：

\n"）;

for（i=0;i<8;i++）

for（j=0;j<5;j++）

he+=valuez（0.1*（i+1）,0.5+0.2*（j+1）,i,j）;

}

voidshuchu（）{//输出数表

inti,j,x,y,k;

doublet;

printf（"数表:

x*,y*,f（x*,y*）,p（x*,y*）:

\n"）;

for（x=0;x<8;x++）{

for（y=0;y<5;y++）{

t=0;

for（i=0;i

for（j=0;j

t=t+aaa_crs[i][j]*pow（0.1*（x+1）,i）*pow（0.5+0.2*（y+1）,j）;

}

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,aaa_val_z[x][y],x+1,y+1,t）;

}

}

}

voidmain（）//主函数

{

inti,j,k,he=0;

for（i=0;i<=10;i++）//计算xi，yi对应的值

for（j=0;j<=20;j++）

he+=valuez（0.08*i,0.5+0.05*j,i,j）;

k=kchoose（）;//计算选择过程中的k值，以及相应的误差值

fvalue（he）;//输出近似表达式的每个值

if（!

he）

shuchu（）;//输出数表

}

三、输出结果：

数表:

xi,yi,f（xi,yi）;

x0=0.000000,y0=0.500000:

f（x0,y0）=4.465040184799e-001

x0=0.000000,y1=0.550000:

f（x0,y1）=3.246832629274e-001

x0=