数值计算上机题答案.docx

数值计算上机题答案.docx

《数值计算上机题答案.docx》由会员分享，可在线阅读，更多相关《数值计算上机题答案.docx（36页珍藏版）》请在冰豆网上搜索。

数值计算上机题答案

数值计算

1．设

，

（1）由递推公式

，从

的几个近似值出发，计算

；

（2）粗糙估计

，用

，计算

；

（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。

（1）

#include

usingnamespacestd;

intmain（）

{

floati=0.1;

floatn=0,ins;

intk;

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

{

n=n+1;

ins=1/n;

i=-5*i+ins;

cout<

}

return0;

}

初值取0.1时，答案为：

0.5

-2

10.3333

-51.4167

257.283

-1286.25

6431.39

-32156.8

160784

-803921

4.01961e+006

-2.0098e+007

1.0049e+008

-5.02451e+008

2.51225e+009

-1.25613e+010

6.28064e+010

-3.14032e+011

1.57016e+012

-7.8508e+012

Pressanykeytocontinue

初值取0.2时，答案为：

0.5

-2.16667

11.0833

-55.2167

276.25

-1381.11

6905.66

-34528.2

172641

-863205

4.31603e+006

-2.15801e+007

1.07901e+008

-5.39503e+008

2.69752e+009

-1.34876e+010

6.74379e+010

-3.3719e+011

1.68595e+012

Pressanykeytocontinue

出现误差的原因在与程序中出现截断误差的原因导致的。

（2）

#include

usingnamespacestd;

intmain（）

{

floati=0.1,n=0,ins;

intk;

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

{

n=n+1;

ins=0.2/n;

i=-0.2*i+ins;

}

cout<

return0;

}

答案如下:

0.0078715

Pressanykeytocontinue

（3）

结果还是相对比较可靠的，结果是收敛的。

产生这种现象的原因是就是初始值取的比较精确，使得不断地递推过程中，答案不断地收敛。

2．求方程

的近似根，要求

，并比较计算量。

（1）在[0，1]上用二分法；

（2）取初值

，并用迭代

；

（3）加速迭代的结果；

（4）取初值

，并用牛顿迭代法；

（5）分析绝对误差。

程序如下：

（1）

#include

#include

#defineE2.71828

usingnamespacestd;

intmain（）

{

floata=0,b=1;

floatf;

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

{

f=pow（E,（（a+b）/2））+10*（a+b）/2-2;

if（f<0）

{

a=（a+b）/2;

b=b;

}

elseif（f>0）

{

a=a;

b=（a+b）/2;

}

else

break;

}

cout<

return0;

}

迭代20次，可以得到满足要求的解0.0905251

（2）

#include

#include

usingnamespacestd;

#defineE2.718282

intmain（）

{

floatx=0;

for（intn=0;n<=20;n++）

{

x=（2-pow（E,x））/10;

}

cout<

return0;

}

答案如下：

0.0905251

Pressanykeytocontinue

此为迭代20次后产生的结果。

（3）

#include

#include

#defineE2.71828

usingnamespacestd;

intmain（）

{

floatx=0.5,xi;

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

{

xi=（2-pow（E,x））/10;

x=1.42857*xi-0.42857*x;

}

cout<

return0;

}

迭代10次得到的结果如下：

0.0925289

Pressanykeytocontinue

（4）

#include

#include

usingnamespacestd;

#defineE2.718282

intmain（）

{

floatx;

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

{

x=x-（pow（E,x）+10*x-2）/（pow（E,x）+10）;

}

cout<

return0;

}

答案如下：

0.0905251

Pressanykeytocontinue

此为迭代5次后产生的结果。

（5）根据二分法误差公式可得，绝对误差小于0.000000042857.、

3．钢水包使用次数多以后，钢包的容积增大，数据如下：

x

2

3

4

5

6

7

8

9

y

6.42

8.2

9.58

9.5

9.7

10

9.93

9.99

10

11

12

13

14

15

16

10.49

10.59

10.60

10.8

10.6

10.9

10.76

试从中找出使用次数和容积之间的关系，计算均方差。

（注：

增速减少，用何种模型）

解题思想：

另拟合曲线的方程式为y=a*exp（bx）;

转换成线性方程有Y=A+bx;（其中，Y=Iny,A=Ina,最后通过最小二乘法计算A,b,转而计算a,b）

计算a,b的程序如下：

#include

#include

#defineN2

#defineE2.71828

usingnamespacestd;

voidchange（float（*）[N+1],float（*）[N+1]）;

intmain（）

{

floaty[15]=

{

6.42,8.2,9.58,

9.5,9.7,10,

9.93,9.99,10.49,

10.59,10.60,10.8,

10.6,10.9,10.76

};

floatx[15],Y[15],a,A,b;

floatarray[2][3]={0};

inti,j,k,m;

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

{

x[i]=i+2;

}

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

{

Y[i]=log（y[i]）;

}

array[0][0]=15;

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

{

array[0][1]+=x[i];

}

array[1][0]=array[0][1];

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

{

array[1][1]+=pow（x[i],2）;

}

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

{

array[0][2]+=Y[i];

}

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

{

array[1][2]+=Y[i]*x[i];

}

float（*p）[N+1]=&array[0];

float*ptr[N];

//使对角元素最大

for（k=0;k

{

for（i=0;i

{

if（abs（*（*（p+k）+k））

change（p+k,p+k+i）;

}

}

cout<<"主对角元素最大的矩阵整理后如下:

"<

for（k=0;k

{

for（i=0;i

{

cout<<*（*（p+k）+i）<<"\t";

}

cout<

}

//开始计算

for（j=0;j

{

ptr[j]=array[j];

}

*ptr[0]=1/（*ptr[0]）;

for（i=1;i

*（ptr[0]+i）=*ptr[0]**（ptr[0]+i）;

for（j=1;j

{

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

{

for（i=1;i

{

*（ptr[j]+i+k）=*（ptr[j]+i+k）-*（ptr[j]+k）**（ptr[k]+i+k）;

}

}

*（ptr[j]+j）=1/（*（ptr[j]+j））;

for（m=1;m

{

*（ptr[j]+j+m）=*（ptr[j]+j+m）**（ptr[j]+j）;

}

}

//此为修正后的方程,高斯消去法的变态形式,以便于以后直接迭代用

cout<<"因子表如下:

"<

for（j=0;j

for（i=0;i

{

cout<<*（ptr[j]+i）<<"";

if（i==N）

cout<

}

cout<

//前代过程如下,矩阵变态形式不变,但是增广项发生了变化

for（j=0;j

{

for（i=0;i

{

*（ptr[j]+N）=*（ptr[j]+N）-*（ptr[j]+i）**（ptr[i]+N）;

}

*（ptr[j]+N）=*（ptr[j]+j）**（ptr[j]+N）;

}

cout<<"前代后的矩阵如下（因子表不变,增广项发生了变化）:

"<

for（j=0;j

for（i=0;i

{

cout<<*（ptr[j]+i）<<"";

if（i==N）

cout<

}

cout<

//前代过程结束后,进行回代过程,可以方便的求出方程的最终解

for（j=N-1;j>=0;j--）

{

for（i=0;i

{

*（ptr[j]+N）=*（ptr[j]+N）-*（ptr[j]+N-1-i）**（ptr[N-1-i]+N）;

}

}

cout<<"得到的解如下:

"<

A=*（ptr[0]+N）;

b=*（ptr[1]+N）;

a=pow（E,A）;

b=b;

cout<

return0;

}

voidchange（float（*a）[N+1],float（*b）[N+1]）

{

intins,i;

for（i=0;i

{

ins=*（*a+i）;

*（*a+i）=*（*b+i）;

*（*b+i）=ins;

}

}

计算结果如下：

主对角元素最大的矩阵整理后如下:

1351495314.78

1513534

因子表如下:

0.0074074111.0741314.78

15-0.032142934

前代后的矩阵如下（因子表不变,增广项发生了变化）:

0.0074074111.07412.33171

15-0.03214290.0313584

得到的解如下:

7.27497

0.0313584

Pressanykeytocontinue

即有a=7.27497;b=0.0313584

则拟合方程为y=7.27497exp0.0313584*x;

4．设

，

，

分析下列迭代法的收敛性，并求

的近似解及相应的迭代次数。

（1）JACOBI迭代；

（2）GAUSS-SEIDEL迭代；

（3）SOR迭代（

）。

（1）

雅克比迭代法的程序如下：

#include

usingnamespacestd;

intmain（）

{

floata[6][6]=

{

{4,-1,0,-1,0,0},

{-1,4,-1,0,-1,0},

{0,-1,4,-1,0,-1},

{-1,0,-1,4,-1,0},

{0,-1,0,-1,4,-1},

{0,0,-1,0,-1,4}

};

floatb[6]={0,5,-2,5,-2,6};

floatx[6]={0};//设x解列向量的初值为0

floatsum;

inti,j,k;

for（k=1;k<=10;k++）//计算10次

{

for（i=0;i<=5;i++）//每循环一次计算每一个x

{

sum=0;

for（j=0;j<=5;j++）//计算sum

{

sum+=a[i][j]*x[j];

}

x[i]=x[i]+（1/a[i][i]）*（b[i]-sum）;

}

}

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

cout<<"x["<

return0;

}

迭代10次得到的解为：

x[0]=0.998445

x[1]=1.99831

x[2]=0.998577

x[3]=1.9987

x[4]=0.998928

x[5]=1.99938

Pressanykeytocontinue

（2）

Gauss—Seidel迭代法程序：

#include

usingnamespacestd;

intmain（）

{

floata[6][6]=

{

{4,-1,0,-1,0,0},

{-1,4,-1,0,-1,0},

{0,-1,4,-1,0,-1},

{-1,0,-1,4,-1,0},

{0,-1,0,-1,4,-1},

{0,0,-1,0,-1,4}

};

floatb[6]={0,5,-2,5,-2,6};

floatx[6]={0};//设x解列向量的初值为0

floatsuma,sumb;

inti,j,k;

for（k=1;k<=10;k++）//计算10次

{

for（i=0;i<=5;i++）//每循环一次计算每一个x

{

suma=0;

for（j=0;j<=i-1;j++）//计算sum

{

suma+=a[i][j]*x[j];

}

sumb=0;

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

{

sumb+=a[i][j]*x[j];

}

x[i]=x[i]+（1/a[i][i]）*（b[i]-suma-sumb）;

}

}

//输出结果

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

cout<<"x["<

return0;

}

程序结果为：

x[0]=0.998445

x[1]=1.99831

x[2]=0.998577

x[3]=1.9987

x[4]=0.998928

x[5]=1.99938

Pressanykeytocontinue

（3）SOR程序将1/a[i][i]前乘上分别乘上加速系数，依题目要求为1.334，1.95，0.95

则求得的解分别为：

w=1.334时：

x[0]=1.00029

x[1]=1.99989

x[2]=1.00015

x[3]=2.00012

x[4]=0.999862

x[5]=2.00009

Pressanykeytocontinue

w=1.95时：

x[0]=2.04027

x[1]=2.40694

x[2]=1.32038

x[3]=1.24714

x[4]=0.429775

x[5]=2.42788

Pressanykeytocontinue

w=0.95时：

x[0]=0.996546

x[1]=1.99615

x[2]=0.996684

x[3]=1.99693

x[4]=0.99741

x[5]=1.99845

Pressanykeytocontinue

5．用逆幂迭代法求

最接近于11的特征值和特征向量，准确到

。

可以对矩阵A进行LU分解（doolittle分解）

LU分解的程序如下：

#include

usingnamespacestd;

#defineN3

intmain（）

{

floata[N][N]=

{

{6,3,1},

{3,2,1},

{1,1,1},

};

floatl[N][N]={0},u[N][N]={0},sum=0;

inti,j,k,r;

for（i=0;i

{

l[i][i]=1;

}

for（k=0;k

{

for（j=k;j

{

sum=0;

for（r=0;r

{

sum+=l[k][r]*u[r][j];

}

u[k][j]=a[k][j]-sum;

}

for（i=k+1;i

{

sum=0;

for（r=0;r

{

sum+=l[i][r]*u[r][k];

}

l[i][k]=（a[i][k]-sum）/u[k][k];

}

}

//输出l,u矩阵结果

cout<<"L矩阵的结果如下:

"<

for（i=0;i

for（j=0;j

{

cout<

if（j==N-1）

cout<

}

cout<<"U矩阵的结果如下:

"<

for（i=0;i

for（j=0;j

{

cout<

if（j==N-1）

cout<

}

return0;

}

输出结果如下:

L矩阵的结果如下:

100

0.510

0.16666711

U矩阵的结果如下:

631

00.50.5

000.333333

Pressanykeytocontinue

然后进行规范化计算，用一下程序计算Y,u的解：

我用直接求解法，不断迭代，程序如下：

#include

#include

usingnamespacestd;

#defineN3

voidchange（double（*）[N+1],double（*）[N+1]）;

intmain（）

{

doublearray[N][N+1]=

{

{6,3,1,1},

{3,2,1,1},

{1,1,1,1}

};

doublemk;

double（*p）[N+1]=&array[0];

inti,j,k,m,n;

double*ptr[N];

for（n=1;n<=10;n++）

{

//使对角元素最大

for（k=0;k

{

for（i=0;i

{

if（abs（*（*（p+k）+k））

change（p+k,p+k+i）;

}

}

//开始计算

for（j=0;j

{

ptr[j]=array[j];

}

*ptr[0]=1/（*ptr[0]）;

for（i=1;i

*（ptr[0]+i）=*ptr[0]**（ptr[0]+i）;

for（j=1;j

{

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

{

for（i=1;i

{

*（ptr[j]+i+k）=*（ptr[j]+i+k）-*（ptr[j]+k）**（ptr[k]+i+k）;

}

}

*（ptr[j]+j）=1/（*（ptr[j]+j））;

for（m=1;m

{

*（ptr[j]+j+m）=*（ptr[j]+j+m）**（ptr[j]+j）;

}

}

//前代过程如下,矩阵变态形式不变,但是增广项发生了变化

for（j=0;j

{

for（i=0;i

{

*（ptr[j]+N）=*（ptr[j]+N）-*（ptr[j]+i）**（ptr[i]+N）;

}

*（ptr[j]+N）=*（ptr[j]+j）**（ptr[j]+N）;

}

//前代过程结束后,进行回代过程,可以方便的求出方程的最终解

for（j=N-1;j>=0;j--）

{

for（i=0;i

{

*（ptr[j]+N）=*