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

振動子の減衰系で、rulinalgの練習も兼ねて実装してみました。この場合のケースならいちいちクロージャ使う必要はなかったですね。

kiito.hatenablog.com

若干ハマった点としては、 rulinalg::vector::Vector<f64> * f64 は実行されるものの、 f64 * rulinalg::vector::Vector<f64> はコンパイルエラーになることです。あまり調べてないですが、組み込み型に他の型引数でtraitを実装するってできないのかな。

#[macro_use]
extern crate rulinalg;

use rulinalg::vector::Vector;

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


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


fn main() {
  let mut vsp = Vector::new(vec![X0, V0]);
  let mut xs = vec![];
  let mut ts = vec![];

  for i in 1..1000 {
    let t = (i as f64) * DELTA_T;
    ts.push(t);
    xs.push(vsp[0]);
    vsp += rk4(&vsp, |vsp| {
      &matrix![0.0, 1.0; - K / M, - A / M] * vsp
    });
  }

  plot(ts, xs);
}

fn rk4<F>(vsp: &Vector<f64>, f: F) -> Vector<f64>
    where F : Fn(&Vector<f64>) -> Vector<f64> {
  let k1 = &(f(vsp) * DELTA_T);
  let k2 = &(f(&(vsp + k1 / 2.0)) * DELTA_T);
  let k3 = &(f(&(vsp + k2 / 2.0)) * DELTA_T);
  let k4 = &(f(&(vsp + k3)) * DELTA_T);
  (k1 + k2 * 2.0 + k3 * 2.0 + k4) / 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();
}

結果です。