現實當中的不少實體問題、工程問題都涉及到微分方程,其中微分方程有常微分方程(Ordinary Differential Equation)和偏微分方程(Partial Differential Equation)之分。一般來說,所謂的常微分方程是指隻有一個自變量的方程,如 u' = 2*u+x+2 。
不少工程問題中涉及的微分方程,我們很難求出方程的解析解,或者說根本不存在精确的解析解。此時,我們需要利用電腦,結合數值分析的方法來近似求出微分方程的相關解,并研究其性質。通過求出多個自變量的值,并求出對應的解,那麼可以繪制出圖形來輔助研究方程的特征。
1 Runge–Kutta算法
根據百度百科的相關介紹,龍格-庫塔(Runge-Kutta)方法是一種在工程上應用廣泛的高精度單步算法。由于此算法精度高,采取措施對誤差進行抑制,是以其實作原理也較複雜。該算法是建構在數學分析的基礎之上的。一般的Runge-Kutta算法的形式如下:

這個公式看起來也比較的抽象,在實際計算中,一般采用4階Runge-Kutta算法進行求值,其計算公式如下:
為了更加直覺,這裡選擇一個示例,假設有如下的一個微分方程需要進行求解:
此時,我們需要計算一下u(0.2)對應的方程解是什麼?那麼此示例的解題過程如下:
實際上,本微分方程的解析解表達式為:
即解析式的解近似為 1.3472599 ,而4階Runge-Kutta算法進行求值的結果為1.3472。由此可知,精确度還是可以的。上述示例來自網址 :
http://www.public.asu.edu/~hhuang38/example_Runge-Kutta.pdf2 F# Runge–Kutta算法實作
下面給出F#的算法實作,示例代碼如下:
let odeInt f (a:float) (b:float) x0 f0 =
let n = 5
let h = (b-a)/float(n)
let x = [| for i in 0 .. n -> a + float(i) * h |]
let u = [| for i in 0 .. (n+1) -> 0. |]
u.[0] <- f0
for i = 0 to n do
let k1 = f(x.[i], u.[i])
let k2 = f(x.[i]+h / 2. , u.[i]+k1*h / 2.)
let k3 = f(x.[i]+h / 2. , u.[i]+k2*h / 2.)
let k4 = f(x.[i]+h , u.[i]+k3*h)
u.[i+1] <- u.[i] + (k1 + 2.*(k2 + k3) + k4) * h / 6.
(x,u)
3 測試
下面給出測試示例,代碼如下:
let testODEInt =
let f1(x,u) = -2.*u + x + 4.
let dx,du = odeInt f1 0. 1.0 0. 1.
printfn "%A" dx
printfn "%A" du
執行示例代碼,結果如下:
[|0.0; 0.2; 0.4; 0.6; 0.8; 1.0|]
[|1.0; 1.3472; 1.61292288; 1.824023499; 1.998505354; 2.148437989; 2.281912828|]
由此可知,在0.2時,u(x)值為1.3472,與示例一緻。