数值计算方法上机实习题答案doc.docx
《数值计算方法上机实习题答案doc.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实习题答案doc.docx（21页珍藏版）》请在冰点文库上搜索。
数值计算方法上机实习题答案doc
数值计算方法上机作业热能工程
1．设In
1
xn
dx，
05
x
（1）
由递推公式In
5In
1
1，从I0的几个近似值出发，计算
I20；
n
解：
易得：
I0
ln6-ln5=0.1823,
程序为：
I=0.182;
forn=1:
20
I=（-5）*I+1/n;
end
I
输出结果为：
I20=-3.0666e+010
（2）
粗糙估计I20，用In1
1In1
1
，计算I0；
5
5n
0.0079
1
x
20
1x
20
0.0095
因为
0
dxI
20
0
dx
6
5
所以取I20
1（0.0079
0.0095）
0.0087
2
程序为：
I=0.0087;
forn=1:
20
I=（-1/5）*I+1/（5*n）;
end
I
I0=0.0083
（3）
分析结果的可靠性及产生此现象的原因
（重点分析原因）。
首先分析两种递推式的误差；
设第一递推式中开始时的误差为
E0
I0I0，递推过程的舍入误差不计。
并记En
In
In，则有En
5En1
（5）nE0。
因为E20
（5）20E0
I20，所此递推式不可靠。
而在第二种递推式中
E0
1E1
（1）nEn，误差在缩小，
5
5
所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中，
考虑误差是否得到控制，
即算法是否数值稳定。
2．求方程ex10x20的近似根，要求xk1xk5104，并比较计算量。
（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）取初值x0
0，并用迭代xk1
2ex
；
10
程序：
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）^2/（z-2*y+x）;
b=x;
end
x
结果：
x=
0.0995
（4）取初值x00，并用牛顿迭代法；
程序：
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）+10*x-2=0'）
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
试从中找出使用次数和容积之间的关系，计算均方差。
（注：
增速减少，用何种模型）
b
设y=f（x）具有指数形式
y
aex（a>0,b<0）。
对此式两边取对数，得lnylnab1
。
记A=lna，B=b，
x
并引入新变量
z=lny，t=1/x。
引入新变量后的数据表如下
x
2
3
4
5
6
7
8
9
t=1/x
0.5000
0.3333
0.2500
0.2000
0.1667
0.1429
0.1250
0.1111
z=lny
1.8594
2.1041
2.2597
2.2513
2.2721
2.3026
2.2956
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
程序：
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=-1.1107
2.4578
由此可得A=2.4578，B=-1.1107，aeA11.6791，b=B=-1.1107
1.1107
方程即为y11.6791ex
数值计算方法上机作业热能工程
计算均方差编程：
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=
0.0594
4
1
0
1
0
0
0
1
4
1
0
1
0
5
0
1
4
1
0
1
2
4．设A
0
1
4
1
，b
，Axb
1
0
5
0
1
0
1
4
1
2
0
0
1
0
1
4
6
分析下列迭代法的收敛性，并求
xk1
xk2
104的近似解及相应的迭代次数。
（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=
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
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=
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）>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=
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．用逆幂迭代法求A321最接近于11的特征值和特征向量，准确到103。
111
程序：
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=
11.0919
特征向量：
my=
0.3845
-1.0000
0.7306
6．用经典R-K方法求解初值问题
y1
2y1
y2
2sinx
，x
y1（0）
2
（1）
y1
2y2
2cosx2sinx
[0,10]，
；
y2
y2（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
y1
2y1
y2
2sinx
[0,10]，
y1（0）
2
（2）
998y1
，x
y2（0）
。
y2
999y2999cosx999sinx
3
和精确解
y1（x）
2ex
sinx
2ex
比较，分析结论。
y2（x）
cosx
程序：
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（'
x
y
（1）
y1
y
（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
运行结果为:
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）：
数值计算方法上机作业热能工程
y（1x2）y0
.
y
（1）y
（1）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=
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）的函数的一种快速方法是将最小二乘法用于下列问题：
≈（a+bx），试用这一方法拟合表4-4给出的中国人口数据。
f（x）（1+cx）
表4-4
次序
年份
人口（亿）
a）
1953
5.82
b）
1964
6.95
c）
1982
10.08
d）
1900
11.34
e）
2000
12.66
解：
把f（x）（1+cx）
≈（a+bx）变成f（x）≈a+bx-cx
f（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=
11.5250
b=
-0.0059
c=
-5.0979e-004
升级会员 