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),

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

所使用公式約略如下。

 

diff03   

 

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

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

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

 

Code Snippet
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. double func(double x) // 欲微分函數
  5. { return sin(x) ;}
  6.  
  7. double FordDiff(double x, double h, double (*fx)(double))
  8. { // 前差微分
  9.     return ( fx(x+h) - fx(x) ) / h;
  10. }
  11.  
  12. double BackDiff(double x, double h, double (*fx)(double))
  13. { // 後差微分
  14.     return ( fx(x) - fx(x-h) ) / h;
  15. }
  16.  
  17. double MidDiff(double x, double h,double (*fx)(double))
  18. { // 中差微分
  19.     return 0.5 * ( fx(x+h) - fx(x-h) ) / h;
  20. }
  21.  
  22. int main()
  23. {
  24.     double x = 1.0, h=0.01, delta;
  25.     double ans = cos(x); // 答案
  26.     double cal;
  27.  
  28.     cal = FordDiff(x, h, func), delta = cal - ans;
  29.     printf("FordDiff : %lf, delta = %lf\n",cal,delta);
  30.     cal = BackDiff(x, h, func), delta = cal - ans;
  31.     printf("BackDiff : %lf, delta = %lf\n",cal,delta);
  32.     cal = MidDiff(x, h, func), delta = cal - ans;
  33.     printf("MidDiff  : %lf, delta = %lf\n",cal,delta);
  34.     return 0;
  35. }

 

輸出結果

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 以 xi+1 代入,令 x-xi = h

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

又 f(xi+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) = h2 * f''(xi) / 2! = h2 * 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)2 * f''(xi) / 2! + (x-xi)3 * f'''(xi) / 3! + ...

f(xi+h) = f(xi) + h * f'(xi) / 1! + h2 * f''(xi) / 2! + h3 * f'''(xi) / 3! + ...

令為 a 式。


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

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

令為 b 式。


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

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

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

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

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

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

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

 

(D) 較高精度之前差法

f(xi + h) = f(xi) + h f'(xi) + h2 f''(xi)/2 + h3 f'''(xi) / 6 + ...  (X 式)

h 用  2h 替代

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

消去 h2 項, (Y) - 4(X) 得到

f(xi+2h) - 4f(xi+h) = -3 f(xi) -2 h f'(xi) + 4h3 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 h2  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 + 4h2 * f''(xi) /4 -  8h3 * f'''(xi) / 6  + .... (W)

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

f(xi+h) = f(xi) + h * f'(xi) / 1 + h2 * f''(xi) / 4 + h3 * f'''(xi) / 6 + ...  (Y)

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

 

找 a(W) + b(X) + c(Y) + d(Z) 可消掉 h2 系數,

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

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

除了消掉 h2 外,還外帶消掉 h3 系數 (必然可找到這組係數無誤) ,

使得最高誤差項壓在 h

 

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

 

進行 -(W) + 8 (X) - 8(Y) + (Z),忽略 h3以後,得到

- 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 * h2 * 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  * h4 * [ f(xi) 六次微分],其他的推導便不再示範。

 

 

 

整理

 

 

diff02

 

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

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

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

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

 

Code Snippet
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. double func(double x) // 欲微分函數
  5. { return sin(x) ;}
  6.  
  7. double FordDiff(double x, double h, double (*fx)(double))
  8. { // 前差微分
  9.     return ( fx(x+h) - fx(x) ) / h;
  10. }
  11.  
  12. double BackDiff(double x, double h, double (*fx)(double))
  13. { // 後差微分
  14.     return ( fx(x) - fx(x-h) ) / h;
  15. }
  16.  
  17. double MidDiff(double x, double h,double (*fx)(double))
  18. { // 中差微分
  19.     return 0.5 * ( fx(x+h) - fx(x-h) ) / h;
  20. }
  21. double Ford3Diff(double x, double h, double (*fx)(double))
  22. { // 三點前差微分
  23.     return 0.5 * ( -fx(x+2*h) + 4*fx(x+h) - 3*fx(x)) / h;
  24. }
  25. double Ford5Diff(double x, double h, double (*fx)(double))
  26. { // 五點前差微分
  27.     return
  28.         ( -25*fx(x)+48*fx(x+h)-36*fx(x+2*h)
  29.                    +16*fx(x+3*h)-3*fx(x+4*h) ) / (12*h) ;
  30. }
  31. double Mid4Diff(double x, double h, double (*fx)(double))
  32. { // 四點中差微分
  33.     return
  34.         (-fx(x+2*h)+fx(x-2*h)+ 8*(fx(x+h)-fx(x-h))) / (12*h);
  35. }
  36.  
  37. int main()
  38. {
  39.     double x = 1.0, h=0.01, delta;
  40.     double ans = cos(x); // 答案
  41.     double cal;
  42.  
  43.     cal = FordDiff(x, h, func), delta = cal - ans;
  44.     printf("FordDiff   : %lf, delta = %+.16e\n",cal,delta);
  45.     cal = BackDiff(x, h, func), delta = cal - ans;
  46.     printf("BackDiff   : %lf, delta = %+.16e\n\n",cal,delta);
  47.  
  48.     cal = MidDiff(x, h, func), delta = cal - ans;
  49.     printf("MidDiff    : %lf, delta = %+.16e\n",cal,delta);
  50.     cal = Ford3Diff(x, h, func), delta = cal - ans;
  51.     printf("Ford3Diff  : %lf, delta = %+.16e\n\n",cal,delta);
  52.  
  53.     cal = Ford5Diff(x, h, func), delta = cal - ans;
  54.     printf("Ford5Diff  : %lf, delta = %+.16e\n",cal,delta);
  55.     cal = Mid4Diff(x, h, func), delta = cal - ans;
  56.     printf("Mid4Diff   : %lf, delta = %+.16e\n\n",cal,delta);
  57.     return 0;
  58. }

 

輸出結果

 

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 該設多少,與挑用的公式、階數、誤差項有絕大的關係。

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()