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

Solve the nth derivative of a given function by using forward 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_₄ + 16f_₃ - 36f_₂ + 48f_₁ - 25f_₀ ) / 12h + O(h^⁴)

f'' = ( 11f_₄ - 56f_₃ + 114f_₂ - 104f_₁ + 35f_₀ ) / 12h^² + O(h^³)

f''' = ( -3f_₄ + 14f_₃ - 24f_₂ + 18f_₁ - 5f_₀ ) / 2h^³ + O(h^²)

f'''' = ( -2f_₅ + 11f_₄ - 24f_₃ + 26f_₂ - 14f_₁ + 3f_₀ ) / h^⁴ + O(h^²)

其中 f_₀=f(x_₀)、 f_₁=f(x_₀+h)、f_₂=f(x_₀+2h)，依此類堆。

副程式：

Option Explicit

Public Function FD1(x As Double, h As Double) As Double
     FD1 = (-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 FD2(x As Double, h As Double) As Double
     FD2 = (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 FD3(x As Double, h As Double) As Double
     FD3 = (-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 FD4(x As Double, h As Double) As Double
     FD4 = (-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 FDn(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
          FDn = FDn + (-1) ^ i * nFa * f(x + (n - i) * h) / (iFa * n_iFa)
          iFa = 1: n_iFa = 1
     Next
     FDn = FDn / (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 "一階導數 = " & FD1(x, h)
     Debug.Print "二階導數 = " & FD2(x, h)
     Debug.Print "三階導數 = " & FD3(x, h)
     Debug.Print "四階導數 = " & FD4(x, h)
     Debug.Print "十階導數 = " & FDn(x, h, 10)
End Sub

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

一階導數 = 7.38888123788608
二階導數 = 7.39628507050393
三階導數 = 7.23958531340329
四階導數 = 7.1385622940312
十階導數 = 12.2372334487863

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