数值微分的计算方法
数值微分的计算方法
内容摘要 求解数值微分问题, 就是通过测量函数在一些离散点上的值, 求得函数的近似导数。本文就所学知识，归纳性地介绍了几种常用的数值微分计算方法。并举例说明计算，实验结果表明了方法的有效性。
关键词 数值微分 Taylor展开式 Lagrange插值 三对角矩阵
引言：数值微分即根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。常见的可以用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值，由此也可导出多点数值微分计算公式。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。
1.Taylor 展开式方法
理论基础：Taylor 展开式
f (x )=f (x 0)+(x -x 0)f '(x 0)
(x -x 0)+
2!
2
f ''(x 0)
(x -x 0)+ +
n !
n
n
f ()(x 0)+
我们借助Taylor 展开式，可以构造函数f (x ) 在点x =x 0的一阶导数和二阶导数的数值微分公式。取步长h >0则
h 2' '
f (x 0+h ) =f (x 0) +hf (x 0) +f (ξ1)
2
'
ξ1∈(x 0, x 0+h ) （1）
所以
f ' (x 0) =
f (x 0+h ) -f (x 0) h ' '
-f (ξ1)
h 2
ξ1∈(x 0, x 0+h ) （2）
同理
h 2' '
f (x 0-h ) =f (x 0) -hf (x 0) +f (ξ2)
2
'
ξ2∈(x 0-h , x 0) （3） ξ2∈(x 0-h , x 0) （4）
f ' (x 0) =
f (x 0) -f (x 0-h ) h ' '
+f (ξ2)
h 2
'
式(2)和式(4)是计算f Taylor 展开式多写几项
(x 0) 的数值微分公式，其截断误差为O (h ) ，为提高精度，将
ξ1∈(x 0, x 0+h ) ξ2∈(x 0-h , x 0)
h 2' ' h 3' ' ' h 4(4)
f (x 0+h ) =f (x 0) +hf (x 0) +f (x 0) +f (x 0) +f (ξ1)
2624
'
h 2' ' h 3' ' ' h 4(4)
f (x 0-h ) =f (x 0) -hf (x 0) +f (x 0) -f (x 0) +f (ξ2)
2624
两式相减得
'
f (x 0+h ) -f (x 0-h ) h 2' ' '
f (x 0) =-f (x 0) +O (h 4) （5）
2h 6
'
上式为计算f ' (x 0) 的微分公式，其截断误差为O(h ), 比式（2）和(4)精度高。
2
两式相加，如果f
(4)
(x ) ∈C [x 0-h , x 0+h ]，则有
f (x 0+h ) -2f (x 0) +f (x 0-h ) h 2(4)
f (x 0) =-f (ξ) ξ∈(x -h , x +h )
00h 212 , (6)
''
式(6)是计算f ' ' (x 0) 的数值微分公式, 其截断误差为O (h 2) 。
例1．设函数f (x ) =ln x , x 0=2, h =0.1，试用数值微分公式计算f ' (2) 的值。 解 由式(2)、式(4)和式(5)分别计算结果为
f ' (2) ≈f ' (2) ≈f ' (2) ≈
ln 2. 1-ln 2
≈0. 4879
0. 1ln 2-ln 1. 9
≈0. 5129
0. 1
ln 2. 1-ln 1. 9
≈0. 5004
0. 2
与真值f
'
(2) =0.5相比，式(5)计算的结果精度较高。
2. 数值微分的Lagrange 插值方法
设函数y =f (x ) 具有n +1个实验数据：(x i , f (x i ) ) (i =0,1, 2, , n ) ，我们希望估计
f ' (x ) 的值，特别x =x i 时，估计f ' (x i ) 的值。
基于插值方法的数值微分做法是，由已知(x i , f (x i ) ) (i =0,1, , n ) 建立Lagrange 插值多项式或Newton 插值多项式（这里以Lagrange 插值方法为例），即
f (x ) =L n (x ) +R n (x )
于是
f ' (x )=L ' n (x )+R ' n (x )
当x =x i 时，有
'
f ' (x i ) =L ' n (x i ) +R n (x i )(i =0,1,2,..., n )
其中
f (n +1) (ξ) '
R n (x i ) =ωn +1(x i )
(n +1)!
'
略去误差项有
f ' (x i ) ≈L ' n (x i )
实际运用中，等距节点更为常见。设
h =
b -a
, x i =a +ih (i =0,1, 2, , n ) , x =a +th ， n
于是有
t t (t -1) 2t (t -1)...(t -n +1) n
f (x ) =y 0+∆y 0+∆y 0+... +∆y 0+R n (x )
1! 2! n !
所以
f ' (x i ) =
d t t (t -1)...(t -n +1) n
{y 0+∆y 0+... +∆y 0}t =i dt 1! n !
h n f n (n +1) (ξ) n
+∏(i -j )
(n +1)! j =0
j ≠i
3. 多点数值微分公式
由于高阶插值的不稳定性，实际应用时多采用n=1,2,4的两点、三点和五点等多点插值型求导公式。
(1) 两点公式(n=1)
1⎧'
f (x ) =(y 1-y 0) 0⎪h ⎨1
⎪f ' (x 1) =(y 1-y 0)
h ⎩
（2）三点公式(n=2)
R ' (x 0) =-
h ' '
f (ξ) 2 (7) h
R ' (x 1) =f ' ' (ξ)
2
⎧' 1h 2(3) '
(-3y 0+4y 1-y 2) R (x 0) =f (ξ) ⎪f (x 0) =2h 3⎪
1h 2(3) ⎪' （8） '
(-y 0+y 2) R (x 1) =-f (ξ) ⎨f (x 1) =2h 6⎪
h 2(3) ' ⎪f ' (x ) =1(y -4y +3y ) R (x 0) =f (ξ) 2012
⎪2h 3⎩
(3)五点公式（n =4）
⎧' 1f (x ) =(-25y 0+48y 1-36y 2+16y 3-3y 4) 0⎪
12h ⎪
⎪f ' (x ) =1(-3y -10y +18y -6y +y )
101234
⎪12h ⎪' 1
(y 0-8y 1+8y 3-y 4) ⎨f (x 2) =12h ⎪
1⎪'
f (x ) =(-y 0+6y 1-18y 2+10y 3+3y 4) 3⎪12h ⎪1⎪f ' (x 4) =(3y 0-16y 1+36y 2-48y 3+25y 4)
12h ⎩h 4(5)
R (x 0) =f (ξ)
5
h 4(5) '
R (x 1) =-f (ξ)
20
（9） h 4(5) '
R (x 2) =-f (ξ)
30h 4(5) '
R (x 3) =-f (ξ)
20h 4(5) '
R (x 0) =f (ξ)
5
'
. 05，分别用三点公式和五点公式计算f ' (2) 的近似值。例2．设f (x ) =ln x ，取h =0
解:由式(8)有
f ' (2) =
1
(-3f (2) +4f (2. 05) -f (2. 10)) =0. 499802861 2⨯0. 05
1
(-f (1. 95) +f (2. 05)) =0. 500104205 2⨯0. 05
f ' (2) =f ' (2) =
1
(f (1. 90) -4f (1. 95) +3f (2)) =0. 499779376 2⨯0. 05
由式(9)有
f ' (2) =
1
(f (1. 90) -8f (1. 95) +6f (2. 05) -f (2. 10)) =0. 499999843
12⨯0. 05
与真值f ' (2) =0.5相比，三点公式已有相当满意精度，而五点公式的结果是十分满意的。 4. 数值微分的隐式格式
前述的Taylor 展式法和数值微分格式均称为显式格式，即直接由已知的
f (x i ) (i =0,1,2, , n ) ，经过适当的算术四则运算，立即可得f ' (x i ) 的近似值。显式格式
优点是计算方便，工作量小，缺点是数值不稳定。为克服后一缺点，隐式格式常常具有数值稳定性。
数值微分的隐式格式建立方法常用通过Taylor 展开式方法或数值积分方法等不同途径。
首先我们用Taylor 展开式方法来推导数值微分的隐式格式。
由式(5)和式(6)，我们用x k 代替x 0得：
f ' (x k ) =
' '
1h 2' ' '
[f (x k +h ) -f (x k -h )]-f (x k ) +O (h 4) 2h 6
1h 2(4)
f (x k ) =2[f (x k +h ) -2f (x k ) +f (x k -h )]-f (ξ)
h 12
所以也有
1' h 2(5) ' '
f (x k ) =2[f (x k +h ) -2f (x k ) +f (x k -h )]-f (ξ)
h 12
' ' '
将最后f
'"
(x k ) 表达式代入f ' (x k ) 表达式可得
f ' (x k ) =
11
[f (x k +h ) -f (x k -h )]-[f ' (x k +h ) -2f ' (x k ) +f ' (x k -h )]+O (h 4) 2h 6
略去误差项O (h 4) ，并且m k 表示f ' (x k ) 的近似值，则有
3
[f (x k +1) -f (x k -1)](k =1, 2,..., n ) h
同样，用数值积分方法也可推导数值微分的隐式格式。
m k -1+4m k +m k +1=
⎰
b
x k +1
x k -1
f (x ) dx =⎰
x k +1
x k -1
f ' (x k ) dx
(k =1, 2,..., n )
根据simpson 公式⎰f (x ) dx =
b -a a +b
[f (a ) +4f () +f (b )]有
a 62
2h
f (x k +1) -f (x k -1) =[m k -1+4m k +m k +1],(k =1, 2,..., n )
63
m k -1+4m k +m k +1=[f (x k +1) -f (x k -1)](k =1, 2,..., n )
h
上式是关于n +1个未知量m 0, m 1, m n 的n -1个方程，如果已知m 0=
3
f ' (x 0) , m n =f ' (x n ) ，记d k =[f (x k +1) -f (x k -1)]，故有方程组：
h
⎡41⎤⎡m 1⎤⎡d 1-m 0⎤
⎥⎢141⎥⎢m ⎥⎢d
22⎥⎢⎥⎢⎥⎢
⎢141⎥⎢m 3⎥⎢d 3⎥
⎥=⎢⎥ （10） ⎢⎥⎢ ⎥⎢⎥⎢⎥⎢
⎢141⎥⎢m n -2⎥⎢d n -2⎥
⎥⎢⎥⎢⎥⎢m d -m 14⎢⎥⎢⎥⎢⎣⎦⎣n -1⎦⎣n -1n ⎥⎦
这就是求m 1, m 2, , m n -1的线性方程组。由于系数矩阵是严格对角占优的三对角矩阵，因此非奇异，解存在且唯一，可由追赶法求解，且数值稳定。
例3．设函数f (x ) =ln x 在节点x =1.5, 1.6, 1.7, 1.8, 1.9, 上的函数值，及
f ' (1.5) , f ' (2.0) 如表1所示。试用数值微分隐式格式式(10)和式(11)求出相应节点上一阶
导数的近似值。
表1 例3节点数据
解：由式(10)有
-0. 666666667⎡41⎤⎡m 1⎤⎡3. 75489429⎤
⎢141⎥⎢m ⎥⎢⎥3. 533491052⎢⎥⎢⎥=⎢⎥ ⎢141⎥⎢m 3⎥⎢⎥3. 33676905⎢⎥⎢⎥⎢⎥
m 143.16081554-0. 500000000⎣⎦⎣4⎦⎣⎦
解方程组得：
m 1=0.[1**********]483
m 2=0.[1**********]067 m 3=0.[1**********]249
m 4=0.[1**********]938
与f
'
(x )=
1
在各相应节点上数值相比，上述结果约有五位有效数字，精度较高。h 越小，x
精度越高。 附追赶法程序： %追赶法
%定义三对角矩阵A 的各组成单元。方程为Ax=d
% a为主对角线下面的次对角线，b 为主对角线，c 为主对角线上面的次对角线，d 为等式右端向量。 % A=[4 1 0 0 % 1 4 1 0 % 0 1 4 1 % 0 0 1 4]
function x=zhuiganfa(a,b,c,d) a=[0 1 1 1];c=[1 1 1];b=[4 4 4 4];
d=[3.75489429-0.666666667,3.53349105,3.33676905,3.16081554-0.5]; n=length(b); u0=0;y0=0;a(1)=0; L(1)=b(1)-a(1)*u0; y(1)=(d(1)-y0*a(1))/L(1); u(1)=c(1)/L(1); for i=2:(n-1)
L(i)=b(i)-a(i)*u(i-1); y(i)=(d(i)-y(i-1)*a(i))/L(i); u(i)=c(i)/L(i); end
L(n)=b(n)-a(n)*u(n-1); y(n)=(d(n)-y(n-1)*a(n))/L(n); x(n)=y(n); for i=(n-1):-1:1
x(i)=y(i)-u(i)*x(i+1); end
运行结果： ans =
0.[1**********]483 0.[1**********]067 0.[1**********]249 0.[1**********]938
>> x=1.5:0.1:2.0; y=log(x);
dy=[1/1.5 0.[1**********]483 0.[1**********]067 0.[1**********]249 0.[1**********]938 1/2.0]; plot(x,y,'b-') hold on
plot(x,1./x,'r-o') hold on
plot(x,dy,'g-*')
作图为：
后序：除了本文介绍的常用方法，还有其它许多比较成熟的方法计算数值微分，可以提高稳定近似导数的收敛率。
参考文献
[1]施吉林, 刘淑珍, 陈桂芝. 计算机数值方法(第三版)[M]. 高等教育出版社.2009 [2]吴天毅. 数值积分中点公式的改进[J]. 天津理工大学学报, 2007, (03) . [3]宫浩. 数值微分算法研究及应用[D]黑龙江大学, 2009
[4]尹秀玲, 闫立梅. 一种数值微分方法[J]德州学院学报, 2007,(02)