（完整版）数值计算方法上机实习题答案

（完整版）数值计算⽅法上机实习题答案
1．设?+=1
05dx x
x I n
n ，（1）由递推公式n
I I n n 1
51+-=-，从0I 的⼏个近似值出发，计算20I ；解：易得：0I =ln6-ln5=0.1823, 程序为：
I=0.182; for n=1:20
I=(-5)*I+1/n; end I
输出结果为：20I = -3.0666e+010 （2）粗糙估计20I ，⽤n
I I n n 51
5111+-
=--，计算0I ；因为 0095.05
6 0079.01020
201
020
≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2
1
20=+=
I 程序为：I=0.0087; for n=1:20
I=(-1/5)*I+1/(5*n); end I
0I = 0.0083
（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。

⾸先分析两种递推式的误差；设第⼀递推式中开始时的误差为000I I E '-=，递推过程的舍⼊误差不计。

并记n
n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。

因为=20E 20020)5(I E >>-，所此递推式不可靠。

⽽在第⼆种递推式中n n E E E )5
1(5110-==-=Λ，误差在缩⼩，
所以此递推式是可靠的。

出现以上运⾏结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，
即算法是否数值稳定。

2．求⽅程0210=-+x e x
的近似根，要求4
1105-+?<-k k x x ，并⽐较计算量。

（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）取初值00=x ，并⽤迭代10
21x k e x -=+；
程序：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
结果：x =
0.0995
（4）取初值00=x ，并⽤⽜顿迭代法；程序：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')
3．钢⽔包使⽤次数多以后，钢包的容积增⼤，数据如下：
试从中找出使⽤次数和容积之间的关系，计算均⽅差。

（注：增速减少，⽤何种模型）设y=f(x)具有指数形式x
b
ae
y =（a>0,b<0）。

对此式两边取对数，得x
b
a y 1
ln ln +=。

记A=lna ，B=b ，
程序：
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 2.3758]; polyfit(t,z,1) 结果：
ans = -1.1107 2.4578
由此可得 A=2.4578，B=-1.1107，6791.11==A
e a ，b=B=-1.1107
⽅程即为x
e
y 1107.16791.11-=
计算均⽅差编程：
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 a=y(i); b=x(i);
c=c+(a-f(b))^2; end
averge=c/15 结果：averge =
0.0594
4．设
----------------=410100141010014101101410010141001014A ，
--=625250b ，b x =A 分析下列迭代法的收敛性，并求42
110-+≤-k
k x x 的近似解及相应的迭代次数。

（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);
运⾏结果为:
y =
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)
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 =
1.0000
2.0000
1.0000
2.0000
1.0000
n =
15
（3） SOR 迭代（95.0,
95.1,334.1=ω）。

程序：
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 =
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
5．⽤逆幂迭代法求
=111123136A 最接近于11的特征值和特征向量，准确到3
10-。

程序：
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 %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
特征向量： my =
0.3845 -1.0000 0.7306
6．⽤经典R-K ⽅法求解初值问题
（1）-+-='++-='x x y y y x y y y sin 2cos 22sin 22212211
，]10,0[∈x ，
==3)0(2
)0(2
1y y ；程序：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 =
0 1 2 3 4 5 6 7 8 9 10 y =
2.0000
3.0000 1.5775 1.2758 1.1802 -0.1457 0.2406 -0.8903 -0.7202 -0.6170 -0.9454 0.2971 -0.2745 0.9652 0.6589 0.7557 0.9901 -0.1449 0.4124 -0.9109 -0.5440 -0.8389 （2）??
-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212211
，]10,0[∈x ，
==3
)0(2
)0(21y y 。

和精确解+=+=--x
e x y x
e x y x
x cos 2)(sin 2)(21⽐较，分析结论。

程序：function ydot=lorenzeq1(x,y)
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]);
fprintf(' x y(1) y1 y(2) y2\n') for j=1:length(y)
fprintf('%4d %.4f %.4f %.4f %.4f\n',x(j),y(j,1),y1(j),y(j,2),y2(j)) end
运⾏结果为:
x y(1) y1 y(2) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5772 1.5772 1.2759 1.2761 2 1.1800 1.1800 -0.1455 -0.1455 3 0.2407 0.2407 -0.8904 -0.8904 4 -0.7202 -0.7202 -0.6169 -0.6170 5 -0.9454 -0.9454 0.2972 0.2971 6 -0.2745 -0.2745 0.9648 0.9651 7 0.6588 0.6588 0.7554 0.7557 8 0.9900 0.9900 -0.1448 -0.1448 9 0.4124 0.4124 -0.9106 -0.9109 10 -0.5439 -0.5439 -
0.8389 -0.8390 结论：R-K ⽅法求解的结果精度较⾼。

7．⽤有限差分法求解边值问题（h=0.1）：
==-=+-''1
)1()1(0
)1(2y y y x y . 程序为： h=0.1;
n=(1-(-1))/h+1; x(1)=-1;x(n-1)=1; y(1)=1;y(n-1)=1; for i=1:n-1
x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)^2); B(i)=2/(h^2)+q(i); end
for i=1:n-2
C(i)=-1/(h^2); end
H=diag(B)+diag(C,1)+diag(C,-1); g(1)=0+1/(h^2); g(n-1)=0+1/(h^2); for i=2:n-2 g(i)=0; 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
0.9072
8．拟合形如f（x）≈（a+bx）/(1+cx)的函数的⼀种快速⽅法是将最⼩⼆乘法⽤于下列问题：f（x）(1+cx)≈（a+bx），试⽤这⼀⽅法拟合表4-4给出的中国⼈⼝数据。

解：把f（x）(1+cx)≈（a+bx）变成f（x）≈a+bx-cx f（x）则近似看成基函数是1，x，-x*f（x）⽽数据是（x i，f（x i））的最⼩⼆乘拟合问题，程序如下：
x=[1953 1964 1982 1900 2000]';
y=[5.82 6.95 10.08 11.34 12.66]';
A=[ones(5,1) x -x.*y];
Z=A\y;
a=Z(1)
b=Z(2)
c=Z(3)
结果：
a =
11.5250
b =
-0.0059
c =
-5.0979e-004。