Skip to content

Instantly share code, notes, and snippets.

@dc1394
Created December 26, 2023 12:02
Show Gist options
  • Save dc1394/0a21860f52507e6c8189a51242a5bdc4 to your computer and use it in GitHub Desktop.
Save dc1394/0a21860f52507e6c8189a51242a5bdc4 to your computer and use it in GitHub Desktop.
Twitterの1D-Newton法の速度比較のコード(Rust版)
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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment