计算方法课程设计报告.docx
《计算方法课程设计报告.docx》由会员分享，可在线阅读，更多相关《计算方法课程设计报告.docx（27页珍藏版）》请在冰点文库上搜索。
计算方法课程设计报告
福建工程学院
数理学院
2014级信息与计算科学
计算方法
课程设计
****************
学号：
**********
班级：
1402
成绩：
实验要求
1．应用自己熟悉的算法语言编写程序，使之尽可能具有通用性。
2．上机前充分准备，复习有关算法，写出计算步骤，反复检查，调试程序。
（注：
在练习本上写，不上交）
3．完成计算后写出实验报告，内容包括:
算法步骤叙述，变量说明，程序清单，输出计算结果，结构分析和小结等。
（注：
具体题目具体分析，并不是所有的题目的实验报告都包含上述内容！
）
4．独立完成，如有雷同，一律判为零分！
5．上机期间不允许做其他任何与课程设计无关的事情，否则被发现一次扣10分，被发现三次判为不及格！
非特殊情况，不能请假。
旷课3个半天及以上者，直接判为不及格。
一、基本技能训练
1、误差分析
实验1.3求一元二次方程的根
实验目的：
研究误差传播的原因与解决对策。
问题提出：
求解一元二次方程
实验内容：
一元二次方程的求根公式为
用求根公式求解下面两个方程：
实验要求：
（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）：
formatlong
a=1;b=-10^10;c=1;
x1=（（-1）*b+sqrt（b^2-4*a*c））/2*a
x2=（（-1）*b-sqrt（b^2-4*a*c））/2*a
运行结果为：
x1=1.000000000000000e+010
x2=0
然后
由符号解求的该方程的真值程序为：
symsx;
y=x^2-10^10*x+1;
s=solve（y,x）;
vpa（s）
运行结果：
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）
运行结果为：
x1=1.000000000000000e+010
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];
cond（A1,1）
cond（A1,2）
cond（A1,inf）
运行结果：
ans=19.2831683168042e+000
ans=8.99393809015365e+000
ans=18.3564356435280e+000
由矩阵条件数判断出矩阵A2是病态矩阵。
（2）高斯列主元消去法程序：
方程组
（1）：
A1=[0.3*10^（-15）59.1431;5.291-6.130-12;11.2952;1211];
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）;
A1（k+1:
n,k+1:
n）=A1（k+1:
n,k+1:
n）-A1（k+1:
n,k）*A1（k,k+1:
n）;
else
break
end
end
L1=tril（A1,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A1,0）
forj=1:
n-1
b1（j）=b1（j）/L（j,j）;
b1（j+1:
n）=b1（j+1:
n）-b1（j）*L（j+1:
n,j）;
end
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）;
end
b1
（1）=b1
（1）/U（1,1）;
x1=b1
运行结果：
方程组
（2）：
A2=[10-701;-32.0999999999962;5-15-1;0102];
b2=[8;5.9000000000001;5;1];
n=4;
fork=1:
n-1
a=max（abs（A2（k:
n,k）））;
[p,k]=find（A2==a）;
B=A2（k,:
）;c=b2（k）;
A2（k,:
）=A2（p,:
）;b2（k）=b2（p）;
A2（p,:
）=B;b2（p）=c;
ifA2（k,k）~=0
A2（k+1:
n,k）=A2（k+1:
n,k）/A2（k,k）;
A2（k+1:
n,k+1:
n）=A2（k+1:
n,k+1:
n）-A2（k+1:
n,k）*A2（k,k+1:
n）;
else
break
end
end
L1=tril（A2,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A2,0）
forj=1:
n-1
b2（j）=b2（j）/L（j,j）;
b2（j+1:
n）=b2（j+1:
n）-b2（j）*L（j+1:
n,j）;
end
b2（n）=b2（n）/L（n,n）;
forj=n:
-1:
2
b2（j）=b2（j）/U（j,j）;
b2（1:
j-1）=b2（1:
j-1）-b2（j）*U（1:
j-1,j）;
end
b2
（1）=b2
（1）/U（1,1）;
x2=b2
运行结果：
（3）不选主元的高斯消去法程序：
方程组
（1）：
clear
formatlong
A1=[0.3e-1559.1431;5.291-6.130-12;11.2952;1211];
b1=[59.17;46.78;1;2];
n=4;
fork=1:
n-1
A1（k+1:
n,k）=A1（k+1:
n,k）/A1（k,k）;
A1（k+1:
n,k+1:
n）=A1（k+1:
n,k+1:
n）-A1（k+1:
n,k）*A1（k,k+1:
n）;
end
L1=tril（A1,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A1,0）
forj=1:
n-1
b1（j）=b1（j）/L（j,j）;
b1（j+1:
n）=b1（j+1:
n）-b1（j）*L（j+1:
n,j）;
end
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）;
end
b1
（1）=b1
（1）/U（1,1）;
x1=b1
运行结果:
方程组
（2）：
clear
formatlong
A2=[10-701;-32.0999999999962;5-15-1;0102];
b2=[8;5.9000000000001;5;1];
n=4;
fork=1:
n-1
A2（k+1:
n,k）=A2（k+1:
n,k）/A2（k,k）;
A2（k+1:
n,k+1:
n）=A2（k+1:
n,k+1:
n）-A2（k+1:
n,k）*A2（k,k+1:
n）;
end
L1=tril（A2,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A2,0）
forj=1:
n-1
b2（j）=b2（j）/L（j,j）;
b2（j+1:
n）=b2（j+1:
n）-b2（j）*L（j+1:
n,j）;
end
b2（n）=b2（n）/L（n,n）;
forj=n:
-1:
2
b2（j）=b2（j）/U（j,j）;
b2（1:
j-1）=b2（1:
j-1）-b2（j）*U（1:
j-1,j）;
end
b2
（1）=b2
（1）/U（1,1）;
x2=b2
运行结果：
（4）分析小元对计算结果的影响
通过观察计算结果，分析可知，小元对计算结果的影响很大，小元的存在会使得到的计算结果有很大的误差。
3、插值
4、数值积分
实验2.4：
复化求和公式计算定积分
实验内容：
计算下列各式右端定积分的近似值
实验要求：
（1）分别用复化梯形公式、复化Simpson公式计算，要求绝对误差限为
；
（2）将计算结果与精确解作比较，并比较各种算法的计算量。
实验步骤：
（1）复化梯形公式和复化Simpson程序：
式
（1）：
a=2;b=3;
symsx;
f=1/（x^2-1）;
n=8;
h=（b-a）/n;
s=0;w=0;
fori=1:
n-1
x=a+h*i;
s=s+subs（f,x）*2;
end
w=（s+subs（f,2）+subs（f,3））*h/2;
f4=2*w
f1=0;f2=0;
fori=1:
n-1
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）;
end
end
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;
symsx;
f=1/（x^2+1）;
n=100;
h=（b-a）/n;
s=0;w=0;
fori=1:
n-1
x=a+h*i;
s=s+subs（f,x）*2;
end
w=（s+subs（f,2）+subs（f,3））*h/2;
f4=4*w
f1=0;f2=0;
fori=1:
n-1
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）;
end
end
f3=（h/3）*（subs（f,2）+subs（f,3）+4*f1+2*f2）*4
运行结果：
f4=3.1176
f3=3.1256
式（3）：
a=0;b=1;
symsx;
f=3^x;
n=100;
h=（b-a）/n;
s=0;w=0;
fori=1:
n-1
x=a+h*i;
s=s+subs（f,x）*2;
end
w=（s+subs（f,2）+subs（f,3））*h/2;
f4=w
f1=0;f2=0;
fori=1:
n-1
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）;
end
end
f3=（h/3）*（subs（f,2）+subs（f,3）+4*f1+2*f2）
运行结果：
f4=1.9805
f3=1.9271
式（4）：
a=1;b=2;
symsx;
f=x*exp（x）;
n=100;
h=（b-a）/n;
s=0;w=0;
fori=1:
n-1
x=a+h*i;
s=s+subs（f,x）*2;
end
w=（s+subs（f,2）+subs（f,3））*h/2;
f4=w
f1=0;f2=0;
fori=1:
n-1
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）;
end
end
f3=（h/3）*（subs（f,2）+subs（f,3）+4*f1+2*f2）
运行结果：
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&&kx0=y;
f=diff（y）;
f1=diff（f）;
if（f1~=0）
y=x0-f/f1;
end
k=k+1;
end
fprintf（'迭代的次数为k=%d\n',k）;
fprintf（'\n'）;
fprintf（'方程的根为x*=%.7f\n',y）;
运行结果：
迭代的次数为k=2
方程的根为x*=32.0000000
（3）调节偏心率查看收敛情况：
调节不用的偏心率的运行结果都趋于32，由此判断出该方程为收敛方程。
实验结论：
调节不同的偏心率程序运行出来的结果没有很大的偏差，所以该方程是收敛的。
2
（1）、一维插值问题的应用
（1）（机翼加工问题）书籍机翼轮廓上的数据如下表所示
x/m
0
3
5
7
9
11
12
13
14
15
y/m
0
1.2
1.7
2.0
2.1
2.0
1.8
1.2
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）
plot（x0,y0,'k+',x,y2,'r'）
grid
title（'spline'）
运行结果：
（2）、二维插值问题的应用
（山区地貌图）在某山区（平面区域
内，单位:
m）测得一些地点的高度（单位:
m）如下表所示。
2400
1430
1450
1470
1320
1280
1200
1080
940
2000
1450
1480
1500
1550
1510
1430
1300
1200
1600
1460
1500
1550
1600
1550
1600
1600
1600
1200
1370
1500
1200
1100
1550
1600
1550
1380
800
1270
1500
1200
1100
1350
1450
1200
1150
400
1230
1390
1500
1500
1400
900
1100
1060
0
1180
1320
1450
1420
1400
1300
700
900
y/x
0
400
800
1200
1600
2000
2400
2800
试作出该山区的地貌图和等值线图。
试验程序：
x=0:
400:
2800;
y=0:
400:
2400;
h=[118013201450142014001300700900;
1230139015001500140090011001060;
12701500120011001350145012001150;
13701500120011001550160015501380;
14601500155016001550160016001600;
14501480150015501510143013001200;
1430145014701320128012001080940];
xi=0:
50:
2800;
xi=xi';
yi=0:
50:
2400;
z=interp2（x,y,h,xi,yi,'cubic'）;
figure
（1）
mesh（xi,yi,z）
figure
（2）
contour（x,y,h）
运行结果：
地貌图及等值线图如下：
三、本课程设计的心得体会（500字左右）
通过本次一周计算方法课程设计的训练，我进一步学习和掌握对程序的设计和编写，从中体会到了程序设计的方便和巧妙，了解到在进行编写一个程序之前，要有明确的目标和整体的设计思想，另外某些具体的细节内容也是相当的重要，这些宝贵的编程思想和从中摸索到的经验都是在编程的过程中获得的宝贵财富，这些经验对我以后的编程会有很大的帮助的，我要好好利用。
虽然这次课程设计是在参考程序的基础上进行的，但是我觉得对自己是一个挑战和锻炼。
我很欣慰自己能在程序中加入自己的想法和有关程序内容，也就是对它的程序改进了一番，并有创新。
但是我感觉自己的创新不够典型，总之还不是很满意。
另外由于时间的紧迫和对知识的了解不够广泛，造成了系统中还存在许多不足，功能上还不够完善。
以后我会继续努力，大胆创新。
这次课程设计让我充分认识到了自己的不足，认识到了动手能力的重要性。
我会在以后的学习中更加努力自己，提高自己，让自己写出更好更完善的程序，为以后的编程打好基础。
同时，一周的课程设计周，使我对计算方法这门课程有了更多的了解和收获，在不断的思考题目和编程中，又温习了一遍老师所教授的知识，对我大有裨益。
通过这次课程设计实习，我从中学到了很多，也真正领悟到了“态度决定一切”这句话的真正含义。
首先这次课程设计使我懂得了理论与实际相结合的重要性。
只有理论知识是远远不够的，只有把所学的理论
升级会员 