Visual Basic 數值分析 - Solve the nth derivative of a given function by using backward difference method.

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

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

f 的 n 階微分：

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

f' = ( 3f_{_-4} - 16f_{_-3} + 36f_{_-2} - 48f_{_-1} + 25f_₀ ) / 12h + O(h^⁴)

f'' = ( 11f_{_-4} - 56f_{_-3} + 114f_{_-2} - 104f_{_-1} + 35f_₀ ) / 12h^² + O(h^³)

f''' = ( 3f_{_-4} - 14f_{_-3} + 24f_{_-2} - 18f_{_-1} + 5f_₀ ) / 2h^³ + O(h^²)

f'''' = ( -2f_{_-5} + 11f_{_-4} - 24f_{_-3} + 26f_{_-2} - 14f_{_-1} + 3f_₀ ) / h^⁴ + O(h^²)

其中 f_₀=f(x_₀)、 f_{_-1}=f(x_₀-h)、f_{_-2}=f(x_₀-2h)，依此類堆。

副程式：

Option Explicit

Public Function BD1(x As Double, h As Double) As Double
     BD1 = (3 * f(x - 4 * h) - 16 * f(x - 3 * h) + 36 * f(x - 2 * h) - _
               48 * f(x - h) + 25 * f(x)) / (12 * h)
End Function

Public Function BD2(x As Double, h As Double) As Double
     BD2 = (11 * f(x - 4 * h) - 56 * f(x - 3 * h) + 114 * f(x - 2 * h) - _
              104 * f(x - h) + 35 * f(x)) / (12 * (h ^ 2))
End Function

Public Function BD3(x As Double, h As Double) As Double
     BD3 = (3 * f(x - 4 * h) - 14 * f(x - 3 * h) + 24 * f(x - 2 * h) - _
               18 * f(x - h) + 5 * f(x)) / (2 * (h ^ 3))
End Function

Public Function BD4(x As Double, h As Double) As Double
     BD4 = (-2 * f(x - 5 * h) + 11 * f(x - 4 * h) - 24 * f(x - 3 * h) + _
               26 * f(x - 2 * h) - 14 * f(x - h) + 3 * f(x)) / (h ^ 4)
End Function

Public Function BDn(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
     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
          BDn = BDn + (-1) ^ i * nFa * f(x - i * h) / (iFa * n_iFa)
          iFa = 1: n_iFa = 1
     Next
     BDn = BDn / (h ^ n)
End Function

範例：

例如要求 f=e^{^x} 在 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 "一階導數 = " & BD1(x, h)
     Debug.Print "二階導數 = " & BD2(x, h)
     Debug.Print "三階導數 = " & BD3(x, h)
     Debug.Print "四階導數 = " & BD4(x, h)
     Debug.Print "十階導數 = " & BDn(x, h, 10)
End Sub

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

一階導數 = 7.38893080689799
二階導數 = 7.38379283523175
三階導數 = 7.27673427084063
四階導數 = 7.21305836139052
十階導數 = 4.50163906151601

此題正確結果應為 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.