F#实现Runge–Kutta算法求解常微分方程

    现实当中的不少物理问题、工程问题都涉及到微分方程，其中微分方程有常微分方程（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.pdf

2 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，与示例一致。