令鋒

（肇慶學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣東 肇慶 526061）

1 問題的提出

數(shù)值微分是指將函數(shù)y=f(x)在給定點(diǎn)a的導(dǎo)數(shù)用點(diǎn)a附近節(jié)點(diǎn)上函數(shù)值的線性組合近似表示[1-2].科學(xué)研究與工程技術(shù)領(lǐng)域中的許多問題都需要計(jì)算函數(shù)f(x)在給定點(diǎn)的導(dǎo)數(shù)，當(dāng)f(x)表達(dá)式復(fù)雜或函數(shù)僅由表格給出時(shí)，就需要進(jìn)行數(shù)值微分，中心差商公式是數(shù)值微分采用的一種重要算法.

記h為步長(zhǎng)，f(x)在節(jié)點(diǎn)a-h，a及a+h處的函數(shù)值分別為f(a-h)，f(a)和f(a+h)，則求f(x)在點(diǎn)a處一階導(dǎo)數(shù)的中點(diǎn)公式及截?cái)嗾`差分別為

求函數(shù)f(x)在點(diǎn)a處二階導(dǎo)數(shù)的三點(diǎn)公式及截?cái)嗾`差分別為

對(duì)公式（1），從截?cái)嗾`差角度分析，h越小計(jì)算結(jié)果越準(zhǔn)確；但從舍入誤差或計(jì)算穩(wěn)定性角度考慮，h越小，f(a-h)與f(a+h)的值將越接近，兩個(gè)相近數(shù)相減會(huì)導(dǎo)致有效位數(shù)嚴(yán)重?fù)p失，因而用中心差商公式（1）計(jì)算一階導(dǎo)數(shù)時(shí)步長(zhǎng)不宜太小或太大，存在最優(yōu)步長(zhǎng)，當(dāng)h取值小于最優(yōu)步長(zhǎng)后，計(jì)算誤差反而會(huì)增大[2-4].同理，采用公式（3）計(jì)算二階導(dǎo)數(shù)時(shí)也存在最優(yōu)步長(zhǎng).

中心差商公式（1）及（3）都具有二階計(jì)算精度，且公式簡(jiǎn)單，編程方便，是數(shù)值微分值得采用的方法.但由于步長(zhǎng)很小時(shí)計(jì)算結(jié)果對(duì)步長(zhǎng)變化敏感，大多數(shù)教材沒有進(jìn)一步討論兩點(diǎn)公式與三點(diǎn)公式的計(jì)算機(jī)實(shí)現(xiàn)，個(gè)別教材中雖然提出宜采用變步長(zhǎng)方法計(jì)算，但沒有給出詳細(xì)的分析、說明和示例，尤其是沒有給出計(jì)算終止條件[5]，因而編程實(shí)現(xiàn)算法存在困難.本文通過分析確定中心差商公式中最優(yōu)步長(zhǎng)面臨的問題，提出以計(jì)算結(jié)果滿足精度要求作為計(jì)算終止條件和以步長(zhǎng)最接近最優(yōu)步長(zhǎng)作為計(jì)算終止條件兩種算法，以期能夠采用中心差商公式的變步長(zhǎng)算法求得函數(shù)在給定點(diǎn)具有較高精度的一階及二階導(dǎo)數(shù).

2 中心差商公式最優(yōu)步長(zhǎng)的確定

設(shè)f(a)，f(a-h)和f(a+h)的舍入誤差分別為ε0，ε1和ε2，ε=max{|ε0|,|ε1|,|ε2|}，則計(jì)算f′(a)與f′′(a)的舍入誤差分別為

根據(jù)求極值的方法，要使誤差 取最小值，h應(yīng)滿足

同理可得，使誤差E2(h)達(dá)到最小值的最優(yōu)步長(zhǎng)為

求出M1與M2的工作較繁瑣，在某些情形下甚至無法實(shí)現(xiàn)，從而確定最優(yōu)步長(zhǎng)面臨諸多困難.為獲得具有較高精確度的計(jì)算結(jié)果，通常采用中心差商公式變步長(zhǎng)算法.以計(jì)算一階導(dǎo)數(shù)為例，首先選擇一個(gè)初始步長(zhǎng)h0，為方便計(jì)算，通常取h0=1，利用公式（1）求出G(h0)，然后將步長(zhǎng)減半得h1，再計(jì)算G(h1)，如此反復(fù)，使步長(zhǎng)在逐步減半的過程中逐步接近最優(yōu)步長(zhǎng)，從而獲得精度較高的計(jì)算結(jié)果.因而，在變步長(zhǎng)中心差商微分法的計(jì)算機(jī)實(shí)現(xiàn)過程中，確定正確的計(jì)算終止條件是算法得以實(shí)現(xiàn)的前提.

3 中心差商公式變步長(zhǎng)算法的計(jì)算終止條件

設(shè)對(duì)初始步長(zhǎng)h0經(jīng)n次減半后得到的步長(zhǎng)為h0，函數(shù)y=f(x)在點(diǎn)a相應(yīng)的導(dǎo)數(shù)f′(a)及f″(a)的近似值分別為G(hn)和T(hn)，以下給出兩種計(jì)算終止條件.

3.1 滿足精度要求時(shí)終止計(jì)算

如果給定了計(jì)算精度要求ε，則前后兩次步長(zhǎng)減半得到的導(dǎo)數(shù)近似值應(yīng)滿足如下條件

此即給定精度要求時(shí)變步長(zhǎng)數(shù)值微分法的計(jì)算終止條件.

上述算法容易編程實(shí)現(xiàn)且避免了直接計(jì)算最優(yōu)步長(zhǎng)，可望獲得理想的計(jì)算結(jié)果.但由于步長(zhǎng)過小時(shí)誤差反而會(huì)更大，當(dāng)要求的精度很高時(shí)，通過將步長(zhǎng)逐步減半的方法可能無法獲得滿足精度要求的計(jì)算結(jié)果.實(shí)際計(jì)算時(shí)，為避免舍入誤差不斷增大出現(xiàn)較大誤差數(shù)值錯(cuò)誤結(jié)果，可預(yù)設(shè)一個(gè)步長(zhǎng)最大減半次數(shù)M，當(dāng)步長(zhǎng)減半次數(shù)超過最大減半次數(shù)時(shí)停止計(jì)算.

若以計(jì)算結(jié)果滿足精度要求作為計(jì)算終止條件，用三點(diǎn)公式（3）計(jì)算函數(shù)f(x)在點(diǎn)a的二階導(dǎo)數(shù)算法可描述如下：

算法1：數(shù)值微分的變步長(zhǎng)三點(diǎn)公式算法

1） 輸入函數(shù)f(x)以及點(diǎn)a的值.

2）輸入計(jì)算精度要求ε及步長(zhǎng)最大減半次數(shù)M.

3） 步長(zhǎng)減半次數(shù)n=0,h=2-n=1，計(jì)算T0(f(a-h)-2f(a)+f(a+h))∕(h2).

3） 步長(zhǎng)減半：n=1,h=2-n=1∕2，計(jì)算T1(f(a-h)-2f(a)+f(a+h))∕(h2).

4）如果 |T1-T0|＞ε，轉(zhuǎn)步驟5）；否則，轉(zhuǎn)步驟6）.

5）如果n≤M,則n=n+1，h=2-n,T0=T1,T1(f(a-h)-2f(a)+f(a+h))∕(h2)，轉(zhuǎn)步驟4）；否則，輸出出錯(cuò)信息，轉(zhuǎn)步驟7）.

6）輸出計(jì)算結(jié)果T1.

7）結(jié)束.

3.2 步長(zhǎng)最接近最優(yōu)步長(zhǎng)時(shí)終止計(jì)算

在步長(zhǎng)逐步減半并接近最優(yōu)步長(zhǎng)hopt過程中，計(jì)算精度將越來越高，而步長(zhǎng)小于hopt后，繼續(xù)對(duì)步長(zhǎng)減半數(shù)值計(jì)算的誤差將會(huì)不斷增大.因此，經(jīng)n-1次步長(zhǎng)減半后得到的步長(zhǎng)hn-1與最優(yōu)步長(zhǎng)hopt若滿足如下關(guān)系

相應(yīng)的導(dǎo)數(shù)近似值則滿足如下條件

此即步長(zhǎng)最接近最優(yōu)步長(zhǎng)時(shí)變步長(zhǎng)數(shù)值微分法的計(jì)算終止條件,G(h)為步長(zhǎng)為hn-1時(shí)f′(a)的近似值.

上述算法既可避免直接確定最優(yōu)步長(zhǎng)，又可獲得具有較高計(jì)算精度的計(jì)算結(jié)果，因而是一種非常實(shí)用的有效算法.

若以步長(zhǎng)最接近最優(yōu)步長(zhǎng)作為計(jì)算終止條件，用變步長(zhǎng)中點(diǎn)公式計(jì)算函數(shù)f(x)在點(diǎn)a的一階導(dǎo)數(shù)的算法可描述如下：

算法2：數(shù)值微分的變步長(zhǎng)中點(diǎn)公式算法

1）輸入函數(shù)f(x)，輸入點(diǎn)a的值.

2）步長(zhǎng)減半次數(shù)n=0，h=2-n=1，計(jì)算G0=[f(a+h)-f(a-h)]∕(2h).

3）步長(zhǎng)減半：n=0，h=2-n=1∕2，計(jì)算G1=[f(a+h)-f(a-h)]∕(2h).

4）n=n+1，h=2-n，計(jì)算Gn=[f(a+h)-f(a-h)]∕(2h).

5）如果 |Gn)-Gn-1|≥ |Gn-1-Gn-2|，轉(zhuǎn)步驟6）；否則，轉(zhuǎn)步驟4）.

6）輸出計(jì)算結(jié)果Gn-1.

7）結(jié)束.

4 說明性算例

算例1用變步長(zhǎng)三點(diǎn)公式計(jì)算函數(shù)f(x)=ex在輸入的點(diǎn)x處的二階導(dǎo)數(shù)近似值，要求誤差不超過0.5×10-6及0.5×10-9，并統(tǒng)計(jì)步長(zhǎng)減半次數(shù).

解 根據(jù)算法1編寫C語言程序生成Threepoints.com文件，計(jì)算結(jié)果如下：

請(qǐng)輸入要計(jì)算二階導(dǎo)數(shù)值的點(diǎn)a:1

請(qǐng)輸入精度要求epsilon:0.0000001

點(diǎn)a=1.000000處的二階導(dǎo)數(shù)為f′′(a)=2.7182818353，步長(zhǎng)共減半11次.

請(qǐng)輸入要計(jì)算二階導(dǎo)數(shù)值的點(diǎn)a:1

請(qǐng)輸入精度要求epsilon:0.000000001

步長(zhǎng)減半次數(shù)已達(dá)25次，仍未達(dá)到精度要求.

算例2取初始步長(zhǎng)h0=1，用變步長(zhǎng)中點(diǎn)公式計(jì)算函數(shù)f(x)=x2e-x在點(diǎn)x=1.0,1.5,…,5.0處的一階導(dǎo)數(shù)及誤差,并統(tǒng)計(jì)對(duì)應(yīng)的最佳步長(zhǎng).

解 根據(jù)算法2編寫C語言程序生成文件Midpoints.com，計(jì)算結(jié)果如下：

5 總結(jié)

本研究提出了中心差商公式變步長(zhǎng)算法計(jì)算機(jī)實(shí)現(xiàn)時(shí)的兩種計(jì)算終止條件，若以計(jì)算結(jié)果滿足精度要求作為計(jì)算終止條件，可避免直接計(jì)算最優(yōu)步長(zhǎng)，方法簡(jiǎn)單，設(shè)計(jì)程序方便，可以獲得滿足一般精度要求的計(jì)算結(jié)果.但當(dāng)要求的精度很高時(shí)，此方案的計(jì)算精度可能不能令人十分滿意.若以步長(zhǎng)最接近最優(yōu)步長(zhǎng)作為計(jì)算終止條件，既可避免直接確定最優(yōu)步長(zhǎng)，又可獲得具有較高計(jì)算精度的計(jì)算結(jié)果，因而是一種非常實(shí)用的有效算法.