3 people like it.

Fourth order Runge-Kutta ODE solver

This is a simple and direct implementation of fourth order runge-kutta ordinary differential equation solver algorithm. In the main function three use cases are shown.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
open System

// -- | ODE, y = e ^ x
let mydydx1 x y = 
    exp(x)
    
// -- | ODE, y = x * e ^ 3 * x - 2 * y
let mydydx2 x y = 
    x * exp (3.0 * x) - 2.0 * y
    
//-- | ODE, y = r * y
let mydydx3 x y = 
    let r = 0.5
    r * y    

// -- | Fourth order runge-kutta algorithm calculations
let rungekutta4' x y h f = 
    let k1 = h * f x y
    let k2 = h * f (x + 0.5 * h) (y + k1 * 0.5)
    let k3 = h * f (x + 0.5 * h) (y + k2 * 0.5)
    let k4 = h * f (x + h) (y + k3)
    y + 1.0/6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

// -- | Fourth order runge-kutta solver
let rec rungekutta4 x y h n f xx yy = 
    let y' = rungekutta4' x y h f
    if (x >= n) then (xx,yy)
    else rungekutta4 (x+h) y' h n f (List.append xx [x+h]) (List.append yy [y'])


// -- | The main entry point.
let main() =
    Console.WriteLine("Welcome to fourth order Runge-Kutta ODE solver!")
    let v = rungekutta4 0.0 1.0 0.01 1.0 mydydx1 [] []
    printfn "%A" ( (snd v).Item ((List.length (snd v)) - 1))
    let v = rungekutta4 0.0 0.0 0.01 1.0 mydydx2 [] []
    printfn "%A" ( (snd v).Item ((List.length (snd v)) - 1))
    let v = rungekutta4 0.0 2.0 0.01 3.0 mydydx3 [] []
    printfn "%A" ( (snd v).Item ((List.length (snd v)) - 1))

main()
namespace System
val mydydx1 : x:float -> y:'a -> float

Full name: Script.mydydx1
val x : float
val y : 'a
val exp : value:'T -> 'T (requires member Exp)

Full name: Microsoft.FSharp.Core.Operators.exp
val mydydx2 : x:float -> y:float -> float

Full name: Script.mydydx2
val y : float
val mydydx3 : x:'a -> y:float -> float

Full name: Script.mydydx3
val x : 'a
val r : float
val rungekutta4' : x:float -> y:float -> h:float -> f:(float -> float -> float) -> float

Full name: Script.rungekutta4'
val h : float
val f : (float -> float -> float)
val k1 : float
val k2 : float
val k3 : float
val k4 : float
val rungekutta4 : x:float -> y:float -> h:float -> n:float -> f:(float -> float -> float) -> xx:float list -> yy:float list -> float list * float list

Full name: Script.rungekutta4
val n : float
val xx : float list
val yy : float list
val y' : float
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val append : list1:'T list -> list2:'T list -> 'T list

Full name: Microsoft.FSharp.Collections.List.append
val main : unit -> unit

Full name: Script.main
type Console =
  static member BackgroundColor : ConsoleColor with get, set
  static member Beep : unit -> unit + 1 overload
  static member BufferHeight : int with get, set
  static member BufferWidth : int with get, set
  static member CapsLock : bool
  static member Clear : unit -> unit
  static member CursorLeft : int with get, set
  static member CursorSize : int with get, set
  static member CursorTop : int with get, set
  static member CursorVisible : bool with get, set
  ...

Full name: System.Console
Console.WriteLine() : unit
   (+0 other overloads)
Console.WriteLine(value: string) : unit
   (+0 other overloads)
Console.WriteLine(value: obj) : unit
   (+0 other overloads)
Console.WriteLine(value: uint64) : unit
   (+0 other overloads)
Console.WriteLine(value: int64) : unit
   (+0 other overloads)
Console.WriteLine(value: uint32) : unit
   (+0 other overloads)
Console.WriteLine(value: int) : unit
   (+0 other overloads)
Console.WriteLine(value: float32) : unit
   (+0 other overloads)
Console.WriteLine(value: float) : unit
   (+0 other overloads)
Console.WriteLine(value: decimal) : unit
   (+0 other overloads)
val v : float list * float list
val printfn : format:Printf.TextWriterFormat<'T> -> 'T

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.printfn
val snd : tuple:('T1 * 'T2) -> 'T2

Full name: Microsoft.FSharp.Core.Operators.snd
val length : list:'T list -> int

Full name: Microsoft.FSharp.Collections.List.length

More information

Link:http://fssnip.net/oY
Posted:9 years ago
Author:Antonio Prestes GarcĂ­a
Tags: mathematics , algorithms , scientific computing