ODE in F#

module ODE
let isBetween (a1, a2) x =
    if a1 < a2 then a1 <= x && x<=a2
    else            a2 <= x && x<=a1

let SolveEuler (f:float*float->float) (h:float) (init:float*float) (T:float) =
    let t0, _ = init
    Seq.unfold (fun ty ->
                    let t_, y_ = fst ty, snd ty
                    if not (isBetween(t0, T) t_) then None
                    else
                        Some (ty, (t_ + h, y_ + h*(f ty)))
                ) init

let SolveRungeKutta (f:float*float->float) (h:float) (init:float*float) (T:float) =
    let t0, _ = init
    Seq.unfold (fun ty ->
                    let t_, y_ = fst ty, snd ty                   
                    if not (isBetween(t0, T) t_) then None
                    else
                        let k1 = f ty
                        let k2 = f (0.5*h+t_, y_+0.5*h*k1)
                        let k3 = f (0.5*h+t_, y_+0.5*h*k2)
                        let k4 = f (    h+t_, y_+h*k3)
                        let t,y = t_+h, y_ + h*(k1+2.0*k2+2.0*k3+k4)/6.0
                        Some (ty, (t, y))
                ) init
Ebook Download
View all
Learn
View all