Edison.X. Blog

1.  若微分用定義方式做，效果有限，若欲求之 eps 過小 (如精度達 1E-9, 1E-12 )，一些情況上特別容易發生 除零之窘境。

2.  數值分析之叢書不知道對微分法說明完不完整，即使跳過不講大概也不會讓人太過意外。但較意外的是，非線性方程式求解中的「牛頓法」，微分函式不是用數值分析方式算出來的，而是事先用筆算算出來的。

從定義做

一般微積分課本裡對於微分之定義大概如此  

f ' (x) = [ f(x+h) -  f(x) ] / h  , h -> 0

當然可微包含的三個條件就不再重述了。

我們假設  f(x) = x * x ，要求 f'(2) 之值為何，

若取 h = 0.01，按上述的定義計算  f'(x) = [ f(2.01) - f(2) ] / 0.01 = 4.01；

再取 h = 0.001，按上述的定義計算 f'(x) = [f(2.001) - f(2) ] / 0.001 = 4.001;

實際上和 f'(2) 之誤差，與 h 成正比。

前差近似、後差近似、央央近似

然而上述所使用的模型，是屬較差的微分模型，所使用的模型為 「二點前差」微分法。

在圖形上取得兩點，分別為 x，得到 f(x)，與 x 往前看一個 h 之 f(x+h)，

而微分所使用的方法絕不只這項，還有「後差除」、「中央差除」等方式，

所使用公式約略如下。

注意到誤差項，那很重要，若對於輸出之微分誤差精度有所要求，

它可拿來做為概估之 h 推算。根據上述之三種不同方法，這次我們以 f(x) = sin(x) ，取 h = 0.01，

以 C 語言計算在 x=1 之微分值為何，並與 cos(1) 做結果比較。

輸出結果

FordDiff : 0.536086, delta = -0.004216
BackDiff : 0.544501, delta = 0.004198
MidDiff : 0.540293, delta = -0.000009

雖然都呼叫了兩次 func，但誤差以中央差除法較小，約與 h 之平方成正比。

除了上述三種方法外，另外對於一階微分上也有「多取幾個點」方式做計算。

推導過程其實與 Taylor series 有密切關係，不會很難，有興趣可借數值分析書回來看，

鑑於筆者相簿容量與撰文時間有限，只挑幾項導導。

Taylor Series 方式推導

(A) 前差法

(1) 先對  f(x) 以 Taylor series 對點 xi 做展開

f(x) = f(xi) + (x-xi) * f'(xi) / 1!  + (x-xi)^2 * f''(xi) / 2! + ...

(2)  x 以 x_i+1 代入，令 x-xi = h

f(x_i+1) = f(xi) + h * f'(xi) + h^2 * f''(xi) / 2! + ....

又 f(x_i+1) 可寫成 f(xi + h)，變成

f(xi+h) = f(xi) + h * f'(xi) + h^2 * f''(xi) / 2! + ....

(3) 只前兩項做近似

f(xi+h) = (fxi) + h * f'(xi) 

f'(xi) = [f(xi+h) - f(xi)] / h

誤差以 Taylor series 第三項代表，

O(h) = h² * f''(xi) / 2! = h² * f''(xi) / 2

(B) 後差法 

和前差法一樣的方式，只是將 x = xi-1、h = xi-1 - x 做替換，不再示範。

(C) 中央差法

(1) 對 x = xi 做展開式

f(x) = f(xi) + (x-xi) * f'(xi) / 1!  + (x-xi)^2 * f''(xi) / 2! + ...

(2) x 以 xi+h 代入，令 x-xi = h，代入 (1)

f(xi+h) = f(xi) + (x-xi) * f'(xi) / 1! + (x-xi)² * f''(xi) / 2! + (x-xi)³ * f'''(xi) / 3! + ...

f(xi+h) = f(xi) + h * f'(xi) / 1! + h² * f''(xi) / 2! + h³ * f'''(xi) / 3! + ...

令為 a 式。

(3) x 以 xi-h 代入，令 xi - x = h，代入 (1) 

f(xi-h) = f(xi) - h * f'(xi) / 1! + h² * f''(xi) / 2! -  h³ * f'''(xi) / 3!  + ....

令為 b 式。

(4) 連立 a, b 式 < 消掉 h² >

f(xi+h) = f(xi) + h * f'(xi) / 1! + h² * f''(xi) / 2! + h³ * f'''(xi) / 3! + ... (a)

f(xi-h) = f(xi) - h * f'(xi) / 1! + h² * f''(xi) / 2! -  h³ * f'''(xi) / 3!  + .... (b)

(a) - (b) 再除 2 得到

[ f(xi+h) - f(xi-h) ] / 2 = h * f'(xi) + h³ * f'''(xi) / 3! + ...

取前兩項得到  f'(xi) = [ f(xi+h) - f(xi-h) ] / 2h，

故最高誤差項約為  O(h) = h³ * f'''(xi) / 3! =  h³ * f'''(xi) / 6

(D) 較高精度之前差法

f(xi + h) = f(xi) + h f'(xi) + h² f''(xi)/2 + h³ f'''(xi) / 6 + ...  (X 式)

h 用  2h 替代

f(xi + 2h) = f(xi) + 2h f'(xi) + 4h² f''(xi)/2 + 8h³ f'''(xi) / 6 + ... (Y 式)

消去 h² 項， (Y) - 4(X) 得到

f(xi+2h) - 4f(xi+h) = -3 f(xi) -2 h f'(xi) + 4h³ f'''(xi) / 6 +...

忽略 h3 後之影響得到

f'(xi) = -[ f(xi+2h) - 4f(xi+h) + 3 f(xi) ] / 2h

f'(xi) = [ -f(xi+2h) +4 f(xi+h) - 3f(xi) ] / 2h

誤差項懶得再推，有興趣自己算，

O(h) = 1/3 h²  f'''(xi)

注意，這種方式可以一直加精度下去 < 但會有極限 >

一開始列出 taylor series , f(x+h) = .... 後，

這裡是用 2h 取代 h 得到第二個點的展開，

若再拿 3h 取代 h 得到第三個點的展開 (三個等式)，再進行消去，得到的精度會更高，

可依此類推。

(E) 較高精度之中央差分 

f(x) = f(xi) + (x-xi) * f'(xi) / 1!  + (x-xi)^2 * f''(xi) / 2! + ...

f(xi-2h) = f(xi) - 2h * f'(xi) / 1 + 4h² * f''(xi) /4 -  8h³ * f'''(xi) / 6  + .... (W)

f(xi-h) = f(xi) - h * f'(xi) / 1 + h² * f''(xi) / 4 -  h³ * f'''(xi) / 6  + ....  (X)

f(xi+h) = f(xi) + h * f'(xi) / 1 + h² * f''(xi) / 4 + h³ * f'''(xi) / 6 + ...  (Y)

f(xi+2h) = f(xi) + 2h * f'(xi) / 1! + 4h² * f''(xi) /4  + 8h³ * f'''(xi) / 6 + ... (Z)

找 a(W) + b(X) + c(Y) + d(Z) 可消掉 h² 系數，

所以挑出來的系數不只一組，但把握一點原則：

消掉後所剩之誤差愈小愈好，妙的是有組系數 [a,b,c,d]，

除了消掉 h² 外，還外帶消掉 h³ 系數 (必然可找到這組係數無誤) ，

使得最高誤差項壓在 h⁴ 

這是一個線性代數問題了。不再細解。

進行 -(W) + 8 (X) - 8(Y) + (Z)，忽略 h³以後，得到

- f(xi - 2h) + 8 f(xi-h) - 8 f(xi+h) + f(xi+2h) 

= (2 -8 -8 +2) * h * f'(xi) + (-4 + 8 - 8 +4) / 4 * h² * f''(xi)

= -12 h * f'(xi) 

其他的稍加整理，得 f'(xi)

f'(xi) = f(xi - 2h) - 8f(xi-h) + 8(xi+h) - f(xi+2h) / 12h

最後最高誤差項算出來為 O(h) = 1/90  * h⁴ * [ f(xi) 六次微分]，其他的推導便不再示範。

整理

在「三點近似」裡注意，雖用了三點前差近似，

但誤差卻比雙點中差差值近似來得大， 這也是筆者強調誤差項原因之一。

若對其它的型式有興趣，可到 wiki 網頁，拉到 Higher Order Method 看述敘。

再以 f(x) = sin(x) 為例，C 語言以不同方式撰寫，觀查其結果。

輸出結果

FordDiff : 0.536086, delta = -4.2163248562707700e-003
BackDiff : 0.544501, delta = +4.1983148694582084e-003

MidDiff : 0.540293, delta = -9.0049934062808035e-006
Ford3Diff : 0.540320, delta = +1.7799082280500755e-005

Ford5Diff : 0.540302, delta = -1.0524227045394241e-009
Mid4Diff : 0.540302, delta = -1.8009915780936581e-010

注意到，(雙點中央差除) 結果比 (三點前差除) 結果來得好；而 (四點中央差除)  比 (五點前差除) 結果來得好。

常見近似整理 

假設 f(i) = f(x)，f(i+1) = f(x+h)，f(i+n) = f(x + n*h) ，誤差 O(h)，

一些常見的公式如下。

(A) 向前差近似

雙點 

fi' =  [ f(i+1) - f(i) ] / h ,  

O(h) = -1/2 * h * f''i < fi 二次微分 >

三點 

fi' =  [ -f(i+2) + 4f(i+1) -3 f(i) ] / 2h ,  

O(h) = 1/3 * h*h * f'''i < fi 三次微分 >

四點 

fi' =  [ 2f(i+3) - 9f(i+2) + 18 f(i+1) -11 f(i) ] / 6h , 

O(h) = -1/4 * h*h * h * f''''i < fi 四次微分 >

(B) 向後差近似

雙點

fi' =  [ f(i) - f(i-h) ] / h ,  

O(h) = -1/2 * h * f''i  < fi 二次微分 >

三點 

fi' =  [ 3f(i) - 4f(i-1) + f(i-2) ] / 2h ,  

O(h) = 1/3 * h*h * f'''i < fi 三次微分 >

四點 

fi' =  [ 11f(i) - 18f(i-1) + 9 f(i-2) - 2f(i-3) ] / 6h ,  

O(h) = 1/4 * h*h * h * f''''i  < fi 四次微分 >

(C) 中央近似

雙點

fi' = [f(i+1) - f(i-1)] / 2h  , 

O(h) = -1/6 * h *h * f'''i < fi 三次微分 >

三點 

 中央近似沒有三點

四點

f' = [ -f(i+2) + 8 f(i+1) - 8 f(i-1) + f(i-2) ] / 12h , 

O(h) = 1 / 30 * h * h * h * h * f'''''i < fi 五分微分 >

誤差估算

我們以 f(x) = sin(x) ，求 f'(1)，取 h = 0.01 ，拿雙點中央近似來講，

fi' = [f(i+1) - f(i-1)] / 2h  , 算出來之 f'(1) 約為 5.402933008747335e-001,

估算差項誤差 O(h) = -1/6 * h *h * f'''i < fi 三次微分 >，

O(h) = 1 /6 * 1e-2 * 1e-2 * cos(1e-2) = 1.6665833340277756e-005, 

故可以推斷，剛剛算出來的 5.402933008747335e-001，與實際的微分值，

誤差在 O(h) = 1.6665833340277756e-005 以內，

而實際上與算出來的誤差是在 -9.0049934062808035e-006 ，

確實在 O(h) 範圍內。

換句話說，若以 IEEE754 而言，忽略 f''i 之影響，要得到 2.2E-16 之精確值，

h 約要設 根號 [ (2.2E-16) * 6 ] ，約為 9E-8 左右；

然而若換作四點的中央誤差，由於 

O(h) = 1 / 30 * h * h * h * h * f'''''i < fi 五分微分 >

反推回去，h 約只要設成 3E-4 即可，而 h=3E-4 代入 O(h) 後，

O(h) = 2.7 e -16

實際上代入後，雖沒辦法得到 2.7E-16，但也得到了 8.5E-14 之誤差 < 因浮點數本身便有誤差>。舉此例主要說明，裡面的 h 該設多少，與挑用的公式、階數、誤差項有絕大的關係。