基于积分的数值微分算法

徐会林，郜星军

(河南理工大学数学与信息科学学院，河南 焦作454000)

0 前言

数值微分，即由给定函数的测量数据近似求其导数，这类问题在科学研究及工程实践中有广泛的应用价值。如在磁共振电阻抗成像技术(MREIT)中，介质的导电率是利用生物组织内部电流产生的磁场信息重建的。基于主磁场方向磁感应强度的调和Bz算法是当前主流的导电率重建算法，而如何由磁感应强度的测量数据计算其二阶导数是关键步骤［1］，数值微分的具体应用还体现在数字图像特征检测问题［2］、期权定价问题［3］、放射性废弃物的存储问题［4］等方面。

数值微分问题的求解，除标准的正则化方法外［5～8］，还有一些其特有的方法，如有限差分法［9］、平滑化方法［10］等。众所周知，微分和积分是互逆的2种基本运算。积分的计算通常要借助微分实现，反过来，也可利用积分求微分，其实质是将微分运算等价转化为第一类积分方程的求解问题［11，12］。

本文将基于微分与积分之间的互逆关系，构造数值微分的实现算法。通过把数值微分问题等价转化为第一类积分方程的求解问题，给出了一元函数前两阶导数的近似计算方法，并通过数值算例说明了解的数值有效性。

1 一阶数值微分算法

设一元函数f(x)∈C1［0，1］，u(x)为其一阶导函数，即u(x)=f'(x)∈C［0，1］，则函数u和f满足第一类的Volterra型积分方程

不失一般性，设f(0)=0，则有

此时，求导数u的问题就转化为积分方程(2)的求解问题［11］。在实际问题中，函数f(x)的表达式一般是未知的，已知的只是其在某些离散点上的取值。此时，导数u或等价的积分方程(2)的求解需引入数值算法。为此，引入数值积分公式将方程(2)左端的积分项进行离散，将积分方程的求解近似为线性方程组的求解，这就是求解积分方程的数值积分法［13］。

首先将区间［0，1］进行n等分，得n+1个离散节点xi=ih，i=0，1，…，n，其中h=1/n为离散步长。为方便记号，将f(x)和u(x)在离散节点上的取值分别记为fi=f(xi)，ui=u(xi)。当x= x1时，引入梯形求积公式可得

当x=xi，2≤i≤n时，引入复化梯形求积公式可得

联立式(3)、式(4)可得如下线性方程组

当u0=u(0)=f'(0)已知时，求解可得

实际计算时，u0往往是未知的，它的求解需借助其它方法，如利用向前差分格式求解可得u0= f(x1)/h。

当x=xi，2≤i≤n时，引入复化中点求积公式可得

联立式(7)、式(8)可得如下线性方程组

求解线性方程组(9)可得一阶导数在中点处的近似值。

2 二阶数值微分算法

设f(x)∈C2［0，1］，u(x)为其二阶导数，即u (x)=f″(x)∈C［0，1］。不失一般性，总是可以假设f(0)=f(1)=0，否则f(x)－f(0)+x(f(0)－f (1))满足该假设且与f(x)有相同的二阶导数。因此有

求解该常微分方程可知函数f(x)及其导数u(x)满足如下关系［14］

方程(11)为第一类的Fredholm型积分方程，引入复化梯形求积公式可得，当x=xi，0≤i≤n时，

当i取0，1，…，n时，可得n+1个方程，联立得如下线性方程组

注意到x0=0，xn=1，所以左端系数矩阵的第一行、第一列、最后一行以及最后一列的元素均为零，将方程组(13)简化可得

求解可得二阶导数在离散节点处的近似值。

3 数值实验

在前两节中，分别给出了基于积分的一阶及二阶数值微分算法。在本节中，将通过算例说明上述算法的数值有效性。

算例1:已知函数f(x)=x4－2x2+x，x∈［0，1］，求其前两阶导数。

图1 n=40时，分别利用式(5)和式(9)求得的一阶导数的近似值uT，uM以及精确的一阶导数f'的图像

表1 算例1中等分数n取值不同时，一阶导数的近似值uM和uT的误差水平

下面，利用式(14)计算二阶导数的近似值。在图2中，给出了当n=40时，利用式(14)求得的二阶导数的近似值u及精确值f″的图像。从图2中可以看出，二阶导数的近似值与精确导数是基本吻合的。进一步，在表2中给出了二阶导数的近似值的误差水平，从中可以看出近似解的计算是二阶收敛的。

算例2:已知函数f(x)=sin(5πx)，x∈［0，1］求其前两阶导数。

图2 n=40时，利用式(14)求得的二阶导数的近似值u及精确的二阶导数f″的图像

首先，分别利用式(5)和式(9)计算一阶导数的近似值。在图3中，给出了当n=40时，分别利用式(5)和式(9)求得的一阶导数的近似值uT，uM以及精确的一阶导数f'的图像。从中可以看出，近似解uM在数值计算中的优越性。进一步，在图4中给出了当n=40时，利用式(14)求得的二阶导数的近似值u及精确的二阶导数f″的图像。从中可以看出，二阶导数的近似值与精确导数是基本吻合的。近似导数的误差分析与算例1是类似的，在此不再赘述。

表2 算例2中等分数n取值不同时，二阶导数近似值u的误差水平

图3 n=40时，分别利用式(5)和式(9)求得的一阶导数的近似值uT，uM以及精确的一阶导数f'的图像

图4 n=40时，利用式(14)求得的二阶导数的近似值u及精确的二阶导数f″的图像

［1］ Seo J K，Woo E J.Magnetic resonance electrical impedance tomography(MREIT)［J］.SIAM Review，2011，53(1):40－68.

［2］ Demigny D.On optimal linear filtering for edge detection［J］.IEEE Transactions on Image Processing，2002，11(7):728－737.

［3］ Hein T，Hofmann B.On the nature of ill-posedness of an inverse problem arising in option pricing［J］.Inverse Problems，2003，19(6):1319－1338.

［4］ Mohankumar N，Auerbach S M.On the use of higherorder formula for numerical derivatives in scientific computing［J］.Communications in Computational Physics，2004，161(3):109－118.

［5］ Wang Y B，Jia X Z，Cheng J.A numerical differentiation method and its application to reconstruction of discontinuity［J］.Inverse Problems，2002，18(6):1461－1476.

［6］ Wei T，Hon Y C，Wang Y B.Reconstruction of numerical derivatives from scattered noisy data［J］.Inverse Problems，2005，21(2):657－672.

［7］ Wang Z Z，Wen R S.Numerical differentiation for high orders by an integration method［J］.Journal of Compu-tational and Applied Mathematics，2010，234(3):941－948.

［8］ 王泽文，温荣生.一阶和二阶数值微分的Lanczos方法［J］.高等学校计算数学学报，2012，34(2): 160－178.

［9］ Lu S，Pereverzev S V.Numerical differentiation from a viewpoint of regularization theory［J］.Mathematics of Computation，2006，75(256):1853－1870.

［10］ Murio D A，Mejia C E，Zhan S.Discrete mollification and automatic numerical differentiation［J］.Computers and Mathematics with Applications，1998，35(5): 1－16.

［11］ Ahn S，Choi U J，Ramm A G.A scheme for stable numerical differentiation［J］.Journal of Computational and Applied Mathematics，2006，186(2):325－334.

［12］ Xu H L，Liu J J.Stable numerical differentiation for the second order derivatives［J］.Advances in Computational Mathematics，2010，33(4):431－447.

［13］ 沈以淡.积分方程(第三版)［M］.北京:清华大学出版社，2012.

［14］ 徐会林，王泽文.二阶数值微分的积分方程方法［J］.江西科学，2009，27(1):108－112.

［15］ 孙志忠，袁慰平，闻震初.数值分析(第二版)［M］.南京:东南大学出版社，2002.