Solve the nth derivative of a given function by using central difference method. (線上執行) (原始碼下載)

一函數在某區間的積分值可以利用 Romberg integration 求得，而一函數在某點的微分值可以利用有限差分(finite difference)求出。有限差分可以分為前向差分(forward difference)、後向差分(backward difference)、中央差分(central difference)，以及結合外插觀念與差分原理所衍生出較準確之理查生外插法(Richardson extrapolation)，在此主要介紹中央差分法，利用 Taylor series expansion 可以導出函數 f 在 x=x0 的各階導數公式( mesh size=h )，詳細的推導過程請參閱市面上的數值分析書籍。

f 的 n 階微分 (n為偶數)：

f的 n 階微分 (n為奇數)：

上述公式的準確度只到O(h2)，對於較常用之 f 的 1 ~ 4 階微分可以用以下較準確之公式：

f' = ( f-2 - 8f-1 + 8f1 - f2 ) / 12h + O(h4)

f'' = ( -f-2 + 16f-1 - 30f0 + 16f1 - f2 ) / 12h2 + O(h4)

f''' = ( f-3 - 8f-2 + 13f-1 - 13f1 + 8f2 - f3 ) / 8h3 + O(h4)

f'''' = ( -f-3 + 12f-2 - 39f-1 + 56f0 - 39f1 + 12f2 - f3 ) / 6h4 + O(h4)

其中 f0=f(x0)、 f-1=f(x0-h)、f-2=f(x0-2h)、 f1=f(x0+h)、f2=f(x0+2h)，依此類堆。

　

副程式：

Option Explicit

Public Function CD1(x As Double, h As Double) As Double
     CD1 = (f(x - 2 * h) - 8 * f(x - h) + 8 * f(x + h) - f(x + 2 * h)) / (12 * h)
End Function

Public Function CD2(x As Double, h As Double) As Double
     CD2 = (-f(x - 2 * h) + 16 * f(x - h) - 30 * f(x) + 16 * f(x + h) _
                - f(x + 2 * h)) / (12 * (h ^ 2))
End Function

Public Function CD3(x As Double, h As Double) As Double
     CD3 = (f(x - 3 * h) - 8 * f(x - 2 * h) + 13 * f(x - h) - 13 * f(x + h) _
                + 8 * f(x + 2 * h) - f(x + 3 * h)) / (8 * (h ^ 3))
End Function

Public Function CD4(x As Double, h As Double) As Double
     CD4 = (-f(x - 3 * h) + 12 * f(x - 2 * h) - 39 * f(x - h) + 56 * f(x) - _
               39 * f(x + h) + 12 * f(x + 2 * h) - f(x + 3 * h)) / (6 * (h ^ 4))
End Function

Public Function CDn(x As Double, h As Double, n As Long) As Double
     Dim i As Long, j As Long, nFa As Long, iFa As Long, n_iFa As Long, even As Boolean
     even = (n Mod 2 = 0)
     nFa = 1: iFa = 1: n_iFa = 1
     For i = n To 2 Step -1
          nFa = nFa * i
     Next
     For i = 0 To n
          For j = i To 2 Step -1
               iFa = iFa * j
          Next
          For j = n - i To 2 Step -1
               n_iFa = n_iFa * j
          Next
          If even Then
               CDn = CDn + (-1) ^ i * nFa * 2 * f(x + (n / 2 - i) * h) / (iFa * n_iFa)
          Else
               CDn = CDn + (-1) ^ i * nFa * (f(x + ((n - 1) / 2 - i) * h) + _
                          f(x + ((n + 1) / 2 - i) * h)) / (iFa * n_iFa)
          End If
          iFa = 1: n_iFa = 1
     Next
     CDn = CDn / (2 * (h ^ n))
End Function

範例：

例如要求 f=ex 在 x=2 的微分值，使用 mesh size=0.1 可以寫出以下程式碼：

'對於不同的函數 f , 只要更改以下副程式即可
Public Function f(x As Double) As Double
     f = Exp(x)
End Function
Private Sub Command1_Click()
     Dim x As Double, h As Double
     x = InputBox("請輸入 x 值", , "2")
     h = InputBox("請輸入間距 h 值", , "0.1")
     Debug.Print "一階導數 = " & CD1(x, h)
     Debug.Print "二階導數 = " & CD2(x, h)
     Debug.Print "三階導數 = " & CD3(x, h)
     Debug.Print "四階導數 = " & CD4(x, h)
     Debug.Print "十階導數 = " & CDn(x, h, 10)
End Sub

執行後在即時運算視窗可得以下結果：

一階導數 = 7.38903143940491
二階導數 = 7.38904788153459
三階導數 = 7.38901291661409
四階導數 = 7.3890345157294
十階導數 = 7.41786188029891

此題正確結果應為 7.389056099，由此可知愈高階導數所得結果誤差愈大，主要誤差來自於 truncation error，少部分誤差來自於 roundoff error。一般來說減小 mesh size 可以增加計算結果的準確度，但是 mesh size 小至某一程度後，會使得 roundoff error 的效應變得很明顯，因此太小的 mesh size 不見得能得到最準確之結果。

以上副程式只適用於解決連續函數之微分問題，對於一不連續函數或是具有奇異點(singularity)之函數將不適用；對於一週期性函數(例如：sin x)，mesh size 不可等於該函數之週期；對於某些較特殊之函數，必須注意該函數之值域(range)，例如 f=log(x)，欲求 x=0.1 時之微分值，使用 mesh size=0.2，此時若用後向差分或中央差分將會產生問題，因為 x< 0 並不在 x 的定義域內(domain)，因此也不在 log(x) 的值域內，此時應該使用前向差分才不致發生這類問題。

　

[ 上一個 | 首頁 | 數值分析 | 下一個 ]

This page was written by Jaric on Jul. 31 , 1999. All rights reserved.

Total pageview since 4/6/1999.