例題3-1 利用 向前差 向後差 中央差 求f(x)=sin(x)的微分

例題3-1 利用　向前差　向後差　中央差　求f(x)=sin(x)的微分

'''
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 (3次微分 (xi) / 3! + ...
f(xi+h) = f(xi) + h * f'(xi) / 1! + h2 * f''(xi) / 2! + h3 * f (3次微分(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(3次微分(xi) / 3!  + .... 令為 b 式。
(4) 連立 a, b 式 < 消掉 h2 >
f(xi+h) = f(xi) + h * f'(xi) / 1! + h2 * f''(xi) / 2! + h3 * f(3次微分(xi) / 3! + ... (a)
f(xi-h) = f(xi) - h * f'(xi) / 1! + h2 * f''(xi) / 2! -  h3 * f(3次微分(xi) / 3!  + .... (b)
(a) - (b) 再除 2 得到 [ f(xi+h) - f(xi-h) ] / 2 = h * f'(xi) + h3 * f(3次微分(xi) / 3! + ...
取前兩項得到  f'(xi) = [ f(xi+h) - f(xi-h) ] / 2h，
故最高誤差項約為  O(h) = h3 * f(3次微分'(xi) / 3! =  h3 * f'(3次微分(xi) / 6
'''
import math

def func( x):    # // 欲微分函數
     return math.sin(x)

def FordDiff( x, h, fx):
    # // 前差微分
    return ( ( fx(x+h) - fx(x) ) / h)

def BackDiff(x,  h,  fx):
    # // 後差微分
    return (( fx(x) - fx(x-h) ) / h)

def MidDiff(x,  h,  fx):
    # // 中差微分
    return (0.5 * ( fx(x+h) - fx(x-h) ) / h)

x = 1.0
h1=[0.01 , 0.005 , 0.01 , 0.05 , 0.1 , 0.5]

ans =math.cos(x)     #// 答案

n=6
for i in range (0 , n):
    h=h1[i]
    cal = FordDiff(x, h, func) + 0.0
    delta = cal - ans + 0.0
    print('h= {%10.6f} '  %(h))
    print("FordDiff : {%10.6f}, delta: {%10.6f}  ,delta_percent : {%10.6f}" %( cal , delta ,  abs(delta/ans)*100 )  )
    cal = BackDiff(x, h, func)
    delta = cal - ans;
    print("BackDiff : {%10.6f}, delta: {%10.6f}  ,delta_percent  : {%10.6f}" %(cal ,  delta ,   abs(delta/ans)*100 )  )
    cal = MidDiff(x, h, func)
    delta = cal - ans;
    print("MidDiff  : {%10.6f}, delta: {%10.6f}  ,delta_percent  : {%10.6f} \n" %(cal , delta , abs(delta/ans)*100 ) )

輸出結果

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX3-1.py =============

h= {  0.010000} 

FordDiff : {  0.536086}, delta: { -0.004216}  ,delta_percent : {  0.780364}

BackDiff : {  0.544501}, delta: {  0.004198}  ,delta_percent  : {  0.777031}

MidDiff  : {  0.540293}, delta: { -0.000009}  ,delta_percent  : {  0.001667} 

h= {  0.005000} 

FordDiff : {  0.538196}, delta: { -0.002106}  ,delta_percent : {  0.389768}

BackDiff : {  0.542404}, delta: {  0.002101}  ,delta_percent  : {  0.388934}

MidDiff  : {  0.540300}, delta: { -0.000002}  ,delta_percent  : {  0.000417} 

h= {  0.010000} 

FordDiff : {  0.536086}, delta: { -0.004216}  ,delta_percent : {  0.780364}

BackDiff : {  0.544501}, delta: {  0.004198}  ,delta_percent  : {  0.777031}

MidDiff  : {  0.540293}, delta: { -0.000009}  ,delta_percent  : {  0.001667} 

h= {  0.050000} 

FordDiff : {  0.519045}, delta: { -0.021257}  ,delta_percent : {  3.934370}

BackDiff : {  0.561110}, delta: {  0.020807}  ,delta_percent  : {  3.851047}

MidDiff  : {  0.540077}, delta: { -0.000225}  ,delta_percent  : {  0.041661} 

h= {  0.100000} 

FordDiff : {  0.497364}, delta: { -0.042939}  ,delta_percent : {  7.947135}

BackDiff : {  0.581441}, delta: {  0.041138}  ,delta_percent  : {  7.613968}

MidDiff  : {  0.539402}, delta: { -0.000900}  ,delta_percent  : {  0.166583} 

h= {  0.500000} 

FordDiff : {  0.312048}, delta: { -0.228254}  ,delta_percent : { 42.245665}

BackDiff : {  0.724091}, delta: {  0.183789}  ,delta_percent  : { 34.015880}

MidDiff  : {  0.518069}, delta: { -0.022233}  ,delta_percent  : {  4.114892} 

>>>