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

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

一函數在某區間的積分值可以利用 Romberg integration 求得，而一函數在某點的微分值可以利用有限差分(finite difference)求出。有限差分可以分為前向差分(forward difference)、後向差分(backward difference)、中央差分(central difference)，以及結合外插觀念與差分原理所衍生出較準確之理查生外插法(Richardson extrapolation)，在此主要介紹理查生外插法，此法主要利用較低準確度的差分法，先求得不同 mesh size 之下的導數值，再利用外插法消去影響誤差最大之數值，而求得較準確之導數值。

以下的副程式先使用準確度為 O(h^²) 之中央差分法，依序求得 mesh size 為 h、h/2、h/4、h/8.....etc之導數值，再利用理查生外插法求得較準確的導數值。使用以下副程式時，若初始 h 設的太小或是欲求較高階之導數，可能會發生溢位或除以零的問題，必須特別注意。

副程式：

Option Explicit

'x=欲求某點之導數
'h=初始 mesh size (不可設的太小)
'n=導數之階數
'tol=容許誤差
Public Function Richardson(x As Double, h As Double, n As Long, tol As Double) As Double
     Dim i As Long, j As Long, T() As Double, L As Long
     ReDim T(1 To 3) As Double
     T(1) = CDn(x, h, n)
     T(2) = CDn(x, h / 2, n)
     T(3) = CDn(x, h / 4, n)
     j = 3
     Do While Abs(T(UBound(T)) - T(UBound(T) - 2)) > tol
          For i = 2 To UBound(T) Step 2
               T(i - 1) = T(i)
          Next
          T(2) = CDn(x, h / (2 ^ j), n)
          ReDim Preserve T(1 To UBound(T) + 2)
          For L = 2 To j
               If L <> j Then
                    T(L * 2) = ((4 ^ (L - 1)) * T(L * 2 - 2) - T(L * 2 - 3)) / ((4 ^ (L - 1)) - 1)
               Else
                    T(UBound(T)) = ((4 ^ (L - 1)) * T(UBound(T) - 1) - _
                                            T(UBound(T) - 2)) / ((4 ^ (L - 1)) - 1)
               End If
          Next
          j = j + 1
     Loop
     Richardson = T(UBound(T))
End Function

其中 CDn 函數請參閱本站的 Solve the nth derivative of a given function by using central difference method.

範例：

例如要求 f=e^{^x} 在 x=2 的微分值，使用初始 mesh size=1.0 可以寫出以下程式碼：

'對於不同的函數 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, tol As Double
     tol = 0.00001
     x = InputBox("請輸入 x 值", , "2")
     h = InputBox("請輸入間距 h 值", , "1")
     Debug.Print "一階導數 = " & Richardson(x, h, 1, tol)
     Debug.Print "二階導數 = " & Richardson(x, h, 2, tol)
     Debug.Print "三階導數 = " & Richardson(x, h, 3, tol)
     Debug.Print "四階導數 = " & Richardson(x, h, 4, tol)
End Sub

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

一階導數 = 7.38905609850807
二階導數 = 7.38905609870055
三階導數 = 7.38905609655925
四階導數 = 7.38905750432817

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