数值计算方法上机实习题答案.docx
《数值计算方法上机实习题答案.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实习题答案.docx（14页珍藏版）》请在冰豆网上搜索。
数值计算方法上机实习题答案
1．设
，
（1）由递推公式
，从
的几个近似值出发，计算
；
解：
易得：
ln6-ln5=0.1823,
程序为：
I=0.182;
forn=1:
20
I=（-5）*I+1/n;
end
I
输出结果为：
=-3.0666e+010
（2）粗糙估计
，用
，计算
；
因为
所以取
程序为：
I=0.0087;
forn=1:
20
I=（-1/5）*I+1/（5*n）;
end
I
（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。
首先分析两种递推式的误差；设第一递推式中开场时的误差为
，递推过程的舍入误差不计。
并记
，那么有
。
因为
，所此递推式不可靠。
而在第二种递推式中
，误差在缩小，所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2．求方程
的近似根，要求
，并比拟计算量。
（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=
（2）取初值
，并用迭代
；
程序：
x=0;
a=1;
whileabs（x-a）>5*1e-4
a=x;
x=（2-exp（x））/10;
end
x
结果：
x=
（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）^2/（z-2*y+x）;
b=x;
end
x
结果：
x=
（4）取初值
，并用牛顿迭代法；
程序：
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=
（5）分析绝对误差。
solve（'exp（x）+10*x-2=0'）
3．钢水包使用次数多以后，钢包的容积增大，数据如下：
x
2
3
4
5
6
7
8
9
y
10
10
11
12
13
14
15
16
试从中找出使用次数和容积之间的关系，计算均方差。
〔注：
增速减少，用何种模型〕
设y=f（x）具有指数形式
〔a>0,b<0〕。
对此式两边取对数，得
。
记A=lna，B=b，并引入新变量z=lny，t=1/x。
引入新变量后的数据表如下
x
2
3
4
5
6
7
8
9
t=1/x
z=lny
10
11
12
13
14
15
16
程序：
t=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625];
z=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758];
polyfit（t,z,1）
结果：
ans=
由此可得A=2.4578，B=-1.1107，
方程即为
计算均方差编程：
x=[2:
16];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
f（x）=11.6791*exp（-1.1107./x）;
c=0;
fori=1:
15
a=y（i）;
b=x（i）;
c=c+（a-f（b））^2;
end
averge=c/15
结果：
averge=
4．设
，
，
分析以下迭代法的收敛性，并求
的近似解及相应的迭代次数。
（1）JACOBI迭代；
程序：
functiony=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;
whilenorm（y-x0）>1e-4
x0=y;
y=B*x0+f;n=n+1;
end
y
n
以文件名jacobi.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,x0）;
运行结果为:
y=
n=
28
（2）GAUSS-SEIDEL迭代；
程序：
functiony=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;
whilenorm（y-x0）>10^（-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,x0）;
运行结果为：
y=
n=
15
（3）SOR迭代〔
〕。
程序：
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）>10^（-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=
n=
13
y=
n=
241
y=
n=
17
5．用逆幂迭代法求
最接近于11的特征值和特征向量，准确到
。
程序：
function[mt,my]=maxtr（A,p,ep）
n=length（A）;
B=A-p*eye（n）;
v0=ones（n,1）;
k=1;
v=B*v0;
whileabs（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=
特征向量：
my=
6．用经典R-K方法求解初值问题
〔1〕
，
，
；
程序：
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〕
，
，
。
和准确解
比拟，分析结论。
程序：
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）];
以文件名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（'xy
（1）y1y
（2）y2\n'）
forj=1:
length（y）
fprintf（'%4d%.4f%.4f%.4f%.4f\n',x（j）,y（j,1）,y1（j）,y（j,2）,y2（j））
end
运行结果为:
xy
（1）y1y
（2）y2
结论：
R-K方法求解的结果精度较高。
7．用有限差分法求解边值问题〔h=0.1〕：
.
程序为：
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）^2）;
B（i）=2/（h^2）+q（i）;
end
fori=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）;
fori=2:
n-2
g（i）=0;
end
y=H\g'
运行结果为：
y=
8．拟合形如f〔x〕≈〔a+bx〕/（1+cx）的函数的一种快速方法是将最小二乘法用于以下问题：
f〔x〕（1+cx）≈〔a+bx〕，试用这一方法拟合表4-4给出的中国人口数据。
表4-4
次序年份人口〔亿〕
a）
b）
c）
d）
e）
解：
把f〔x〕（1+cx）≈〔a+bx〕变成f〔x〕≈a+bx-cxf〔x〕那么近似看成基函数是1，x，-x*f〔x〕而数据是〔xi，f〔xi〕〕的最小二乘拟合问题，程序如下：
x=[19531964198219002000]';
y=[5.826.9510.0811.3412.66]';
A=[ones（5,1）x-x.*y];
Z=A\y;
a=Z
（1）
b=Z
（2）
c=Z（3）
结果：
a=
b=
c=
升级会员 