読者です 読者をやめる 読者になる 読者になる

チラシの裏

プログラミングとか色々

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

学校の課題で改良オイラー法とルンゲクッタ法を実装しました。今回も言語の指定は無かったので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は連立常微分方程式 \frac{d\mathbf{x}}{dt}=\mathbf{f}\left(t,\mathbf{x}\right) の初期値 t=t_0, \mathbf{x}=\mathbf{x_0} における数値解をt=t_0からh刻みでn個返します。

例えば、一階微分方程式\frac{di}{dt}=100-100iの初期値t=0, i=0における数値解を、刻み幅を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で適当にプロットすると下の図のようになります。 f:id:fetburner:20140126203017p:plain