徐會(huì)林，郜星軍

(河南理工大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，河南 焦作454000)

0 前言

數(shù)值微分，即由給定函數(shù)的測(cè)量數(shù)據(jù)近似求其導(dǎo)數(shù)，這類問(wèn)題在科學(xué)研究及工程實(shí)踐中有廣泛的應(yīng)用價(jià)值。如在磁共振電阻抗成像技術(shù)(MREIT)中，介質(zhì)的導(dǎo)電率是利用生物組織內(nèi)部電流產(chǎn)生的磁場(chǎng)信息重建的?；谥鞔艌?chǎng)方向磁感應(yīng)強(qiáng)度的調(diào)和Bz算法是當(dāng)前主流的導(dǎo)電率重建算法，而如何由磁感應(yīng)強(qiáng)度的測(cè)量數(shù)據(jù)計(jì)算其二階導(dǎo)數(shù)是關(guān)鍵步驟［1］，數(shù)值微分的具體應(yīng)用還體現(xiàn)在數(shù)字圖像特征檢測(cè)問(wèn)題［2］、期權(quán)定價(jià)問(wèn)題［3］、放射性廢棄物的存儲(chǔ)問(wèn)題［4］等方面。

數(shù)值微分問(wèn)題的求解，除標(biāo)準(zhǔn)的正則化方法外［5～8］，還有一些其特有的方法，如有限差分法［9］、平滑化方法［10］等。眾所周知，微分和積分是互逆的2種基本運(yùn)算。積分的計(jì)算通常要借助微分實(shí)現(xiàn)，反過(guò)來(lái)，也可利用積分求微分，其實(shí)質(zhì)是將微分運(yùn)算等價(jià)轉(zhuǎn)化為第一類積分方程的求解問(wèn)題［11，12］。

本文將基于微分與積分之間的互逆關(guān)系，構(gòu)造數(shù)值微分的實(shí)現(xiàn)算法。通過(guò)把數(shù)值微分問(wèn)題等價(jià)轉(zhuǎn)化為第一類積分方程的求解問(wèn)題，給出了一元函數(shù)前兩階導(dǎo)數(shù)的近似計(jì)算方法，并通過(guò)數(shù)值算例說(shuō)明了解的數(shù)值有效性。

1 一階數(shù)值微分算法

設(shè)一元函數(shù)f(x)∈C1［0，1］，u(x)為其一階導(dǎo)函數(shù)，即u(x)=f'(x)∈C［0，1］，則函數(shù)u和f滿足第一類的Volterra型積分方程

不失一般性，設(shè)f(0)=0，則有

此時(shí)，求導(dǎo)數(shù)u的問(wèn)題就轉(zhuǎn)化為積分方程(2)的求解問(wèn)題［11］。在實(shí)際問(wèn)題中，函數(shù)f(x)的表達(dá)式一般是未知的，已知的只是其在某些離散點(diǎn)上的取值。此時(shí)，導(dǎo)數(shù)u或等價(jià)的積分方程(2)的求解需引入數(shù)值算法。為此，引入數(shù)值積分公式將方程(2)左端的積分項(xiàng)進(jìn)行離散，將積分方程的求解近似為線性方程組的求解，這就是求解積分方程的數(shù)值積分法［13］。

首先將區(qū)間［0，1］進(jìn)行n等分，得n+1個(gè)離散節(jié)點(diǎn)xi=ih，i=0，1，…，n，其中h=1/n為離散步長(zhǎng)。為方便記號(hào)，將f(x)和u(x)在離散節(jié)點(diǎn)上的取值分別記為fi=f(xi)，ui=u(xi)。當(dāng)x= x1時(shí)，引入梯形求積公式可得

當(dāng)x=xi，2≤i≤n時(shí)，引入復(fù)化梯形求積公式可得

聯(lián)立式(3)、式(4)可得如下線性方程組

當(dāng)u0=u(0)=f'(0)已知時(shí)，求解可得

實(shí)際計(jì)算時(shí)，u0往往是未知的，它的求解需借助其它方法，如利用向前差分格式求解可得u0= f(x1)/h。

當(dāng)x=xi，2≤i≤n時(shí)，引入復(fù)化中點(diǎn)求積公式可得

聯(lián)立式(7)、式(8)可得如下線性方程組

求解線性方程組(9)可得一階導(dǎo)數(shù)在中點(diǎn)處的近似值。

2 二階數(shù)值微分算法

設(shè)f(x)∈C2［0，1］，u(x)為其二階導(dǎo)數(shù)，即u (x)=f″(x)∈C［0，1］。不失一般性，總是可以假設(shè)f(0)=f(1)=0，否則f(x)－f(0)+x(f(0)－f (1))滿足該假設(shè)且與f(x)有相同的二階導(dǎo)數(shù)。因此有

求解該常微分方程可知函數(shù)f(x)及其導(dǎo)數(shù)u(x)滿足如下關(guān)系［14］

方程(11)為第一類的Fredholm型積分方程，引入復(fù)化梯形求積公式可得，當(dāng)x=xi，0≤i≤n時(shí)，

當(dāng)i取0，1，…，n時(shí)，可得n+1個(gè)方程，聯(lián)立得如下線性方程組

注意到x0=0，xn=1，所以左端系數(shù)矩陣的第一行、第一列、最后一行以及最后一列的元素均為零，將方程組(13)簡(jiǎn)化可得

求解可得二階導(dǎo)數(shù)在離散節(jié)點(diǎn)處的近似值。

3 數(shù)值實(shí)驗(yàn)

在前兩節(jié)中，分別給出了基于積分的一階及二階數(shù)值微分算法。在本節(jié)中，將通過(guò)算例說(shuō)明上述算法的數(shù)值有效性。

算例1:已知函數(shù)f(x)=x4－2x2+x，x∈［0，1］，求其前兩階導(dǎo)數(shù)。

圖1 n=40時(shí)，分別利用式(5)和式(9)求得的一階導(dǎo)數(shù)的近似值uT，uM以及精確的一階導(dǎo)數(shù)f'的圖像

表1 算例1中等分?jǐn)?shù)n取值不同時(shí)，一階導(dǎo)數(shù)的近似值uM和uT的誤差水平

下面，利用式(14)計(jì)算二階導(dǎo)數(shù)的近似值。在圖2中，給出了當(dāng)n=40時(shí)，利用式(14)求得的二階導(dǎo)數(shù)的近似值u及精確值f″的圖像。從圖2中可以看出，二階導(dǎo)數(shù)的近似值與精確導(dǎo)數(shù)是基本吻合的。進(jìn)一步，在表2中給出了二階導(dǎo)數(shù)的近似值的誤差水平，從中可以看出近似解的計(jì)算是二階收斂的。

算例2:已知函數(shù)f(x)=sin(5πx)，x∈［0，1］求其前兩階導(dǎo)數(shù)。

圖2 n=40時(shí)，利用式(14)求得的二階導(dǎo)數(shù)的近似值u及精確的二階導(dǎo)數(shù)f″的圖像

首先，分別利用式(5)和式(9)計(jì)算一階導(dǎo)數(shù)的近似值。在圖3中，給出了當(dāng)n=40時(shí)，分別利用式(5)和式(9)求得的一階導(dǎo)數(shù)的近似值uT，uM以及精確的一階導(dǎo)數(shù)f'的圖像。從中可以看出，近似解uM在數(shù)值計(jì)算中的優(yōu)越性。進(jìn)一步，在圖4中給出了當(dāng)n=40時(shí)，利用式(14)求得的二階導(dǎo)數(shù)的近似值u及精確的二階導(dǎo)數(shù)f″的圖像。從中可以看出，二階導(dǎo)數(shù)的近似值與精確導(dǎo)數(shù)是基本吻合的。近似導(dǎo)數(shù)的誤差分析與算例1是類似的，在此不再贅述。

表2 算例2中等分?jǐn)?shù)n取值不同時(shí)，二階導(dǎo)數(shù)近似值u的誤差水平

圖3 n=40時(shí)，分別利用式(5)和式(9)求得的一階導(dǎo)數(shù)的近似值uT，uM以及精確的一階導(dǎo)數(shù)f'的圖像

圖4 n=40時(shí)，利用式(14)求得的二階導(dǎo)數(shù)的近似值u及精確的二階導(dǎo)數(shù)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］ 王澤文，溫榮生.一階和二階數(shù)值微分的Lanczos方法［J］.高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，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］.北京:清華大學(xué)出版社，2012.

［14］ 徐會(huì)林，王澤文.二階數(shù)值微分的積分方程方法［J］.江西科學(xué)，2009，27(1):108－112.

［15］ 孫志忠，袁慰平，聞?wù)鸪?數(shù)值分析(第二版)［M］.南京:東南大學(xué)出版社，2002.