数值分析B定.docx

数值分析B定.docx

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

数值分析B定

《数值分析B》计算实习题目三

一、算法设计方案

1.确定数表

，其中i=0,1,…,10；j=0,1,…,20。

a）将

代入非线性方程组解出

。

b）利用所给的

二维数表，用分片二次代数插值将上步得出的

值代入计算

的值。

c）由

得出数表

2.求

在区域

上的一个近似表达式

，其中r=0,1,…,k；s=0,1,…k。

a）以函数组

为基函数，构成以{

}为参数的曲面族

进行曲面拟合。

b）当曲面对数据的拟合精度

时停止拟合。

3.求

的值,其中i=1,2,…,8；j=1,2,…,5，以观察

逼近

的效果。

a）用第1步中确定数表

的方法将

代入求

。

b）将

代入第2步求出的

求

。

二、全部源程序

#include"stdio.h"

#include"math.h"

#include"stdlib.h"

#defineE11e-12//---------------------------------------------------------------迭代求解精度

#defineE21e-7//-----------------------------------------------------------------------拟合精度

intGauss（doublea[4][4],doubleb[4],doublec[4]）;//---------------------------函数声明

voidNewton（doublex,doubley,doubler[]）;

doublechazhi（doublet,doubleu）;

doublepxy（doublex,doubley,intz）;

intnihe（）;

intDoolittle（doublea[11][11],intr）;

voidhuidai1（doublea[11][11],doubleb[11],doublex[11],intr）;

//------------------------------------------------------------------------------------------变量声明

doubleX[11],Y[21];//-------------------------------------------------插值节点横,纵坐标

doublexyBiao[11][21];//------------------------------------------关于z、x、y的二维数表

doubleC[11][11];//---------------------------函数z=f（x,y）的近似表达式的系数矩阵

intmain（）//---------------------------------------------------------------------------------主函数

{

doublex,y,t,u;

doublex1[8],y1[5],fxy[8][5],r[4],q[4];

inti,j;

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

X[i]=0.08*i;

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

Y[j]=0.5+0.05*j;

printf（"\n数表:

xi,yi,f（xi,yi）：

\n"）;//-----------------建立关于z、x、y的二维数表

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

{

x=X[i];

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

{

y=Y[j];

q[0]=0.5;

q[1]=1.5;

q[2]=-0.3;

q[3]=1.5;

Newton（x,y,q）;//----------Newton法解非线性方程组从x，y得到t，u

t=q[0];

u=q[1];

xyBiao[i][j]=chazhi（t,u）;

printf（"x=%1.2f",X[i]）;

printf（",y=%1.2f\n",Y[j]）;

printf（"f（x,y）=%1.11e\n",xyBiao[i][j]）;

}

}

intk=nihe（）;

//-------------计算f（x[i],y[j]）,p（x[i],y[j]）的值，以观察p（x,y）逼近f（x,y）的效果

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

x1[i]=0.1*（i+1）;

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

y1[i]=0.5+0.2*（i+1）;

printf（"\n数表x*,y*,f（x*,y*）,p（x*,y*）：

\n"）;

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

{

x=x1[i];

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

{

y=y1[j];

r[0]=0.5;//-------------------------------------------------------迭代初始向量

r[1]=1.5;

r[2]=-0.3;

r[3]=1.5;

Newton（x,y,r）;//-----------Newton法解非线性方程组从x，y得到t，u

t=r[0];

u=r[1];

fxy[i][j]=chazhi（t,u）;

printf（"x*=%1.1f",x1[i]）;

printf（",y*=%1.1f\n",y1[j]）;

printf（"f（x*,y*）=%1.11e\n",fxy[i][j]）;

printf（"p（x*,y*）=%1.11e\n",pxy（x,y,k））;

}

}

return0;

}

//子函数一：

----------------------------------------------------------列主元素Gauss消去法

intGauss（doublea[4][4],doubleb[4],doublec[4]）

{

inti,j,k,r;

doublem[4],w[4][5];

doublet;

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

{

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

w[i][j]=a[i][j];

w[i][4]=b[i];

}

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

{

r=k;

t=fabs（w[k][k]）;

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

{

if（fabs（w[i][k]）>t）

{

r=i;

t=fabs（w[i][k]）;

}

}

if（r!

=k）

{

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

{

t=w[k][j];

w[k][j]=w[r][j];

w[r][j]=t;

}

}

if（0==w[k][k]）

{

return0;

}

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

{

m[i]=w[i][k]/w[k][k];

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

{

w[i][j]=w[i][j]-m[i]*w[k][j];

}

}

}

for（i=4-1;i>=0;i--）

{

c[i]=w[i][4];

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

c[i]-=w[i][j]*c[j];

c[i]/=w[i][i];

}

return1;

}

//子函数二：

-------------------------------Newton法解非线性方程组从x，y得到t，u

voidNewton（doublex,doubley,doubler[]）

{

doublec[4],s[4],a[4][4];

doublee,max;

inti,j,d=0,t=0,q=0;

do

{

d++;

for（i=0;i<4;i++）//---------------------------------------------初始化系数矩阵

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

a[i][j]=1.0;

a[2][0]=a[3][1]=0.5;

a[0][0]=-0.5*sin（r[0]）;

a[1][1]=0.5*cos（r[1]）;

a[2][2]=-sin（r[2]）;

a[3][3]=cos（r[3]）;

c[0]=2.67-0.5*cos（r[0]）-r[1]-r[2]-r[3]+x;

c[1]=1.07-r[0]-0.5*sin（r[1]）-r[2]-r[3]+y;

c[2]=3.74-0.5*r[0]-r[1]-cos（r[2]）-r[3]+x;

c[3]=0.79-r[0]-0.5*r[1]-r[2]-sin（r[3]）+y;

if（!

Gauss（a,c,s））

{

exit（-1）;

}

max=fabs（s[0]）;

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

{

if（fabs（s[i]）>max）

{

max=fabs（s[i]）;

t=i;

}

}

max=fabs（r[0]）;

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

{

if（fabs（r[i]）>max）

{

max=fabs（r[i]）;

q=i;

}

}

e=fabs（s[t]）/fabs（r[q]）;

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

r[i]+=s[i];

}while（e>E1）;

}

//子函数三：

---------------------------------------------------分片二次代数插值计算h（t,u）

doublechazhi（doublet,doubleu）

{

doubletuBiao[6][6]={-0.50,-0.34,0.14,0.94,2.06,3.50,-0.42,-0.50,-0.26,0.30,1.18,2.38,-0.18,-0.50,-0.50,-0.18,0.46,1.42,0.22,-0.34,-0.58,-0.50,-0.10,0.62,0.78,-0.02,-0.50,-0.66,-0.50,-0.02,1.50,0.46,-0.26,-0.66,-0.74,-0.50};//----------------------------------------------关于z、t、u的二维数表

doublet2[6]={0.0,0.2,0.4,0.6,0.8,1.0};

doubleu2[6]={0.0,0.4,0.8,1.2,1.6,2.0};

doublet1[3],u1[3];

intm,n,g,h;

m=int（10*t）%2==0?

5*t:

5*t+1;

n=int（5*t）%2==0?

2.5*t:

2.5*t+1;

m=0==m?

m+1:

m;

m=5==m?

m-1:

m;

n=0==n?

n+1:

n;

n=5==n?

n-1:

n;

for（g=m-1;g

{

t1[g-m+1]=1.0;

for（h=m-1;h

{

if（g!

=h）

{

t1[g-m+1]*=t-t2[h];

t1[g-m+1]/=t2[g]-t2[h];

}

}

}

for（g=n-1;g

{

u1[g-n+1]=1.0;

for（h=n-1;h

{

if（g!

=h）

{

u1[g-n+1]*=u-u2[h];

u1[g-n+1]/=u2[g]-u2[h];

}

}

}

doublez=0;

for（g=m-1;g

for（h=n-1;h

z+=t1[g-m+1]*u1[h-n+1]*tuBiao[g][h];

returnz;

}

//子函数四：

------------------------------------------------计算f（x,y）的近似表达式p（x,y）

doublepxy（doublex,doubley,intz）

{

doublep=0.0;

for（intr=0;r<=z;r++）

for（ints=0;s<=z;s++）

p+=C[r][s]*pow（x,r）*pow（y,s）;

returnp;

}

//子函数五：

----------------------------------------------------------------------------函数拟合

intnihe（）

{

inti,j,k;

doublewucha;

doubleA[11][21],B[11][11],V[11][11],D[21][11],G[21][11],Q[11][11];

printf（"\n选择过程的k和wucha值：

\n"）;

for（k=1;k<11;k++）

{

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

{

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

{

B[i][j]=1;

for（intt=0;t

B[i][j]*=X[i];

}

}

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

{

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

{

V[i][j]=0;

for（intt=0;t<11;t++）

V[i][j]+=B[t][i]*B[t][j];

}

}

if（!

Doolittle（V,k+1））

{

exit（-1）;

}

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

{

doublex[11],b[11];

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

{

b[i]=0;

for（intt=0;t<11;t++）

b[i]+=B[t][i]*xyBiao[t][j];

}

huidai1（V,b,x,k+1）;

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

{

A[i][j]=x[i];

}

}

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

{

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

{

G[i][j]=1;

for（intt=0;t

{

G[i][j]*=Y[i];

}

}

}

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

{

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

{

Q[i][j]=0;

for（intt=0;t<21;t++）

Q[i][j]+=G[t][i]*G[t][j];

}

}

if（!

Doolittle（Q,k+1））

{

exit（-1）;

}

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

{

doublex[11],b[11];

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

b[i]=G[j][i];

huidai1（Q,b,x,k+1）;

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

D[j][i]=x[i];

}

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

{

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

{

C[i][j]=0;

for（intt=0;t<21;t++）

C[i][j]+=A[i][t]*D[t][j];

}

}

wucha=0;

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

{

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

{

doublev=xyBiao[i][j]-pxy（X[i],Y[j],k）;

wucha+=v*v;

}

}

printf（"k=%d\n",k）;

printf（"wucha=%1.11e\n",wucha）;

if（wucha

{

printf（"\n达到精度要求时k和wucha值：

\n"）;

printf（"k=%d\n",k）;

printf（"wucha=%1.11e\n",wucha）;

break;

}

}

printf（"\n达到精度要求时p（x,y）中的系数C[i][j]：

\n"）;

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

{

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

printf（"C[%d][%d]=%1.11e\n",i,j,C[i][j]）;

}

returnk;

}

//子函数六：

----------------------------------------------------------------------Doolittle分解

intDoolittle（doublea[11][11],intr）

{

inti,j,k,t;

for（k=0;k

{

for（j=k;j

{

for（t=0;t

{

a[k][j]-=a[k][t]*a[t][j];

}

}

if（0==a[k][k]）

{

return0;

}

for（i=k+1;i

{

for（t=0;t

{

a[i][k]-=a[i][t]*a[t][k];

}

a[i][k]/=a[k][k];

}

}

return1;

}

//子函数七：

--------------------------------------------求Ly=b得y，进而由Ux=y得x

voidhuidai1（doublea[11][11],doubleb[11],doublex[11],intr）

{

inti,j;

doubley[11];

for（i=0;i

{

y[i]=b[i];

for（j=0;j

y[i]-=a[i][j]*y[j];

}

for（i=r-1;i>=0;i--）

{

x[i]=y[i];

for（j=i+1;j

{

x[i]-=a[i][j]*x[j];

}

x[i]/=a[i][i];

}

}

三、运算结果