数值计算方法上机实习题答案(最新整理)

fprintf(' x y(1) y1
%norm(v-v0)>ep k=k+1;
数值计算方法上机作业
q=v; u=v/norm(v,inf) v=B*u; v0=q; end mt=1/norm(v,inf)+p my=u 主界面中输入：A=[1 -2 -3]; maxtr(A,11,0.001) 结果为： 特征值： mt =
11.0919
（2）
y1 y2
2 y1 y2 2 sin 998 y1 999 y2
x 999
cos
x
999
sin
x
，
x
[0,10]
，
y1 (0) y2 (0)
2 3
。
和精确解
y1
(
x)
y2 (x)
2e x 2e x
sin x cos x
比较，分析结论。
程序：function ydot=lorenzeq1(x,y)
（1） JACOBI 迭代； 程序： function y=jacobi(a,b,x0) D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1); B=D\(L+U); f=D\b; y=B*x0+f;n=1; while norm (y-x0)>1e-4
x0=y; y=B*x0+f;n=n+1; end y n 以文件名 jacobi.m 保存。 程序： a=[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]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]'; jacobi(a,b,x0);
(5) 20
E0
I 20
，所此递推式不可靠。而在第二种递推式中
E0
1 5
E1
(
1)n 5
En
，误差在缩小，
所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，
即算法是否数值稳定。
2． 求方程 e x 10x 2 0 的近似根，要求 xk1 xk 5 104 ，并比较计算量。
1.0000 2.0000 1.0000 2.0000 1.0000 2.0000
数值计算方法上机作业
热能工程
n=
15
（3） SOR 迭代（ 1.334, 1.95, 0.95 ）。
程序： function y=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; while norm(y-x0)>10^(-4)
x0=y; y=lw*x0+f;n=n+1; end y n 以文件名 sor.m 保存。 程序： a=[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]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]'; c=[1.334 1.95 0.95]; for i=1:3 w=c(i); sor(a,b,w,x0); end 运行结果分别为： y=
热能工程
特征向量： my =
0.3845 -1.0000 0.7306 6．用经典 R-K 方法求解初值问题
（1）
y1 y2
2 y1
y1
y2 2sin x 2 y2 2 cos x
wenku.baidu.com
2
sin
x
，
x
[0,10]
，
y1 y2
(0) (0)
2 3
；
程序：function ydot=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=
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
for n=1:20
I=(-1/5)*I+1/(5*n);
end
I
I 0 = 0.0083
（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
首先分析两种递推式的误差；设第一递推式中开始时的误差为
E0
I0
I
0
，递推过程的舍入误差不计。
并记 En I n I n ，则有 En 5En1 (5)n E0 。因为 E20
a=y(i); b=x(i); c=c+(a-f(b))^2;
end averge=c/15 结果：averge =
0.0594
4 1 0 1 0 0
0
1 4 1 0 1 0
5
4．设
A
0
1
4
1
0
1
，
b
2
，
Ax
b
1 0 1 4 1 0
5
0 0
1 0
0 1
1 0
4 1
41
2 6
分析下列迭代法的收敛性，并求 xk1 xk 2 104 的近似解及相应的迭代次数。
（1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4
c=(b+a)/2;
数值计算方法上机作业
if exp(c)+10*c-2>0 b=c;
else a=c; end end c 结果：c =
0.0903
（2）
取初值 x0
0 ，并用迭代 xk1
2 ex 10
数值计算方法上机作业
热能工程
1． 设 I n
1 xn dx ，
0 5 x
（1）
由递推公式 I n
5I n1
1 n ，从 I 0 的几个近似值出发，计算 I20 ；
解：易得： I0 ln6-ln5=0.1823,
程序为： I=0.182; for n=1:20
I=(-5)*I+1/n; end I
ydot=[-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)];
以文件名 lorenzeq1.m 保存。
程序：x=0:10;
y1=2*exp(-x)+sin(x);
y2=2*exp(-x)+cos(x);
[x,y]=ode45('lorenzeq1',[0:10],[2;3]);
b
设 y=f(x)具有指数形式 y ae x （a>0,b<0）。对此式两边取对数，得 ln y ln a b 1 。记 A=lna，B=b，
x
并引入新变量 z=lny，t=1/x。引入新变量后的数据表如下
x
2
3
4
5
6
7
t=1/x 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429
输出结果为： I20 = -3.0666e+010
（2）
粗糙估计 I 20
，用 I n1
1 5
I n1
1 5n
，计算 I 0 ；
因为
0.0079
1 0
x 20 dx
6
I 20
1 0
x 20 dx 0.0095 5
所以取 I 20
1 (0.0079 0.0095) 2
0.0087
程序为：I=0.0087;
程序：
t=[0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667
0.0625];
z=[1.8594 2.1041 2.2597 2.2513 2.2721 2.3026 2.2956 2.3016 2.3504 2.3599 2.3609 2.3795 2.3609 2.3888
z=lny 1.8594 2.1041 2.2597 2.2513 2.2721 2.3026
8 0.1250 2.2956
9 0.1111 2.3016
10
11
12
13
14
15
16
0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625
2.3504 2.3599 2.3609 2.3795 2.3609 2.3888 2.3758
2.3758];
polyfit(t,z,1)
结果：
ans = -1.1107 2.4578
由此可得 A=2.4578，B=-1.1107， a e A 11.6791 ，b=B=-1.1107
1.1107
方程即为 y 11.6791e x
数值计算方法上机作业
热能工程
计算均方差编程： x=[2:16]; y=[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]; f(x)=11.6791*exp( -1.1107./x); c=0; for i=1:15
n=
17
6 5．用逆幂迭代法求 A 3
3 2
1 1
最接近于
11
的特征值和特征向量，准确到10
3
。
1 1 1
程序： function [mt,my]=maxtr(A,p,ep) n=length(A); B=A-p*eye(n); v0=ones(n,1); k=1; v=B*v0; while abs(norm(v,inf)-norm(v0,inf))>ep
结果：x =
0.0995
（4） 取初值 x0 0 ，并用牛顿迭代法；
程序：x=0; a=0;b=1; while abs(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)+10*x-2=0')
数值计算方法上机作业
运行结果为: y=
1.0000 2.0000 1.0000 2.0000 1.0000 2.0000
热能工程
n=
28
（2） GAUSS-SEIDEL 迭代； 程序： function y=seidel(a,b,x0) 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; while norm(y-x0)>10^(-4)
3．钢水包使用次数多以后，钢包的容积增大，数据如下：
x 2 3 4 567 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 试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）
x0=y; y=G*x0+f;n=n+1; end y n 以文件名 deisel.m 保存。 程序： a=[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]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]'; jacobi(a,b,x0); 运行结果为： y=
0 1 2 3 4 5 6 7 8 9 10 y=
数值计算方法上机作业
热能工程
2.0000 1.5775 1.1802 0.2406 -0.7202 -0.9454 -0.2745 0.6589 0.9901 0.4124 -0.5440
3.0000 1.2758 -0.1457 -0.8903 -0.6170 0.2971 0.9652 0.7557 -0.1449 -0.9109 -0.8389
；
程序：x=0; a=1; while abs(x-a)>5*1e-4
a=x; x=(2-exp(x))/10; end x 结果：x =
0.0905
（3） 加速迭代的结果； 程序：x=0; a=0;b=1; while abs(b-a)>5*1e-4
a=x; y=exp(x)+10*x-2; z=exp(y)+10*y-2; x=x-(y-x)^2/(z-2*y+x); b=x; end x