计算方法与计算 实验一误差分析

《计算方法与分析》实验报告

实验一、误差分析

一、上机前的准备工作

1、复习和掌握与本次实验有关的教学内容。

2、根据本次实验要求，在纸上编写算法及上机的程序，并经过人工模拟运行检验，减少不必要的错误，提高上机效率。切忌不编程序、不作人工检查就进行程序输入，这只能使上机调试的难度增加，甚至可能带来学习自信心的下降，影响后续课程的学习。

二、上机实验步骤

1、启动开发环境；

2、建立源程序文件，输入源程序；

3、编译产生目标程序，连接生成可执行程序，运行程序，输出结果；

4、对数值计算结果进行误差分析，讨论数值算法的收敛性与稳定性；

5、整理实验报告。

三、实验报告

实验报告是记录实验工作全过程的技术文档，实验报告的撰写是科学技术工作的一个组成部分。《数值分析》实验报告包括下列要求：

1、实验原理；

2、实验内容和要求；

3、数值算法描述，包括数据输入、数据处理和数据输出；

4、算法的实现

（1）给出具体的计算实例，

（2）经调试正确的源程序清单，

（3）对具体的数值例子给出数值结果；

5、计算结果的误差分析，算法的收敛性与稳定性的讨论；

6、实验心得。

误差分析

误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时，由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法的好坏会影响到数值结果的精度。 一、实验目的

1、通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；

2、通过上机计算，了解误差、绝对误差、误差界、相对误差界的有关概念；

3、 通过上机计算，了解舍入误差所引起的数值不稳定性。

二、算法实例 例 1.1 用差商h

a f h a f a f )

()()(-+≈

'求x x f ln )(=在3=x 处导数的近似值。取

1.0=h ，1000.0=h ，h =0.000 000 000 000 001和h =0.000 000 000 000 000 1分别用MATLAB 软件计算，取十五位数字计算。

解： 在MATLAB 工作窗口输入下面程序

>>a=3;h=0.1;y=log(a+h)-log(a);yx=y/h 运行后得

yx = 0.32789822822991. 将此程序中h 改为0.000 1，运行后得

yx = 0.33332777790385.

后者比前者好。再取h = 0.000 000 000 000 001，运行后得

yx = 0.44408920985006，

不如前者好。取h = 0.000 000 000 000 000 1，运行后得

yx = 0，

算出的结果反而毫无价值。

例1.2 分别求方程组b AX =在下列情况时的解，其中A ???

?

??=01.11

11

. （1）???? ??=22b ; （2）???

? ??=01.22b . 解： （1） 首先将方程组b AX =化为同解方程b A X 1-=，然后在MATLAB 工

作窗口输入程序

>> b=[2,2]';A=[1,1;1,1.01]; X=A\b

运行后输出当???? ??=22b 时，b AX =的解为???? ??=02X ； （2）同理可得，当???? ??=01.22b 时，b AX =的解为???? ??=11X .

例1.3 计算e 的近似值。 解：泰勒级数

e +++++++=!!4!3!21432　

　n x x x x x n

x

)(∞<<-∞x ,

取1=x ，得

e ++++++

+=!1

!41!31!2111　

　n . （1.1） 这是一个无限过程，计算机无法求到精确值。只能在（1.1）取有限项时计算，再估计误差。如果取有限项

!!!!)( n s n 1

413121111+

+++++=

作为e 的值必然会有误差，根据泰勒余项定理可知其截断误差为

e !)1()1(　

+=-n e s n θ

)10(<<θ.

如果取（1.1）的前九项，输入程序

>> n=8; s=1;S=1; for k=1:n s=s*k; S=S+1/s, end s,S,

R=3/(s*(n+1))

>>n=8; s=1;S=1; for k=1:n s=s*k; S=S+1/s; end s,S,

R=3/(s*(n+1))

运行结果：

s =40320 S =2.7183

R =8.2672e-006

或

>>S1=1+1+1/2+1/(1*2*3)+1/(1*2*3*4)+1/(1*2*3*4*5)+1/

(1*2*3*4*5*6)+1/(1*2*3*4*5*6*7)+1/(1*2*3*4*5*6*7*8),

R1=3/(1*2*3*4*5*6*7*8*9)

运行后结果

S =8.267195767195768e-006 R =2.71827876984127 颠倒！！！！！！！！！

因为截断误差为

e ),10(101968.267!93

!)18()1(6-8<

-θθ　

　e s 所以e 的近似值e ≈≈++++++++=!81!71!61!51!41!31!2111)1(8

s 2.718 28. 例1.4 取282.718作为e 的四舍五入近似值时，求其绝对误差和相对误差。 解：在MATLAB 工作窗口输入程序

>>juewu=exp(1)-2.71828

运行后输出结果为

juewu = 1.828 459 045 505 326e-006

例1.5 计算?2

sin π

x

x

d x 的近似值，并确定其绝对误差和相对误差。 解 因为被积函数

x

x

sin 的原函数不是初等函数，故用泰勒级数求之。 ++-+-=!

!!!1sin 42９　７５　３８

６x x x x x x )(∞<<-∞x , （1.5） 这是一个无限过程，计算机无法求到精确值。可用（1.5）的前四项

!

!!142７５　３６x x x -

+-代替被积函数x x sin ，得

?=2

sin π

x x y d ?≈20(πx !

!!142７５　３６x x x -+-)d x =!7)2(!5)2(!3)2(27

5

3

７５　３?-?+?-π

π

π

π=y ?. 根据泰勒余项定理和交错级数收敛性的判别定理，得到绝对误差

!

99)

2(?9

?<-=πy

y R = WU ， 在MATLAB 命令窗口输入计算程序如下：

syms x

f=1-x^2/(1*2*3)+x^4/(1*2*3*4*5)-x^6/(1*2*3*4*5*6*7) y=int(f,x,0,pi/2),y1=double(y)

y11=pi/2-(pi/2)^3/(3*3*2)+(pi/2)^5/(5*5*4*3*2)-(pi/

2)^7/(7*7*6*5*4*3*2)

inf=int(sin(x)/x,x,0,pi/2) ,infd=double(inf) WU =(pi/2)^9/(9*9*8*7*6*5*4*3*2), R =infd-y11

因为运行后输出结果为： =y 1.370 762 168 154 49，y

?=1.370 744 664 189 38，=R 1.750 396 510 491 47e-005， WU= 1.782 679 830 970 664e-005410-<. 所

以，y ?的绝对误差为=ε410-，故?=20sin π

x

x

y d 7 1.370≈x 。y

?的相对误差为 =r ε7

1.37010?4

-=

y ε<0.007 3%. 例 1.6 设计三种算法求方程01522

=-+x x 在)3,2(的一个正根*x 的近似

值，并研究每种算法的误差传播情况.

解：为解已知方程，我们可以设计如下三种算法，然后将计算结果与此方程的精确解5.2*=x 比较，观察误差的传播.

算法1 将已知方程化为同解方程=x 2

215x -.取初值20=x ，按迭代公式

2

1

215k k x x -=+ 依次计算 ,,,,21n x x x ，结果列入表1–3中。

算法2 将已知方程化为同解方程1

215

+=

x x .取初值20=x ，按迭代公式

12151+=

+k k x x

依次计算 ,,,,21n x x x ，结果列入表1–1中。

算法3 将已知方程化为同解方程1

415

22+-+-=x x x x x .取初值20=x ，按迭代

公式为 1

415

22

1

+-+-=+k k k

k k x x x x x 依次计算 ,,,,21n x x x ，结果列入表1–1中。

我们为这三种算法的计算编写两套MATLAB 程序如下：

（1）MATLAB 主程序

function [k,juecha,xiangcha,xk]= liti112(x0,x1,limax) % 输入的量--x0是初值, limax 是迭代次数和精确值x;

% 输出的量--每次迭代次数k 和迭代值xk,

% --每次迭代的绝对误差juecha 和相对误差xiangcha ， x(1)=x0;

for i=1:limax

x(i+1)=fl(x(i));%程序中调用的fl.m juecha = abs(x(i)-x1);

xiangcha = juecha /( abs(x(i))+eps);

xk=x(i);k=i-1;[(i-1),juecha,xiangcha,xk] end

（2）MATLAB 调用函数程序及其计算结果 ① 算法2的MATLAB 调用函数程序

function y1=fl(x) y1=15/(2*x+1);

② 将MATLAB 主程序和调用函数程序分别命名liti112.m 和fl.m ，分别保存为M 文件，然后在MATLAB 工作窗口输入命令

>> [k,juecha,xiangcha,xk]= liti112(2,2.5,100)

③ 运行后输出计算结果列入表1–1和表 1-2中。

④ 将算法2的MATLAB 调用函数程序的函数分别用y1=15-2*x^2和y1=x-(2*x^2+x-15)/(4*x+1)代替，得到算法1和算法3的调用函数程序，将其保存，运行后将三种算法的前8个迭代值821,,,x x x 列在一起（见表 1-1），进行比较.将三种算法的821,,,x x x 对应的绝对误差和相对误差的值列在一起（见表 1-2），进行比较。

表 1-1 例1.6中三种算法的计算结果

算 法

迭代次数

算法1的迭代结果

算法2的迭代结果

算法3的迭代结果

2 2.000 000 00 2.000 000 00 1 7 3.000 000 00 2.555 555 56 2 -8

3 2.142 857 1

4 2.500 550 06 3 -13 763 2.837 837 84 2.500 000 06 4 -378 840 323 2.246 963 56 2.500 000 00

5 -2.870 41710? 2.24

6 963 56 2.500 000 00 6 -1.64

7 83510? 2.321 774 84 2.500 000 00 7

-5.430 77010?

2.657 901 65

2.500 000 00

99

-Inf

2.500 000 01

2.500 000 00

表 1-2 例1.6中三种算法计算结果的误差

算法 迭 代 次 数

算法1的误差 算法2的误差 算法3的误差 绝对误差

相对误差

绝对误差

相对误差

绝对误差

相对误差

0 0.500 000 00 0.250 000 00 0.500 000 00

0.250 000 00 0.500 000 00 0.250 000 00 1 4.500 000 00 0.642 857 14 0.500 000 00 0.166 666 67 0.055 555 60 0.021 739 13 2 85.500 000 00 1.030 120 48 0.357 142 86 0.1666 666 70 0.000 550 10 0.000 219 97 3 13 765.500 00 0.000 100 02 0.337 837 84 0.119 047 62 0.000 000 06 0.000 000 02 4 378 840 326 1.000 000 01

0.253 036 44 0.112 612 62 0.000 000 00

0.000 000 00

5 2.870 399 811710?

1 0.230 287 04 0.084 345 48 0 0 6 1.647 839 013510?

1 0.178 225 16 0.076 76

2 47 0 0 7 5.430 746 807010?

1 0.157 901 65 0.059 408 39 0 0

99

Inf

NaN

0.000 000 01

0.000 000 00

例1.7 求数）

181(71915-+?=-x 的近似值。 解 （1）直接用MATLAB 命令

>> x=(7^15)*(sqrt(1+8^(-19))-1)

运行后输出结果

x = 0.

问题出现在两个相近的数19

81-+与1相减时，计算机运行程序

>>sqrt(1+8^(-19))-1

运行后输出结果

ans = 0.

由于计算机硬件只支持有限位机器数的运算，因此在计算中可能引入和传播舍入误差.因为有效数字的严重损失，导致输出18119-+-的结果为0，计算机不能再与数157继续进行真实的计算，所以，最后输出的结果与x 的精确值不符。

（2）如果化为

1

8

187)18

1(719

191519

15++?=

-+?=---x ，

再用MATLAB 命令

>> x=(7^15)*( (8^(-19))/(sqrt(1+8^(-19))+1))

运行后输出结果

x = 1.6471e-005

这是因为18

119

-+-化为

1

8

1819

19++--后，计算机运行程序

>> x= (8^(-19))/(sqrt(1+8^(-19))+1)

运行后的结果为

x =3.4694e-018

由于有效数字的损失甚少，所以运算的结果-18103.4694?再与157继续计算，最后输出的结果与x 的精确值相差无几。

例1.8 求数)13030ln(2--=y 的近似值。 解 （1）直接用MATLAB 程序

>> x=30;x1= sqrt(x^2-1)

运行后输出结果

x1 = 29.9833

输入MATLAB 程序

>> x=30; x1=29.9833;y=log(x-x1)

运行后输出结果

y = -4.0923

（2）因为)13030ln(2--中的30=x 很大，如果采用倒数变换法

1

112

21-+=

--=x x x x z ，

即

1

30301ln

)13030ln(2

2-+=--)190030ln(-+-=.

输入MATLAB 程序

>> x=30;y=-log(x+sqrt(x^2-1))

运行后输出结果

y = -4.0941

（3）输入MATLAB 程序

>> x=30; y=log(x-sqrt(x^2-1))

运行后输出结果

y = -4.0941

可见，（2）计算的近似值比（1）的误差小。

参加计算的数，有时数量级相差很大.如果不注意采取相应的措施，在它们的加减法运算中，绝对值很小的那个数经常会被绝对值较大的那个数“吃掉”，不能发挥其作用，造成计算结果失真。

例1.9 请在16位十进制数值精度计算机上利用软件MA TLAB 计算下面的两个数

0.30.1111111111111111*++=x 和0.30.11111111111111111*++=y

将计算结果与准确值比较，解释计算结果。 解 在MA TLAB 工作窗口输入下面程序

>> x=111111*********+0.1+0.3, y=1111111111111111+0.1+0.3

运行后输出结果

x = 1.111111*********e+014，y =1.111111*********e+015

从输出的结果可以看出，x *x =，而y *

y ≠.为什么*

y 仅仅比*x 多一位1，而y *

y ≠呢？

这是因为计算机进行运算时，首先要把参加运算的数写成绝对值小于1而“阶码”相同的数，这一过程称为数的“对阶”。在16位十进制数值精度计算机上利用软件MATLAB 计算这两个数，把运算的数*

x 写成浮点规格化形式为

，

1515

15*103000**********.0001010000000000000.000100111111111111111.0?+?+?=x

在16位十进制数值精度计算机上，三项的数都表示为小数点后面16位数字的数与15

10

之积，所以，计算机没有对数进行截断，而是按原来的三个数进行计算。因此，计算的结果x *

x =。而

16

16

16*10030000000000000.00010010000000000000.000101111111111111111.0?+?+?=y

三项的数都表示写成绝对值小于1而“阶码”都为16

10的数以后，第一项的纯小数的

小数点后面有16位数字.但是，后两项数的纯小数的小数点后面有17位数字，超过了16位十进制数值精度计算机的存储量，计算机对后两项的数都进行截断最后一位，即后两项的数都是16位机上的零，再进行计算，所以计算结果与实际不符。

三、实验任务

对20,,2,1,0 =n ，计算定积分

?+=1

5dx x x y n

n .

算法1：利用递推公式

151

--=

n n y n

y , 20,,2,1 =n , 取 ?

≈-=+=1

0182322.05ln 6ln 5

1

dx x y . >> a=0.182322;n=1;b=(1/n)-(5*a)

算法2：利用递推公式

n n y n y 5

1

511-=

- 1,,19,20 =n . 注意到

???=

≤+≤=1

1

020201

0201051

515611261dx x dx x x dx x , 取 008730.0)126

11051(20120≈+≈

y . 思考：从计算结果看，哪个算法是不稳定的，哪个算法是稳定的。

>> a=0.008730;n=20;b=(1/5)*1/n-(1/5)*a

解:

10:,(6),:

n y 1

n n n

-≈=

由递推公式可得当 =1,2,,20 时以及0.182322=0.182322 小数点后面取位使用Matlab 计算得出下列各值 -5y y

Ni

由以上计算结果可以看出，当n ≥8时，出现负值，完全失真，但与原题中：

编号

计算结果

1 0.0884

2 0.0580

3 0.0433

4 0.033

5 5 0.0325

6 0.0042

7 0.1219

8 -0.4845

9 2.5336 10 -12.5680 11 62.9309 12 -314.5712 13 1.5729e+003 14 -7.8644e+003 15 3.9322e+004 16 -1.9661e+005 17 9.8305e+005 18 4.9153e+006 19 -2.4576e+007 20

1.2288e+008

01111

01000

1

15555min max x n n n

x x x x x x dx dx dx

x x ≤≤≤≤++<<++???116(1)

5(1)

n n n y ++<<

完全不符，

~~

11|||||5()|n n n n n E y y y y --=-=--

02020||5||n n E E E =误差被放大5

因此第一种算法是不稳定的。

我们来看第2种方法：

-120

:20191,(6)

,:

n 1

1

55

n n y y n =≈由递推公式可得

当 =，，

， 时以及0.008730=0.008730 小数点后面取位使用Matlab 计算得出下列各值 -

y

编号

计算结果

1

0.1823 2 0.0884 3 0.0580 4 0.0431 5 0.0343 6 0.0285 7 0.0243 8 0.0212 9 0.0188 10 0.0169 11 0.0154 12 0.0141 13 0.0130 14 0.0120 15 0.0112 16 0.0105 17 0.0099 18 0.0093 19 0.0089 20 0.0083

由上看来当n 减少时,

也随之降低,并没有波动,因此看来也较为稳定。

2020

1

||||

5

1

|5m n n m

E E E -=

|误差被缩小了

倍

四、实验心得

对于这门课程来说，是第一次接触，从专升本上来的我，对于Matlab 软件操作简直是一窍不通，论是带公式或者是编写代码。这是我学习这门课程的最大障碍！第一次上课，听老师讲解Matlab 软件的操作，我学会了一点，但也局限于把文档里的代码，往Matlab 里一粘贴，看看它是否与文档中，所记录的答案一致。有的错，就用红色的标注，我也只能做到这里了，用Matlab 计算一些，可能按按计算器也能得出的结果，可谓是心有余而力不足。 这是我第一次上机课所学会的，比如一些基本的格式： 一般MATLAB 命令格式为：

[输出参数1，输出参数2，……]=（命令名）（输入参数1，输入参数2，……） 输出参数用方括号，输入参数用圆括号如果输出参数只有一个可不使用 括号。

对于Matlab 软件来说，能操作好它，对于眼前来说，能将自己的毕业论文写好，往远了说，在以后的工作中，碰到此类问题，也有能力解决。

对于课程，本人还有一些想法，既然是上机课，我觉得学会Matlab 是很重要的，能不能针对Matlab 来多做些习题，毕竟上机时间比较充裕，因而曾加习题，从易到难，提高对Matlab 的操作水平！！

~

~111

|||||()|

5

n n n n n E y y y y --=-=--~

n y