数值计算方法上机实习题答案

数值计算方法上机作业热能工程

1.设I=f1X，'dx,

nn5+X

(1)由递推公式I=—5I+1，从入的几个近似值出发，计算L；

nn—1fi020

解：易得：I=ln6-ln5=0.1823,

o

程序为：

1=0.182;

forn=1:20

l=(-5)*l+1/n;

end

I

输出结果为：I=-3.0666e+010

20

(2)粗糙估计I”，用ii=—计算J

2055/70

因为0.0079xf1-^-dx</<A«0.0095

n620c5

所以取、=1(0.0079+0.0095；=0.0087

程序为:1=0.0087;

forn=1:20

l=(-1/5)*l+1/(5*n);

end

I

I=0.0083

o

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

首先分析两种递推式的误差；设第一递推式中开始时的误差为E=l—I',递推过程的舍入误差不计。

000

并记E=I—I',则有E=-5E=•••=(-5)nEo因为E=

nnnnn—1020

(-5)20E»l,所此递推式不可靠。而在第二种递推式中E=-1E=-=(-l)nE,误差在缩小,

020051sn

所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，

即算法是否数值稳定。

2.求方程ex+10x—2=0的近似根，要求|x-x|<5x10-4,并比较计算量。

(1)在［0，1］上用二分法;

程序：a=0;b=1.0;

whileabs(b-a)>5*1e-4

c=(b+a)/2;

数值计算方法上机作业热能工程

ifexp(c)+10*c-2>0

b=c;

elsea=c;

end

end

c

结果:c=

0.0903

2-gx

(2)取初值X=0,并用迭代X=-----------

0k+110

程序：x=0;

a=1；

whileabs(x-a)>5*1e-4

a=x;

x=(2-exp(x))/10;

end

x

结果：x=

0.0905

(3)加速迭代的结果；

程序:x=0;

a=0;b=1;

whileabs(b-a)>5*1e-4

a=x;

y=exp(x)+10*x-2;

z=exp(y)+10*y-2;

x=x-(y-x)A2/(z-2*y+x);

b=x;

end

结果：x=

0.0995

(4)取初值x=0,并用牛顿迭代法;

o

程序:x=0;

a=0;b=1;

whileabs(b-a)>5*1e-4

a=x;

数值计算方法上机作业热能工程

x=x-(exp(x)+10*x-2)/(exp(x)+10);

b=x;

end

x

结果:

x=

0.0905

(5)分析绝对误差。

solve('exp(x)+l0Ax-2=0')

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

X23456789

y6.428.29.589.59.7109.939.99

10111213141516

10.4910.5910.60J0.810.610.9J0.76

试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）

。1

设y=f（x）具有指数形式y=aex（a>0,b<0）。对此式两边取对数，得lny=lna+b—。记A=lna，B=b,

并引入新变量z=lny,t=1/xo引入新变量后的数据表如下

X23456789

t=1/x0.50000.33330.25000.20000.16670.14290.12500.1111

z=lny1.85942.10412.25972.25132.27212.30262.29562.3016

10111213141516

0.10000.09090.08330.07690.07140.06670.0625

2.3758I

2.35042.35992.36092.37952.36092.3888

程序：

t=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.0667

0.0625];

z=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.3888

2.3758];

polyfit(t,z,1)

结果：

ans=-1.11072.4578

由此可得A=2.4578,B=-1.1107,a=eA=11.6791,b=B=-1.1107

1.1107

方程即为y=11.6791ex

数值计算方法上机作业热能工程

运行结果为:

y=

1.0000

2.0000

1.0000

2.0000

1.0000

2.0000

n=

28

(2)GAUSS-SEIDEL迭代；

程序：

functiony=seidel(a,b,xO)

D=diag(diag(a));

U=-triu(a,1);

L=-tril(a.-1)；

G=(D-L)\U;

f=(D-L)\b;

y=G*x0+f;n=1;

whilenorm(y-x0)>10A(-4)

x0=y;

y=G*x0+f;n=n+1;

end

y

n

以文件名deisel.m保存。

程序：

a=[4-10-100；-14-10-10；0-14-10-1；-10-14-10；0-10-14-1;00-10-14];

b=[05-25-26],;

x0=[000000],;

jacobi(a.b.xO);

运行结果为：

y=

1.0000

2.0000

1.0000

2.0000

1.0000

2.0000

数值计算方法上机作业热能工程

n=

15

(3)SOR迭代(①=1.334,1.95,0.95)。

程序：

functiony=sor(a,b,w,x0)

D=diag(diag(a));

U=-triu(a,1);

L=-tril(a,-1)；

lw=(D-w*L)\((1-w)*D+w*U);

f=(D-w*L)\b*w;

y=lw*x0+f；n=1;

whilenorm(y-x0)>10A(-4)

x0=y;

y=lw*x0+f;n=n+1;

end

y

n

以文件名sor.m保存。

程序：

a=[4-10-100；-14-10-10；0-14-10-1；-10-14-10；0-10-14-1;00-10-14]；

b=[05-25-26],;

x0=[000000],;

c=[1.3341.950.95];

fori=1:3

w=c(i);

sor(a,b,w,x0);

end

运行结果分别为：

y=

1.0000

2.0000

1.0000

2.0000

1.0000

2.0000

n=

数值计算方法上机作业热能工程

13

y=

1.0000

2.0000

1.0000

2.0000

1.0000

2.0000

n=

241

y=

1.0000

2.0000

1.0000

2.0000

1.0000

2.0000

n=

17

(631)

5.用逆幕迭代法求人=321最接近于11的特征值和特征向量，准确到10-3。

Uli

程序：

function[mt,my]=maxtr(A,p,ep)

n=length(A);

B=A-p*eye(n);

v0=ones(n,1);

k=1；

v=B*vO;

whileabs(norm(v,inf)-norm(vO,inf))>ep

%norm(v-v0)>ep

k-k+1;

数值计算方法上机作业热能工程

q=v;

u=v/norm(v,inf)

v=B*u;

vO=q;

end

mt=1/norm(v,inf)+p

my=u

主界面中输入:A=[1-2-3];maxtr(A,11,0.001)

结果为：

特征值：

mt=

11.0919

特征向量：

my=

0.3845

-1.0000

0.7306

6.用经典R-K方法求解初值词题

y"=-2y+y+2sinAfy(0)=2

⑴1二乙,xe[0,10].f*/c

y?=-2y,+2cosx-2sinx[y^(0)=3

程序：functionydot=lorenzeq(x,y)

ydot=[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)]

以文祥民lorenzeq.m保存。

主窗口输入：[x,y]=ode45('lorenzeq',[0:10],[2;3])

运行结果为：

x=

0

1

2

3

4

5

6

7

8

9

10

y=

数值计算方法上机作业热能工程

2.00003.0000

1.57751.2758

1.1802-0.1457

0.2406-0.8903

-0.7202-0.6170

-0.94540.2971

-0.27450.9652

0.65890.7557

0.9901-0.1449

0.4124-0.9109

-0.5440-0.8389

yf=-2y+y+2sin.xy(0)=2

(2)'；।2"99cosx-999sin,xG[0,10],“

y：=998y「999y,+9xy(0)=3

2

和精确解〈乙y(x)=)2e-x+sinx比较，分析结论。

y(X)=2e-x+cosA

程序：functionydot=lorenzeq1(x,y)

ydot=[-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)];

以文件名lorenzeql.m保存。

程序:x=0:10;

y1=2*exp(-x)+sin(x);

y2=2*exp(-x)+cos(x);

[x,y]=ode45('lorenzeq11,[0:10],|2;3]);

fprintffxy(i)y!y(2)y2\n')

forj=1:length(y)

fprintf('%4d%.4f1%.4f%.4f%.4f\n：x(j),y(j,1),y1(j),y(j,2),y2(j))

end

运行结果为：

Xy(i)yiy(2)y2

02.00002.00003.00003.0000

11.57721.57721.27591.2761

21.18001.1800-0.1455-0.1455

30.24070.2407-0.8904-0.8904

4-0.7202-0.7202-0.6169-0.6170

5-0.9454-0.94540.29720.2971

6-0.2745-0.27450.96480.9651

70.65880.65880.75540.7557

80.99000.9900-0.1448-0.1448

90.41240.4124-0.9106-0.9109

10-0.5439-0.5439-G.8389-0.8390

结论：R-K方法求解的结果精度较高。

7.用有限差分法求解边值问题(h=0.1)：

数值计算方法上机作业热能工程

y〃_(l+x2))，=0

y(-i)=y(i)=i

程序为：

h=0.1;

n=(1-(-1))/h+1;

x(1)=-1;x(n-1)=1;

y(1)=1;y(n-1)=1;

fori=1:n-1

x(i)=x(1)+(i-1)*h;

q(i)=(1+x(i)A2);

B(i)=2/(M2)+q(i);

end

fori=1:n-2

C(i)=-1/(hA2);

end

H=diag(B)+diag(C,1)+diag(C,-i)；

g⑴=0+"(M2);

g(n-1)=0+1/(M2);

fori=2:n-2

g(i)=o；

end

y=H\g'

运行结果为：

y=

0.9027

0.8235

0.7592

0.7074

0.6661

0.6338

0.6095

0.5922

0.5814

0.5767

0.5778

0.5846

0.5974

0.6163

0.6420

0.6752

0.7167

0.7680

0.8308

数值计算方法上机作业