以前のコードを改造して、 x
と v
それぞれで微分方程式を指定できるようにしました。クロージャ便利ですね。
use gnuplot::{Figure, Caption, Color}; const K: f64 = 1.0; const M: 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 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, |x, v| { ( v, (- K * x - A * v) / M ) }); x += dx; v += dv; ts.push(t); xs.push(x); } plot(ts, xs); } 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(); }