Create binary to output energy data for symplectic vs RK4

src/bin/comp_rk4.rs

@@ -0,0 +1,41 @@
+use dpsim::*;
+
+use std::ops::{Add, Mul};
+
+fn rk4<T>(y: T, dt: f64, f: impl Fn(T) -> T) -> T
+where T: Add<Output = T> + Mul<f64, Output = T> + Copy
+{
+	let k1 = f(y);
+	let k2 = f(y + k1 * dt * 0.5);
+	let k3 = f(y + k2 * dt * 0.5);
+	let k4 = f(y + k3 * dt);
+	y + (k1 + k2 * 2.0 + k3 * 2.0 + k4) * (dt / 6.0)
+}
+
+fn rk4_evolve<T>(state: &mut T, h: f64)
+where T: Hamiltonian + Add<Output = T> + Mul<f64, Output = T> + Copy
+{
+	*state = rk4(*state, h, |s| s.dots());
+}
+
+fn main() {
+	let omega = 500.0;
+	let h = 1e-3;
+
+	let mut rk4 = State {
+		q: DoublePendulum { theta: std::f64::consts::PI - 0.01, phi: 0.0 },
+		p: DoublePendulum { theta: 0.0, phi: 0.0 },
+	};
+	let mut aug = Augmented::from(rk4, omega);
+
+	println!("time,rk4,tao");
+	for i in 0..36000000u32 {
+		rk4_evolve(&mut rk4, h);
+		aug.evolve(h);
+		if i % 1024 == 0 {
+			println!(
+				"{},{},{}",
+				h * (i as f64), rk4.energy(), aug.average().energy());
+		}
+	}
+}