使用插值法公式组成数字电路进行计算的计算机
使用插值法公式组成数字电路进行计算的计算机是一种可以使用插值法计算积分值，导数值，函数值的数字计算机，它由按键，液晶显示器，中央处理器组成。按键输入的程序保存在磁带上面，中央处理器在程序的作用下按照牛顿向前插值公式，牛顿向后插值公式，斯特灵插值公式，贝塞尔插值公式，拉格朗日外插法公式进行计算，按照数学计算公式进行编程。
中央处理器由程序语句判断执行电路，程序语句判断控制电路，端口，加法器，减法器，乘法器，除法器，n次方计算器，对数计算器，三角函数计算器构成。
键盘输入的程序按每行保存在磁带中，程序语句判断电路根据键盘输入的程序的关键字判断电路执行相应的操作，例如输入ADD，电路执行加法操作，程序语句判断控制电路根据键盘输入的程序的关键字控制电路的工作，例如输入NIUDUN DIEDAI，电路将上面计算电路执行多次，进行牛顿迭代计算。
它的相关资料下载网址为：
链接: https://pan.baidu.com/s/1426ypZw_bd9NvfazElx3MQ?pwd=pmjn 提取码: pmjn 复制这段内容后打开百度网盘手机App，操作更方便哦
https://photo.baidu.com/photo/wap/albumShare/invite/nsQhGUEoOS?from=linkShare
https://photo.baidu.com/photo/wap/albumShare/invite/xBzDCVctvX?from=linkShare
https://photo.baidu.com/photo/wap/albumShare/invite/KoDqrAArsD?from=linkShare
https://photo.baidu.com/photo/wap/albumShare/invite/wjcLIUupA?from=linkShare
链接:https://pan.baidu.com/s/1giJmPpB2iaq2fJY259XVBw?pwd=9ug8
提取码:9ug8
链接:https://pan.baidu.com/s/1l7d4UX3hNFhFlUYvLIjtAQ?pwd=7ocu
提取码:7ocu
https://www.aliyundrive.com/s/5Uxz1AsL7G5
微云文件分享:插值法数字计算机下载地址:https://share.weiyun.com/5nwXzVag
https://115.com/s/sw6qrr833u5?password=b453#
插值法数字计算机
访问码：b453
https://kdocs.cn/join/ge5wqfb?f=101
第一部分使用插值法公式组成数字电路进行计算的计算机
该计算器首先通过晶振产生32768HZ的谐振方波信号，再经过分频电路将这个方波信号的频率降低为100HZ,，即周期为0.01秒，再将这个100HZ的信号接入到按键的公共端，按键共有60个，它们的一端接到一起，另外一端分别接到倍频器上。相当于这些按键并联在一起，当某个按键被按下时，100HZ的信号就会接入到倍频器上，经过倍频后，频率变为1HZ，
为什么按键上面的频率是100HZ，这是因为100HZ的频率，周期是1毫秒，通常使用者按下按键的时间在1毫秒左右，所以，只有这个频率的信号才会在按下按键时输入到后级电路中。键值计算电路由十进制转二进制电路组成，当有数字键按下时，对应的数字按键输出端输出对应的数值。数值按键的输出端接上或门，或门两两相接，最后输出一个或门，当有任何计算符号按键按下时，或门输出高电平，或门后面接上计数器，计数器记录按键按下的次数，当有按键按下时，计数器将对应的次数输入到加法器，加法器给键值乘以10，100，1000，等倍数。当连续按2次按键时，需要用乘法器给键值乘以10，连续按下3次按键时，需要用乘法器给键值乘以100，依次类推。所有数值按键的输出端连接到一起，输出到计算符号电路，进行计算。计算符号编码电路产生对应计算符号的编码，输送给计算符号按键电路。用计算符号按键输入计算符号±×÷,cos,sin,ln,log,等，
当RS触发器的输入端R,S都是1时，触发器保持输出端没有变化。利用这个特点，当按键输入高电平1时，电路输出高电平1给存储器，当按键断开输入低电平0时，RS触发器仍然给存储器输入1，当清零键按下时，RS触发器的S端输入0，触发器给存储器输入0，存储器清零。
当有按键按下时RS触发器Q输出1， Q 输出0，按下清零键以后，RS触发器Q端输出0， Q 端输出0
按键编码器产生二进制编码，每个编码对应一个按键。
当数字键1，按下时，这个与门输出0000001给后面计算电路，所有按键存储器后面两两之间接上或门，或门后面再接上或门，最后接上计数器，当按键按下时，计数器变为1，对应的存储器输出对应键值。当按键按下第二次时，计数器输出2，输出两位数字，当按键按下第三次时，计数器输出3，输出三位数字。
经过两个异或门和一个或门以后输出高电平111111111，这使后面的与门输出按键的数值到寄存器1，
当开始输入时，按清零键，计算机按键输入为0.此时，开始输入字符，将字符输入到寄存器1，
按键输入的程序存储在磁带A上面，超强磁性磁带的基材由50%醋酸酯DAC，50%醋酸酯TAC构成，超强磁性磁带的磁性粉末粘合剂有1%氯乙烯，1%醋酸乙烯共聚体，1%苯乙烯-丁二烯共聚体，1%硝化纤维素。1%纤维素，1%丁腈橡胶，1%丙烯酸酯橡胶，1%无定形聚酯，1%氨酯橡胶，1%聚氨基甲酸乙酯树脂，环氧树脂，密胺树脂，1%醋酸乙烯，1%丙烯酸酯丁基系的软质树脂，超强磁性磁带的磁性粉末分散剂由10ml乙醇，20g尿素，10ml双氧水，10g蔗糖，20g聚乙二醇4000，油酸钾皂试剂20g，黄色色素10g,司盘80试剂10ml，氧化铝10g，氨水50g，大豆油10g，α-烯基磺酸钠5g,十二烷基苯磺酸钠5g,烯丙基磺酸钠5g,二甲苯磺酸钠5g，椰子油脂肪酸渗透二乙醇酰胺6501日化，1%卵磷脂组成，磁性粉末稳定剂有对氯乙烯系粘合剂，使用硬脂酸钡等金属无机盐。磁性粉末防带静电剂是在磁性层内渗入炭黑或石墨等固体导电粉末。超强磁性磁带的磁性粉由二氧化铬，三氧化二铁，铬化铁，氧化镍，氧化钴，氧化钇，镝，二氧化锰。把磁性粉末，粘合剂，增塑剂，稳定剂，分散剂，加入水中，使各个磁性粉末相互溶解到水里，再球磨机混合均匀，最后用刮片涂覆到基材上面。
注意：收音机磁带使用涂着四氧化三铁的硝酸纤维素条，铁芯（铁氧体/羟基铁芯），0.32-0.45mm变压器钢片，线圈（0.08mm漆包线1200-1500匝），放音头间隙0.02mm，工作间隙0.5mm，磷铜萡/黄铜箔，
磁带录音机电路如下：
按键电路如下：
计算机中央处理器CPU电路原理图
程序语句判断电路
程序关键字判断电路，程序关键字判断电路，查询到关键字，并执行该关键字所要求的功能。
程序计算符号判断电路，程序计算符号判断电路，查询到计算符号，并执行该计算符号所要求的功能。
数据判断电路，程序数据判断电路，查询到数据符号，并执行该数据符号所形成的数据。
字符判断电路，程序字符判断电路，查询到字符，并执行该字符的功能。
磁带程序判断执行电路原理图。
语句执行电路，按照语句判断的输出，执行这条语句，输出到CPU端口并执行。
关于数字电路加法器，计数器，分频器的电路可参见《中国集成电路大全》丛书，《中国集成电路大全编写委员会编，国防工业出版社1987年出版.
该计算器首先通过晶振产生32768HZ的谐振方波信号，再经过分频电路将这个方波信号的频率降低为100HZ,，即周期为0.01秒，再将这个100HZ的信号接入到按键的公共端，按键共有60个，它们的一端接到一起，另外一端分别接到倍频器上。相当于这些按键并联在一起，当某个按键被按下时，100HZ的信号就会接入到倍频器上，经过倍频后，频率变为1HZ为什么按键上面的频率是100HZ，这是因为100HZ的频率，周期是1毫秒，通常使用者按下按键的时间在1毫秒左右，所以，只有这个频率的信号才会在按下按键时输入到后级电路中。键值编码电路由二进制编码电路组成，当有按键按下时，对应的按键输出端输出对应的按键编码。每个按键的输出端接上或门，或门两两相接，最后输出一个或门，当有任何计算按键按下时，或门输出高电平，这个或门在和每个按键的输出端接上与门，这些与门在两两之间接上或门，最后一个或门接上按键寄存器。按键寄存器将输入的按键输出保存到磁带寄存器A中，计算机CPU通过算法语言关键字判断语句，计算符号判断电路，中断判断电路，定时器判断电路，数据判断电路，选择性的判断执行磁带存储器A中的按键输入程序。计算机CPU通过执行电路执行上面语句判断电路输出的内容。最后将执行结果通过IO端口输出，并用液晶显示器显示出来。
如果出现PROGRAM BEGIN说明程序开始，与门导通，如果出现空格说明前面是一个关键字，或字符或数据，与门导通。如果出现回车说明前面是一个程序段，需要执行这段程序，与门导通。
关键字比较电路，和每个关键字的代码相互比较，如果代码相同·，执行该关键字的功能。
字符比较电路，和每个字符的代码相互比较，如果代码相同·，执行该字符的功能。
数据比较电路，和每个数据的代码相互比较，如果代码相同·，产生该数据的二进制编码。
磁带程序判断执行电路原理图
磁带程序执行控制电路原理图
出现BSEJS Ⅴ#GONGSHI 时，程序将执行贝塞尔插值法Ⅴ#公式计算电路。
出现NDXQJS LG FANGCHENG B时，程序将执行牛顿向前插值法计算电路算电路。
出现LGLRJS Ⅸ#GONGSHI时，程序将执行拉格朗日反插法公式Ⅸ导数计算电路。
计算机原理图如下：
第二部分 插值法计算电路
下面的资料可参见《数值分析》，美国J.B.斯卡勃罗著，陈盖民，余介石译，科学出版社1958年出版。
例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,故
0
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）+
2
=0.4970679364+0.000819310+0.0000000053=0.4971498727
此结果准确到它的末位数字。
牛顿向前插值法计算电路程序：
PROGRAM BEGIN 程序开始
JISUAN2# 进入JISUAN2#选择电路
DISPLAY ON 计算机液晶显示器显示已知数值
lg3.141=0.4970679364,
lg3.142=0.4972061807,
lg3.143=0.4973443810,
lg3.144=0.4974825374,
lg3.145=0.4976206498,
DISPLAY OFF
NDXQJS LG FANGCHENG B 执行牛顿向前插值法计算电路
MAKEA X0=3.141,X1=3.142,X2=3.143,X3=3.144,X4=3.145,
将方程式的系数用X0,X1,X2,X3,X4存储器存储起来
MAKEB LGX0=0.4970679364,LGX1=0.4972061807,LG2=0.4973443810,LG4=0.4974825374,LG4=0.4976206498,
MAKEC X=3.1415926536计算方程组u,△y的公式是固定公式，对应电路也是固定电路
MAKED H=0.001可以将键盘输入的X0,X1,X2,X3,X4分别代入电路中的寄存器，计算函数值
RUN NDXQJS 与门导通，接通牛顿向前插值法计算电路，开始计算
SAVE JISUAN2# ZU 保存计算中所有JISUAN2#电路寄存器的数据
OUTPUT JISUAN2# 1#IO 将计算电路JISUAN1#寄存器中的数据输送到计算机1#端口
DISPLAY ON 计算机显示器显示计算得到的数值lgx
LGX
DISPLAY OFF 停止显示
PROGRAM END 程序结束
用数字电路加法器，乘法器，除法器，各种数值表示上面的函数，
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
表示上面公式的电路图如下：
例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
由于表中的差分是完全有规则的并且很快的减小，所以这个结果或许可以准确到末一位数字。
用贝赛尔半数插值Ⅴ#公式计算积分的程序：
PROGRAM BEGIN 程序开始
JISUAN3# 计算机液晶显示器显示积分公式
DISPLAY ON
dφ
F(φ)= ∫
2
1-sin φ/2
DISPLAY OFF
BSEJS Ⅴ#GONGSHI 与门导通，选择贝赛尔半数插值Ⅴ#公式计算电路
MAKEA F(21°)=0.370634373,F(22°)=0.388705151,F(23°)=0.406834931,F(24°)=0.425026420,F(25°)=0.443282329,F(26°)=0.461605362,
用键盘输入积分公式的初始值F(21°)=0.370634373,F(22°)=0.388705151,F(23°)=0.406834931,F(24°)=0.425026420,F(25°)=0.443282329,F(26°)=0.461605362,
MAKEB A=21,B=22,C=,23,D=24,E=25,F=26
计算△y的公式是固定公式，对应电路也是固定电路
MAKEC M=23.5
RUN BSEJS Ⅴ#GONGSHI 与门导通，接通贝赛尔半数插值Ⅴ#公式计算电路，开始计算
SAVE JISUAN2# ZU 保存计算中所有JISUAN2#电路寄存器的数据
OUTPUT JISUAN2# 1#IO 将计算电路JISUAN1#寄存器中的数据输送到计算机1#端口
DISPLAY ON 计算机显示器显示计算得到的数值F(23.5°)
DISPLAY OFF 停止显示
PROGRAM END 程序结束
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。
两数相减，得到的数据进行平方计算，如果是正数，平方后得到正数，开方后得到正数，加上原来的正数，输出正数。
如果加法器输出0，或门输出1，与门输出1，状态寄存器A保存1，证明两数相减得到的是负数，反之，则是正数
当23.5处于23，24之间时，两个与门输出不同的电平，当两个异或门输入不相同电平时，异或门输出高电平，也就是说，当23.5处于23和24之间时，23和24末端的与门输出高电平，
第二部分 拉格朗日公式，反插法计算电路
下面的资料可参见《数值分析》，美国J.B.斯卡勃罗著，陈盖民，余介石译，科学出版社1958年出版。
例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
-
1809439004689185 27012817923189185*0.48+ *0.49
=-0.0207787+0.157737+0.369928-0.0299495=0.476937
准确到六位小数的真值为0.476936, 注：这个问题的计算除非有计算机可用，否则应该用对数来做。
(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
用拉格朗日反插法公式Ⅸ计算导数的程序：
PROGRAM BEGIN 程序开始
JISUAN3#
选择保存计算数据到JISUAN3#寄存器，同时标志寄存器保存3，说明今天进行了3次计算
DISPLAY ON 计算机液晶显示器显示积分公式
2
2 x -x
∫ e dx x
√π 0
0.4846555 0.46
0.4937452 0.47
0.5027498 0.48
0.5116683 0.49
DISPLAY OFF
LGLRJS Ⅸ#GONGSHI 与门导通，选择拉格朗日反插法公式Ⅸ导数计算电路
MAKEA x0=0.46,x1=0.47,x2=0.48,x3=0.49
用键盘输入积分公式的自变量x的各个初始值 x0=0.46,x1=0.47,x2=0.48,x3=0.49
MAKEB y=0.5,y0=0.4937455,y1=0.4937452,y2=0.5027498,y3=0.5116683
用键盘输入积分公式的y的各个初始值 y=0.5,y0=0.4937455,y1=0.4937452,y2=0.5027498,y3=0.5116683
RUN LGLRJS Ⅸ#GONGSHI
计算x的公式是固定公式，对应电路也是固定电路，与门导通，接通拉格朗日反插法公式Ⅸ导数计算电路，开始计算。
SAVE JISUAN2# ZU 保存计算中所有JISUAN2#电路寄存器的数据
OUTPUT JISUAN2# 1#IO 将计算电路JISUAN1#寄存器中的数据输送到计算机1#端口
DISPLAY ON 计算机显示器显示计算得到的数值x
X
DISPLAY OFF 停止显示
PROGRAM END 程序结束
(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
用数字电路加法器，乘法器，除法器，各种数值表示上面的函数，
第三部分 克莱姆法则计算电路
下面的资料可参见《工程数学》，林益编，高等教育出版社2003年出版
用克莱姆法则计算线性方程组的程序：
PROGRAM BEGIN 程序开始
JISUAN2#
选择保存计算数据到JISUAN2#寄存器，同时标志寄存器保存2，说明今天进行了2次计算
DISPLAY ON 计算机液晶显示器显示线性方程组
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
DISPLAY OFF
KLMJS SIYUAN FANGCHENG B
与门导通，选择线性方程组的克莱姆法则计算电路
MAKEB a11=2,a12=1.a13=-5,a14=1,b1=8
将方程式的系数用a11,a12,a13,a14,b1存储器存储起来
MAKEB a21=1,a22=-3.a23=0,a24=-6,b2=9
MAKEB a31=0,a32=2.a33=-1,a34=2,b3=-5
MAKEB a41=1,a42=4.a43=-7,a44=6,b4=0
计算方程组D,D1,D2,D3,D4的公式是固定公式，对应电路也是固定电路, 可以将键盘输入的a11=2,a12=1.a13=-5,a14=1,b1=8分别代入电路中的寄存器，计算函数值。
RUN KLMJS与门导通，接通克莱姆法则计算电路，开始计算
SAVE JISUAN2# ZU 保存计算中所有JISUAN2#电路寄存器的数据
OUTPUT JISUAN2# 1#IO将计算电路JISUAN2#寄存器中的数据输送到计算机1#端口
DISPLAY ON计算机显示器液晶显示计算得到的方程组的根X1,X2,X3,X4
X1,X2,X3,X4
DISPLAY OFF 停止显示
PROGRAM END 程序结束
计算过程如下所示：
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
a a a a
11 12 13 14
a a a a
D= 21 22 23 24
a a a a
31 32 33 34
a a a a
41 42 43 44
b a a a
1 12 13 14
b a a a
D = 2 22 23 24
1
b a a a
3 32 33 34
b a a a
4 42 43 44
a b a a
11 1 13 14
a b a a
D = 21 2 23 24
2
a b a a
31 3 33 34
a b a a
41 4 43 44
a a b a
11 12 1 14
a a b a
D = 21 22 2 24
3
a a b a
31 32 3 34
a a b a
41 42 4 44
a a a b
11 12 13 1
a a a b
D = 21 22 23 2
4
a a a b
31 32 33 3
a a a b
41 42 43 4
于是方程组唯一的解为
D
1
x = =3
1 D
D
2
x = =-4
2 D
D
3
x = =1
3 D
D
4
x = =1
4 D
用数字电路加法器，乘法器，除法器，各种数值表示上面的函数，
第五部分 线性方程组的误差计算电路
下面的资料可参见《数值分析》，美国J.B.斯卡勃罗著，陈盖民，余介石译，科学出版社1958年出版
例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。
计算线性方程组的解的误差程序：
PROGRAM BEGIN 程序开始
JISUAN2#
选择保存计算数据到JISUAN2#寄存器，同时标志寄存器保存2，说明今天进行了2次计算
DISPLAY ON 计算机液晶显示器显示线性方程组
3.21x-4.86y=5.73
2.13x+8.63y=12.65
DISPLAY OFF
HLSWCJS ERYUAN FANGCHENG B
与门导通，选择线性方程组的克莱姆法则计算电路
MAKEB a11=3.21,a12=-4.86,b1=5.73, 将方程式的系数用a11,a12,b1存储器存储起来
MAKEC a21=2.13,a22=-8.63,b2=12.65，
可以将键盘输入的误差系数M=0.005,代入公式进行计算
GET WUCHA M=0.005
RUN HLSWCJS 与门导通，接通行列式误差计算电路，开始计算
SAVE JISUAN2# ZU保存计算中所有JISUAN2#电路寄存器的数据
OUTPUT JISUAN2# 1#IO 将计算电路JISUAN1#寄存器中的数据输送到计算机1#端口
DISPLAY ON 计算机显示器显示计算得到的方程组的根X，Y,△X,△Y,
X，Y,△X,△Y,
DISPLAY OFF 停止显示
PROGRAM END 程序结束
用数字电路加法器，乘法器，除法器，各种数值表示上面的函数，
第六部分 插值法
下面的资料可参见《数值分析》，美国J.B.斯卡勃罗著，陈盖民，余介石译，科学出版社1958年出版。
差分、插值法的牛顿公式
. 15.引言
插值法曾经被定义为：在表格的各行之间读数的一种技术；在初等数学中，这个名词常常用来表示一种过程，从函数的一组给定的或列表的数值当中，来计算函数的中间值。然而插值法的一般问题，无论如何要比这个广泛得多。在高等数学中，我们常常会处理一些函数，它们的解析形式，或者是完全不知道的，或者是具有某种性质（复杂性或其他）以至难以照所需的运算进行。不论那一种情况，都需要把给定函数换成另一个很容易处理的函数。这种把一个给定函数用另一个简单的函数来代替或表达的运算，便构成了插值法这个名词的广义含义。
注：作者认为这种用简单函数去代替复杂者的过程，可以叫做“解析置换法”（analytical replacement)或“函数置换法”（functional substitution)原理。因此，插值法的一般问题乃是：
借助于函数自变量取各个定值时所得的一组函数给定值，去表达一个已知的或未知的函数，使之成为预先选定的形式。
譬如，令y=f(x)表示这样一个函数，它在自变量x取x ,x ,x ,…,x 等值时，
0 1 2 n
具有y ,y ,y ,…,y 等各函数值；
0 1 2 n
又作一个较简单的任意函数φ(x)，使它在x ,x ,x ,…,x 各点，
0 1 2 n
与f(x)具有相同的值，这时若在给定的全区间上，把f(x)用φ(x)来代替，那么这个过程便构成了插值法，而函数φ(x)便是一个插值法公式。函数φ(x)可以有各种各样的形式。当φ(x)是一个多项式时，用φ(x)表达f(x)的过程，称为抛物线的或多项式的插值法，当φ(x)是一个有限项三角级数时，这过程便是三角的插值法。同样，φ(x)可以是指数函数，勒让得多项式，贝塞尔函数等等级数。在实用问题中，我们常常选择一个最简单的函数作为φ(x)，来表达在所考虑区间上的给定函数。由于多项式是最简单的函数，所以我们常常采用多项式作为φ(x)，并且几乎所有的插值标准公式都是多项式的公式。如果给定函数具有周期性，那么最好用三角级数去表达它。要证明用多项式或三角级数去代替一个给定函数是正确的这件事，就有赖于维尔斯德拉斯在1885年所证明的两个定理，注：一个实变数的所谓任意函数的解析表达法：这两个定理的证明，见菲赫金哥茨(Фвггевгольд)微积分学教程第二十章，709页（1955年中译本三卷三分册，597-600页）。
这两个定理可以叙述如下：
Ⅰ.任何一个在区间(a,b)内连续的函数，可以用一个多项式在该区间内表达到任何需要的准确程度；即能够找到一个多项式P(x)，对区间(a,b)内任何x值，都能令│f(x)-P(x)│<ε成立，这里ε是任何预先制定的正量；
Ⅱ.以2π为周期的连续函数，都能用一个具有
g(x)=a +a sinx+a sin2x+…+a sinnx+b cosx+b cos2x+…+b cosnx
0 1 2 n 1 2 n
形式的有限项三角级数把它表达出来；即对于所考虑的区间内的一切x值都有│f(x)-g(x)│<δ的关系，这里δ是任何预先制定的正量。
这两个定理的几何意义是：在画出y=f(x),y=f(x)+ε及y=f(x)-ε的图形之后，便可找出一个多项式或有限项三角级数，不论ε多么小，对所有在a与b之间的x值，都能使它的图形始终保持在y=f(x)+ε和y=f(x)-ε所画成的区域之内（见图1）。所以这些定理的意义就是说：一个给定函数可以被多项式或有限项三角级数所代替而且能达到任何需要的准确程度。
16.差分
倘若用y ,y ,y ,…,y 表示任何y=f(x)的一组函数值，
0 1 2 n
那么y -y ,y -y ,y -y ,…,y -y 就叫做函数y的一阶差分，
1 0 2 1 3 2 n n-1
将这些差分用△y ,△y ,△y 等等表示，则得
0 1 2
△y =y -y ， △y =y -y ，…,△y =y -y ，△y =y -y ，
0 1 2 1 2 1 n-1 n n-1 n n+1 n
2 2
这些一阶差分之间的差，叫作二阶差分，它们用△ y ,△ y 等等表示，即
0 1
2
△ y =△y -△y =y -2y +y
0 1 0 2 1 0
2
△ y =△y -△y =y -2y +y
1 2 1 2 2 0
…………………………………
下面的差分表表明所有各阶差分是如何形成的：
表1.对角差分表
x y
△y 2
△ y 3
△ y 4
△ y 5
△ y 6
△ y 7
△ y 8
△ y
x
0 y
0
△y
0
x
1 y
1 2
△ y
0
△y
1 3
△ y
0
x
2 y
2 2
△ y
1 4
△ y
0
△y
2 3
△ y
1 5
△ y
0
x
3 y
3 2
△ y
2 4
△ y
1 6
△ y
0
△y
3 3
△ y
2 5
△ y
1 7
△ y
0
x
4 y
4 2
△ y
3 4
△ y
2 6
△ y
1 8
△ y
0
△y
4 3
△ y
3 5
△ y
2 7
△ y
1
x
5 y
5 2
△ y
4 4
△ y
3 6
△ y
2
△y
5 3
△ y
4 5
△ y
3
x
6 y
6 2
△ y
5 4
△ y
4
△y
6 3
△ y
5
x
7 y
7 2
△ y
6
x
8 y
8 △y
7
这个表叫作对角差分表，大部分的差分表都属于这一种，但是对很多用途来说，还是用一种更紧凑的表，名叫水平差分表的比较合适。在水平差分表中，不同阶的差分是以下标去代替指数而表示的。
m m
注：也有用记号▽来表示的，即以▽ y 代替本书的△ y
n n
我们可以把前面的差分表，使用水平差分表的记号重新写成水平形式如下：
表2.水平差分表
x y △y △ y
2 △ y
3 △ y
4 △ y
5 △ y
6 △ y
7 △ y
8
x
0 y
0
x
1 y
1 △ y
1 1 △ y
2 0
x
2 y
2 △ y
1 2 △ y
2 2
x
3 y
3 △ y
1 3 △ y
2 3 △ y
3 3
x
4 y
4 △ y
1 4 △ y
2 4 △ y
3 4 △ y
4 4
x
5 y
5 △ y
1 5 △ y
2 5 △ y
3 5 △ y
4 5 △ y
5 5
x
6 y
6 △ y
1 6 △ y
2 6 △ y
3 6 △ y
4 6 △ y
5 6 △ y
6 6
x
7 y
7 △ y
1 7 △ y
2 7 △ y
3 7 △ y
4 7 △ y
5 7 △ y
6 7 △ y
7 7
x
8 y
8 △ y
1 8 △ y
2 8 △ y
3 8 △ y
4 8 △ y
5 8 △ y
6 8 △ y
7 8 △ y
8 8
为了要看出在水平差分表与对角差分表之间，同阶差分的关系，我们在表3和表4中，给出以y表示的两类差分。由这两个表中可以看出，在顶上的对角线两者是相同的。而在表3中，底部向上斜的对角线是和表4底部的水平线相同的。还可以举个例来说，由表3得
3
△ y =y -3y +3y -y
1 4 3 2 1
同样从表4可得
△ y =y -3y +3y -y
3 4 4 3 2 1
因此,
3
△ y =△ y
1 3 4
浏览一下表3和表4便可看出：附上指数和下标的△之间的一般关系为
m
△ y =△ y （从y 起向前数）
k m k+m k
或 m
△ y =△ y （从y 起向前数）
m n n-m n
其中m表示差分的阶，而k和n则表示列表数值的次序。表3.对角差分表
y
△y 2
△ y 3
△ y 4
△ y 5
△ y
y
0
y -y
1 0
y
1 y -2y +y
2 1 0
y -y
2 1 y -3y +3y -y
3 2 1 0
y
2 y -2y +y
3 2 1 y -4y +6y -4y +y
4 3 2 1 0
y -y
3 2 y -3y +3y -y
4 3 2 1 y -5y +10y -10y +5y -y
5 4 3 2 1 0
y
3 y -2y +y
4 3 2 y -4y +6y -4y +y
5 4 3 2 1
y -y
4 3 y -3y +3y -y
5 4 3 2 y -5y +10y -10y +5y -y
6 5 4 3 2 1
y
4 y -2y +y
5 4 3 y -4y +6y -4y +y
6 5 4 3 2
y -y
5 4 y -3y +3y -y
6 5 4 3 y -5y +10y -10y +5y -y
7 6 5 4 3 2
y
5 y -2y +y
6 5 4 y -4y +6y -4y +y
7 6 5 4 3
y -y
6 5 y -3y +3y -y
7 6 5 4 y -5y +10y -10y +5y -y
8 7 6 5 4 3
y
6 y -2y +y
7 6 5 y -4y +6y -4y +y
8 7 6 5 4
y -y
7 6 y -3y +3y -y
8 7 6 5
y
7 y -2y +y
8 7 6
y -y
8 7
y
8
表4.水平差分表
y △y △ y
2 △ y
3 △ y
4 △ y
5
y
0
y
1 y -y
1 0
y
2 y -y
2 1 y -2y +y
2 1 0
y
3 y -y
3 2 y -2y +y
3 2 1 y -3y +3y -y
3 2 1 0
y
4 y -y
4 3 y -2y +y
4 3 2 y -3y +3y -y
4 3 2 1 y -4y +6y -4y +y
4 3 2 1 0
y
5 y -y
5 4 y -2y +y
5 4 3 y -3y +3y -y
5 4 3 2 y -4y +6y -4y +y
5 4 3 2 1 y -5y +10y -10y +5y -y
5 4 3 2 1 0
y
6 y -y
6 5 y -2y +y
6 5 4 y -3y +3y -y
6 5 4 3 y -4y +6y -4y +y
6 5 4 3 2 y -5y +10y -10y +5y -y
6 5 4 3 2 1
y
7 y -y
7 6 y -2y +y
7 6 5 y -3y +3y -y
7 6 5 4 y -4y +6y -4y +y
7 6 5 4 3 y -5y +10y -10y +5y -y
7 6 5 4 3 2
y
8 y -y
8 7 y -2y +y
8 7 6 y -3y +3y -y
8 7 6 5 y -4y +6y -4y +y
8 7 6 5 4 y -5y +10y -10y +5y -y
8 7 6 5 4 3
对角差分与水平差分之间的关系，可以用一个数值的例子作进一步的说明。下表表示，当函数y=sinhx取一组等间隔值时，它的两类差分值。
x y
△y 2
△ y 3
△ y 4
△ y 5
△ y 6
△ y
1.5 2.12928
0.24629
1.6 2…37557 0.002377
0.27006 0.000271
1.7 2.64563 0.002648 0.000026
0.29654 0.000297 0.000003
1.8 2.94217 0.002945 0.000029 0.000001
0.32599 0.000326 0.000004
1.9 3.26816 0.003271 0.000033
0.35870 0.000359
2.0 3.62686 0.003630
0.39500
2.1 4.02186
x y △y △ y
2 △ y
3 △ y
4 △ y
5 △ y
6
1.5 2.12928
1.6 2…37557 0.24629
1.7 2.64563 0.27006 0.002377
1.8 2.94217 0.29654 0.002648 0.000271
1.9 3.26816 0.32599 0.002945 0.000297 0.000026
2.0 3.62686 0.35870 0.003271 0.000326 0.000029 0.000003
2.1 4.02186 0.39500 0.003630 0.000359 0.000033 0.000004 0.000001
可以看出：不论写成对角差分还是水平差分，七个函数值的差分值是相同的。在某些工作中水平差分具有显著的优点。例如在微分方程的数值解中，在后面的函数值常常是已知的，而在前面的常常是未知的。此时水平差分表中，在同一条线上所列的各阶差分，就好像是函数的最后已知值，而这些差分正是求函数的次一计算值时所要用的。对这种情况来说，水平型差分表要比对角型来得更紧凑更方便一些。在另一方面来说，我们在下一章开始研究中心差分插值法时，就会觉得对角差分表对那个用途来说又是特别好的。
17.在列表数值中一个误差的影响
令y ,y ,y ,…,y 表示函数真值，并假定y 的值受到误差ε的影响，
0 1 2 n 5
因而它的错误值为y +ε，因此，这些y的逐次差分，有如下所示：
5
表5.在列表数值中，一个误差所引起的影响
y
△y 2
△ y 3
△ y 4
△ y
y
0
△y
0
y
1 2
△ y
0
△y
1 3
△ y
0
y
2 2
△ y
1 4
△ y
0
△y
2 3
△ y
1
y
3 2
△ y
2 4
△ y +ε
1
△y
3 3
△ y +ε
2
y
4 2
△ y +ε
3 4
△ y -4ε
2
△y +ε
4 3
△ y -3ε
2
y +ε
5 2
△ y -2ε
4 4
△ y +6ε
3
△y -ε
5 3
△ y +3ε
4
y
6 3
△ y +ε
5 4
△ y -4ε
4
△y
6 3
△ y -ε
5
y
7 2
△ y
6 4
△ y +ε
5
△y
7 3
△ y
6
y
8 2
△ y
7 4
△ y -4ε
6
△y
8 3
△ y
7
y
9 2
△ y
8
△y
9
y
10
此表说明：一个误差的影响随着逐次差分而增大，而ε的系数是具有正负相同的二项式系数，并且在任何一纵列的差分中，误差的代数和是0. 它还表明：在差分中误差最大误差是和错误列表值同在一条水平线上。下表系表明水平差分表中，一个误差的影响：这里误差的影响与前表相同，但在此表中任何阶的第一个错误差分都和错误列表值同在一条水平线上。
表6
y △ y
1 △ y
2 △ y
3 △ y
4
y
0
y
1 △ y
1 1
y
2 △ y
1 2 △ y
2 2
y
3 △ y
1 3 △ y
2 3 △ y
3 3
y
4 △ y
1 4 △ y
2 4 △ y
3 4 △ y
4 4
y +ε
5 △ y +ε
1 5 △ y +ε
2 5 △ y +ε
3 5 △ y +ε
4 5
y
6 △ y -ε
1 6 △ y -2ε
2 6 △ y -3ε
3 6 △ y -4ε
4 6
y
7 △ y
1 7 △ y +ε
2 7 △ y +3ε
3 7 △ y +6ε
4 7
y
8 △ y
1 8 △ y +ε
2 8 △ y -ε
3 8 △ y -4ε
4 8
y
9 △ y
1 9 △ y
2 9 △ y
3 9 △ y +ε
4 9
y
10 △ y
1 10 △ y
2 10 △ y
3 10 △ y +ε
4 10
这里误差的影响与前表相同，但在此表中任何阶的第一个错误差分都和错误列表值同在一条水平线上。按照差分表中一个误差蔓延的规律，便使得我们可以循着误差的踪迹而追究它的来源，并将它改正。让我们来考虑下面的表，作为在列表函数中检查并改正误差过程的一个实例。
x y △ y
1 △ y
2 △ y
3 △ y
4 ε
0.10 0.09983
0.15 0.14944 0.04961
0.20 0.19867 0.04923 -0.00038
0.25 0.24740 0.04873 -0.00050 -0.00012
0.30 0.29552 0.04812 -0.00061 -0.00011 0.00001
0.40 0.38945 0.04655 -0.00083 -0.00009 0.00004
0.45 0.43497 0.04552 -0.00103 -0.00020 -0.00011 ε
0.50 0.47943 0.04446 -0.00106 -0.00003 0.00017 -4ε
0.55 0.52269 0.04326 -0.00120 -0.00014 -0.00011 6ε
0.60 0.56464 0.04195 -0.00131 -0.00011 0.000043 -4ε
0.65 0.60519 0.04055 -0.00140 -0.00009 0.00002
0.70 0.64422 0.03903 -0.00150 -0.00012 -0.00003
这里，三阶差分在靠近直列中央的部分是很不规则的，而四阶差分更加不规则。在各直列中，这种不规则性，都从对应于x=0.10的水平线上开始。由于四阶差分的代数和为1，故四阶差分的平均值仅仅约为第五位小数的0.1个单位。因此这个例题中所求的四阶差分大部分是积聚的误差。现在参照表6我们可得
-4ε=-11，6ε=17等等.
因为ε接近于3个单位。
由于(y +ε)-ε=y ,所以与x=0.40相对应的y的真值为0.38945-0.00003=0.38942。
k k
现在便可以校正各个差分的直列了，并且将会发现三阶差分实际上是常数。倘若函数的几个列表值都受到误差的影响，则函数的逐次差分将变成不规则的，但此时要确定各个误差的来源和大小，却不是一件容易的事。当表中每个y值受到大小为ε的误差的影响时，
注：这句话的意思是指y 受有错误ε 的影响。
k k
由表3和表4可以明显的看出，各个三阶差分所受到误差的影响为ε -3ε +3 -ε ，
k k-1 k-2 k-3
各个四阶差分受到误差的影响为ε -4ε +6ε -4ε +ε ，
k k-1 k-2 k-3 k-4
以下可依次类推。在实用问题中，函数y的列表数值时由量测或计算而获得的，因之它们就容易受到量测误差的影响，或由于将计算结果抹尾至给定位数而产生误差的影响。不论那一种情况，都会依取差分的过程而增大，并且只是它们，已足够引起高阶差分变成不规则的。
注：对于函数列表值中误差的详细讨论，见莱斯：插值法的理论和实践（Rice:Theory and Practice of Interpolation),第7-15页及46-47页。又见皮尔曼：数学差近方法讲义(O.Biermann:Voriesungun uber Matematiuche Nahe rungsmethoden),138页等。
因为这个缘故，我们通常不宜使用四阶以上的差分。
18.差分和导数的关系
差分和导数的关系有时是需要知道的。它们之间的基本关系是
n n (n)
△ f(x)=(△x) f (x+θn△x),0<θ<1 (18.1)
和
n
△ f(x) n
lim = f (x)
△x→0 n
(△x)
这些关系的推导，见瓦莱-浦逊的“无穷小解析教程”（Vallee-Poussins Cours dAnalyse Infinitesimale),Ⅰ（第四版，1921），第72-73页，现在加以证明如次：
注：这里补充的证明是译自Franklin`s A Treatise on Advanced Calculus,第93节。
设函数f(x)在全部区间x,x+n△x上有n阶导数，而它的（n-1)阶导数在这区间的端点是连续的。对于n=1，关系式（18.1）成立，它即是拉格朗日（Lagrange)中值定理。如此我们只能从1，2，…，n-1的情形来证明它对n为真，那就是用数学归纳法证实这个结果。如f(x)的n阶和n-1阶导数满足所需要的条件，那么函数
△f(x)=f(x+△x)-f(x) (1)
也会如此，由此，它的(n-1)阶导数在区间内存在，而它的(n-2)阶导数在端点连续，照归纳法上的假设，可用n-1代替n，来对函数△f(x)应用方程(18.1)这样得
n n-1 n-1
△ f(x)=△ f(x+△x)-△ f(x)
n-1 (n-1) (n-1)
=△x {f [x+△x+θ`(n-1)△x]-f [x+θ`(n-1)△x]} (2)
这里θ是介于0和1之间的一个适当数值。由于f(x)的导数方面的条件，表明方程(2)的右端在开区间内有一阶导数，而函数本身在闭区间内连续。如此便可对括号的差，应用拉格朗日中值定理，得到下面的结果： n (n) △x f [x+θ△x+θ(n-1)△x] (3) 这里的θ``也和θ一样，是介于0与1之间的一个适当数值，因此有
θ``△x+θ`(n-1)△x=θn△x (4)
而θ也是介于0与1之间的一个适当数值。
n
把这方程和从△ f(x)得到的（3）式结合，我们求得
n n (n)
△ f(x)=△x f (x+θn△x) (5)
这是对于n的结果，因而用归纳法完成了证明。特别指出，如果函数f(x)在闭区间内有连续的n阶导数，就可满足条件，因而（18.1）成立。其次，假设函数f(x)对于一个特别值有一个有限的n阶导数，那么，在这值的某邻域内。它有(n-1)阶导数，因而有连续的(n-2)阶导数。由此，对于△x充分小的值，可以推出方程（2）。
(n-1)
但是，因为函数f(x)在x有n阶导数，故f (x)在x是可微的，故有
(n-1) (n-1)
f [x+△x+θ`(n-1)△x]-f (x)
(n)
=[1+θ`(n-1)]△x[f (x)+α] (6)
与
(n-1) (n-1)
f [x+θ`(n-1)△x]-f (x)
(n)
=θ`(n-1)△x[f (x)+α`] (7)
其中α和α`都是无穷小，跟着△x同趋于0. 由这两个方程和方程（3），就有
n n (n)
△ f(x)=△x [f (x)+β],
即,
n
△ f(x) (n)
lim = f (x)
△x→0 n
(△x)
这里β是一个无穷小。
19.多项式的差分
现在让我们来计算n次多项式的逐次差分，我们有
n n-1 n-2
y=f(x)=ax +bx +cx +…+kx+L (1)
n n-1 n-2
y+△y=a(x+h) +b(x+h) +c(x+h) +…+k(x+h)+L (2)
这里h=△x, 从(2)减去(1)，得到
n n n-1 n-1 n-2 n-2
△y=a[(x+h) -x ]+b[(x+h) -x ]+c[(x+h) -x ]+…+kh
n n-1
按二项式定理展开(x+h) ,(x+h) 等各量，我们有
n n-1 n(n-1) 2 n-2 n(n-1)(n-2) 3 n-3 n
△y=a[x +nhx + h x + h x +…-x ]+
2 3!
n-1 n-2 (n-1)(n-2) 2 n-3 n-1
+b[x +(n-1)hx + h x +…-x ]+
2
n-1 n-2 (n-2)(n-3) 2 n-4 n-2
+c[x +(n-2)hx + h x +...-x ]+...+kh
2
即
n-1 2 n(n-1) n-2
△y=anhx +[ah +b(n-1)h]x +
2
n(n-1)(n-2) 2 (n-1)(n-2) n-3
+[ah +bh +ch(n-2)]x +...
3！ 2
n-2 n-3
现在如果△x(=h)是常数，各个括号内x ,x 等的系数便是常数，可以用单独的常系数b,c来代替它们。因此我们有
n-1 n-2 n-3
△y=anhx +bx +cx +…+kx+L (3)
于是，n次多项式的一阶差分是另一个n-1次的多项式。求二阶差分时，可在(3)中给予x以增量△x=h而有
n-1 n-2 n-3
△y+△(△y)=anh(x+h) +b(x+h) +c(x+h) +…+k(x+h)+L (4)
从(4)减去(3)，我们得
2 n-1 n-1 n-2 n-2 n-3 n-3
△(△y)=△ y=anh[(x+h) -x ]+b[(x+h) -x ]+c[(x+h) -x ]+…+k`h
n-1 n-2 n-3 n-4
用二项式定理展开(x+h) ,(x+h) 等等各量且照旧用一个字母来代替x ,x 等量的常数系数，我们便有
2 2 n-2 n-3 n-4
△ y=an(n-1)h x +bx +cx +…+kx+L
于是，一个n-2次多项式的二阶差分是一个n-2次多项式。照此继续计算，得到的n阶差分是一个零次多项式，即
n n n-n n 0 n
△ y=a[n(n-1)(n-2)…1]h x =an!h x =ah n!
所以n次多项式的n阶差分是常数，而更高阶的差分都是零。读者必须记住：这个结果仅当h(即△x)为常数时，即仅当x值取算术级数时方始成立。我们刚才证明的命题，可以叙述如下：
当n次多项式的自变量按算术级数取值时，即取等距值时，它的n阶差分是常数。它的逆命题也是真的，即
当自变量的值取算术级数时，如果列表函数的n阶差分为常数，则该函数为n次多项式。
注：这一命题的证法见莱斯：插值法的理论和实践，24页。
这里的第二个命题使得我们能够把任何一个某阶差分为常数或接近于常数的函数，可以用一个多项式取代替它。譬如在17节中的列表函数就可以用三次多项式去表达它，因为经过校正后的三阶差分近似的是常数。以上的两个命题的较简的论证如下。
注：这个证法根据别席考维奇(ЯСБези]ович)的近似计算法(1949年版),35节，见中译本107-108页。
取n次多项式
n n-1 n-k
f(x)=a x +a x +…+a x +…+a
0 1 k n
(k)
那么，它的k阶导数f (x)在k≤n时是n-k次，在k>n时都是0.由泰勒公式有
2 n
△x △x n
△f(x)=f(x+△x)-f(x)=△xf(x)+ f``(x)+...+ f (x) 2! n! 因此，一个n次多项式的一阶差分是一个n-1次多项式。这差分的最高次项即是△xf(x)中的最高次项，因而等于
n-1
△x*na x
0
这就是说，n次多项式的一阶差分中最高次项，等于用△x去乘多项式最高次项的导数。由此便可写出n次多项式△f(x)各阶差分的最高次项如下：
2 n-1
△x *n(n-1)a x 二阶的
0
3 n-2
△x *n(n-1)(n-2)a x 三阶的
0
……………………………………………………
k n-2
△x *n(n-1)(n-k+1)a x k阶的
0
………………………………………………………
k 0 n
△x n(n-1)…1a x =△x n!a n阶的
0 0
这样便证明了上面的正命题。至于逆命题可证明如次：
由于关系式（18.1），可知
n n (n)
△ f(x)=△x f [x+θn△x]=c
即
(n) c
f [x+θn△x]=
n
△x
与x的值无关，而是一个常数，因此f(x)是一个n次多项式。
20.牛顿向前插值公式
我们的下一个问题是找出合适的多项式，它能在给定区间上代替任何的给定函数。令y=f(x)表示一个函数它在自变量x取等距的值x ,x ,x ,…,x 时
0 1 2 n
具有y ,y ,y ,…,y 诸值，
0 1 2 n
再令φ(x)表示一个n次多项式。该多项式可写成下面的形式。
φ(x)=a +a (x-x )+a (x-x )(x-x )+a (x-x )(x-x )(x-x )+a (x-x )(x-x )(x-x )(x-x )
0 1 0 2 0 1 3 0 1 2 4 0 1 2 3
+…+a (x-x )(x-x )(x-x )…(x-x ) (1)
n 0 1 2 n-1
现在我们来确定系数a ,a ,a ,…,a 等，以便能使
0 1 2 n
φ(x )=y ,φ(x )=y ,…,φ(x )=y ,
1 1 2 2 n n
将x的各值x ,x ,x ,…,x 逐次代入(1),同时命φ(x )=y
0 1 2 n 0 0
φ(x )=y ,…,并且记住x -x =h,x -x =2h,…,则得
1 1 0 1 2 0
y =a 或a =y
0 0 0 0
y =a +a (x -x )=y +a h
1 0 1 1 0 0 1
y -y △y
1 0 0
a = =
1 h h
y =a +a (x -x )+a (x -x )(x -x )
2 0 1 2 0 2 2 0 2 1
y -y
1 0
=y + +a (2h)(h)
0 h 2
2
y -2y +y △ y
2 1 0 0
a = =
2 2 2
2h 2h
y =a +a (x -x )+a (x -x )(x -x ) +a (x -x )(x -x )(x -x )
3 0 1 3 0 2 3 0 3 1 2 3 0 3 1 3 2
y +y y -2y +y
1 0 2 1 0
=y + (3h)+ (3h)(2h)+a (3h)(2h)h
0 h 2 3
2h
2
y -3y +3y -y △ y
3 2 1 0 0
a = =
3 3 3
6h 3!h
y =a +a (x -x )+a (x -x )(x -x ) +a (x -x )(x -x )(x -x )
4 0 1 4 0 2 4 0 4 1 2 4 0 4 1 4 2
+a (x -x )(x -x )(x -x )(x -x )
2 4 0 4 1 4 2 4 3
y +y y -2y +y
1 0 2 1 0
=y + (4h)+ (4h)(3h)+
0 h 2
2h
y -3y +3y -y
3 2 1 0
(4h)(3h)(2h)h+a (4h)(3h)(2h)h
3 4
6h
4
y -4y +6y -4y +y △ y
4 3 2 1 0 0
a = =
4 4 4
4!h 4!h
继续进行这种计算系数的方法，我们便可求出
5
△ y
0
a = ,
5 5
5!h
6
△ y
0
a = ,…,
6 6
6!h
n
△ y
0
a = ,
n n
n!h
将这些a ,a ,a ,…,a 各值代入(1)，则得
0 1 2 n
2
△y △ y
0 0
φ(x)=y + (x-x )+ (x-x )(x-x )+
0 h 2 0 1
2h
3
△ y
0
+ (x-x )(x-x )(x-x )+,
3 0 1 2
3!h
4
△ y
0
+ (x-x )(x-x )(x-x ) (x-x )+…+
4 0 1 2 3
4!h
n
△ y
0
+ (x-x )(x-x )(x-x ) (x-x )…(x-x ) (2)
n 0 1 2 3 n-1
n!h
这就是用x表示的牛顿向前插值公式。此公式可以使用变数更换法加以简化，首先我们把（2）写成下列的等价形式：
2
(x-x ) △ y (x-x ) (x-x )
0 0 0 1
φ(x)=y +△y + +
0 0 h 2 h h
3
△ y (x-x ) (x-x ) (x-x )
0 0 1 2
+ +
3 h h h
4
△ y (x-x ) (x-x ) (x-x ) (x-x )
0 0 1 2 3
+ +… (3)
3 h h h h
现在令
x-x
0
=u
h
即, x=x +hu
0
于是，由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值之用。
0
始点y 可以是任何一个列表值，而y仅含有选作始点的后继y值。
0
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
2h
继续进行这种计算系数的方法，我们便可得到
△ 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
2h
△ 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
n
u=
h
或
x=x +hu
u
于是由
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 hx-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 0u(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值之用。
n
现在我们举一些例题来说明牛顿公式的用法。
例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,故
0
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）+
2
=0.4970679364+0.000819310+0.0000000053=0.4971498727
此结果准确到它的末位数字。
注：像在上面的例题中，这种在给定诸值的范围以外计算函数值的过程称为外插法。使用它时应该慎重从事，不过假设已经知道函数在靠近给定诸值范围的端点是平滑的变动的，并且假如h取得尽量的小，那么我们在给定诸值范围以外相距h的地方进行外插，通常是很稳妥的。
例3.在1918年1月1日亮每小时的赤纬已在下表给出，试求在北京时间3时35分15秒时的赤纬。
小时 赤纬 △ y
1 △ y
2 △ y
3
0 8°2953.7`` 1 8°1819.4 -11`34.2
2 8°643.5`` -1135.9 -1.6
3 7°556.1`` -1137.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，故
n
t-t 8
n -0时24分45秒 -1485
u= = = =-0.3569
h h 8
3600
u+1=0.6431
将这些值代入（Ⅱ）并用δ表示所需求的赤纬，则得
(0.6431)(-0.3569）
δ=7°4327``.2+(-0.3569)(-1138.9）+ (-1.5)
h
=7°4327``.2+49.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°4327``.2+(1)(-1138.9)+ (-1.5)=7°3146``.8 2 在美国天文历与航海历(American Ephemeris and Nautical Almanac)中所给的真值为7°3146.9, 因此外插值的误差只不过0.1而已。
第三章 插值法 中心差分公式
22.引言
牛顿公式（Ⅰ）和（Ⅱ）是基本的，并且对于几乎所有的插值情况来说，都是可以应用的，不过就一般而论它不如另一类公式所谓中心差分公式收敛得快。后一类公式所用的差分，都尽可能取在通过对角差分表所作水平线的附近。浏览表3，可以看出：这些差分所含的各函数值，都在通过该函数值所在水平线的上下。所以中心差分公式特别适合于内插靠近列表值中央部分的函数值之用。最重要的中心差分公式是两个已知的所谓斯特灵(Stirling)公式和贝赛尔公式。
注：中心差分公式还有高斯的一种。但是柯姆利(J.L.Comrie)曾指出，它实际上很少用到。皮尔逊(K.Pearson)也指摘这个公式，杰佛里(Jeffreys)认为这公式确是无用，但它与牛顿公式有直接关系，最易推导，由它就很容易推出其他公式。关于高斯公式以及由它证明斯特灵公式和贝赛尔公式的过程，可参看别席考维奇同书第43到45节（中译本125-129页）。推导它们的方法有好几种，不过最简单的一种是由牛顿公式（Ⅰ）用代数变换法来推导的。
23.斯特灵插值公式
要推导斯特灵公式，我们首先写出对角差分表，
然后把列表值y ,以及跟通过y 所作的水平线紧挨着的那些差分，
0 0
都标记出来，以便加以特别注意。这些量在表7中，系用黑体字印出。
以y 为始点的牛顿公式（Ⅰ）写为
0
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
△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
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
2
各种记法间的关系显见是：
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
©m = , (d)m =
5 2 7 2
…………………………….
所以这些m就表示：
跟通过y 所作水平线紧挨着的一些奇阶差分的算术平均值。
0
2 3
我们目前的目标，是把△y ,△ y ,△ y ,等等，
0 0 0
要用各个m以及位于通过y 水平线上的各个偶阶差分表示出来，
0
2
这件事可采用消去法来做，即从△y ,△ y 等等起顺着对角线向上的进行到右边去，
0 0
直至到达水平线上的各量时为止。为了对这件事的进行有所帮助，
2 4 6 8
对于偶阶差分△ y ,△ y ,△ y ,△ y 等等，
-1 -2 -3 -4
不论在以下任何代数运算的地方出现，都在它们的下面画线，画线的目的是为了提醒我们下面画线的量是不消去的。从差分的定义，我们可得
2
△ y =△y -△y
-1 0 -1
所以,
2
△y =△ y -△y
0 -1 -1
但从(a)可知
△y =2m -△y
-1 1 0
所以,
2
△y =△ y +2m -△y
0 -1 1 0
所以,
2
△y =m +△ y /2 (1)
0 1 -1
2
要想求出所需量表示的△ y 的值，则有
0
注：读者可能认为这些解法过程相当麻烦，
2
也不容易看出与上文“从△y ,△ y 等等起顺着对角线向上的进行到右边去，
0 0
直至到达水平线上的各量时为止”一句话的关系。
2 3
现在用图解来说明△ y 与△ y 的求法如下：
0 0
（到水平线上为止）
3
△ y
-2
2 4
△ y △ y
-1 -2
3
△ y
-1
2
△ y
0
沿着对角线向上
3 5
△ y △ y
-2 -3
4 6
△ y △ y
-2 -3
3 5
△ y △ y
-1 -2
4
△ y
-1
3
△ y 3
0 从△ y 沿着对角线向上
0
这里有两个三数关系，即
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
3
便得到了△ y 。这种图解说明，不难使它一般化。
0
2n 2
在求△ y 时，有n +n个三数关系n个二数关系。
0
2n 2
图解中除去△ y 外，有n +3n个量，但有n+1个是消不去的。
0
2 2 3 2
因此，恰好从 n +n+n=n +2n个关系中，消去n +3n-(n+1）=n +2n-1个量，
2n+1 2
在求△ y 时，有(n+1) 个三数关系，n+1个二数关系。
0
2n+1 2
除去△ y 外，n +4n+2个量，但n+1个是不消去的。
0
2 2 2 2
因此，恰好从(n+1) +(n+1)=n +3n+2个关系中，消去n +4n+2-(n+1)=n +3n+1个量。
注释结束。
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
并且自©已知
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
§ △ 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
® △ y =△ y -△ y
-4 -3 -4
并且，自(d)可知
7 7
(s) △ y =2m -△ y
-3 7 -4
7
自(s)减去®，并解出△ 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)代入§，并解出△ 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
©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
然而，从©有
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
-
2 2 △ y +△ y+ △ y + 2 2 5! -2 4 4
(v -1/4)(v -9/4)(v -25/4) -2 -1 -
2 2 2 (2n-1) 2n 2n+…+ 6! 2 2
(v -1/4)(v -9/4)…[v - ] △ y +△ y
2 -n -n+1 -
2 2 2 (2n-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
-
1809439004689185 27012817923189185*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!
©若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
©e=1+1+ + + +…+ +
2! 3! 4! (n-1)! n!
这里,
θ
e
R = ;
n n!
但由于e<3而θ≤1，所以显然有
θ
3
(d) R = ;
n n!
关于R 更确定的公式可求之如下：
n
把级数©写到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准确到小数任意多少位的值，在展开式©中所必需取的项数。譬如，我们打算用级数©求准到小数十位的e值，那么我们可由方程
1
=0.00000000005
(n-1)(n-1)!
求出n。
靠阶乘倒数表的帮助，便可求出n-1=13或n=14。所以我们应该在级数©中取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
©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不会很大。
