注释结束。
2
要想求出所需量表示的△ y 的值，则有
0
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
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
3
要找出△ y 我们从下式着手（见上页注）： 0
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
5
自(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
4
要找出△ y 我们从下式着手：
0
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
7
自(s)减去(r)，并解出△ y
-3
7 8
(t) △ y =m -△ y /2
-3 7 -4
6
将(t)代入(q)，并解出△ y
-2
6 6 8
(u) △ y =△ y +m -△ y /2
-2 -3 7 -4
5
将(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
或
2
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
0
u=
h
在这个公式当中一共有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
△y 2
△ y 3
△ y 4
△ y 5
△ y 6
△ y 7
△ y 8
△ y
y
-4
△y
-4
y
-3 2
△ y
-4
△y
-3 3
△ y
-4
y
-2 2
△ y
-3 4
△ y
-4
△y
-2 3
△ y
-3 5
△ y
-4
y
-1 2
△ y
-2 4
△ y
-3 6
△ y
-4
△y
-1 3
△ y
-2 5
△ y
-3 7
△ y
-4
y
0 2
△ y
-1 4
△ y
-2 6
△ y
-3 8
△ y
-4
△y
0 3
△ y
-1 5
△ y
-2 7
△ y
-3
y
1 2
△ y
0 4
△ y
-1 6
△ y
-2 8
△ y
-3
△y
1 3
△ y
0 5
△ y
-1 7
△ y
-2
y
2 2
△ y
1 4
△ y
0 6
△ y
-1
△y
2 3
△ y
1 5
△ y
0
y
3 2
△ y
2 4
△ y
1
△y
3 3
△ y
2
y
4 2
△ y
3
△y
4
y
5
现在我们令
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 为始点的牛顿公式(Ⅰ)。
0
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
2
欲求得△ y 我们从下式着手：
0
由定义可知
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
2
欲求△ y 则自定义得 0
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
4
自(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
4
欲求△ y 则从下式着手： 0
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
6
自(m)减去(L)并解出△ y
-2
6 7
(n) △ y =m +△ y /2 -2 6 -3
5
令(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 +△ yu(u-1)(u+1)(u-2)(u+2)(u-3) -3 -3
-
+…+ 6! 2 2n 2n △ y +△ yu(u-1)(u+1)(u-2)(u+2)...(u-n)(u+n-1) -n -n
-
+ (2n)! 2u(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
[135...(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
(v -1/4)(v -9/4)(v -25/4) -2 -1 -
+…+ 6! 2 2
(v -1/4)(v -9/4)...[v - ] △ y +△ y
2 -n -n+1 -
+ 6! 2 2
(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
2 x -x
f(x)= ∫ e dx
√π 0
的各值。求积分在x=0.5437时的值。
解：这里我们取x =0.54,x=0.5437。由于h=0.01,故
0
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)使用斯特灵公式（Ⅲ）则得
2
(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)+
2
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,则
0
1.7489-1.75 0.0011
u= =- =-0.11
0.01 0.01
表B
x -x
e
△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 +
2
0.0121 (0.0121-1) (0.0000001749-0.0000001727)
-
(0.0000173774)-0.11 + 2 6 2 (0.0121-1)
+0.0121 (22)
24
=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)+
2
0.39*0.39-0.25 0.0000175523+0.0000173774
-
( )( ) +
2 2(0.39*0.39-0.25)
+(0.39) (0.0000001749)+ 6
(0.390.39-0.25)(0.390.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,这时可得
0
这个数值也准确到十位小数，不过这个级数要比前面的情况收敛得稍慢一些；从贝赛尔公式所得出的两个级数都比斯特灵公式所得出的这个级数要收敛得稍慢一些。
注意，在这里自然而然的会产生这样一个问题：
答案是：它们的准确度不相上下，对于给定的差分表来说，
收敛的迅速与否同公式(Ⅲ)中u的大小，以及公式(Ⅵ)中v的大小有关。u和v的值愈小，则级数收敛得愈快。所以我们应该经常选择始点x 以便能促使u和v尽可能地小。
0
在多数情况下，选择起始点使得-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(φ)= ∫
2
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
所以,
y
1
A =
1 (x -x )(x -x )...(x -x )
1 0 1 2 1 n
同样我们可以求得
y
2
A =
2 (x -x ) (x -x )(x -x )...(x -x )
2 0 2 1 2 3 2 n
……………………………
y
2
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)，
n
但分母显然是∏` (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
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+
90897180943270128 9089790046179231
153445*62548*116683 153445*62548*27498
-
*0.48+ *0.49
=-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 逼除，则得
0
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
0
[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
4
[1] [1] [1] [1] △ y
u (u -1)(u -2)(u -3) 0
- (2)
4! △y
0
其三次近似值为
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
4
[2] [2] [2] [2] △ y
u (u -1)(u -2)(u -3) 0
- (3)
4! △y
0
依次类推以至更高次的近似值。现在我们举例说明这个方法。
例1.现在已经给出概率积分
2
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并使用贝塞尔公式。所以我们可得
0
x =0.47,h=0.01,y=1/2=0.5
0
将这个y值以及表中查到的适应量代入贝塞尔公式（Ⅵ），则得
2
(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.194623(0.194623) -0.25
=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,故得
0
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代入斯特灵公式（Ⅲ），则得
2
62=61.9785+0.6199u+0.0031u
或,
2
31u +6199u=215
所以,
2
-6199+ (6199) +431215
u=
62
-6199+6201.15 2.15
= = =0.0347
62 62
因为h=0.01,而x=x +hu,故得
0
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
D=
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 所在的第一行，
11
第一列后剩下的元素按原来的位置构成的二阶行列式，称为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）行列式和方阵是两个完全不同的概念，
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，则该方程组有唯一解
D
j
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
D=
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
于是方程组唯一的解为
D
1
x = =3
1 D
D
2
x = =-4
2 D
D
3
x = =1
3 D
D
4
x = =1
4 D
如果线性方程组右端的常数b (i=1,2,...,n)全为0，
i
即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 2 f(x)=sec x
2
△x=cos x△y弧度 (5)
或
2
(△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/sinxcosx=2M/sin2x
△x=sin2x△y/2M=1.1513sin2x△y
或，
△x<1.16sin2x△y弧度 (9)
而，
(△x)``<238000sin2x△y秒 (10)
3.指数函数，
x
f(x)=e
x
f`(x)=e
x
△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
0
维尔斯德拉斯(Weierstrass)的p函数, 即勒让得(Legredre)多项式P (x)等的导数，
n
因此，我们应用这些表就能定出自变量和它的误差。椭圆积分是两个自变量的函数，每一个自变量的误差显然不能唯一的决定，但是应用公式(6.1)及等效原理，便能求出计算自变量误差的公式。譬如，令I表示椭圆积分，而用F(θ,φ)表示两个自变量的函数，则得
I=F(θ,φ),
故,
ӘF ӘF
△I=( )△θ+( )△φ
Әθ Әφ
假定右端两项相等，则得
ӘI ӘI
△θ= ,△φ=
ӘF ӘF
2 2
Әθ Әφ
若知道了微分的误差△I，我们就可以从这两个公式求出θ及φ的对应误差。注意，把公式(3)和(5)比较一下，就会看出：由角的正切去求角，常常比由角的正弦去求角所产生的
2
误差要小，这是因为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弧度（约数）=41, 所以，由x的logsin来求x，不可避免的误差就可能大到4秒。在另一方面来说，倘若由正切对数去求x，则从（9）可得 △<1.16*(√3/2)*0.00005=0.000005弧度=1
因此这个误差仅为上述误差的四分之一。
上面的公式简单的证实了计算工作者久已知道的事：由角的正切或余切去求角，要比由角的正弦或余弦去求角来得更准确。注：在依靠表面而求得的结果中，定出最大可能误差的问题是颇有几分复杂的。读者在鲁路西著“数值计算讲义”(J.Luroths 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
n项以后的余项，则可得下列的有用形式：
1.泰勒公式（1）：
2
(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）：
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.马克劳林公式：
n
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)的积分形式，在微积分教科书中通常是不给出的，
n
所以我们要举例说明如何去应用它。
例.求ln(x+h)的展开式中，n项以后的余项。
解：这里f(x)=lnx,
2
f`(x)=-(1/x );
3
f``(x)=2/x ;
Ⅳ 4
f (x)=-(6/x );
…………………………………….
n-1
(n) (-1) (n-1)!
f (x)=
n
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) 的最大值。
n
注：下面所讨论的都是指R (x)与有关各式的绝对值，但是没有写出绝对值的记号
n
n-1
于是把绝不会大于1的因子(-1) 略去，则得
n
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.
n
1 n
( )(0.01) =0.00000005
n
显然可以看出，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
这里,
4
(0.01)
R< =0.0000000025
4
所以只要在展开式中取三项，我们就能得到准确到八位小数的结果。
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，
n
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)
n
如果m是正或负的分数，或是负的整数，那么二项式展开式只在│x│<1时才成立。
m
并且，除m是正整数外，像(a+b) 的二项式，在展开之前必须要把它写成下列形式：
m b m
a (1+ ) 若a>b ;
a
或,
m b m
b (1+ ) 若b>a ;
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!
这里,
θ
e
R = ;
n n!
但由于e<3而θ≤1，所以显然有
θ
3
(d) R = ;
n n!
关于R 更确定的公式可求之如下：
n
把级数(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
n
之和，此几何级数之和是
1 n
=
1 n-1
1-
n
因此,
1 n
(e)R <
n n! n-1
或,
1
R <
n (n-1)(n-1)!
用公式(e)我们便能求出：要得到e准确到小数任意多少位的值，在展开式(c)中所必需取的项数。譬如，我们打算用级数(c)求准到小数十位的e值，那么我们可由方程
1
=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 的上限，则有
n
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)
1
中的对应项，而后一个级数是以 为公比的几何级数，
2
(2m+1)
其和为,
1
1
1-
2
(2m+1)
或,
2
(2m+1)
4m(m+1)
因此,
2
1 (2m+1)
R <2( )
n 2m+1
(2n+1)(2m+1) 4m(m+1)
1
=
2n-1
2m(m+1)(2n+1)(2m+1)
所以,
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
10
因为由(b)式得
1 1 1 1
=
2 45(2n+1)(9) 2 10
10
或
2n-1 3
(2n+1)(9) =510 =500000000
所以我们用试探方法求得n的约为4.1，当n=5时，这个对数将准确到十一位小数。
12d)某些n阶导数。
要计算级数的余项，就必须求出给定函数的n阶导数，为了使R 的计算便利起见，
n
所以我们给出了下面的某些简单函数的n阶导数一览表，记号D表示对x进行微分，即
d
D =
dx
n x x n
(a)D a =a (lna)
n
(b)D sinx=sin[x+n(π/2)]
n
(c)D cosx=cos[x+n(π/2)]
n n
n 1 (-1) n!b
(d)D ( )=
a+bx n+1
(a+bx)
n
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)=
n
(a+bx)
n
n lnx (-1) n
(g)D ( )= [lnx-(1+1/2+1/3+...+1/n)]
x n+1
x
n-1 1
(-1) 2(n-1)!cos[n*arcsin( )]
2
n 2 1+x
(h)D ln(1+ x)=
2 n/2
(1+x )
n-1 1
(-1) 2(n-1)!sin[n*arcsin( )]
2
n 1+x
(i)D arctanx=
2 n/2
(1+x )
n 1
(-1) n!sin[(n+1)*arcsin( )]
2
n 1 1+x
(j)D =
2 2 (n+1)/2
1+x (1+x )
n
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 ;
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）可得
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
x=
a b -a b
1 2 2 1
c a -c a
2 1 1 2
y=
a b -a b
1 2 2 1
所以很明显的，如果令c =c =k即可看出，
1 2
x及y的大小是和c 及c 的大小成正比的，这时便有
1 2
k( b -b )
1 2
x=
a b -a b
1 2 2 1
k( a -a )
1 2
y=
a b -a b
1 2 2 1
因此，要想得到误差△x及△y的上限，我们就应该使得（4）式右端尽可能大些，要做到这一点，应该在（4）式右端取各值绝对值之和。这应该注意的是：方程（8）只不过是（1）的微分，所以只要把给定方程微分一下，就能立刻把它写出来。
线性联立方程组的解，其固有误差的上限现在可按下面的程序把它求出：
1.依通常的方法，把给定方程组的未知数解出；
2.求出给定的方程组各方程的微分，并把结果写成（4）式的形式；
3.把误差的假定值△a 等等，以及在第一步所求得的x和y值，
1
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
D=
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 3z dz z
1 2 3x x dx 1 2 3 -
y y dy (3) 1 2 3z 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不会很大。
