使用插值法公式组成数字电路进行计算的计算机-转让饼屋-中国面包师贴吧/中国烘焙网贴吧/烘焙师招聘求职贴吧

于是，由x =x +h,x =x +2h,...便可得出

1 0 2 0

x-x x-(x +h) x-x -h x-x

1 0 0 1 h

= = = - =u-1

h h h h h

x-x x-(x +2h) x-x

2 0 0 2h

= = - =u-1

h h h h

……………………………………..

x-x x-[(x +(n-1)h] x-x

n-1 0 0 (n-1)h

= = - =u-(n-1)=u-n+1

h h h h

(x-x ) (x-x )

0 1

把这些 , 等值代入（3），则得

h h

u(u-1) 2 u(u-1)(u-2) 2

φ(x)=φ(x +hu)=g(u)=y +u△y + △ y + △ y +

0 0 0 2! 0 3! 0

u(u-1)(u-2)(u-3) 4 u(u-1)(u-2)...(u-n+1) n

△ y +…+ △ y (1)

4! 0 n! 0

这就是通常所写的牛顿向前插值公式的形式，此后我们将说它是牛顿公式（Ⅰ）。

△的系数是二项式的系数。

注：晚近的历史研究证明，这个公式实在是詹姆士*格列高利(James Gregory)早在1670年首先发明的。按我国九章算术刘徽辑注本（公元263年）第七章“盈不足”术，即是一种最简单的一次插值法。后来在历法应用上，这法得到推广。隋刘焯的皇极历（公元600年著）开始用等距二次插值公式，唐李淳风的麟德历（公元664年著）也用这法，但均未提出这种算法的名称。唐张遂（即僧一行）的大衍历（公元724年著），有“步日漉术”，是不等距的二次插值法，又有“步月漉术”和“步五星术”，则是等距的情形。这些称谓，皆是就历法应用上而命名的。宋秦九韶的数学九章（1247），始应用这类方法到历法以外的问题，而称为“招法”。先此有宋沈括（1031-1095）的“圆法”和“妥法”等名称，惜法均不详。元郭守敬的授时历（1281）有“平立定三差法”，把插值法推到等距三次的情形。至元朱世杰的“四元玉监”（1303），把插值法叫做“招差术”，演算方法与现在插值法相同。可见祖国从六世纪到十三世纪，对内插法已广泛应用，比西洋早一千多年。可见李严著“中算家的内插法研究”一卷（1957年，科学出版社本）。

它之所以叫“向前”插值公式，

是因为公式中所包含的函数列表值都是自y 起向上进到它的右边去的（自y 前向），

0 0

而没有一个是该值进到左边去的缘故。由于这个缘故，这个公式主要用在一组列表值中内插靠近开端部分的y值，以及外插自y 后退（向左）一段短距离的y值之用。

始点y 可以是任何一个列表值，而y仅含有选作始点的后继y值。

21.牛顿向后插值公式

前节诸公式不能用来内插列表值中靠近末端的y值。要对这种情况导出一个公式，我们就写出下列形式的多项式φ(x)：

φ(x)=a +a (x-x )+a (x-x )(x-x )+a (x-x )(x-x )*(x-x )

0 1 n 2 n n-1 3 n n-1 n-2

+a (x-x )(x-x )(x-x )(x-x ) +...+a (x-x )(x-x )(x-x )...(x-x ) (1)

4 n n-1 n-2 n-3 n n n-1 n-2 1

然后我们确定系数a ,a ,a ,...,a ，以便使得

0 1 2 n

φ(x )=y ,φ(x )=y ,...,

n n n-1 n-1

将x的各值x ,x 等等代入（1），同时命

n n-1

φ(x )=y ,φ(x )=y ,...,

n n n-1 n-1

则得

y =a 或a =y

n 0 0 n

y =a +a (x -x )=y +a (-h)

n-1 0 1 n-1 n n 1

y -y △ y

n n-1 1 n

a = =

1 h h

y =a +a (x -x )+a (x -x )(x -x )

n-1 0 1 n-1 n 2 n-2 n n-2 n-1

y -y

n n-1

=y + (-2h)+a (-2h)(-h)

n h 2

y -2y +y △ y

n n-1 n-2 2 n

a = =

2 h 2

继续进行这种计算系数的方法，我们便可得到

△ y

3 0

a = ,

3 3

3!h

△ y

4 n

a = ,…,

4 4

4!h

△ y

n 0

a = ,

n n

n!h

将这些a ,a ,a 等等各值代入(1)则得

0 1 2

△ y △ y

1 n 2 n

φ(x)=y + (x-x )+ (x-x )(x-x )+

h n 2 n n-1

△ y

3 n

+ (x-x )(x-x )(x-x )+

3 n n-1 n-2

3!h

△ y

4 n

+ (x-x )(x-x )(x-x ) (x-x )+…+

3 n n-1 n-2 n-3

4!h

△ y

n n

+ (x-x )(x-x )(x-x )(x-x )...(x-x ) (2)

n n n-1 n-2 n-3 1

n!h

这就是用x表示的牛顿向后插值公式。它可以像第20节那样，使用变数更换法来加以简化。首先我们写出（2）的等价形式

(x-x ) △ y (x-x ) (x-x )

n 2 n n n-1

φ(x)=y +△ y +

n 1 n h 2 h h

△ y (x-x ) (x-x ) (x-x )

3 n n n-1 n-2

+ +

3! h h h

△ y (x-x ) (x-x ) (x-x ) (x-x )

4 n n n-1 n-2 n-3

+ +…+

4! h h h h

△ y (x-x ) (x-x ) (x-x )

n n n n-1 1

+ …… (3)

3! h h h

现在令

x-x

或

x=x +hu

于是由

x =x +h,x =x +2h,...,

n-1 n n-2 n

则得

x-x x-(x -h) x-x +h x-x

n-1 n n n h

= = = + =u+1

h h h h h

x-x x-(x -2h) x-x

n-2 n n 2h

= = + =u+2

h h h h

……………………..

x-x x-[x -(n-1)h] x-x

1 n n (n-1)h

= = + =u+n-1

h h h h

(x-x ) (x-x )

n n-1

把这些 , 等值代入（3），则得

h h

u(u-1) u(u-1)(u-2)

φ(x)=φ(x +hu)= φ(u)=y +u△y + △ y + △ y +

0 0 0 2! 2 0 3! 3 0

u(u-1)(u-2)(u-3) u(u-1)(u-2)...(u-n+1)

△ y +…+ △ y (Ⅱ)

4! 4 0 n! n 0

这就是通常所写的牛顿向后插值公式的形式。此后我们将说它是牛顿公式（Ⅱ）。可以看出：这个公式使用水平差分，而向前插值公式使用对角差分。（Ⅱ）之所以被称为“向后”插值公式，这是因为它所含的列表值，

都是自y 后退到左边去的而没有一个是自y 进到右边去的缘故，

n n

这个公式主要用来插值：

一组列表值中靠近末端的y值，以及外插自y 向前（向右）一段短距离的y值之用。

现在我们举一些例题来说明牛顿公式的用法。

例1.试求lgπ已知

lg3.141=0.4970679364, lg3.142=0.4972061807, lg3.143=0.4973443810,

lg3.144=0.4974825374, lg3.145=0.4976206498,

解：我们先做出如下的差分表：

2 3

x y=lgx △y △ y △ y

3.141 0.4970679864

0.0001382443

3.142 0.4972061807 -0.0000000440

0.0001382003 0.0000000001

3.143 0.4972063810 -0.0000000439

0.0001381564 -0.0000000001

3.144 0.4974825374 -0.0000000440

0.0001381124

3.145 0.4976206498

这里x=π=3.1415926536，x =3.141,h=0.001,故

x-x 3.1415926536-3.141

u= 0 = =0.5926536

h 0.001

u-1=-0.4073464,...依次类推

以这些值代入第20节中的{Ⅰ)式，即得, 第20节中的{Ⅰ)式如下所示

u(u-1) 2 u(u-1)(u-2) 2

φ(x)=φ(x +hu)=g(u)=y +u△y + △ y + △ y +

0 0 0 2! 0 2! 0

u(u-1)(u-2)(u-3) 4 u(u-1)(u-2)...(u-n+1) n

△ y +...+ △ y (Ⅰ)

4! 0 n! 0

0.5926536(-0.4073464)(-440)

lgπ=0.4970679364+0.5926536（1382443）+

=0.4970679364+0.000819310+0.0000000053=0.4971498727

此结果准确到它的末位数字。

注：像在上面的例题中，这种在给定诸值的范围以外计算函数值的过程称为外插法。使用它时应该慎重从事，不过假设已经知道函数在靠近给定诸值范围的端点是平滑的变动的，并且假如h取得尽量的小，那么我们在给定诸值范围以外相距h的地方进行外插，通常是很稳妥的。

例3.在1918年1月1日亮每小时的赤纬已在下表给出，试求在北京时间3时35分15秒时的赤纬。

小时赤纬 △ y

1 △ y

2 △ y

0 8°29`53.7``

1 8°18`19.4`` -11`34.2``

2 8°6`43.5`` -11`35.9`` -1.6``

3 7°55`6.1`` -11`37.4`` -1.5`` 0.1``

4 7°43`27.2`` -11`38.9`` -1.5`` 0.0``

解：因为所求赤纬接近于给定诸值的末端，所以我们使用牛顿公式（Ⅱ），同时我们作出上列水平差分表。用t来表示以小时计算的时间则，

t =4,t=3时25分15秒，h=1，故

t-t 8

n -0时24分45秒 -1485

u= = = =-0.3569

h h 8

3600

u+1=0.6431

将这些值代入（Ⅱ）并用δ表示所需求的赤纬，则得

(0.6431)(-0.3569）

δ=7°43`27``.2+(-0.3569)(-11`38``.9）+ (-1``.5)

=7°43`27``.2+4`9``.4+0``.2

=7°43`36``.8

例4.使用上题的数据，试求t=5时时月亮的赤纬。

解：这里，t=t =5,t =4

n+1 n

t -t

n+1 n h

u= = =1,u+1=2

h h

代入（Ⅱ）可得

(1)(2)

δ=7°43`27``.2+(1)(-11`38``.9)+ (-1``.5)=7°31`46``.8

在美国天文历与航海历(American Ephemeris and Nautical Almanac)中所给的真值为7°31`46``.9, 因此外插值的误差只不过0``.1而已。

第三章插值法中心差分公式

22.引言

牛顿公式（Ⅰ）和（Ⅱ）是基本的，并且对于几乎所有的插值情况来说，都是可以应用的，不过就一般而论它不如另一类公式所谓中心差分公式收敛得快。后一类公式所用的差分，都尽可能取在通过对角差分表所作水平线的附近。浏览表3，可以看出：这些差分所含的各函数值，都在通过该函数值所在水平线的上下。所以中心差分公式特别适合于内插靠近列表值中央部分的函数值之用。最重要的中心差分公式是两个已知的所谓斯特灵(Stirling)公式和贝赛尔公式。

注：中心差分公式还有高斯的一种。但是柯姆利(J.L.Comrie)曾指出，它实际上很少用到。皮尔逊(K.Pearson)也指摘这个公式，杰佛里(Jeffreys)认为这公式确是无用，但它与牛顿公式有直接关系，最易推导，由它就很容易推出其他公式。关于高斯公式以及由它证明斯特灵公式和贝赛尔公式的过程，可参看别席考维奇同书第43到45节（中译本125-129页）。推导它们的方法有好几种，不过最简单的一种是由牛顿公式（Ⅰ）用代数变换法来推导的。

23.斯特灵插值公式

要推导斯特灵公式，我们首先写出对角差分表，

然后把列表值y ,以及跟通过y 所作的水平线紧挨着的那些差分，

0 0

都标记出来，以便加以特别注意。这些量在表7中，系用黑体字印出。

以y 为始点的牛顿公式（Ⅰ）写为

u(u-1) 2 u(u-1)(u-2) 2

y=y +u△y + △ y + △ y +

0 0 2！ 0 3！ 0

u(u-1)(u-2)(u-3) 4 u(u-1)(u-2)(u-3)(u-4) 5

+ △ y + △ y+... (Ⅰ)

4！ 0 5！

表7

△y 2

△ y 3

△ y 4

△ y 5

△ y 6

△ y 7

△ y 8

△ y

-4

△y

-4

-3 2

△ y

-4

△y

-3 3

△ y

-4

-2 2

△ y

-3 4

△ y

-4

△y

-2 3

△ y

-3 5

△ y

-4

-1 2

△ y

-2 4

△ y

-3 6

△ y

-4

△y

-1 3

△ y

-2 5

△ y

-3 7

△ y

-4

0 2

△ y

-1 4

△ y

-2 6

△ y

-3 8

△ y

-4

△y

0 3

△ y

-1 5

△ y

-2 7

△ y

-3

1 2

△ y

0 4

△ y

-1 6

△ y

-2 8

△ y

-3

△y

1 3

△ y

0 5

△ y

-1 7

△ y

-2

2 2

△ y

1 4

△ y

0 6

△ y

-1

△y

2 3

△ y

1 5

△ y

3 2

△ y

2 4

△ y

△y

3 3

△ y

4 2

△ y

△y

2 3 4 5

y=y +C △y +C △ y +C △ y +C △ y +C △ y +... (B)

0 1 0 2 0 3 0 4 0 5 0

其中各个C表示各个二项式系数。

注：有用符号δ，δ，...来记中心差分的，为了便于比较起见，列出下表：

阶

下标 1 2 3 4

-2 y

-2

△y =δy =△ y

-2 -3/2 1 1

-1 y

-1 2 2

△y =δ y =△ y

-2 -3/2 3 1

△y =δy =△ y

-2 -1/2 1 0 3 3

△ y =δ y =△ y

-2 -1/2 2 1

0 y

0 2 2

△y =δ y =△ y

-1 0 2 1 4 4

△ y =δ y =△ y

-2 -1/2 2 1

△y =δy =△ y

0 1/2 1 1 3 3

△ y =δ y =△ y

-1 1/2 3 2

1 y

1 2 2

△y =δ y =△ y

0 1 2 2

△y =δy =△ y

0 3/2 1 1

2 y

各种记法间的关系显见是：

m m

△ y =δ y =△ y

k k+m/2 m k+m

即三种差分记法的下标的关系是，对角差分下标加阶数一半，得中心差分下标，再加阶数一半，得水平差分下标。现在我们命

3 3

△y +△y △ y +△ y

-1 0 -2 -1

(a)m = , (b)m =

1 2 3 2

5 5 7 7

△ y +△ y △ y +△ y

-3 -2 -4 -3

(c)m = , (d)m =

5 2 7 2

…………………………….

所以这些m就表示：

跟通过y 所作水平线紧挨着的一些奇阶差分的算术平均值。

2 3

我们目前的目标，是把△y ,△ y ,△ y ,等等，

0 0 0

要用各个m以及位于通过y 水平线上的各个偶阶差分表示出来，

这件事可采用消去法来做，即从△y ,△ y 等等起顺着对角线向上的进行到右边去，

0 0

直至到达水平线上的各量时为止。为了对这件事的进行有所帮助，

2 4 6 8

对于偶阶差分△ y ,△ y ,△ y ,△ y 等等，

-1 -2 -3 -4

不论在以下任何代数运算的地方出现，都在它们的下面画线，画线的目的是为了提醒我们下面画线的量是不消去的。从差分的定义，我们可得

△ y =△y -△y

-1 0 -1

所以,

△y =△ y -△y

0 -1 -1

但从(a)可知

△y =2m -△y

-1 1 0

所以,

△y =△ y +2m -△y

0 -1 1 0

所以,

△y =m +△ y /2 (1)

0 1 -1

要想求出所需量表示的△ y 的值，则有

注：读者可能认为这些解法过程相当麻烦，

也不容易看出与上文“从△y ,△ y 等等起顺着对角线向上的进行到右边去，

0 0

直至到达水平线上的各量时为止”一句话的关系。

2 3

现在用图解来说明△ y 与△ y 的求法如下：

0 0

（到水平线上为止）

△ y

-2

2 4

△ y △ y

-1 -2

△ y

-1

△ y

沿着对角线向上

3 5

△ y △ y

-2 -3

4 6

△ y △ y

-2 -3

3 5

△ y △ y

-1 -2

△ y

-1

△ y 3

0 从△ y 沿着对角线向上

这里有两个三数关系，即

3 2 3 3 3 4

△ y -△ y =△ y ，△ y -△ y =△ y ，

0 -1 -1 -1 -2 -2

一个二数关系，即

3 3

2m =△ y +△ y

3 -1 -1

3 3 3

从这三个方程中，消去△ y 和△ y ，便得到了△ y

-1 -2 0

与上相同，可知道这里有四个三数关系，

3 3 4 5 5

两个二数关系从这六个方程中消去△ y ,△ y ,△ y ,△ y ,△ y ,五个量，

-2 -1 -1 -3 -2

便得到了△ y 。这种图解说明，不难使它一般化。

2n 2

在求△ y 时，有n +n个三数关系n个二数关系。

2n 2

图解中除去△ y 外，有n +3n个量，但有n+1个是消不去的。

2 2 3 2

因此，恰好从 n +n+n=n +2n个关系中，消去n +3n-(n+1）=n +2n-1个量，

2n+1 2

在求△ y 时，有(n+1) 个三数关系，n+1个二数关系。

2n+1 2

除去△ y 外，n +4n+2个量，但n+1个是不消去的。

2 2 2 2

因此，恰好从(n+1) +(n+1)=n +3n+2个关系中，消去n +4n+2-(n+1)=n +3n+1个量。

注释结束。

要想求出所需量表示的△ y 的值，则有

3 2 2

△ y =△ y -△ y

-1 0 -1

或

2 2 2

(e)△ y =△ y -△ y

0 -1 -1

然而从定义已知

4 3 3

(f)△ y =△ y -△ y

-2 -1 -2

并且由(b)已知

3 3

(g)△ y =2m -△ y

-1 3 -2

故自(g)减去(f)：并解出△ y 可得

-1

3 4

(h)△ y =m +△ y /2

-1 3 -2

将(h)代入(e)则得

2 2 4

△ y =△ y +m +△ y /2 （2）

0 -1 3 -2

要找出△ y 我们从下式着手（见上页注）：

4 3 3

△ y =△ y -△ y

-1 0 -1

或

3 3 4

(i)△ y =△ y -△ y

0 -1 -1

4 4

=m +△ y /2+△ y

3 -2 -1

(自(h)得出

然而

5 4 4

△ y =△ y -△ y

-2 -1 -2

4 4 5

(j)△ y =△ y -△ y

-1 -2 -2

又

6 5 5

(k)△ y =△ y -△ y

-3 -2 -2

并且自(c)已知

3 4

(L)△ y =m +△ y

-1 3 -2

自(L)减去(k)并解出△ y .

-2

5 6

(m)△ y =m +△ y /2

-2 2 -3

将(m)代入(j)得

4 4 6

(n)△ y =△ y +m +△ y /2

-1 -2 5 -3

将(n)代入(i)得

3 4 6

△ y =m +3△ y /2+m +△ y /2 (3)

0 3 -2 5 -3

要找出△ y 我们从下式着手：

5 4 4

△ y =△ y -△ y

-1 0 -1

或

4 4 5

(o) △ y =△ y +△ y

0 -1 -1

4 6 5

=△ y +m +△ y /2+△ y

-2 5 -3 -1

(自(n)得出

然而

6 5 5

(p) △ y =△ y -△ y

-2 -1 -2

5 6

=△ y +m +△ y /2

-1 5 -3

得自(m), 并且

7 6 6

(q) △ y =△ y -△ y

-3 -2 -3

又

8 7 7

(r) △ y =△ y -△ y

-4 -3 -4

并且，自(d)可知

7 7

(s) △ y =2m -△ y

-3 7 -4

自(s)减去(r)，并解出△ y

-3

7 8

(t) △ y =m -△ y /2

-3 7 -4

将(t)代入(q)，并解出△ y

-2

6 6 8

(u) △ y =△ y +m -△ y /2

-2 -3 7 -4

将(u)代入(p)，并解出△ y

-1

5 6 8

(v) △ y =m -3△ y /2+m -△ y /2

-1 5 -3 7 -4

将(v)代入(o)得

4 4 5 8

△ y =△ y +2m +2△ y +m -△ y /2

0 -2 5 -3 7 -4

现在把(1),(2),(3),(4)代入(B)，我们便得

2 2 4

y=y +C (m +△ y /2)+C (m +△ y +△ y /2)+

0 1 1 -1 2 2 -1 2

4 6

+C (m +m +3△ y /2+△ y /2)+

3 3 5 -2 -3

4 6 8

+C (2m +m +△ y +2△ y +△ y /2)

4 5 7 -2 -3 -4

或

y=y +C m +(C /2+C )△ y +(C +C )m +

0 1 1 1 2 -1 2 3 3

4 5 6

+(C /2+3C /2+C )△ y +m , △ y 等等各项,

2 3 4 -2 -3

将各个C及m用它们的值置换之，即得

△y +△y

-1 0 u u(u-1) 2

y=y +u +( + )△ y +

0 2 2 2 -1

3 3

△y +△y

u(u-1) u(u-1)(u-2) -2 -1

+( + ) +

2 6 2

u(u-1) 3u(u-1)(u-2) 3u(u-1)(u-2)(u-3) 4

+( + + )△ y + ……

2 12 2 -2

或

3 3

△y +△y 2 2 △ y +△ y

-1 0 u 2 u(u -1) -2 -1

y=y +u + △ y + +

0 2 2 -1 3! 2

2 2

u (u -1) 4

△ y +...

4! -2

以上面所举的计算纲领，继续进行，我们就得到斯特灵公式，即

3 3

△y +△y 2 2 △ y +△ y

-1 0 u 2 u(u -1) -2 -1

y=y +u + △ y + +

0 2 2 -1 3! 2

5 5

2 2 2 2 2 2 △ y +△ y

u (u -1) 4 u(u -1 )(u -2 ) -3 -2

△ y + +

4! -2 5! 2

2n-1 2n-1

2 2 2 2 2 2 2 2 △ y +△ y

u(u -1 )(u -2 )(u -3 )...[u -(n-1) ] -n -(n-1)

(2n-1)! 2

2 2 2 2 2 2 2 2

u(u -1 )(u -2 )(u -3 )...[u -(n-1) ] 2n

△ y (Ⅲ)

(2n)! -n

其中,

x-x

在这个公式当中一共有2n+1项，并且多项式与给定函数在下列2n+1个点上彼此相重合：

u=-n,-(n-1),-(m-2),...,-2,-1,0,1,2,...,n-2,n-1,n;

或

x=x -nh,x -(n-1)h,...,x -h,x ,x +h,...,x +(n-1)h x +nh

0 0 0 0 0 0 i 0

24.贝赛尔插值公式

贝赛尔公式的推导与斯特灵公式的推导相似。首先我们和前面一样写出对角差分表，然后，在y 和y 中央作一水平线把紧挨着水平线的各个量都标记出来，以便加以特别注意。

0 1

这些量在下表以黑体字印出。

表8

△y 2

△ y 3

△ y 4

△ y 5

△ y 6

△ y 7

△ y 8

△ y

-4

△y

-4

-3 2

△ y

-4

△y

-3 3

△ y

-4

-2 2

△ y

-3 4

△ y

-4

△y

-2 3

△ y

-3 5

△ y

-4

-1 2

△ y

-2 4

△ y

-3 6

△ y

-4

△y

-1 3

△ y

-2 5

△ y

-3 7

△ y

-4

0 2

△ y

-1 4

△ y

-2 6

△ y

-3 8

△ y

-4

△y

0 3

△ y

-1 5

△ y

-2 7

△ y

-3

1 2

△ y

0 4

△ y

-1 6

△ y

-2 8

△ y

-3

△y

1 3

△ y

0 5

△ y

-1 7

△ y

-2

2 2

△ y

1 4

△ y

0 6

△ y

-1

△y

2 3

△ y

1 5

△ y

3 2

△ y

2 4

△ y

△y

3 3

△ y

4 2

△ y

△y

现在我们令

3 3

△y +△y △ y +△ y

0 1 -1 0

(a)m = (b)m =

0 2 2 2

4 4 6 6

△ y +△ y △ y +△ y

-3 -1 -3 -2

(c)m = (d)m =

4 2 6 2

8 8

△ y +△ y

-4 -3

(e)m = +…

7 2

在此情况下，各个n就分别表示：

纵标y 和y 的算术平均值以及紧挨着普通y 的水平线上下的各个偶阶差分的算术平

0 i

均值。其次，我们像在23节做过的那样，写出y 为始点的牛顿公式(Ⅰ)。

2 n

我们的问题是要把y ,△y ,△ y ,...,△ y 等用各个m以及位于通道y 的水平线

0 0 0 0 1/2

上的各个奇阶差分表示出来。这件事可以用消去法的过程来做，

2 3

即自△ y ,△ y 等等起顺着对角线向上进行到右边去，

0 0

直至达到水平线的各量为止。在以下的运算中，对于水平线上的各个奇阶差分都在下面画线，以表示它们是不消去的。由定义可得

△y=y -y

1 0

所以,

y =y -△y

0 1 0

然而自(a)

y =2m -△y

所以,

y =2m -y -△y

0 0 0 0

y =m -△y /2 (1)

0 0 0

欲求得△ y 我们从下式着手：

由定义可知

3 2 2

△ y =△ y -△ y

-1 0 -1

所以,

2 3 2

△ y =△ y -△ y

-1 -1 -1

然而由(b),

2 2

△ y =2m -△ y

-1 2 0

所以,

2 3 3

△ y =△ y +2m -△ y

0 -1 2 0

或

2 2

△ y =m +△ y /2 (2)

0 2 -1

欲求△ y 则自定义得

4 3 3

△ y =△ y -△ y

-1 0 -1

所以,

3 3 4

(f) △ y =△ y +△ y

0 -1 -1

然而，从(c)有

4 4

(g) △ y =2m +△ y

-1 3 -2

并且由定义

5 4 4

(h) △ y =△ y -△ y

-2 -1 -2

自(g)减去(h)并解出△ y

-1

4 5

(i) △ y =m +△ y /2

-1 4 -2

将(i)代入(f)则得

3 3 5

△ y =△ y +m +△ y /2 (3)

0 -1 4 -2

欲求△ y 则从下式着手：

5 4 4

△ y =△ y -△ y （得自定义）

-1 0 -1

所以,

4 4 5

(j) △ y =△ y +△ y

0 -1 -1

5 5

=m +△ y /2+△ y （得自(i)）

4 -2 -1

现在由定义,

6 5 5

(k) △ y =△ y -△ y

-2 -1 -2

又由定义,

7 6 6

(L) △ y =△ y -△ y

-3 -2 -3

并且自(d)得

6 6

(m) △ y =2m -△ y

-2 6 -3

自(m)减去(L)并解出△ y

-2

6 7

(n) △ y =m +△ y /2

-2 6 -3

令(k)等于(n)并解出△ y

-1

5 5 7

(o) △ y =△ y +m +△ y /2

-1 -2 6 -3

将(o)代入(j)可得

4 5 7

△ y =m +△ y /2+m +△ y /2 (4)

0 4 -2 6 -3

2 3

现在把这些y ,△ y ,△ y ,等等的值代入23节(A)式中，即得

0 0 0

u(u-1) u(u-1) u(u-1)(u-2) 3

y=m -△y /2+u△y + m +[ + ]△ y +

0 0 0 2 2 4 6 -1

u(u-1)(u-2) u(u-1)(u-2)(u-3) u(u-1)(u-2) u(u-1)(u-2)(u-3) 5

+[ + ]m×[ + ]△ y +

6 24 12 16 -2

5 7

+△ y,m 及△ y 各项。

6 -3

简化之，并将各个m用它们的值代换则得

2 2

y +y △ y +△ y

0 1 u(u-1) -1 0

y= +(u-1/2)△y + +

2 0 2 2

2 2

△ y +△ y

u(u-1)(u-1/2) 3 u(u-1)(u+1)(u-2) -2 -1

△ y + +...

2 -1 4! 2

现在因为△y =y -y ,故两项可以化作y +u△y ，因此得

0 1 0 0 0

2 2

△ y +△ y

u(u-1) -1 0 u(u-1)(u-1/2) 3

y=y +u△y + + △ y +

0 0 2 2 3! -1

4 4

△ y +△ y

u(u-1)(u+1)(u-2) -2 -1

+...

4! 2

继续进行像上面所作的那样计算，我们就会获得贝赛尔插值公式：

2 2

△ y +△ y

u(u-1) -1 0 u(u-1)(u-1/2) 3

y=y +u△y + + △ y +

0 0 2 2 3! -1

4 4

△ y +△ y

u(u-1)(u+1)(u-2) -2 -1 u(u-1)(u-1/2)(u+1)(u-2) 5

+ △ y +

4! 2 5! -2

u(u-1)(u-1/2)(u+1)(u-2) 5

+ △ y +

5! -2

6 6

△ y +△ y

u(u-1)(u+1)(u-2)(u+2)(u-3) -3 -3

+ +…+

6! 2

2n 2n

△ y +△ y

u(u-1)(u+1)(u-2)(u+2)...(u-n)(u+n-1) -n -n

+ +

(2n)! 2

u(u-1)(u-1/2)(u+1)(u-2)(u+2)...(u-n)(u+n-1) 2n+1

+ △ y (Ⅳ)

(2n+1)! -n

如果在这个公式中，我们命u=1/2,则可得简单公式如下：

2 2 4 4

y +y △ y +△ y △ y +△ y

0 1 1 -1 0 3 -2 -1

y= - +

2 8 2 128 2

6 6

△ y +△ y

5 -3 -2

- +…+

1024 2

2n 2n

2 △ y +△ y

[1*3*5...(2n-1)] -n -n+1

+(-1) (Ⅴ)

2 2

2 (2n)！

贝赛尔公式的这个重要特殊情况，叫做半数插值公式。它可用来计算函数在任何两个给定值中间的值。若令u-1/2=v或u=v+1/2便可获得在形式上更对称更方便的贝赛尔公式。在公示（Ⅳ）进行这种代换则得

2 2

y +y 2 △ y +△ y 2

0 1 (v -1/4) -1 0 v(v -1/4) 3

y= +v△y + + △ y +

2 0 2 2 3! -1

4 4 2 2

2 2 △ y +△ y v(v -1/4)(v -9/4)

(v -1/4)(v -9/4) -2 -1 5

+ + △ y +

2 2 5! -2

4 4

2 2 △ y +△ y

(v -1/4)(v -9/4)(v -25/4) -2 -1

+ +…+

6! 2

2 2 2 (2n-1) 2n 2n

(v -1/4)(v -9/4)...[v - ] △ y +△ y

2 -n -n+1

+ +

6! 2

2 2 2 (2n-1)

(v -1/4)(v -9/4)...[v - ]

2 2n+1

+ △ y (Ⅵ)

6! -n

在公式(Ⅳ)与(Ⅵ)中，共有2n+2项，而由它们表达出来的多项式和给定函数在2n+2个点上，即在点

u=-n,-n+1,-n+2,...,-1,0,1,2,...,n,n+1;,

2n+1 2n+1 2n+1 3 1 1 3 2n-1 2n-1

v=- ,- ,... ,- ,... ,- ,- , , ,…, ,

2 2 2 2 2 2 2 2 2

x=x -nh,x -(n-1)h,...,x -h,x ,x +h,...,x +nh,x +(n+1)h

0 0 0 0 0 0 0

上彼此重合一致。其中v的零点是x +h/2，而u的是x 。

0 0

我们现在将把斯特灵和贝赛尔公式应用到一些数值实例上去。

例1.表A给出当x取某些等距值时，概率积分

2 x -x

f(x)= ∫ e dx

√π 0

的各值。求积分在x=0.5437时的值。

解：这里我们取x =0.54,x=0.5437。由于h=0.01,故

x-x

0 0.5437-0.54 0.0037

u= = = =0.37

h 0.01 0.01

表A

x f(x)

△f(x) 2

△ f(x) 3

△ f(x) 4

△ f(x)

0.51 0.5292437

0.008655

0.52 0.5378987 -0.0000896

0.0085654 -0.0000007

0.53 0.5464641 -0.0000902 0

0.0084751 -0.0000007

0.54 0.5549392 -0.0000910 0

0.0083841 -0.0000007

0.55 0.5633233 -0.0000917 1

0.0082924 -0.0000006

0.56 0.5716157 -0.0000923

0.0082001

0.57 0.5798158

(a)使用斯特灵公式（Ⅲ）则得

(0.008655+0.0083841） (0.37)

f(0.5437)=0.5549392+0.37 + (-910)+

2 2

0.37(0.37*0.37-1) (-7-7)

2 6

=0.5549392+0.00311895-0.00000623+0.00000004=0.5580520

(b)要用贝赛尔公式求f(0.5437),则使用(Ⅵ)要方便一些。这里

v=u-1/2=0.37-0.50=-0.13

代入(Ⅵ)，则得

0.5549392+0.5633233

f(0.5437)= +(-0.13)(0.0083841)+

0.0169-0.25 (-0.000091-0.0000917) -0.13(0.0169-0.25)(-7)

+ +

2 2 6

=0.55913125-0.08108993+0.00001065=0.5580520

-x

例2.在x取等距值时，e 的各值给出如表B，

-x

求x=1.7489时e 之值。

解：(a)使用斯特灵公式, 这里我们取x=1.7489,x =1.76,h=0.01,则

1.7489-1.75 0.0011

u= =- =-0.11

0.01 0.01

表B

x -x

△f(x) 2

△ f(x) 3

△ f(x) 4

△ f(x)

1.72 0.1790661479

-0.0017817379

1.73 0.1772844100 0.0000177285

-0.0017640094 -0.0000001762

1.74 0.1755204006 0.0000175523 0.0000000013

-0.0017464571 -0.0000001749

1.75 0.1737739435 0.0000173774 0.0000000022

-0.0017290792 -0.0000001727

1.76 0.1720448638 0.0000172047 0.0000000015

-0.0017118750 -0.0000001712

1.77 0.1703329888 0.0000170335

-0.0016948415

1.78 0.1686381473

代入（Ⅲ）便得

(-0.0017464571-0.0017290797)

f(1.7489)=0.1737739435-0.11 +

0.0121 (0.0121-1) (0.0000001749-0.0000001727)

+ (0.0000173774)-0.11 +

2 6 2

(0.0121-1)

+0.0121 (22)

=0.1737739435+0.00019115452+0.00000010513-0.00000000315

或,

-1.7489

f(1.7489)=e =0.1739652000

该值准确到十位小数。

(b)使用贝赛尔公式。由于数值1.7489在区间1.74-1.75当中要比它在区间1.75-1.76当中更靠近中央，所以我们选取x=1.74以便能使v值可能的小。因此可得

1.7489-1.74

u= =0.89

0.01

v=u-1/2=0.89-0.50=0.39,

所以,

0.1755204006+0.1737739435

f(1.7489)= +0.39(-0.0017464571)+

0.39*0.39-0.25 0.0000175523+0.0000173774

+ ( )( ) +

2 2

(0.39*0.39-0.25)

+(0.39) (0.0000001749)+

(0.39*0.39-0.25)(0.39*0.39-2.25) (13+22)

24 2

=0.17464717205-0.00068111827-0.00000085490+0.00000000111+0.00000000001

或, f(1.7489)=0.1739652000,

与上面所得相同。我们也可以取x =1.75,此时v=-0.61,这时可得

这个数值也准确到十位小数，不过这个级数要比前面的情况收敛得稍慢一些；从贝赛尔公式所得出的两个级数都比斯特灵公式所得出的这个级数要收敛得稍慢一些。

注意，在这里自然而然的会产生这样一个问题：

答案是：它们的准确度不相上下，对于给定的差分表来说，

收敛的迅速与否同公式(Ⅲ)中u的大小，以及公式(Ⅵ)中v的大小有关。u和v的值愈小，则级数收敛得愈快。所以我们应该经常选择始点x 以便能促使u和v尽可能地小。

在多数情况下，选择起始点使得-0.5≤u≤0.5以及-0.5≤v≤0.5是可能的，譬如在例题1中起始点是如此选定的：u=0.37,v=-0.13, 而在例题2中u=-0.11,v=0.39。可注意的是，在第一个问题中贝赛尔公式收敛得较快而在第二例题中斯特灵公式收敛的较快；它的原因是在第一种情形中v比u要小，而在第二情形中u比v要小。一般的规则可以这样说：在靠近区间的中央部分进行插值时，比如说u=0.25到0.75（v=-0.25到0.25), 贝赛尔公式能得出较准的结果：但在靠近区间的开端和末端部分进行插值时，比如说u=-0.25到0.25，斯特灵公式却能给出较好的结果。这个问题的其它方面可见第五章。

例3.在φ取某些等距值时，椭圆积分

dφ

F(φ)= ∫

1-sin φ/2

表C

2 3 4

φ F(φ) △F △ F △ F △ F

21° 0.370634373

0.018070778

22° 0.388705151 0.000059002

0.018129780 0.000002707

23° 0.406834931 0.000061709 0.000000004

0.018191489 0.000002711

24° 0.425026420 0.000064420 -0.000000007

0.018255909 0.000002704

25° 0.443282329 0.000067124

0.018323033

26° 0.461605362

的值给出如下表，求F(22.5°)的值。

解：因为我们所求的函数值介乎两个给定列表值正中央，所以我们使用半数值公式（Ⅴ）。因此可得

2 2 4 4

y +y 1 △ y +△ y 3 △ y +△ y

y= 0 1 - -1 0 + -2 -1

2 8 2 128 2

2 2n 2n

5 △ y +△ y [1*3*5...(2n-1)] △ y +△ y

- -3 -2 +…+(-1) -n -n+1 (Ⅴ)

1024 2 2n 2

2 (2n)！

贝赛尔公式的这个重要特殊情况，叫做半数插值公式。

0.406834931+0.425026420 1 0.000061709+0.000064420

F(22.5°)= -

2 8 2

3 0.00000000004-0.00000000007

128 2

=0.1459306755-0.0000078831=0.415922792

由于表中的差分是完全有规则的并且很快的减小，所以这个结果或许可以准确到末一位数字。

第七部分拉格朗日公式，反插值法

下面的资料可参见《数值分析》，美国J.B.斯卡勃罗著，陈盖民，余介石译，科学出版社1958年出版。

拉格朗日公式，反插值法

Ⅰ.插值法的拉格朗日公式

25.引言。

前一章推导出的插值公式，只在自变量取等间距值时才可应用。但有时要从自变量取等间距值来求得函数值，这是不方便的甚至是不可能的。遇到这种情况时，就希望有一个插值公式，式中仅包含在手边能有的那些数据。现在我们就要推导这样的一个公式。

26.拉格朗日公式

令(x ,y ),(x ,y ),(x ,y ),...,(x ,y )表示y=f(x)中，】

0 0 1 1 2 2 n n

变量x和y相对应的n+1对值。我们用下列形式的n次多项式

φ(x)=A (x-x )(x-x )(x-x )...(x-x )+

0 1 2 3 n

+A (x-x )(x-x )(x-x )...(x-x )+

1 0 2 3 n

+A (x-x )(x-x )(x-x )...(x-x )+

2 0 1 3 n

+…………………………+

+A (x-x )(x-x )(x-x )...(x-x ) (1)

n 0 1 2 n-1

来代换给定函数。多项式中一共有n+1项，而每一项中有n个因子。

其次，我们定出A ,A ,A ,...,A 等n+1个常数，

0 1 2 n

以便能使φ(x )=y ,φ(x )=y ,...,φ(x )=y 。

0 0 1 1 2 2

在(1)式中令x=x ,φ(x )=y 则得

0 0 0

y =A (x -x )(x -x )...(x -x )

0 0 1 0 1 2 1 n

所以,

A =

1 (x -x )(x -x )...(x -x )

1 0 1 2 1 n

同样我们可以求得

A =

2 (x -x ) (x -x )(x -x )...(x -x )

2 0 2 1 2 3 2 n

……………………………

A =

n (x -x )(x -x )...(x -x )

n 0 n 2 n n-1

把这些A值代入(1)，我们可得

(x-x )(x-x )...(x-x )

0 2 n

φ(x)= y +

(x -x )(x -x )...(x -x ) 0

0 1 0 2 0 n

(x-x )(x-x )...(x-x )

0 2 n

+ y +

(x -x )(x -x )...(x -x ) 1

1 0 1 2 1 n

(x-x ) (x-x )(x-x )...(x-x )

0 1 3 n

+ y +…+

(x -x ) (x -x )(x -x )...(x -x ) 2

2 0 2 1 2 3 2 n

(x-x )(x-x )...(x-x )

0 1 n-1

+ y (Ⅶ)

(x -x )(x -x )...(x -x ) n

n 0 n 1 n n-1

这个公式也可以写成下列形式

∏ (x)

n n

φ(x)= ∑ y ,i=0,1,...,n (Ⅷ)

i=0 (x-x )∏` (x) i

i n

注：原书误作∏` (x)，

但分母显然是∏` (x )=(x -x )(x -x )...(x -x )(x -x )...(x -x )

n i i 0 i 1 I i-1 I i+1 I n

其中,

∏ (x)=(x-x )(x-x )...(x-x )

n 0 1 n

∏` (x)=d∏ (x)/dx

n n

公式（Ⅶ）和（Ⅷ）称为拉格朗日插值公式，其自变量取等间距值或不取等间距值都可以。这里可注意的是，拉格朗日公式并不含有有关函数的逐次差分，并且其中没有什么可供我们去估计所得结果可靠性的东西。由于拉格朗日公式只是两个变量间的关系，其中任一个均可当作自变量，所以很显然，若把y看作自变量，则我们便可写出x是y的函数的公式来。因此在（Ⅶ）的右端将x同y互换一下，便得

(y-y )(y-y )...(y-y )

1 2 n

φ(y)= x +

(y -y )(y -y )...(y -y ) 0

0 1 0 2 0 n

(y-y )(y-y )...(y-y )

0 2 n

+ x +

(y -y )(y -y )...(y -y ) 1

1 0 1 2 1 n

(y-y )(y-y )...(y-y )

0 1 n

+ x +…+

(y -y )(y -y )...(y -y ) 2

2 0 2 1 2 n

(y-y )(y-y )...(y-y )

0 1 n-1

+ x (Ⅸ)

(y -y )(y -y )...(y -y ) n

n 0 n 1 n n-1

拉格朗日公式的主要用途有二：

（1）当函数自变量的给定值不是等间距时，求出函数的任何值。

（2）求出与给定函数值相对应的自变量的值。第二个问题可用公式（Ⅸ）解出。

现在我们举两个例题来说明这些用途。

例1.下表给出x和lgx的某些对应值。

试计算lg323.5的值。

x 321.0 322.8 324.2 325.0

lgx 2.50651 2.50893 2.51081 2.51188

解：这里x=323.5,x =321.0,x =322.8,x =324.2,x =325.0

0 1 2 3

把这些值代入（Ⅶ）则得

(323.5-322.8)(323.5-324.2)(323.5-325.0)

lg323.5= *2.50651+

(321-322.8)(321-324.2)(321-325)

(323.5-321)(323.5-324.2)(323.5-325)

*2.50893+

(321-321)(321-324.2)(322.8-325)

(323.5-321)(323.5-322.8)(323.5-325)

*2.51081+

(324.2-321)(324.2-322.8)(324.2-325)

(323.5-321)(323.5-322.8)(323.5-324.2)

*2.51188

(325-321)(325-322.8)(325-324.2)

=-0.07996+1.18794+1.83897-0.43708=2.50987

这个结果准确到末尾数字。

例2.下表给出与某些x值相对应的概率积分的值，问x为何值时等于该积分等于1/2?

2 x -x

∫ e dx x

√π 0

0.4846555 0.46

0.4937452 0.47

0.5027498 0.48

0.5116683 0.49

解：命概率积分的值为y，则

y=1/2=0.5,x =0.46,x =0.47,x =0.48,x =0.49

0 1 2 3

将这些值代入（Ⅸ），则得

(0.5-0.4937452)(0.5-0.5027498)(0.5-0.5116683)

x= *0.46+

(0.4846555-0.4937452)(0.4846555-0.5027498)(0.4846555-0.5116683)

(0.5-0.4846555)(0.5-0.5027498)(0.5-0.5116683)

*0.47+

(0.4937452-0.4846555)(0.4937452-0.5027498)(0.4937452-0.5116683)

(0.5-0.4846555)(0.5-0.4937452)(0.5-0.5116683)

*0.48+

(0.4937452-0.4846555)(0.5027498-0.4937452)(0.5027498-0.5116683)

(0.5-0.4846555)(0.5-0.4937452)(0.5-0.5027498)

*0.49

(0.5116683-0.4846555)(0.5116683-0.4937452)(0.5116683-0.5027498)

62548*27498*116683 153445*27498*116683

=- *0.46+ *0.47+

90897*180943*270128 90897*90046*179231

153445*62548*116683 153445*62548*27498

+ *0.48+ *0.49

180943*90046*89185 270128*179231*89185

=-0.0207787+0.157737+0.369928-0.0299495=0.476937

准确到六位小数的真值为0.476936, 注：这个问题的计算除非有计算机可用，否则应该用对数来做。

注：这个问题的计算除非有计算机可用，否则应该用对数来做。

注意：读者通过前面二个例题的计算，将会觉得拉格朗日公式应用起来很麻烦，并且含有很多计算。使用它还必须十分小心与慎重，因为假如自变量的值取得不密集的话，则易使其结果很不准确。由于这个缘故，除了牛顿、斯特灵和贝塞尔公式不能应用的情况就不应该使用拉格朗日公式。

Ⅱ.反插值法

27.定义

如果一个给定的函数值介乎两个列表值之间，那么，要去求出它的对应自变量值的这种过程就是反插值法。反插值法问题有好几种解法，但在本书中我们只讲三种。

28.拉格朗日公式法

处理这个问题的方法之一，是使用具有像拉格朗日公式（Ⅸ）那种的式子，而式中的x，系当作y的函数来表出的。前节的例2实在就是反插值法的问题，所以我们不再解释这个方法。

29.逐步求近法

第二种方法是逐步求近法或叠代法。要明白如何用这个方法，我们可考虑一下牛顿公式（Ⅰ），即

u(u-1) 2 u(u-1)(u-2) 2

y=y +u△y + △ y + △ y +

0 0 2 0 3! 0

u(u-1)(u-2)(u-3) 4 u(u-1)(u-2)...(u-n+1) n

△ y +…+ △ y (1)

4! 0 n! 0

移项并用△y 逼除，则得

2 3 4

y-y u(u-1)△ y u(u-1)(u-2) △ y u(u-1)(u-2)(u-3) △ y

0 0 0 0

u= - - - (1)

△y 2△y 3! △y 4! △y

0 0 0 0

要得到u的一次近似值，我们略去所有高于一阶的差分，即得

y-y

[1] 0

u =

△y

[1]

将u 代入（1）的右端则得二次近似值，于是可得

2 2

y-y [1] [1] △ y [1] [1] [1] △ y

[2] 0 u (u -1) 0 u (u -1)(u -2) 0

u = - -

△y 2 △y 3! △y

0 0 0

[1] [1] [1] [1] △ y

u (u -1)(u -2)(u -3) 0

- (2)

4! △y

其三次近似值为

2 3

y-y [2] [2] △ y [2] [2] [2] △ y

[3] 0 u (u -1) 0 u (u -1)(u -2) 0

u = - -

△y 2 △y 3! △y

0 0 0

[2] [2] [2] [2] △ y

u (u -1)(u -2)(u -3) 0

- (3)

4! △y

依次类推以至更高次的近似值。现在我们举例说明这个方法。

例1.现在已经给出概率积分

2 x -x

∫ e dx x

√π 0

的列表值（表见下页），

解：这里最好使用中心差分公式。

2 3 4

x y △y △ y △ y △ y

0.45 0.4754818

0.0091737

0.46 0.4846555 -0.000084

0.0090897 -0.0000011

0.47 0.4937452 -0.0000851 0.0000001

0.0090046 -0.000001

0.48 0.5027498 -0.0000861 0.0000002

0.0089185 -0.0000008

0.49 0.5116683 -0.0000869

0.0088316

0.50 0.524999

由观察看出，所求之x值，在0.47与0.48之间，由粗略的线性插值法可知它约为0.47(2/3).. 因此我们取x =0.47并使用贝塞尔公式。所以我们可得

x =0.47,h=0.01,y=1/2=0.5

将这个y值以及表中查到的适应量代入贝塞尔公式（Ⅵ），则得

(v -0.25) v(v -0.25)

0.5=0.4982475+0.0090046v+ (-0.0000856)+ (-0.0000010)

2 6

移项并用0.0090046逼除，则得

2 2

v=0.194623-(v -0.25)(-0.004753)-v(v -0.25)(-0.00000185) (4)

将(4)式右端所有高于一阶的各项略去，则得v的一次近似值，即

[1]

v =0.194623,

将此v值代入(4)式右端，我们便求出二次近似值为

[2] 2 2

v =0.194623-[(0.194623) -0.25](-0.004753)-0.194623[(0.194623) -0.25](-0.0000185)

=0.194623-0.001008-0.000001=0.193612

现在将此v值代入（4）式右端，我们可得

[3]

v =0.194623-0.0010101-0.000001=0.193612

此值与上值相差不大，所以我们不必再求更高次的近似值。由于,

u=v+1/2,x=x +hu,故得

u=0.693612,

x=0.47+0.01(0.693612)=0.47693612

此值准确到六位小数。

注：在此例题中，v值不可能获得五位以上的可靠值，因为(4)式右端系使用近似数0.0090046除得之结果，它的第五位数字是不可靠的。事实上v值仅在头四位数字准确。假如所有高于二阶的差分均被略去，那么反插法问题只是一个解二次方程问题而已，下例就说明这一点。

例2.给定sinhx=62,求x

解：作出差分表如下，我们发现所有高于二阶的差分均为零。我们还可以看出，所求之x值比4.82要稍微大一些，因此我们取x=0.482,并使用斯特灵公式。

x y=sinhx

△y 2

△ y 3

△ y

4.80 60.7511

0.6106

4.81 0.0062

0.6168 0

4.82 61.9785 0.0062

0.623 0

4.83 62.6015 0.0062

0.6292

4.84 63.2307

把y=62代入斯特灵公式（Ⅲ），则得

62=61.9785+0.6199u+0.0031u

或,

31u +6199u=215

所以,

-6199+ (6199) +4*31*215

-6199+6201.15 2.15

= = =0.0347

62 62

因为h=0.01,而x=x +hu,故得

x=4.82+0.01(0.0347)=4.8203

第八部分克莱姆法则

下面的资料可参见《工程数学》，林益编，高等教育出版社2003年出版

1.3.1行列式的概念

a a

设A 11 12 是二阶方阵，

a a

21 22

称a a -a a 为A的行列式或二阶行列式，记为│A│或

11 12 21 22

a a

11 12

a a

21 22

a a

11 12

│A│= =a a -a a

a a 11 12 21 22

21 22

利用二阶行列式的概念，二元线性方程组

a x +a x =b

11 1 12 2 1

{

a x +a x =b

21 1 22 2 2

的解可以简化的表示成

D D

1 2

x = .x =

1 D 2 D

其中,

a a

11 12

a a

21 22

b a

1 12

D =

1 b a

2 22

a b

11 1

D =

2 a b

21 2

这里，自然要求D≠0，D称为方程组系数矩阵的行列式，简称为系数行列式。

D ,D 分别是将系数行列式D中的第1，2列对应的换为方程组的常数项而得的二阶行

1 2

列式。若A是一个三阶行列式，即

a a a

11 12 13

a a a

AD= 21 22 23

a a a

31 32 33

a b a a a a

11 1 21 23 21 22

│A│=a -a +a

11 a b 12 a a 13 a a

21 2 31 33 31 32

=a (a a -a a )-a (a a -a a )+a (a a -a a )

11 22 33 23 32 12 21 33 23 31 13 21 32 22 31

=a a a +a a a +a a a -a a a -a a a -a a a

11 22 33 12 23 31 13 32 21 13 22 31 12 21 33 11 32 23

对三阶行列式划去a 所在的第一行，

第一列后剩下的元素按原来的位置构成的二阶行列式，称为a 的余子式，记为M ，即

11 11

a a

22 23

M =

11 a a

32 33

类似的，

a a

21 23

M =

12 a a

31 33

a a

21 22

M =

13 a a

31 32

i+j

A =(-1) M (i,j=1,2,3)，称A 为元素a 的代数余子式，如：

ij ij ij ij

1+1

A =(-1) M =M

11 11 11

1+2

A =(-1) M =M

12 12 12

1+3

A =(-1) M =M

13 13 13

于是三阶行列式也可以定义为

a a a

11 12 13

a a a =a A +a A +a A

21 22 23 11 11 12 12 13 13

a a a =a M +a M +a M

31 32 33 11 11 12 12 13 13

上式右边也称为三阶行列式按第一行的展开式，类似的可以定义n阶行列式。

a a ….. a

11 12 1n

a a ….. a

定义,对n阶行列式A= 21 22 2n

a a …..a

31 32 3n

其行列式,

│A│=a A +a A +...+a A

i1 i1 i2 i2 in in

其中A 为a 的代数余子式（j=1,2,....,n)

ij ij

注意：（1）上述定义是对行列式按第i行展开（i=1,2,,...,n)，其实按行列式的任何一列展开也可作为行列式的定义。

（2）行列式和方阵是两个完全不同的概念，

n阶方阵是 n 个数按一定方式排列成的数表，而n阶行列式是这个数表按一定的运算法则所确定的一个数。

（3）行列式还有其他的定义方法，读者可自行阅读其他相关资料。

例1，用定义计算行列式

4 0 2

3 1 -1

2 2 3

解，按三阶行列式的定义，得

4 0 2

1 -1 3 -1 3 1

3 1 -1 =4 -0 +2 =28

2 3 2 3 2 2

2 2 3

1.3.3.克莱姆法则

在1.3.1中，我们借助于二阶行列式表示了二元线性方程组的解。下面把这一结论推广到n元线性方程组

定理1（克莱姆法则）含n个未知量、n个方程式的线性方程组

AX=B

其中A=(a ) ,X=(x ) ,B=(b ) ，

ij n×n j n×1 i n×1

若其系数行列式D=│A│≠0，则该方程组有唯一解

X = ,j=1,2,3,...n

j D

其中,

a ……a b a …..a

11 1j-1 1 1j+1 1n

D = a …..a b a …..a

j 21 2j-1 2 2j+1 2n

….. … …

a …..a b a …..a

n1 nj-1 n nj+1 nn

证明略.

例6，解线性方程组：

2x +x -5x +x =8

1 2 3 4

x -3x -6x =9

1 2 4

{

2x -x +2x =-5

2 3 4

x +4x -7x +6x =0

1 2 3 4

解：系数行列式为

2 1 -5 1

1 -3 0 -6

0 2 -1 2

1 4 -7 6

-3 0 -6 1 0 -6 1 -3 -6 1 -3 0

=2 2 -1 2 -1 0 -1 2 -5 0 2 2 -1 0 2 -1

4 -7 6 1 -7 6 1 4 6 1 4 -7

=2(18+84-24-42)-1(-6-6+14)-5(12-6+12-8)-1(-14+3+4)

=72-50-2+7=27

8 1 -5 1

9 -3 0 -6

D = =81

1 -5 2 -1 2

0 4 -7 6

2 8 -5 1

1 9 0 -6

D = =-108

2 0 -5 -1 2

1 0 -7 6

2 1 8 1

1 -3 9 -6

D = =-27

3 0 2 -5 2

1 4 0 6

2 1 -5 8

1 -3 0 9

D = =27

4 0 2 -1 -5

1 4 -7 0

于是方程组唯一的解为

x = =3

1 D

x = =-4

2 D

x = =1

3 D

x = =1

4 D

如果线性方程组右端的常数b (i=1,2,...,n)全为0，

即B为零矩阵，则称AX=0为齐次线性方程组。

显然，若D≠0，由于b =0(i=1,2,...,n)，所以D =0，

j j

从而当系数行列式D≠0时，齐次线性方程组AX=0有唯一解；

x =x =...=x =0

1 2 n

每个未知数的值都为零的解称为零解；反之，至少有一个未知数的值不等于零的解称为非零解。注意到任何齐次线性方程组总有解（因为它至少有零解），那么，在什么条件下齐次线性方程组有非零解呢？

定理2，齐次线性方程组AX=0为零解的充要条件是稀疏行列式D=0,

其中A=(a ) ,X=(x )

ij n×n j n×1

例7，判断齐次线性方程组

x +x +2x +3x =0

1 2 3 4

x +2x +3x -x =0

1 2 3 4

{

3x -x -x -2x =0

1 2 3 4

2x +3x -x -x =0

1 2 3 4

是否有非零解。

解，因该方程组的系数行列式

1 1 2 3

1 2 3 -1

D= =-153≠0

3 -1 -1 -2

2 3 -1 -1

所以方程组只有非零解

例8，k为何值时，齐次线性方程组

(5-k)x+2y+2z=0

{ 2x+(1-k)y=0

2x+(4-k)z=0

有非零解？

解：因

5-K 2 2

D= 2 6-K 0 =(5-k)(6-k)(4-k)-4(6-k)-4(4-k)

=(5-k)(2-k)(8-k)

2 0 4-K

所以，当D=0即k=2,5或8时，方程组有零解。最后，请读者注意，克莱姆法则的优点是不仅指出了解的存在性，而且还具体给出了解的表达式，但用克莱姆法则解线性方程组时有两个前提条件，一是方程个数与未知数个数相同，二是系数行列式不等于0，这使得克莱姆法则在应用时有很大的局限性。我们将在6.6中进一步研究线性方程组的理论和解法。

第九部分近似计算的准确度

下面的资料可参见《数值分析》，美国J.B.斯卡勃罗著，陈盖民，余介石译，科学出版社1958年出版。

11.由列表函数决定自变量的准确度

在许多问题中，必须计算含有一个未知量的某函数，于是就要从函数的列表值来定出这个未知量。像由对数表求出真数，及由三角函数表求出角度，都是这一类的例子。如果被计算的函数带有误差，那么，由这个函数定出的自变量就必然也多少有点出入。本节的目的便是研究所求自变量值的准确度。在单项式函数表中所列载的都是一个自变量的函数，

设自变量为x,列表函数为y，则

y=f(x)

由此可以近似的得出关系式

△y=f`(x)△x,

由此可得，

△y

△x= (11.1)

f`(x)

这是由函数表计算自变量的误差的基本方程，其中△y表示被计算函数的表值的误差，而△x表示自变量的对应误差。应该注意：△x的大小是依赖于函数的误差、函数的性质和自变量本身的大小三者而定的。现在我们把(11.1)应用到几个列表函数上去。

1.对数

f`(x)=1/x

故由(11.1)得

△x=x△y (1)

(b)f(x)=lgx

f`(x)=M/x.其中M=0.43429，△x=x△y/M=2.3026x△y,

因此，△x<2.31x△y (2)

2.三角函数。

(a)f(x)=sinx，f`(x)=cosx，

△x=△y/cosx=secx△y弧度 (3)

或，

(△x)``=206264.8secx△y秒 (4)

(b)f(x)=tanx

f`(x)=sec x

△x=cos x△y弧度 (5)

或

(△x)``=206264.8cos x△y秒 (6)

(c)f(x)=lg sinx

f`(x)=Mcosx/sinx=Mcotx

△x=△y/Mcotx=2.3026tanx*△y弧度 (7)

或

(△x)``<475000tanx*△y秒 (8)

(d)f(x)=lg tanx

2 2

f`(x)=Msec x/tan x=M/sinx*cosx=2M/sin2x

△x=sin2x*△y/2M=1.1513sin2x*△y

或，

△x<1.16sin2x*△y弧度 (9)

而，

(△x)``<238000sin2x*△y秒 (10)

3.指数函数，

f(x)=e

f`(x)=e

△x=△y/e (11)

4.其他列表函数。

当给定函数的导数已知或容易求得时，我们便可以用基本方程（11.1）算出任何自变量的误差。例如，在剑克和爱姆德的函数表中(Jahnke and Emde`s Funktionentafeln)便用表列出了：logΓ(x+1),[注：logΓ(x+1)=log(x+1)!]

x 2

误差函数∫ e dx

维尔斯德拉斯(Weierstrass)的p函数, 即勒让得(Legredre)多项式P (x)等的导数，

因此，我们应用这些表就能定出自变量和它的误差。椭圆积分是两个自变量的函数，每一个自变量的误差显然不能唯一的决定，但是应用公式(6.1)及等效原理，便能求出计算自变量误差的公式。譬如，令I表示椭圆积分，而用F(θ,φ)表示两个自变量的函数，则得

I=F(θ,φ),

故,

ӘF ӘF

△I=( )△θ+( )△φ

Әθ Әφ

假定右端两项相等，则得

ӘI ӘI

△θ= ,△φ=

ӘF ӘF

2 2

Әθ Әφ

若知道了微分的误差△I，我们就可以从这两个公式求出θ及φ的对应误差。注意，把公式(3)和(5)比较一下，就会看出：由角的正切去求角，常常比由角的正弦去求角所产生的

误差要小，这是因为cos x小于secx的缘故；后者可以具有1到∞之间的任何值，而前者其值不会超过1. 公式(7)和(9)更清楚的说明了由角的正切去求角的优越性。从(9)可以明显的看出：由于sin2x不会超于1，故x的误差很少会超过y的误差；可是由(7)可以看出：当角度从它的正弦对数来定出时，x的误差可能是y的误差的若干倍。让我们考虑一个数值的例子。

假定我们从五位的正弦对数表来求x，由于表中所列的值都是已末尾的数，所以，由表本身的固有误差所引起的△y的值，就可能会大到0.000005，取x=60°并代入（7），则得△x=2.3.26√3*0.000005=0.00002弧度（约数）=4``1, 所以，由x的logsin来求x，不可避免的误差就可能大到4秒。在另一方面来说，倘若由正切对数去求x，则从（9）可得

△<1.16*(√3/2)*0.00005=0.000005弧度=1``

因此这个误差仅为上述误差的四分之一。

上面的公式简单的证实了计算工作者久已知道的事：由角的正切或余切去求角，要比由角的正弦或余弦去求角来得更准确。注：在依靠表面而求得的结果中，定出最大可能误差的问题是颇有几分复杂的。读者在鲁路西著“数值计算讲义”(J.Luroth`s Vorie-sungen uber numerisches Rschen ,Leipzig,1900)一书中将会找到这问题的巧妙处理法。然而，这个问题的实用重要性是不大的，因为这种计算的误差极少会组合起来，造成它们的最大聚合影响的，它们在计算进行中都彼此相抵消了。

12.级数近似法的准确度。

将函数展成幂级数，算出头几项的值，这样去求得一个函数的数值，常比用其他方法求值容易些。事实上，这个方法，有时是计算函数数值唯一可能的方法。把函数展成幂级数的一般方法是应用泰勒公式进行的。这个公式的两种标准形式如下：

2 n-1 n

(x-a) (x-a) (n-1) (x-a) (n-1)

f(x)=f(a)+(x-a)f`(x)+ f``(a)+...+ f (a)+ f [a+θ(x-a)] ,0<θ<1 (1)

2! (n-1)! n!

2 n-1 n

h h (n-1) h (n)

f(x+h)=f(x)+hf`(x)+ f``(x)+...+ f (x)+ f (x+θh) ,0<θ<1 (2)

2! (n-1)! n!

在(1)式中令a=0，则得马克劳林(Maclaurin)公式：

2 n-1 n

x x (n-1) x (n)

f(x)=f(0)+xf`(0)+ f``(0)+...+ f (0)+ f (θx) ,0<θ<1

2! (n-1)! n!

这三个公式中，每一个末项都是n项以后的余项。这个余项是我们在本节中将要关心的量，余项的形式并不只限于上面的形式，其他有用的形式将在下面给出。

12a)在泰勒和马克劳林级数中的余项，令R 表示泰勒和马克劳林展开式中，

n项以后的余项，则可得下列的有用形式：

1.泰勒公式（1）：

(x-a) (n)

(a)R (x)= f [a+θ(x-a)], 0<θ<1

n n!

1 x-a (n) n-1

(b)R (x)= ∫ f (x-t)t dt

n (n-1)! 0

2.泰勒公式（2）：

h (n)

(a)R (x)= f (x+θh), 0<θ<1

n n!

1 h (n) n-1

(b)R (x)= ∫ f (x+h-t)t dt

n (n-1)! 0

3.马克劳林公式：

x (n)

(a)R (x)= f (θx), 0<θ<1

n n!

1 x (n) n-1

(b)R (x)= ∫ f (x-t)t dt

n (n-1)! 0

可以看出：第二种形式（积分形式）是完全确定的，并且不含有未定因素θ。可是无论用哪一种形式都必须先求出f(x)的n阶导数。

由于R (x)的积分形式，在微积分教科书中通常是不给出的，

所以我们要举例说明如何去应用它。

例.求ln(x+h)的展开式中，n项以后的余项。

解：这里f(x)=lnx,

f`(x)=-(1/x );

f``(x)=2/x ;

Ⅳ 4

f (x)=-(6/x );

…………………………………….

n-1

(n) (-1) (n-1)!

f (x)=

n-1 (n-1)! x 1 n-1

R (x)=(-1) ∫ t dt

n (n-1)! 0 n

(x+h-t)

现在由于t由0变到h，所以命被积函数中t=h便可求出R (x) 的最大值。

注：下面所讨论的都是指R (x)与有关各式的绝对值，但是没有写出绝对值的记号

n-1

于是把绝不会大于1的因子(-1) 略去，则得

h t dt 1 h n-1 1 h 1 h n

R (x)< ∫ = ∫ t dt= = ( )

n 0 n n 0 n n n x

x x x

假定x=1,h=0.01,于是h/x=0.01，所以如果我们要知道在ln1.01的开展式中，需要取多少项才能使所得结果准确到七位小数，可以取R ≤0.00000005.

1 n

( )(0.01) =0.00000005

显然可以看出，n=4便会使余项比能容许的误差小很多。因此我们在ln(x+h)的展开式中取四项。读者可以很容易证明：由余项的第一种形式也会给出像刚才求得的相同结果。

12b)交错级数，所谓交错级数就是指：各项是正负相同的无穷级数。这种级数只要：

(a)每项的绝对值小于前一项，并且(b)当n趋于无穷大时，第n项的极限是0，那么它就收敛。交错级数在应用数学中常会碰到，并且它最能满足计算的目的。因为要定出计算结果的误差经常是很容易的。决定误差的规则可简述如下：在收敛的交错级数中，计算到任一项为止所引起的误差，经常小于所舍弃部分的头一项。譬如，由于

2 3 4 5

x x x x

ln(1+x)=x- + - + -….

2 2 2 2

故有,

2 3

(0.01) (0.01)

ln(1.01)=0.01- + +R

2 3

这里,

(0.01)

R< =0.0000000025

所以只要在展开式中取三项，我们就能得到准确到八位小数的结果。

12c)一些重要的级数和它们的余项。

下面是一些最有用的级数和它们的余项，不过交错级数并不包括在内，因为它的余项可以根据上述规则来计算。

1.二项式级数

m m(m-1) 2 m(m-1)(m-2) 3 m(m-1)(m-2)...(m-n+2) n-1

(1+x) =1+mx+ x + x +...+ x +R

2! 3! (n-1)! n

其中，

(a)在一切情形下

m(m-1)(m-2)...(m-n+1) n m-n

R = x (1+θx) ,0<θ<1

n n!

(b)若x>0,n>m

m(m-1)(m-2)...(m-n+1) n

R < x

n n!

(c)若x<0及n>m，

m(m-1)(m-2)...(m-n+1) x

R <

n n! n-m

(1+x)

(d)若-1<m<0，

n m

R <│x │(1+x)

如果m是正或负的分数，或是负的整数，那么二项式展开式只在│x│<1时才成立。

并且，除m是正整数外，像(a+b) 的二项式，在展开之前必须要把它写成下列形式：

m b m

a (1+ ) 若a>b ;

或,

m b m

b (1+ ) 若b>a ;

注：都应就绝对值来说，则│a│<│b│与│b│<│a│

2.指数级数：

2 3 n-1 n

x x x x x θx

e =1+x+ + +…+ + e

2! 3! (n-1)! n!

2 n-1 n

x (xlga) (xlga) (xlga) θx

a =1+lga+ +…+ + a

2! (n-1)! n!

如果在(a)式中令x=1,则得计算e的级数如下：

1 1 1 1 e

(c)e=1+1+ + + +…+ +

2! 3! 4! (n-1)! n!

这里,

R = ;

n n!

但由于e<3而θ≤1，所以显然有

(d) R = ;

n n!

关于R 更确定的公式可求之如下：

把级数(c)写到n项以上，则得

1 1 1 1 1 1 1

e=[1+1+ + + +…+ ]+ + + +…

2! 3! 4! (n-1)! n! (n+1)! (n+2)!

这里n项以后的余项是

1 1 1

R = + + +…

n n! (n+1)! (n+2)!

1 1 1

= (1 + + +…)

n! n+1 (n+1)(n+2)

在右端括弧内的量，显然小于几何级数

1 1

1 + + +…

n 2

之和，此几何级数之和是

1 n

1 n-1

因此,

1 n

(e)R <

n n! n-1

或,

R <

n (n-1)(n-1)!

用公式(e)我们便能求出：要得到e准确到小数任意多少位的值，在展开式(c)中所必需取的项数。譬如，我们打算用级数(c)求准到小数十位的e值，那么我们可由方程

=0.00000000005

(n-1)(n-1)!

求出n。

靠阶乘倒数表的帮助，便可求出n-1=13或n=14。所以我们应该在级数(c)中取14项。如果我们要计算准确到小时100位的e值，那么我们可用同样的方法求出：我们应该在级数(e)中取71项。

3.对数级数

1 1 1 1

(a)ln(m+1)=lnm+2[ + + +…+ ]+R

2m+1 3 5 2m-1 n

3(2m+1) 5(2m+1) (2n-1)(2m+1)

要找出R 的上限，则有

1 1 1

R =2[ + + +...]

n 2m+1 2m+3 2m+5

(2n+1)(2m+1) (2n+3)(2m+1) (2n+5)(2m+1)

在括号内的级数，在第一项以后，它每一项都小于级数

1 1 1

+ + +...

2m+1 2m+3 2m+5

(2n+1)(2m+1) (2n+3)(2m+1) (2n+5)(2m+1)

或,

1 1 1

[1+ + +... ]

2m+1 2 4

(2n+1)(2m+1) (2m+1) (2m+1)

中的对应项，而后一个级数是以为公比的几何级数，

(2m+1)

其和为,

(2m+1)

或,

(2m+1)

4m(m+1)

因此,

1 (2m+1)

R <2( )

n 2m+1

(2n+1)(2m+1) 4m(m+1)

2n-1

2m(m+1)(2n+1)(2m+1)

所以,

(b)R <

n 2n-1

2m(m+1)(2n+1)(2m+1)

例1，试由(a)式取三项来计算ln2。由于m=1,n=3，故

1 1 1

ln2=2[ + + ]=0.693004

3 3 5

3(3) 5(3)

再由(b)式可得

1 1

R < =0.000147

n 2 5

2(7)(3)

它影响小数第四位。由于ln2取到八位小数的真值是0.69314718，所以上面求得值的误差是0.000143，它小于0.000147.

例2.求准确到小数的ln5之值。

这里m=4，而

1 1

R =

n 2 10

因为由(b)式得

1 1 1 1

2 4*5(2n+1)(9) 2 10

或

2n-1 3

(2n+1)(9) =5*10 =500000000

所以我们用试探方法求得n的约为4.1，当n=5时，这个对数将准确到十一位小数。

12d)某些n阶导数。

要计算级数的余项，就必须求出给定函数的n阶导数，为了使R 的计算便利起见，

所以我们给出了下面的某些简单函数的n阶导数一览表，记号D表示对x进行微分，即

D =

n x x n

(a)D a =a (lna)

(b)D sinx=sin[x+n(π/2)]

(c)D cosx=cos[x+n(π/2)]

n n

n 1 (-1) n!b

(d)D ( )=

a+bx n+1

(a+bx)

n 1 (-1) 1*3*5...(2n-1) n

(e)D = b

a+bx n (2n+1)/2

2 (a+bx)

n n

n 1 (-1) (n-1)!b

(f)D ln(a+bx)=

(a+bx)

n lnx (-1) n

(g)D ( )= [lnx-(1+1/2+1/3+...+1/n)]

x n+1

n-1 1

(-1) 2(n-1)!cos[n*arcsin( )]

n 2 1+x

(h)D ln(1+ x)=

2 n/2

(1+x )

n-1 1

(-1) 2(n-1)!sin[n*arcsin( )]

n 1+x

(i)D arctanx=

2 n/2

(1+x )

n 1

(-1) n!sin[(n+1)*arcsin( )]

n 1 1+x

(j)D =

2 2 (n+1)/2

1+x (1+x )

n α+βx (-1) n!

(k)D ( )= [βbcos(n+1)θ+(α+βa)sin(n+1)θ]

2 2 n+1

(x-a) +b bρ

其中，

2 2

ρ= (x-a) +b ;

θ=arctan

x-a

关于n阶导数的广泛研究，读者可参考斯替芬生的插值法(Stef-fensen:Interpolation),231-241页。

13.线性联立方程的解的准确度。

在应用数学的领域中，会碰到一些线性方程组，它的系数及常数项均带有少许的误差。这种误差可能是由于实验数值的不准确后末尾凑整而引起的，所以这样的方程组的解就不免有几分不准确，我们很需要有一种方法能估计这样解的固有误差。考虑简单的方程组

a x+b y=c

1 1 1

{ (1)

a x+b y=c

2 2 2

如果常数的正确值为a +△a ,b +△b 等等。

1 1 1 1

而x及y的对应真值是x+△x及y+△y，那么方程组(1)就成为

(a +△a )(x+△x)+(b +△b )(y+△y)=c +△c

1 1 1 1 1

{ (2)

(a +△a )(x+△x)+(b +△b )(y+△y)=c +△c

2 2 2 2 2

乘出之后，把含有两个误差相乘积的各项都去掉

（例如△a ,△x等等，因为它们是可忽略的），在应用（1）可得

a △x+x△a +b △y+y△b =△c

1 1 1 1

{ (3)

a △x+x△a +b △y+y△b =△c

2 2 2 2

现在，因为△a ,△b 等等都是假定已知的，而x及y可由（1）式求出，故得

1 1

a △x+b △y=△c -(x△a +y△b )

1 1 1 1 1

{ (4)

a △x+b △y=△c -(x△a +y△b )

2 2 2 2 2

既然这两个方程的右端是已知的，所以由（4）式就能很容易地求出△x及△y。应该注意：△x及△y在(4)式中的系数与x及y在（1）式中的系数相同。因此，如果用行列式解这些方程式（克拉茂规则，Cramer`s Ruie)，那么行列式

a b

1 1

a b

2 2

对两组方程均可应用。由于在（1）式中x及y的值是

c b -c b

1 2 2 1

a b -a b

1 2 2 1

c a -c a

2 1 1 2

a b -a b

1 2 2 1

所以很明显的，如果令c =c =k即可看出，

1 2

x及y的大小是和c 及c 的大小成正比的，这时便有

1 2

k( b -b )

1 2

a b -a b

1 2 2 1

k( a -a )

1 2

a b -a b

1 2 2 1

因此，要想得到误差△x及△y的上限，我们就应该使得（4）式右端尽可能大些，要做到这一点，应该在（4）式右端取各值绝对值之和。这应该注意的是：方程（8）只不过是（1）的微分，所以只要把给定方程微分一下，就能立刻把它写出来。

线性联立方程组的解，其固有误差的上限现在可按下面的程序把它求出：

1.依通常的方法，把给定方程组的未知数解出；

2.求出给定的方程组各方程的微分，并把结果写成（4）式的形式；

3.把误差的假定值△a 等等，以及在第一步所求得的x和y值，

4.在第三步所求得方程的右端，取各量的算术和（绝对值之和），并把方程组中的△x,△y等等解出。

例1，解方程组

3.21x-4.86y=5.73

2.13x+8.63y=12.65

假定式中所有数值都是已经末尾的，并且只准确到两位小数。

解，用行列式解这些方程

5.73 -4.36

12.65 8.63 49.4499+55.1540

x= = =2.828

3.21 -4.36 27.7023+9.2868

2.13 8.63

3.21 5.73

2.13 12.65 40.6065-12.2049

y= = =0.768

36.9891 36.9891

在继续进行以前，我们先把x和y的值代入给定的方程组中，看看它们能否满足这些方程，可以看出：它们都能满足。因为给定方程只准到小数两位，

所以误差△a ,△b 等等就不会大于0.005，根据上面的公式(4)得

1 1

a △x+b △y=△c -(x△a +y△b )

1 1 1 1 1

{ (4)

a △x+b △y=△c -(x△a +y△b )

2 2 2 2 2

因此，求x与y的可能误差的方程是

3.21△x-4.36△y=0.005+(2.828+0.768)(0.005)= 0.005(1+2.828+0.768)=0.0230,

2.13△x-8.63△y=0.005+(2.828+0.768)(0.005)=0.0230,

于是

0.023 -4.36

0.023 8.63 0.0230 1 -4.36

△x= = =0.0081

37(注：36.9891≈37) 37 1 8.63

3.21 0.023

2.13 0.023 0.0230 3.21 1

△y= = =0.0067

37(注：36.9891≈37) 37 2.13 1

这些误差影响x的二位小数及y的第三位小数。因此我们可取x=2.83及y=0.768，并理会到两数的末位数字都稍微有些出入。如果在给定的方程中，把第一个方程的y的系数变号，然后再解这个方程，则可求得△x=0.0033及△y=0.00084。

例2.解方程组

1.22x-1.32y+3.96z=2.12

{ 2.12x-3.52y+1.62z=-1.26

4.23x-1.21y+1.09z=3.22

所有各数均已末尾，并只准确到所给出的数字位数。

解.在这里,

1.22 -1.32 3.96

D= 2.12 -3.52 1.62 =64.404516-23.884520=40.520

4.23 -1.21 1.09

55.077264-16.832552

x= =0.9438

40.520

52.555064-12.938452

y= =1.227

40.520

47.612135-21.126204

z= =0.65865

40.520

如上所写的解，可以说明在减法中并没有丧失有效数字。可以看出，当把以上各值代入给定方程，除去第二个方程有0.001的差异外，其余均能满足。因为在系数及常数项中的可能误差不得超过0.005，所以误差方程是

1.22△x-1.32△y+3.96△z=0.005+(0.944+1.227+0.654)(0.005)=0.01912,

2.12△x-3.62△y+1.62△z=0.01912,

4.32△x-1.21△y+1.09△z=0.01912,

于是,

0.01912 -1.32 3.96

0.01912 -3.52 1.62

1 -1.32 3.96

0.01912 -1.21 1.09 0.01912

△x= = 1 -3.52 1.62 =0.0031

40.520 40.520

1 -1.21 1.09

同理△y=0.0020,△z=0.0031，因此x,y及z都准确到小数第2位，我们可取它们为：

x=0.94,y=1.23,z=0.65,

注：所求得的△x,△y及△z的值，都是x,y,z的最大可能误差，这些量的真实误差或许比最大误差要小的多，而且很可能是如此。

例3.解方程组

47.11x+13.72y=40.44

{

13.72x+4y=11.78

除去第二个方程中y的系数是完全正确数外，所有其他的数都是被抹尾的数，并且只准确到所给数字的位数。

解.用行列式来解，即得

40.44 13.72

11.78 4 161.76-161.6216 0.1334

x= = = =0.6865

47.11 18.72 188.44-188.2384 0.2016

13.72 4

3.21 5.73

2.13 12.65 554.9588-554.8368 0.1190

y= = = =0.5903

36.9891 0.2016 0.2016

可注意的是，在这个例中，前三个有效数字虽然已在减法中丧失。然而所得的值仍然能满足给定方程。更可注意的是，在上面的计算中，没有一个数会被我们抹尾，在联立方程组的解法中，不达到最后结果，就没有一个数应该被抹尾，或用任何方法加以省略，否则，就会引入计算的误差，而所得结果也不会满足给定方程。为了计算x和y的误差，便有

△a =△b =△c =△a =△c ≤0.005，及△b =0

1 1 1 2 2

因此误差方程是

47.11△x+13.72△y=0.005+(0.6865+0.5903)(0.005)=0.01138

{

13.72△x+4△y=0.005+0.6865*0.005=0.00843

故,

0.01138 13.72

0.00843 4 0.04552-0.11566

△x= = =-0.35

0.2016 0.2016

47.11 0.01138

13.72 0.01138 0.3971-0.1561

△y= = =1.20

0.2016 0.2016

x及y的可能误差如此的大，竟超过了这些量的本身，这就意味着所求得的x及y值都可能是无价值的。但是从另外的说明可知，x的值还大致准确，而y的值却毫无用处。

注：1.这个例题是从145节例3稍加改变而得的。

译者按：解这题时，要把一切计算求到八位有效数字，得到x=0.6691,y=0.65725才会合用。几何的研究帮组说明这个例题中的某些疑难的，看一下给定方程便可知这两方程的图形是几乎平行的直线。因此，由于方程系数的更改，而引起直线位置的少许变化，便会使得直线交点发生很大的变化。要想得到较可靠的结果，唯一的方法就是要有较准确的系数。在线性方程组的解中，固有误差的实在来源是由于减法中，首几位有效数字的丧失。

任何解方程的方法，都不免要使数量级相同的数相减，而当这两个这样的数几乎相等时，将发生首几位数字的丧失，这就会引起解的固有误差。

注：除了迭代法以外，参看164节。如果用行列式来解方程，而行列式又用子行列式的方法来展开，那么这种首几位数的丧失就能一望而知。在含有许多方程的方程组中，这种解法不合实用，所以首几位数字的丧失就无从找出，因此必然要有一种计算方法来求出解的最大固有误差。

例3.可以用来说明一件事情，即线性方程组的解，可能会比系数及常数项不可靠的多。因此，解线性方程组时，不达到终了，就不应该把外表似乎多余的数字去掉。然后应该将求得的结果代入给定方程内，看它是否能满足这些方程。并且在最后还应当算出解的误差的上限。给出最后结果的数字位数时，不应该超过由误差所能容许的位数。

13a.行列式中的误差

在上节的例3中，可以看出，如果一个行列式中的元素，由于抹尾凑整或其它原因使之成为不是完全正确数时，那么它在展开或求值的过程中，就会丧失最重要的有效数字，致使这行列式的值，受了严重的影响。这中丧失的大小，事先是不能确定的。但是，当一行列式中各个元素有了给定的可能误差时，我们就可以定出这行列式中误差的上限，现在我们取一个三阶行列式来作为说明。

设有

x x x

1 2 3

y y y (1)

1 2 3

z z z

1 2 3

现在，如果各个元素具有误差，其正负号未知，

而大小是△x ,△y 等，它们与x ,y 等比较起来是很小的，

1 1 1 1

于是由此，行列式D将误差△D，使

x +△x x +△x x +△x

1 1 2 2 3 3

D+△D= y +△y y +△y y +△y (2)

1 1 2 2 3 3

z +△z z +△z z +△z

1 1 2 2 3 3

按行列式的加法定理，（2）式右边可以展成八个行列式的和，其中第一个便是原来的行列式D。其余的行列式中，有三个各自含有误差元素一直列，另有三个各自含有误差元素两直列，而剩下一个行列式有三直列的误差元素。含有误差元素一列以上的行列式可以略去，因为在展开时，所得的各项皆将含有这些误差的二次方和三次方，故它们与只含误差一次方各项比较时，是可以略去的。由此可知△D的值是三个行列式的和，每一个只含单独一直列的误差元素。但是这些行列式只是D的微分，故有

dx x x

1 2 3

dD= dy y y +

1 2 3

dz z z

1 2 3

x dx x

1 2 3

+ y dy y +

1 2 3

z dz z

1 2 3

x x dx

1 2 3

+ y y dy (3)

1 2 3

z z dz

1 2 3

即

dD=(y z -y z )dx -(x z -x z )dy +(x y -x y )dz -

2 3 3 2 1 2 3 3 2 1 2 3 3 2 1

-(y z -y z )dx +(x z -x z )dy -(x y -x y )dz +

1 3 3 1 2 1 3 3 3 1 2 1 3 3 1

+(y z -y z )dx -(x z -x z )dy +(x y -x y )dz

1 2 2 1 3 1 2 2 1 3 1 2 2 1 3

只有在各元素的号与误差的号能使(4)式右边的十八项皆成相同符号时，才会产生最大可能的误差——这样的可能性是极少的。方程（4）表明，含有不准确元素的行列式中，误差可能是从零到相当大的一个数中间的某一个数，然而必须记住，（4）中各项，大多数彼此相消，故在一般情况中，dD不会很大。

招烘焙师十名	生日刀叉盘组	香美人皮草店	马卡龙包装盒	布丁瓶奶瓶泡茶瓶	进口食用墨水
数码蛋糕样品	绵羊皮真皮衣	麻花机油炸锅	广告位10元/月	广告位10元/月	数码蛋糕打印机

中国面包师贴吧-楼主(阅：1608/回：0)使用插值法公式组成数字电路进行计算的计算机