改良オイラー法とルンゲクッタ法

学校の課題で改良オイラー法とルンゲクッタ法を実装しました。今回も言語の指定は無かったのでOCamlです。

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)))

関数revised_euler f h (t0, x0) n及びrunge_kutta f h (t0, x0) nは連立常微分方程式 の初期値 における数値解をからh刻みでn個返します。

例えば、一階微分方程式の初期値における数値解を、刻み幅を0.021, 0.020, 0.019と変化させながオイラー法とルンゲクッタ法で求めて出力するプログラムは次のようになります。

let f t [i] = [100. *. (1. -. i)]
let init = (0., [0.])
let ( <*> ) fs xs = List.flatten @@ List.map (fun f ->
  List.map (fun x -> f x) xs) fs

let () =
  [Diffequ.revised_euler; Diffequ.runge_kutta] <*> [f]
    <*> [0.021; 0.02; 0.019] <*> [init] <*> [10]
  |> List.iter (fun xys ->
      List.iter (fun (x, [y]) -> Printf.printf "%f %f\n" x y) xys;
      print_string "\n\n")

この出力と解析解をgnuplotで適当にプロットすると下の図のようになります。