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

中間目標であったLotka-Volterraの方程式がシミュレーションで解けました。

shimaphoto03.com

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


const A: f64 = 0.01;
const B: f64 = 0.05;
const C: f64 = 0.0001;
const X0: f64 = 1000.0;
const Y0: f64 = 100.0;
const DELTA_T: f64 = 0.001;


fn main() {
  let (mut x, mut y) = (X0, Y0);
  let mut xs = vec![];
  let mut ys = vec![];
  let mut ts = vec![];

  for i in 1..500000 {
    let t = (i as f64) * DELTA_T;
    let (dx, dy) = runge_kutta(x, y, |x, y| {
      (
        A * x - C * x * y,
        - B * y + C * x * y
      )
    });
    x += dx;
    y += dy;
    ts.push(t);
    ys.push(y);
    xs.push(x);
  }

  plot(xs, ys);
}

fn runge_kutta<F>(x: f64, v: f64, f: F) -> (f64, f64) 
    where F : Fn(f64, f64) -> (f64, f64) {
  let df = |x, v| {
    let (rx, rv) = f(x, v);
    (rx * DELTA_T, rv * DELTA_T)
  };
  let (k1x, k1v) = df(x, v);
  let (k2x, k2v) = df(x * k1x / 2.0, v + k1v / 2.0); 
  let (k3x, k3v) = df(x + k2x / 2.0, v + k2v / 2.0);
  let (k4x, k4v) = df(x + k3x, v + k3v);
  (
    (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0,
    (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0
  )
}

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();
}

x ウサギ（被食者）の個体数、y キツネ（捕食者）の数が、相平面上でどんな軌道を描いているのか分かります。