数值分析第三次作业.docx

数值分析第三次作业.docx

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

数值分析第三次作业

数值分析第三次上机作业

算法设计

1、求解非线性方程组

采用牛顿法解非线性方程组,将题目中给出的

当作已知量代入题目给定的非线性方程组，求出与

相对应的数组t[i][j]，u[i][j]

2、分片二次代数插值

对所求出的数组t[i][j]，u[i][j]，通过分片二次代数插值运算，得到与数组t[11][21]，u[11][21]对应的数组z[11][21]，得到二元函数z=

，二次插值采用教材给出的分片二次代数插值。

3、曲面插值

利用x[11]，y[21]，z[11][21]建立二维函数表，进行曲面插值计算，逐步提高k值，计算其精度，看其是否满足要求，条件满足则循环结束，并得到曲面拟合的系数矩阵C[r][s]，算法采用教材给出的曲面拟合算法，求出所需矩阵给出，然后按公式进行计算。

4、比较

逼近

的效果

观察和

逼近

的效果时，只需要利用新的点列

重复计算二次代数插值，得到与新的插值节点

对应的

，再与对应的

比较即可，求解

事可以直接使用（3）中的C[r][s]和k值。

源程序（在CodeBlockC/C++集成开发环境下编译通过）

#include

#include

#include

#include

usingnamespacestd;

#defineEQUANUM4

#defineTHTA1E-12

#defineSIGMA1E-7

#defineITIME1000

voidprintMatrix（double*A,intr,ints）

{

for（inti=0;i

{

for（intj=0;j

{

cout<

}

cout<

}

}

voidgetFValue（doubleb[],doubleoffset[],doublex[]）

{

b[1]=-（0.5*cos（x[1]）+x[2]+x[3]+x[4]-offset[1]）;

b[2]=-（x[1]+0.5*sin（x[2]）+x[3]+x[4]-offset[2]）;

b[3]=-（0.5*x[1]+x[2]+cos（x[3]）+x[4]-offset[3]）;

b[4]=-（x[1]+0.5*x[2]+x[3]+sin（x[4]）-offset[4]）;

}

voidgetFDeri（doublecoef[][EQUANUM+1],doublex[]）

{

coef[1][1]=-0.5*sin（x[1]）;

coef[1][2]=1;

coef[1][3]=1;

coef[1][4]=1;

coef[2][1]=1;

coef[2][2]=0.5*cos（x[2]）;

coef[2][3]=1;

coef[2][4]=1;

coef[3][1]=0.5;

coef[3][2]=1;

coef[3][3]=-sin（x[3]）;

coef[3][4]=1;

coef[4][1]=1;

coef[4][2]=0.5;

coef[4][3]=1;

coef[4][4]=cos（x[4]）;

}

//解线性方程组

voidsolveEqua（double（*mat）[EQUANUM+1],doubleb[],doublex[]）

{

for（intk=1;k

{

inti=k;

for（intj=k;j<=EQUANUM;j++）

{

if（mat[j][k]>mat[i][k]）

i=j;

}

for（intj=k;j<=EQUANUM;j++）

{

doubletemp=mat[k][j];

mat[k][j]=mat[i][j];

mat[i][j]=temp;

}

doubletemp=b[k];

b[k]=b[i];

b[i]=temp;

for（inti=k+1;i<=EQUANUM;i++）

{

doublemik=mat[i][k]/mat[k][k];

for（intj=k+1;j<=EQUANUM;j++）

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

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

}

}

x[EQUANUM]=b[EQUANUM]/mat[EQUANUM][EQUANUM];

for（intk=EQUANUM-1;k>=1;k--）

{

x[k]=0;

for（intj=k+1;j<=EQUANUM;j++）

x[k]-=mat[k][j]*x[j];

x[k]+=b[k];

x[k]=x[k]/mat[k][k];

}

}

doublegetVecLen（doublevec[]）

{

doubletemp=0;

for（inti=1;i<=EQUANUM;i++）

if（temp

temp=fabs（vec[i]）;

returntemp;

}

voidNewton（doublex[],doubleoffset[]）

{

chara;

doubledetx[EQUANUM+1];

doubletempx[EQUANUM+1];

for（inti=1;i<=EQUANUM;i++）

x[i]=1;

doubleFDeri[EQUANUM+1][EQUANUM+1];

doubleFValue[EQUANUM+1];

for（inti=1;i

{

getFDeri（FDeri,x）;

getFValue（FValue,offset,x）;

solveEqua（FDeri,FValue,detx）;

if（getVecLen（detx）/getVecLen（x）<=THTA）

break;

else

for（inti=1;i<=EQUANUM;i++）

x[i]+=detx[i];

}

}

voidbiInter（doublet[][21],doubleu[][21],doublez[][21],intxs,intys）

{

doublett[6]={0,0.2,0.4,0.6,0.8,1};

doubleuu[6]={0,0.4,0.8,1.2,1.6,2};

doublezz[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}};

doubleh=0.2,tao=0.4;

intn=5,m=5;

for（inti=0;i

{

for（intj=0;j

{

intxp=0,yp=0;

if（t[i][j]<=tt[1]+h/2）

xp=1;

elseif（t[i][j]>tt[n-1]-h/2）

xp=n-1;

else

{

for（intq=2;q<=n-2;q++）

if（（t[i][j]>tt[q]-h/2）&&（t[i][j]<=tt[q]+h/2））

xp=q;

}

if（u[i][j]<=uu[1]+tao/2）

yp=1;

elseif（u[i][j]>uu[m-1]-tao/2）

yp=n-1;

else

{

for（intq=2;q<=m-2;q++）

if（（u[i][j]>uu[q]-tao/2）&&（u[i][j]<=uu[q]+tao/2））

yp=q;

}

z[i][j]=0;

for（intk=xp-1;k<=xp+1;k++）

{

doublelk=1;

for（intti=xp-1;ti<=xp+1;ti++）

{

if（ti!

=k）

lk*=（t[i][j]-tt[ti]）/（tt[k]-tt[ti]）;

}

for（intr=yp-1;r<=yp+1;r++）

{

doublelr=1;

for（intti=yp-1;ti<=yp+1;ti++）

{

if（ti!

=r）

lr*=（u[i][j]-uu[ti]）/（uu[r]-uu[ti]）;

}

z[i][j]+=lk*lr*zz[k][r];

}

}

}

}

}

voidmatrixMut（double*A,double*B,intr,ints,intt,double*result）

{

for（inti=0;i

{

for（intj=0;j

{

result[i*t+j]=0;

for（intk=0;k

result[i*t+j]+=A[i*s+k]*B[k*t+j];

}

}

}

voidtranspose（double*A,intr,ints,double*result）

{

for（inti=0;i

for（intj=0;j

result[j*r+i]=A[i*s+j];

}

voidcopyMatrix（double*A,intr,ints,intwd,double*result）

{

for（inti=0;i

for（intj=0;j

result[i*s+j]=A[i*wd+j];

}

voidinversion（double*A,intr,double*result）

{

doublecopyA[r*r];

copyMatrix（A,r,r,r,copyA）;

for（inti=0;i

for（intj=0;j

{

if（i!

=j）

result[i*r+j]=0;

else

result[i*r+j]=1;

}

for（inti=0;i

{

intindex=i;

for（intj=i+1;j

if（fabs（copyA[j*r+i]）>fabs（copyA[index*r+i]））

index=j;

if（index!

=i）

{

for（intj=0;j

{

doubletemp=copyA[index*r+j];

copyA[index*r+j]=copyA[i*r+j];

copyA[i*r+j]=temp;

temp=result[index*r+j];

result[index*r+j]=result[i*r+j];

result[i*r+j]=temp;

}

}

doubletemp=copyA[i*r+i];

if（temp!

=1）

{

for（intj=0;j

{

copyA[i*r+j]/=temp;

result[i*r+j]/=temp;

}

}

for（intj=0;j

{

if（（copyA[j*r+i]!

=0）&&（i!

=j））

{

temp=copyA[j*r+i];

for（intk=0;k

{

copyA[j*r+k]-=temp*copyA[i*r+k];

result[j*r+k]-=temp*result[i*r+k];

}

}

}

}

}

double*surFit（double*z,int&kvalue）

{

intxs=11,ys=21,num=9;

doubleB[xs*num];

doubleG[ys*num];

doubleP[xs*ys];

doubleTB[xs*ys];

doubleTG[xs*ys];

doubleBT[xs*ys];

doubleGT[xs*ys];

doubleBB[xs*ys];

double*GG=newdouble[xs*ys];

for（inti=0;i

{

for（intj=0;j

B[j*num+i]=pow（0.08*j,i）;

for（intj=0;j

G[j*num+i]=pow（0.5+0.05*j,i）;

}

doublesig=0;

for（inti=0;i

{

sig=0;

copyMatrix（B,xs,i+1,num,TB）;

transpose（TB,xs,i+1,BT）;

matrixMut（BT,TB,i+1,xs,i+1,BB）;

inversion（BB,i+1,TB）;

copyMatrix（G,ys,i+1,num,TG）;

transpose（TG,ys,i+1,GT）;

matrixMut（GT,TG,i+1,ys,i+1,GG）;

inversion（GG,i+1,GT）;

matrixMut（TB,BT,i+1,i+1,xs,BB）;

matrixMut（BB,z,i+1,xs,ys,GG）;

matrixMut（GG,TG,i+1,ys,i+1,BB）;

matrixMut（BB,GT,i+1,i+1,i+1,GG）;

for（intj=0;j

{

for（intk=0;k

{

doubletemp=0;

for（intp=0;p

for（intq=0;q

temp+=GG[p*（i+1）+q]*B[j*num+p]*G[k*num+q];

P[j*ys+k]=temp;

sig+=（z[j*ys+k]-temp）*（z[j*ys+k]-temp）;

}

}

cout<<"k="<

:

scientific|ios:

:

uppercase）<<"sigma="<

ofstreamout;

out.open（"shubiao1.txt",ios:

:

app）;

out<<"k="<

:

scientific|ios:

:

uppercase）<<"sigma="<

out.close（）;

if（sig

{

kvalue=i;

out.open（"shubiao2.txt",ios:

:

app）;

out<<"k="<

:

scientific|ios:

:

uppercase）<<"sigma="<

for（intk=0;k<=i;k++）

{

for（intj=0;j<=i;j++）

{

out<<"C["<

cout<<"C["<

}

}

out.close（）;

returnGG;

}

}

returnNULL;

}

intmain（）

{

doublet[11][21],u[11][21],z[11][21],offset[EQUANUM+1],x[EQUANUM+1];

intkvalue=0;

for（inti=0;i<=10;i++）

{

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

{

offset[1]=2.67+0.08*i;

offset[2]=1.07+0.5+0.05*j;

offset[3]=3.74+0.08*i;

offset[4]=0.79+0.5+0.05*j;

Newton（x,offset）;

t[i][j]=x[1];

u[i][j]=x[2];

}

}

biInter（t,u,z,11,21）;

ofstreamout;

out.open（"shubiao.txt"）;

for（inti=0;i<=10;i++）

{

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

{

cout<

:

scientific|ios:

:

uppercase）<<"x="<<0.08*i<<"y="<<0.5+0.05*j;

cout<

:

scientific|ios:

:

uppercase）<<"f（x,y）="<

out<

:

scientific|ios:

:

uppercase）<<"x="<<0.08*i<<"y="<<0.5+0.05*j;

out<

:

scientific|ios:

:

uppercase）<<"f（x,y）="<

}

}

out.close（）;

double*GG=surFit（z[0],kvalue）;

for（inti=1;i<=8;i++）

{

for（intj=1;j<=5;j++）

{

offset[1]=2.67+0.1*i;

offset[2]=1.07+0.5+0.2*j;

offset[3]=3.74+0.1*i;

offset[4]=0.79+0.5+0.2*j;

Newton（x,offset）;

t[i-1][j-1]=x[1];

u[i-1][j-1]=x[2];

}

}

biInter（t,u,z,8,5）;

doublep[8][5];

for（inti=1;i<=8;i++）

for（intj=1;j<=5;j++）

{

doubletemp=0;

for（intii=0;ii<=kvalue;ii++）

for（intjj=0;jj<=kvalue;jj++）

temp+=GG[ii*（kvalue+1）+jj]*pow（0.1*i,ii）*pow（0.5+0.2*j,jj）;

p[i-1][j-1]=temp;

}

out.open（"shubiao3.txt"）;

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

{

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

{

cout<

:

scientific|ios:

:

uppercase）;

cout<<"x["<

cout<

:

scientific|ios:

:

uppercase）;

cout<<"p（x,y）="<

cout<<"deta="<

out<

:

scientific|ios:

:

uppercase）;

out<<"x["<

out<

:

scientific|ios:

:

uppercase）;

out<<"p（x,y）="<

out<<"deta="<

}

}

out.close（）;

return0;

}

程序输出

1.数表（xi,yj）和f（xi,yj）

x=0y=0.5f（x,y）=4.46504018481E-001

x=0y=0.55f（x,y）=3.24683262928E-001

x=0y=0.6f（x,y）=2.10159686683E-001

x=0y=0.65f（x,y）=1.03043608316E-001

x=0y=0.7f（x,y）=3.40189556268E-003

x=0y=0.75f（x,y）=-8.87358136380E-002

x=0y=0.8f（x,y）=-1.73371632750E-001

x=0y=0.85f（x,y）=-2.50534611467E-001

x=0y=0.9f（x,y）=-3.20276506388E-001

x=0y=0.95f（x,y）=-3.82668066110E-001

x=0y=1f（x,y）=-4.37795766738E-001

x=0y=1.05f（x,y）=-4.85758941444E-001

x=0y=1.1f（x,y）=-5.26667254884E-001

x=0y=1.15f（x,y）=-5.60638479797E-001

x=0y=1.2f（x,y）=-5.87796538768E-001

x=0y=1.25f（x,y）=-6.08269779090E-001

x=0y=1.3f（x,y）=-6.22189452876E-001

x=0y=1.35f（x,y）=-6.29688378186E-001

x=0y=1.4f（x,y）=-6.30899760003E-001

x=0y=1.45f（x,y）=-6.25956152545E-001

x=0y=1.5f（x,y）=-6.14988546609E-001

x=0.08y=0.5f（x,y）=6.38015226511E-001

x=0.08y=0.55f（x,y）