计算方法课程设计.doc

计算方法课程设计

32
2020年4月19日
文档仅供参考，不当之处，请联系改正。

数理学院

级信息与计算科学

课 程 设 计

姓名： 刘金玉
学号：
班级： 1402
成绩：

实验要求
1． 应用自己熟悉的算法语言编写程序，使之尽可能具有通用性。
2． 上机前充分准备，复习有关算法，写出计算步骤，重复检查，调试程序。（注：在练习本上写，不上交）
3． 完成计算后写出实验报告，内容包括:算法步骤叙述，变量说明，程序清单，输出计算结果，结构分析和小结等。（注：具体题目具体分析，并不是所有的题目的实验报告都包含上述内容！）
4． 独立完成，如有雷同，一律判为零分！
5． 上机期间不允许做其它任何与课程设计无关的事情，否则被发现一次扣10分，被发现三次判为不及格！非特殊情况，不能请假。旷课3个半天及以上者，直接判为不及格。

目 录

一、基本技能训练 4
1、误差分析 4
2、求解非线性方程 6
3、插值 12
4、数值积分 12
二、提高技能训练 16
1、 16
2、 18
三、本课程设计的心得体会（500字左右） 21

一、基本技能训练
1、误差分析
实验1.3 求一元二次方程的根
实验目的：
研究误差传播的原因与解决对策。
问题提出：求解一元二次方程
实验内容：
一元二次方程的求根公式为

用求根公式求解下面两个方程：

实验要求：
（1） 考察单精度计算结果（与真解对比）；
（2） 若计算结果与真解相差很大，分析其原因，提出新的算法（如先求再根据根与系数关系求）以改进计算结果。
实验步骤：
方程（1）：
根据求根公式，写出程序：
format long
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.5630
x2 = -3.5630
然后
由符号解求的该方程的真值程序为：
syms x;
y=x^2+3*x-2;
s=solve(y,x);
vpa(s)
运行结果为：
X1= 0.
X2= -3.4927987
方程（2）：
format long
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*ac))/2*a
运行结果为：
x1 = 1.000e+010
x2 = 0
然后
由符号解求的该方程的真值程序为：
syms x;
y=x^2-10^10*x+1;
s=solve(y,x);
vpa(s)
运行结果：
X1= .0
X2= 0.
由此可知，对于方程（1），使用求根公式求得的结果等于精确值，故求根公式法可靠。而对于方程（2），计算值与真值相差很大，故算法不可靠。
改进算法，对于方程（2）:先用迭代法求x1，再用，求x2
程序为：
syms k
a=[ ];
for i=1:100
a(1)=4;
a(i+1)=(10^10*a(i)-1)^(1/2);
end
x1=a(100)
x2=1/(x1)
运行结果为：
x1 = 1.000e+010
x2 = 1.000e-010
实验结论：
对于方程（1），两种方法在精确到小数点后15位时相同，说明两种算法的结果都是精确的。
对于方程（2），两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算法是不精确的。
原因：方程（2）用求根公式计算时，公式中，b是大数，出现了大数吃掉小数的误差，也出现了两个相近的数相减的误差，因此出现x2=0这样大的误差。改进的结果会比较准确。
2、求解非线性方程
实验2.1 Gauss消去法的数值稳定性实验
实验目的：观察和理解高斯消元过程中出现小主元即 很小时，引起方程组解的数值不稳定性．
实验内容：求解线性方程组 其中
（1），
（2）
实验要求：
（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的
（2）用高斯列主元消去法求得L和U及解向量
（3）用不选主元的高斯消去法求得L和U及解向量
（4）观察小主元并分析对计算结果的影响
实验步骤：
（1） 计算矩阵的条件数程序：
矩阵A1：
A1=[0.3*10^(-15) 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];
cond(A1,1)
cond(A1,2)
cond(A1,inf)
运行结果：
ans = 136.e+000
ans = 68.1e+000
ans = 84.31e+000
由矩阵条件数判断出矩阵A1是病态矩阵。
矩阵A2：
A2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2];
cond(A1,1)
cond(A1,2)
cond(A1,inf)
运行结果：
ans = 19.2e+000
ans = 8.65e+000
ans = 18.0e+000
由矩阵条件数判断出矩阵A2是病态矩阵。
（2） 高斯列主元消去法程序：
方程组（1）：

A1=[0.3*10^(-15) 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];
b1=[59.17;46.78;1;2];
n=4;
for k=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;
if A1(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);
for i=1:n
L1(i,i)=1;
end
L=L1
U=triu(A1,0)
for j=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);
for j=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 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2];
b2=[8;5.1;5;1];
n=4;
for k=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;
if A2(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);
for i=1:n
L1(i,i)=1;
end
L=L1
U=triu(A2,0)
for j=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);
for j=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
format long
A1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 211.2 9 5 2;1 2 1 1];
b1=[59.17;46.78;1;2];
n=4;
for k=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);
for i=1:n
L1(i,i)=1;
end
L=L1
U=triu(A1,0)
for j=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);
for j=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
format long
A2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2];
b2=[8;5.1;5;1];
n=4;
for k=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);
for i=1:n
L1(i,i)=1;
end
L=L1
U=triu(A2,0)
for j=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);
for j=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;
syms x;
f=1/(x^2-1);
n=8;
h=(b-a)/n;
s=0;w=0;
for i=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;
for i=1:n-1
if rem(i,2)~=0
x1=a+i*h;
f1=f1+subs(f,x1);
else rem(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;
syms x;
f=1/(x^2+1);
n=100;
h=(b-a)/n;
s=0;w=0;
for i=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;
for i=1:n-1
if rem(i,2)~=0
x1=a+i*h;
f1=f1+subs(f,x1);
else rem(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;
syms x;
f=3^x;
n=100;
h=(b-a)/n;
s=0;w=0;
for i=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;
for i=1:n-1
if rem(i,2)~=0
x1=a+i*h;
f1=f1+subs(f,x1);
else rem(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;
syms x;
f=x*exp(x);
n=100;
h=(b-a)/n;
s=0;w=0;
for i=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;
for i=1:n-1
if rem(i,2)~=0
x1=a+i*h;

f1=f1+subs(f,x1);
else rem(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&&k
方程的根为 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=[0 3 5 7 9 11 12 13 14 15];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.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('piecewise linear') %图片命名
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

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

2400
2800
试作出该山区的地貌图和等值线图。
试验程序：
x=0:400:2800;
y=0:400:2400;
h=[1180 1320 1450 1420 1400 1300 700 900 ;
1230 1390 1500 1500 1400 900 1100 1060;
1270 1500 1200 1100 1350 1450 1200 1150;
1370 1500 1200 1100 1550 1600 1550 1380;
1460 1500 1550 1600 1550 1600 1600 1600;
1450 1480 1500 1550 1510 1430 1300 1200;
1430 1450 1470 1320 1280 1200 1080 940];
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字左右）
经过本次一周计算方法课程设计的训练，我进一步学习和掌握对程序的设计和编写，从中体会到了程序设计的方便和巧妙，了解到在进行编写一个程序之前，要有明确的目标和整体的设计思想，另外某些具体的细节内容也是相当的重要，这些宝贵的编程思想和从中摸索到的经验都是在编程的过程中获得的宝贵财富，这些经验对我以后的编程会有很大的帮助的，我要好好利用。
虽然这次课程设计是在参考程序的基础上进行的，可是我觉得对自己是一个挑战和锻炼。我很欣慰自己能在程序中加入自己的想法和有关程序内容，也就是对它的程序改进了一番，并有创新。可是我感觉自己的创新不够典型，总之还不是很满意。另外由于时间的紧迫和对知识的了解不够广泛，造成了系统中还存在许多不足，功能上还不够完善。以后我会继续努力，大胆创新。
这次课程设计让我充分认识到了自己的不足，认识到了动手能力的重要性。我会在以后的学习中更加努力自己，提高自己，让自己写出更好更完善的程序，为以后的编程打好基础。同时，一周的课程设计周，使我对计算方法这门课程有了更多的了解和收获，在不断的思考题目和编程中，又温习了一遍老师所教授的知识，对我大有裨益。
经过这次课程设计实习，我从中学到了很多，也真正领悟到了“态度决定一切”这句话的真正含义。首先这次课程设计使我懂得了理论与实际相结合的重要性。只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，把理论付诸实际行动，从实践中得出结论，才能更加深刻的理解书本所学知识，从而提高自己的实际动手能力和独立思考的能力。