数值计算上机报告.docx
《数值计算上机报告.docx》由会员分享，可在线阅读，更多相关《数值计算上机报告.docx（33页珍藏版）》请在冰点文库上搜索。
数值计算上机报告
上海电力学院
现代数值计算上机实验
报告
院 系：
能源与机械工程学院
专业年级：
2014级动力机械140101班
学生姓名：
张焱儒
学号：
14101058
指导老师：
黄建雄
2015年1月4日
1.设
，
（1）由递推公式
，从
的几个近似值出发，计算
；
解：
I0=
=0.1823
计算I20编辑matlab命令如下：
I=0.1823
forn=1:
1:
20,
I=-5*I+1/n;
fprintf（'%.1d%.4f\n',n,I）;
end
结果：
（2）粗糙估计
，用
，计算
；
解：
I20=
使用复合中点公式进行积分，相应的matlab程序如下：
I=0;
forh=0:
0.001:
1,
m=h+0.0005;
I=I+0.001*m^20/（5+m）;
fprintf（'%.1d%.4f\n',m,I）;
end
disp（I）;
fork=1:
20,
n=21-k;
I=0.2*（1/n-I）;
fprintf（'%.1d%.4f\n',n,I）;
end
disp（I）
结果：
程序结束时输出两个I值，第一个表示I20，第二个表示I0；
分别为I20=0.0082
I0=0.1823
（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。
首先分析两种递推式的误差：
设第一递推式中开始时的误差为
，递推过程的舍入误差不计。
并记
，
则有
所此递推式不可靠。
而在第二种递推式中
，误差在缩小，所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中，考虑舍入误差是否得到有效控制或是舍入误差的增长是否影响产生可靠的结果，即算法是否数值稳定。
2.求方程
的近似根，要求
，并比较计算量。
（1）在[0，1]上用二分法；
Matlab程序如下：
a=0;
b=1;
c=b-a;
n=0
whilec>0.0005,
x=（a+b）/2;
f=exp（x）+10*x-2;
iff>0,
b=x;
c=b-a;
elseiff<0,
a=x;
c=b-a;
else
x=x;
c=0;
end
n=n+1;
fprintf（'%.1d%.4f%.4f\n',n,x,c）;
end
结果如下：
解得到;x=0.0903
（2）取初值
，并用迭代
；
采用matlab进行迭代的程序如下：
x=0;
c=1;
n=0;
whilec>0.0005,
m=x;
m=（2-exp（m））/10;
c=abs（m-x）;
x=m;
n=n+1;
fprintf（'%.1d%.4f%.4f\n',n,x,c）;
end
结果：
解得x=0.0905
（3）加速迭代的结果；
采用matlab进行迭代的程序如下：
x=0;
n=0;
a=0;
b=1;
whileabs（a-b）>0.0005,
n=n+1;
a=x;
y=（2-exp（x））/10;
z=（2-exp（y））/10;
x=x-（y-x）^2/（z-2*y+x）;
b=x;
fprintf（'%.1d%.4f%.4f\n',n,x,abs（a-b））;
end
结果如下：
（4）取初值
，并用牛顿迭代法；
Matlab程序如下：
x=0;
a=1;
n=0;
whileabs（a）>0.0005,
n=n+1;
a=（exp（x）+10*x-2）/（exp（x）+10）;
x=x-a;
fprintf（'%.1d%.4f%.4f\n',n,x,abs（a））;
end
运行结果：
（5）分析绝对误差：
迭代次数
二分法
代数式迭代
加速迭代
牛顿迭代
X（k）
Erroe
X（k）
Erroe
X（k）
Erroe
X（k）
Erroe
1
0.5000
0.5000
0.1000
0.1000
0.0905
0.0905
0.0909
0.0909
2
0.2500
0.2500
0.0895
0.0105
0.0905
0．0000
0.0905
0.0004
3
0.1250
0.1250
0.0906
0.0012
4
0.0625
0.0625
0.0905
0.0001
5
0.0938
0.0313
6
0.0781
0.0156
7
0.0859
0.0078
8
0.0898
0.0039
9
0.0918
0.0020
10
0.0908
0.0010
11
0.0903
0.0005
我们可以看到，在运算要求到同一精度的情况下，采用
（1）的二分法运算了11次，采用
（2）的方法运算了4次，采用（3）的加速迭代法运算了2次，采用（4）的牛顿迭代法也需运算2次。
也就是说牛顿的迭代的收敛速度与加速迭代速度都是超线性收敛的，而简单迭代法是线性收敛的。
而二分法收敛速度较慢,所以在工程上不经常使用。
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
试从中找出使用次数和容积之间的关系，计算均方差。
（注：
增速减少，用何种模型）
设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
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
解：
将使用次数x与体积y的关系用matlab采用如下程序绘制在二维坐标系:
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
plot（x,y,'b*-'）;
结果如下：
由数据点分布图可知，拟合曲线y=f（x）随着x的增加而上升，但上升速度由快到慢，当x趋于无穷大时，y趋于某个常数，故曲线有一水平渐进线。
根据上述特征很容易想到用Logistic模型来拟合该曲线。
设y=f（x）的形式为y=aeb/x（a>0,b<0），两边取对数，得lny=lna+b/x。
记A=lna，B=b，并引入新变量m=lny，n=1/x。
引入新变量后的数据表如下：
x
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
n=1/x
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
m=lny
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
Matlab拟合程序如下：
n=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625];
m=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758];
polyfit（n,m,1）
运行的结果：
由此可得A=2.4578B=-1.1107
a=eA=11.6791b=B=-1.1107
由此得到使用次数与容积的函数为
将统计表和函数用matlab绘制在同一坐标系内程序如下：
x1=[2345678910111213141516];
y1=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
x=2:
1:
16;
y=11.679*exp（-1.1107*x.^（-1））;
holdon;
plot（x,y,'ro-'）;
plot（x1,y1,'b*-'）;
结果如图：
计算均方差s，matlab程序如下：
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
s=0;
forn=2:
1:
16,
a=abs（11.679*exp（-1.1107*n.^（-1））-y（n-1））;
s=s+a.^2;
end
s=（s/15）.^（1/2）;
disp（s）;
运算结果均方差S=0.2438
小结：
根据已给的条件计算函数是十分困难的，但通过对离散点的分析及变化规律找出其中的规律，并通过计算来得到实际的函数是十分有用的方法。
本题就是这样做的一个典型，在n=1/x和m=lny的基础上找到了它们之间的关系并通过这种关系来拟合原函数，并最终验证计算结果。
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=
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]';
seidel（a,b,x0）;
运行结果为：
y=
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
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=
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
小结：
三种迭代法都是可以计算出结果，从收敛的角度来看雅可比迭代法，高斯--赛德尔迭代法迭代法都比较好,而超松弛迭代法（SOR）的迭代效果由于w的取值而有很大的变化。
迭代收敛的条件是迭代的迭代矩阵谱半径小于1。
5．用逆幂迭代法求
最接近于11的特征值和特征向量，准确到
。
解：
建立eigIPower函数的m文件。
function[t,y]=eigIPower（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
t=1/norm（v,inf）+p
y=u
命令窗口中输入：
>>A=[631;321;111];eigIPower（A,11,0.001）;
结果为：
特征值：
t=
11.0919
特征向量：
y=
0.3845
-1.0000
0.7306
小结：
逆幂迭代用来寻找最小的特征值，从以上程序可以看出迭代7次即达到了精度要求，证明了逆幂迭代是一种比较好的算法。
6．用经典R-K方法求解初值问题
（1）
，
，
；
解:
采用经典R-K公式计算的MATLAB程序如下：
y1=2;
y2=3;
forh=0:
0.01:
10,
k1=0.01*（-2*y1+y2+2*sin（h））;
l1=0.01*（y1-2*y2+2*cos（h）-2*sin（h））;
k2=0.01*（-2*（y1+0.5*k1）+（y2+0.5*l1）+2*sin（h+0.005））;
l2=0.01*（（y1+0.5*k1）-2*（y2+0.5*l1）+2*cos（h+0.005）-2*sin（h+0.005））;
k3=0.01*（-2*（y1+0.5*k2）+（y2+0.5*l2）+2*sin（h+0.005））;
l3=0.01*（（y1+0.5*k2）-2*（y2+0.5*l2）+2*cos（h+0.005）-2*sin（h+0.005））;
k4=0.01*（-2*（y1+k3）+（y2+l3）+2*sin（h+0.01））;
l4=0.01*（（y1+k3）-2*（y2+l3）+2*cos（h+0.01）-2*sin（h+0.01））;
y1=y1+1/6*（k1+2*k2+2*k3+k4）;
y2=y2+1/6*（l1+2*l2+2*l3+l4）;
ifh==fix（h）;
fprintf（'%.1d%.4f%.4f\n',h,y1,y2）;
else
end
end
结果如下所示：
（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
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-0.8389-0.8390
小结：
经典R-K法是一种四阶的算法，具有很高的精度，这一点通过比较
（1）中的数值解和解析解就可以看得出来，在计算
（2）题时我们可以从数值稳定的相关理论判断此题运用经典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=
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
《方法二》
用有限差分法matlab程序如下：
h=0.1;
n=2/0.1-1;
g
（1）=1/（h.^2）;
g（n）=1/（h^2）;
fori=2:
1:
18,
g（i）=0;
end
g=[g
（1）;g
（2）;g（3）;g（4）;g（5）;g（6）;g（7）;g（8）;g（9）;g（10）;g（11）;g（12）;g（13）;g（14）;g（15）;g（16）;g（17）;g（18）;g（19）];
disp（g）;
fori=1:
1:
19,
forj=1:
1:
19,
ifi==1,
H（1,1）=2/（h.^2）+（1+（-1+0.1*i）.^2）;
H（1,2）=-1/（h.^2）;
elseifi==19,
H（19,18）=-1/（h.^2）;
H（19,19）=2/（h.^2）+（1+（-1+0.1*i）.^2）;
else
ifj==i,
H（i,j）=2/（h.^2）+（1+（-1+0.1*i）.^2）;
elseifj==i-1,
H（i,i-1）=-1/（h.^2）;
elseifj==i+1;
H（i,i+1）=-1/（h.^2）;
else
H（i,j）=0;
end
end
end
end
disp（H）;
y=H\g;
fori=1:
1:
19,
fprintf（'%.4f%.4f\n',-1+0.1*i,y（i,1））
end
运算结果为：
H矩阵为：
Y在各点的近似值为：
X
Y
小结：
有限差分法是将微分方程离散化为差分方程进行求解，而求解差分方程就是解一个三对角矩阵线性方程组，从以上程序可以看出编程很容易实现用追赶法求解三对角矩阵线性方程组，且具有较高的精度。
8．拟合形如f（x）≈（a+bx）/（1+cx）的函数的一种快速方法是将最小二乘法用于下列问题：
f（x）（1+cx）≈（a+bx），试用这一方法拟合表4-4给出的中国人口数据。
（数值实验四）
表4-4
次序年份人口（亿）
a）19535.82
b）19646.95
c）198210.08
d）199011.34
e）200012.66
升级会员 