计算方法课程设计

数理学院2014级信息与计算科学

课程设计

姓名：刘金玉

学号： 3141301240

班级： 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 求一元二次方程的根 实验目的：

研究误差传播的原因与解决对策。 问题提出：求解一元二次方程20ax bx c ++=

实验内容：

一元二次方程的求根公式为

1,22b x a

-+=

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

2210(1)320(2)1010

x x x x +-=-+=

实验要求：

（1） 考察单精度计算结果（与真解对比）；

（2） 若计算结果与真解相差很大，分析其原因，提出新的算法（如先求1x 再

根据根与系数关系求2x ）以改进计算结果。 实验步骤： 方程（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.5632798704

X2= -3.56327987

方程（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*a*c))/2*a

运行结果为：

x1 = 1.0000e+010

x2 = 0

然后

由符号解求的该方程的真值程序为：

syms x;

y=x^2-10^10*x+1;

s=solve(y,x);

vpa(s)

运行结果：

X1= 10000000000.0

X2= 0.6415321827

由此可知，对于方程（1），使用求根公式求得的结果等于精确值，故求根公式法可靠。而对于方程（2），计算值与真值相差很大，故算法不可靠。

改进算法，对于方程（2）:先用迭代法求x1，再用11a

x x

，求x2

c

程序为：

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.0000e+010 x2 = 1.0000e-010 实验结论：

对于方程（1），两种方法在精确到小数点后15位时相同，说明两种算法的结果都是精确的。

对于方程（2），两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算法是不精确的。

原因：方程（2

）用求根公式计算时，b -b 是大数，出现了大数吃掉小数的误差，也出现了两个相近的数相减的误差，所以出现x2=0这样大的误差。改进的结果会比较准确。

2、求解非线性方程

实验2.1 Gauss 消去法的数值稳定性实验

实验目的：观察和理解高斯消元过程中出现小主元即()

k kk a 很小时，引起方程组

解的数值不稳定性．

实验内容：求解线性方程组(1,2)i i A x b i == 其中

（1）15

10.31059.14

3

15.291

6.1301211.29521211A -???

?--

?= ?

? ???，159.1746.7812b ??

?

?= ? ?

??

（2）

210701 3 2.09999999999962 5151 0102

A

-

??

?

- ?

=

?

--

?

??

1

8

5.9000000000001

5

1

b

??

?

?

=

?

?

??

实验要求：

（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.294470872045e+000

ans = 68.4295577125341e+000

ans = 84.3114602051800e+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.2831683168042e+000

ans = 8.99393809015365e+000

ans = 18.3564356435280e+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

运行结果：

18100026.791010

00.47240.1755100.08930.02020.49291L ?

?????

?=??-?

?--?? 11.2

952059.143100 2.835 1.2310

000.801U ??????=??

-??

?? 1[18.9882;3.3378;34.747;33.9865]x =--

方程组（2）：

A2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2]; b2=[8;5.91;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(1:j-1)=b2(1:j-1)-b2(j)*U(1:j-1,j); end

b2(1)=b2(1)/U(1,1); x2=b2 运行结果：

121

00

00.51000.30.410100

0.40.3331L ?????

?=??

--??

?

-?? 1070

10 2.55 1.5006 2.300

00 3.3667U -????-??=?????? 112[0.44410; 1.0;1.0;1.0]x -=?-

（3） 不选主元的高斯消去法程序： 方程组（1）： clear format long

A1=[0.3e-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

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(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); end

b1(1)=b1(1)/U(1,1); x1=b1 运行结果:

15

151510

0017.63710100'37.3310 2.1168103.333100.18900L ?

?

????

?=?????

??

?

151815150.31059.14310 1.0431052.911017.63710'001680000.5U ?????-?-?-???=??--??

??

'1[23.6848;1.0005;0;0]x =

方程组（2）： clear format long

A2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2]; b2=[8;5.91;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 运行结果：

12121

0000.310

0'0.5 2.499810100

0.9999100.4001L ??

??-?

?=??

-??

?

-???

12

1212

10

701

0 1.0106

2.3'0014.998710 5.7495100

00.400 3.3667U -????-??

?=??

???

?

??

'

52[0.36210; 1.0;1.0;1.0]x -=?-

（4） 分析小元对计算结果的影响

通过观察计算结果，分析可知，小元对计算结果的影响很大，小元的存在会使得到的计算结果有很大的误差。

3、插值

4、数值积分

实验2.4：复化求和公式计算定积分 实验内容：

计算下列各式右端定积分的近似值

3122

2011

(1) ln 22 (2) =411dx dx x x π=--+??

122

012(3) 3 (4) e ln 3x x dx xe dx ==?? 实验要求：

（1） 分别用复化梯形公式、复化Simpson 公式计算，要求绝对误差限为

70.510ε-=?；

（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;

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

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 个积分常数为(,,,,,)T a e i M ωΩ，其中a 表示轨道半长径，e 表示轨道偏心率，i 表示轨道倾角， Ω是升交点赤经， ω 是近地点辐角，M 是平近点角。其中第六个积分常数，通常可以用其他的一些量替换，如过近地点时刻p t ，真近点角度 f 与偏近点角度E 。这几个是相互等价的关系，我们这里要讲的就是平近点角 M 与偏近点角E 的一个关系，有如下方程

sin E e E M -=

这就是著名的 Kelper 方程，是在人造卫星运动理论或者行星运动理论中最基本的方程之一。 因为如果我们已经知道卫星轨道量去计算卫星的星历时， 第 6 个根数通常是使用平近点角。这时候需要把平近点角化为偏近点角，也就是需要计算上述的 Kepler 方程。具体理论这里不再阐述，而 Kepler 方程本身是我们这里感兴趣的。可采用 Newton 迭代法计算 Kepler 方程，在实际工程中，一般也是这样做的，因为 Newton 迭代法收敛速度快而且精度比较高。也可用其它算法计算。 实验目的与要求 （1）了解 Kepler 方程；

（2）能够用牛顿迭代方法计算 Kepler 方程； （3）通过调节轨道偏心率 e 查看运算收敛情况。

实验内容及数据来源

编写解 kepler 方程的方法，给定偏心率为 0.01、平近点角为 32 度时计算出偏近点角E 。

实验步骤：

（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

x0=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）、一维插值问题的应用

加工时需要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') 运行结果：

piecewise linear

spline

（2）、二维插值问题的应用

(山区地貌图)在某山区（平面区域()()0,28000,2400?内，单位:m ）测得一些

试验程序：

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)

运行结果：

地貌图及等值线图如下：