改良オイラー法およびルンゲ・クッタ法による常微分方程式の解法

diffequ.ml

let rec iterate_aux acc f init = function
  | 0 -> List.rev acc
  | n -> iterate_aux (init :: acc) f (f init) (n - 1)
let iterate f init = iterate_aux [] f init

(* スカラーとベクトルの積 *)
let ( *^ ) x ys = List.map (( *. ) x) ys
(* ベクトル同士の和 *)
let ( +^ ) xs ys = List.map2 ( +. ) xs ys

(* 改良オイラー法 *)
let revised_euler f h = iterate (fun (t, xs) ->
  let k1 = h *^ f t xs in
  let k2 = h *^ f (t +. h) (xs +^ k1) in
  (t +. h, xs +^ 0.5 *^ (k1 +^ k2)))

(* ルンゲ・クッタ法 *)
let runge_kutta f h = iterate (fun (t, xs) ->
  let k1 = h *^ f t xs in
  let k2 = h *^ f (t +. 0.5 *. h) (xs +^ 0.5 *^ k1) in
  let k3 = h *^ f (t +. 0.5 *. h) (xs +^ 0.5 *^ k2) in
  let k4 = h *^ f (t +. h) (xs +^ k3) in
  (t +. h, xs +^ 1. /. 6. *^ (k1 +^ 2. *^ k2 +^ 2. *^ k3 +^ k4)))