Twitterの1D-Newton法の速度比較のコード（Rust版）

newton.rs

use std::fs::File;
use std::io::Write;

fn main() {
    let mass = 1.0;
    let k = 1.0;
    let dt = 1e-2;
    let nt = 100000000;

    let mut xt = vec![0.0; nt + 1];
    let mut vt = vec![0.0; nt + 1];

    let mut x = 0.0;
    let mut v = 1.0;

    for it in 0..=nt {
        xt[it] = x;
        vt[it] = v;
        let result = runge_kutta_4th(x, v, dt, mass, k);
        x = result.0;
        v = result.1;
    }

    let mut ofs = File::create("result_rust.out").expect("file not found.");
    for it in nt - 1000..=nt {
        write!(
            ofs,
            "{:.7e} {:.7e} {:.7e}\n",
            (it as f64) * dt,
            xt[it],
            vt[it]
        ).expect("cannot write.");
    }
}

fn runge_kutta_4th(x: f64, v: f64, dt: f64, mass: f64, k: f64) -> (f64, f64) {
    let mut xtmp = x;
    let mut vtmp = v;

    let x1 = v;
    let v1 = force(x, mass, k);

    let x2 = v + 0.5 * dt * v1;
    let v2 = force(x + 0.5 * x1 * dt, mass, k);

    let x3 = v + 0.5 * dt * v2;
    let v3 = force(x + 0.5 * x2 * dt, mass, k);

    let x4 = v + dt * v3;
    let v4 = force(x + x3 * dt, mass, k);

    xtmp += (x1 + 2.0 * x2 + 2.0 * x3 + x4) * dt / 6.0;
    vtmp += (v1 + 2.0 * v2 + 2.0 * v3 + v4) * dt / 6.0;

    (xtmp, vtmp)
}

fn force(x: f64, mass: f64, k: f64) -> f64 {
    -x * k / mass
}