计算方法课程设计Word文档格式.doc
《计算方法课程设计Word文档格式.doc》由会员分享，可在线阅读，更多相关《计算方法课程设计Word文档格式.doc（22页珍藏版）》请在冰点文库上搜索。
（1）考察单精度计算结果（与真解对比）；
（2）若计算结果与真解相差很大，分析其原因，提出新的算法（如先求再根据根与系数关系求）以改进计算结果。
实验步骤：
方程
（1）：
根据求根公式，写出程序：
formatlong
a=1;
b=3;
c=-2;
x1=（（-1）*b+sqrt（b^2-4*a*c））/2*a
x2=（（-1）*b-sqrt（b^2-4*a*c））/2*a
运行结果：
x1=0.561552812808830
x2=-3.561552812808830
然后
由符号解求的该方程的真值程序为：
symsx;
y=x^2+3*x-2;
s=solve（y,x）;
vpa（s）
运行结果为：
X1=0.56155281280883027491070492798704
X2=-3.561552812808830274910704927987
方程
（2）：
b=-10^10;
c=1;
x1=1.000000000000000e+010
x2=0
y=x^2-10^10*x+1;
X1=10000000000.0
X2=0.000000000116415321827
由此可知，对于方程
（1），使用求根公式求得的结果等于精确值，故求根公式法可靠。
而对于方程
（2），计算值与真值相差很大，故算法不可靠。
改进算法，对于方程
（2）:
先用迭代法求x1，再用，求x2
程序为：
symsk
a=[];
fori=1:
100
a
（1）=4;
a（i+1）=（10^10*a（i）-1）^（1/2）;
end
x1=a（100）
x2=1/（x1）
x2=1.000000000000000e-010
实验结论：
对于方程
（1），两种方法在精确到小数点后15位时相同，说明两种算法的结果都是精确的。
对于方程
（2），两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算法是不精确的。
原因：
方程
（2）用求根公式计算时，公式中，b是大数，出现了大数吃掉小数的误差，也出现了两个相近的数相减的误差，所以出现x2=0这样大的误差。
改进的结果会比较准确。
2、求解非线性方程
实验2.1Gauss消去法的数值稳定性实验
观察和理解高斯消元过程中出现小主元即很小时，引起方程组解的数值不稳定性．
求解线性方程组其中
（1），
（2）
（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的
（2）用高斯列主元消去法求得L和U及解向量
（3）用不选主元的高斯消去法求得L和U及解向量
（4）观察小主元并分析对计算结果的影响
（1）计算矩阵的条件数程序：
矩阵A1：
A1=[0.3*10^（-15）59.1431;
5.291-6.130-12;
11.2952;
1211];
cond（A1,1）
cond（A1,2）
cond（A1,inf）
ans=136.294470872045e+000
ans=68.4295577125341e+000
ans=84.3114602051800e+000
由矩阵条件数判断出矩阵A1是病态矩阵。
矩阵A2：
A2=[10-701;
-32.0999999999962;
5-15-1;
0102];
ans=19.2831683168042e+000
ans=8.99393809015365e+000
ans=18.3564356435280e+000
由矩阵条件数判断出矩阵A2是病态矩阵。
（2）高斯列主元消去法程序：
方程组
（1）：
b1=[59.17;
46.78;
1;
2];
n=4;
fork=1:
n-1
a=max（abs（A1（k:
n,k）））;
[p,k]=find（A1==a）;
B=A1（k,:
）;
c=b1（k）;
A1（k,:
）=A1（p,:
b1（k）=b1（p）;
A1（p,:
）=B;
b1（p）=c;
ifA1（k,k）~=0
A1（k+1:
n,k）=A1（k+1:
n,k）/A1（k,k）;
n,k+1:
n）=A1（k+1:
n）-A1（k+1:
n,k）*A1（k,k+1:
n）;
else
break
end
L1=tril（A1,0）;
n
L1（i,i）=1;
L=L1
U=triu（A1,0）
forj=1:
b1（j）=b1（j）/L（j,j）;
b1（j+1:
n）=b1（j+1:
n）-b1（j）*L（j+1:
n,j）;
b1（n）=b1（n）/L（n,n）;
forj=n:
-1:
2
b1（j）=b1（j）/U（j,j）;
b1（1:
j-1）=b1（1:
j-1）-b1（j）*U（1:
j-1,j）;
b1
（1）=b1
（1）/U（1,1）;
x1=b1
方程组
（2）：
b2=[8;
5.9000000000001;
5;
1];
a=max（abs（A2（k:
[p,k]=find（A2==a）;
B=A2（k,:
c=b2（k）;
A2（k,:
）=A2（p,:
b2（k）=b2（p）;
A2（p,:
b2（p）=c;
ifA2（k,k）~=0
A2（k+1:
n,k）=A2（k+1:
n,k）/A2（k,k）;
n）=A2（k+1:
n）-A2（k+1:
n,k）*A2（k,k+1:
L1=tril（A2,0）;
U=triu（A2,0）
b2（j）=b2（j）/L（j,j）;
b2（j+1:
n）=b2（j+1:
n）-b2（j）*L（j+1:
b2（n）=b2（n）/L（n,n）;
b2（j）=b2（j）/U（j,j）;
b2（1:
j-1）=b2（1:
j-1）-b2（j）*U（1:
b2
（1）=b2
（1）/U（1,1）;
x2=b2
（3）不选主元的高斯消去法程序：
方程组
（1）：
clear
A1=[0.3e-1559.1431;
A1（k+1:
L1（i,i）=1;
b1（j）=b1（j）/L（j,j）;
b1（j+1:
b1（j）=b1（j）/U（j,j）;
b1（1:
运行结果:
A2（k+1:
（4）分析小元对计算结果的影响
通过观察计算结果，分析可知，小元对计算结果的影响很大，小元的存在会使得到的计算结果有很大的误差。
3、插值
4、数值积分
实验2.4：
复化求和公式计算定积分
计算下列各式右端定积分的近似值
（1）分别用复化梯形公式、复化Simpson公式计算，要求绝对误差限为；
（2）将计算结果与精确解作比较，并比较各种算法的计算量。
（1）复化梯形公式和复化Simpson程序：
式
（1）：
a=2;
f=1/（x^2-1）;
n=8;
h=（b-a）/n;
s=0;
w=0;
x=a+h*i;
s=s+subs（f,x）*2;
w=（s+subs（f,2）+subs（f,3））*h/2;
f4=2*w
f1=0;
f2=0;
ifrem（i,2）~=0
x1=a+i*h;
f1=f1+subs（f,x1）;
elserem（i,2）==0
x2=a+i*h;
f2=f2+subs（f,x2）;
f3=（h/3）*（subs（f,2）+subs（f,3）+4*f1+2*f2）*2
f4=0.4064
f3=0.4055
式
（2）:
a=0;
b=1;
f=1/（x^2+1）;
n=100;
f4=4*w
f3=（h/3）*（subs（f,2）+subs（f,3）+4*f1+2*f2）*4
f4=3.1176
f3=3.1256
式（3）：
f=3^x;
f4=w
f3=（h/3）*（subs（f,2）+subs（f,3）+4*f1+2*f2）
f4=1.9805
f3=1.9271
式（4）：
b=2;
f=x*exp（x）;
f4=7.6769
f3=7.5809
（2）求精确解程序：
y1=log
（2）
y2=pi
y3=2/log（3）
y4=exp
（2）
y1=0.6931
y2=3.1416
y3=1.8205
y4=7.3891
（3）分析
使用复化梯形公式和复化Simpson公式计算所得的结果与真实值比较有很大的误差。
两种公式算法中，复化梯形公式计算量比复化Simpson公式的计算量少，但是两种公式算法精确度都不是很高，有待使用其他精确度高的算法替代。
二、提高技能训练
1、kepler方程的计算
在人造卫星轨道理论中对应的二体问题是可积系统，其中有6个积分常数为，其中表示轨道半长径，表示轨道偏心率，表示轨道倾角，是升交点赤经，ω是近地点辐角，是平近点角。
其中第六个积分常数，通常可以用其他的一些量替换，如过近地点时刻，真近点角度f与偏近点角度。
这几个是相互等价的关系，我们这里要讲的就是平近点角M与偏近点角的一个关系，有如下方程
这就是著名的Kelper方程，是在人造卫星运动理论或者行星运动理论中最基本的方程之一。
因为如果我们已经知道卫星轨道量去计算卫星的星历时，第6个根数通常是使用平近点角。
这时候需要把平近点角化为偏近点角，也就是需要计算上述的Kepler方程。
具体理论这里不再阐述，而Kepler方程本身是我们这里感兴趣的。
可采用Newton迭代法计算Kepler方程，在实际工程中，一般也是这样做的，因为Newton迭代法收敛速度快而且精度比较高。
也可用其它算法计算。
实验目的与要求
（1）了解Kepler方程；
（2）能够用牛顿迭代方法计算Kepler方程；
（3）通过调节轨道偏心率e查看运算收敛情况。
实验内容及数据来源
编写解kepler方程的方法，给定偏心率为0.01、平近点角为32度时计算出偏近点角。
（1）对Kepler方程的了解：
Kepler方程又称作开普勒方程，是二体问题运动方程的一个积分。
它反映天体在其轨道上的位置与时间t的函数关系。
对于椭圆轨道，开普勒方程可以表示为E－esinE=M，式中E为偏近点角，M为平近点角，都是从椭圆轨道的近地点开始起算，沿逆时针方向为正，E和M都是确定天体在椭圆轨道上的运动和位置的基本量。
（2）牛顿迭代法计算开普勒方程程序：
N=100000;
x0=0;
y=32-x0+0.01*sin（x0）;
E=1.0e-6;
k=1;
while（norm（y-x0）>
=E&
&
k<
N）
x0=y;
f=diff（y）;
f1=diff（f）;
if（f1~=0）
y=x0-f/f1;
k=k+1;
fprintf（'
迭代的次数为k=%d\n'
k）;
\n'
方程的根为x*=%.7f\n'
y）;
迭代的次数为k=2
方程的根为x*=32.0000000
（3）调节偏心率查看收敛情况：
调节不用的偏心率的运行结果都趋于32，由此判断出该方程为收敛方程。
调节不同的偏心率程序运行出来的结果没有很大的偏差，所以该方程是收敛的。
2
（1）、一维插值问题的应用
（1）（机翼加工问题）书籍机翼轮廓上的数据如下表所示
x/m
3
5
7
9
11
12
13
14
15
y/m
1.2
1.7
2.0
2.1
1.8
1.0
1.6
加工时需要x每改变0.1m时y值，并画出相应的轮廓曲线。
试验程序：
x0=[035791112131415];
y0=[01.21.72.02.12.01.81.21.01.6];
x=0:
0.1:
15;
y1=interp1（x0,y0,x）;
%分段线性插值
y2=interp1（x0,y0,x,'
spline'
%三次样条插值
subplot（2,1,1）
plot（x0,y0,'
k+'
x,y1,'
r'
grid
title（'
piecewiselinear'
）%图片命名
subplot（2,1,2）
x,y2,'
（2）、二维插值问题的应用
（山区地貌图）在某山区（平面区域内，单位:
m）测得一些地点的高度（单位:
m）如下表所示。
2400
1430
1450
1470
1320
1280
1200
1080
940
2000
1480
1500
1550
1510
1300
1600
1460
1370
1100
1380
800
1270
1350
1150
400
1230
1390
1400
900
1060
1180
1420
700
y/x
2800
试作出该山区的地貌图和等值线图。
400:
2800;
y=0:
2400;
h=[118013201450142014001300700900;
1230139015001500140090011001060;
12701500120011001350145012001150;
13701500120011001550160015501380;
14601500155016001550160016001600;
14501480150015501510143013001200;
1430145014701320128012001080940];
xi=0:
50:
xi=xi'
;
yi=0:
z=interp2（x,y,h,xi,yi,'
cubic'
figure
（1）
mesh（xi,yi,z）
figure
（2）
contour（x,y,h）
地貌图及等值线图如下：
三、本课程设计的心得体会（500字
升级会员 