【Rust】Rustで数値計算プロジェクトを試す その6

ひとまず、Rust習得の第二の中間目標として「Lotka-Volterra方程式を4次のRunge-Kutta法で近似計算する」ことを置いてみます。いきなり2重振り子を計算するのは、極座標とかラグランジアンみたいな概念が必要っぽくて、「大学時代に触れた覚えがあるけど急に思い出せない」ので後にします。

ひとまずRunge-Kutta法の実装です。

shimaphoto03.com

こちらの調和振動子の実装をパクりました。

qiita.com

所有権周りで苦労すると思ってましたが、Copy Traitが実装されているおかげで、関数に渡した後でも使えるので詰まることはありませんでした。

doc.rust-jp.rs

汎用性をもたせるために、関数ポインタとかいろいろ試してみたいですね。

use gnuplot::{Figure, Caption, Color};


const K: f64 = 1.0;
const M: f64 = 1.0;
const X0: f64 = 0.0;
const V0: f64 = 1.0;
const DELTA_T: f64 = 0.01;


fn main() {
  let (mut x, mut v) = (X0, V0);
  let mut xs = vec![];
  let mut ts = vec![];

  for i in 1..1000 {
    let t = (i as f64) * DELTA_T;
    let (dx, dv) = runge_kutta(x, v, t);
    x += dx;
    v += dv;
    ts.push(t);
    xs.push(x);
  }

  plot(ts, xs);
}

fn runge_kutta(x: f64, v: f64, t: f64) -> (f64, f64) {
  let k1v = f(x, v, t) * DELTA_T / M;
  let k1x = v * DELTA_T;
  let k2x = (v + k1v / 2.0) * DELTA_T; 
  let k2v = f(x + k1x / 2.0, v + k1v / 2.0, t + DELTA_T / 2.0) * DELTA_T / M;
  let k3v = f(x + k2x / 2.0, v + k2v / 2.0, t + DELTA_T / 2.0) * DELTA_T / M;
  let k3x = (v + k2v / 2.0) * DELTA_T;
  let k4v = f(x + k3x, v + k3v, t + DELTA_T) * DELTA_T / M;
  let k4x = (v + k3v) * DELTA_T;
  let v = (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0;
  let x = (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
  (x, v)
}

fn f(x: f64, _v: f64, _t: f64) -> f64 {
  - K * x
}

fn plot(x: Vec<f64>, y: Vec<f64>) {
  let mut fg = Figure::new();
  fg.set_terminal("pngcairo", "output.png");
  fg.axes2d()
      .lines(&x, &y, &[Caption("A line"), Color("black")]);
  fg.show();
}