Julia語言 例題3-1 求sin(x)的微分

Julia語言 例題3-1 求sin(x)的微分

向前差  f'(x) =  ( f(x+h) - f(x)) / h
向後差  f'(x) =  ( f(x) - f(x-h)) / h

中央差  f'(x) =  0.5 * ( f(x+h) - f(x-h)) / h

取 h = 0.01

源自於
http://edisonx.pixnet.net/blog/post/90838030-%5Bc%E8%AA%9E%E8%A8%80%E6%95%B8%E5%80%BC%E5%88%86%E6%9E%90%5D-%E4%B8%80%E9%9A%8E%E5%BE%AE%E5%88%86%E6%B3%95


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，

using Printf

#=====================================================
Trigonometric and hyperbolic functions
All the standard trigonometric and hyperbolic functions are also defined:
sin    cos    tan    cot    sec    csc
sinh   cosh   tanh   coth   sech   csch
asin   acos   atan   acot   asec   acsc
asinh  acosh  atanh  acoth  asech  acsch
sinc   cosc   atan2
=====================================================#

function fx(x::Float64) #// 欲微分函數
    return sin(x)
end

function FordDiff(x::Float64 , h::Float64, fx::typeof(func)) #//前差微分
    a=fx(x+h)
    b=fx(x)
    return ((a-b) / h)
end

function BackDiff(x::Float64 , h::Float64, fx::typeof(func)) #//後差微分
    return ((fx(x) - fx(x-h) ) / h )
end

function MidDiff(x::Float64 , h::Float64, fx::typeof(func)) #//中差微分
    return (0.5 * ( fx(x+h) - fx(x-h) ) / h)
end

#=====================================================#
x=1.0
h=0.01
answer = cos(x)  # // 答案
cal = FordDiff(x, h, func)
delta = cal - answer
s=@sprintf("FordDiff : %0.5lf, delta = %0.5lf\n",cal,delta)
println(s)
             
cal = BackDiff(x, h, func)
delta = cal - answer
s=@sprintf("BackDiff : %0.5lf, delta = %0.5lf\n",cal,delta)
println(s)
             
cal = MidDiff(x, h, func)
delta = cal - answer
s=@sprintf("MidDiff : %0.5lf, delta = %0.5lf\n",cal,delta)
println(s)


輸出畫面

FordDiff : 0.53609, delta = -0.00422

BackDiff : 0.54450, delta = 0.00420

MidDiff : 0.54029, delta = -0.00001