qualia_core_db/solvers/calculus/
mechanics.rs

1//! Hamiltonian mechanics primitives with invariant diagnostics.
2
3use super::analysis::{AnalysisError, Vector};
4
5pub fn canonical_poisson_bracket<const N: usize>(
6    df_dq: Vector<N>,
7    df_dp: Vector<N>,
8    dg_dq: Vector<N>,
9    dg_dp: Vector<N>,
10) -> Result<f64, AnalysisError> {
11    Ok(df_dq.dot(dg_dp)? - df_dp.dot(dg_dq)?)
12}
13
14#[repr(C)]
15#[derive(Debug, Clone, Copy, PartialEq)]
16pub struct PhaseState<const N: usize> {
17    pub q: Vector<N>,
18    pub p: Vector<N>,
19}
20
21pub fn stormer_verlet_step<const N: usize, F, G>(
22    state: &mut PhaseState<N>,
23    step: f64,
24    potential_gradient: F,
25    kinetic_gradient: G,
26) -> Result<(), AnalysisError>
27where
28    F: Fn(Vector<N>) -> Vector<N>,
29    G: Fn(Vector<N>) -> Vector<N>,
30{
31    if !step.is_finite() || step == 0.0 {
32        return Err(AnalysisError::InvalidDomain);
33    }
34    state.q.validate()?;
35    state.p.validate()?;
36    let half_force = potential_gradient(state.q).scale(-0.5 * step)?;
37    state.p = state.p.add(half_force);
38    state.q = state.q.add(kinetic_gradient(state.p).scale(step)?);
39    state.p = state.p.add(potential_gradient(state.q).scale(-0.5 * step)?);
40    state.q.validate()?;
41    state.p.validate()
42}
43
44#[repr(C)]
45#[derive(Debug, Clone, Copy, PartialEq)]
46pub struct InvariantDrift {
47    pub initial: f64,
48    pub final_value: f64,
49    pub absolute_drift: f64,
50    pub relative_drift: f64,
51}
52
53pub fn invariant_drift(initial: f64, final_value: f64) -> Result<InvariantDrift, AnalysisError> {
54    if !initial.is_finite() || !final_value.is_finite() {
55        return Err(AnalysisError::NonFinite);
56    }
57    let absolute_drift = (final_value - initial).abs();
58    Ok(InvariantDrift {
59        initial,
60        final_value,
61        absolute_drift,
62        relative_drift: absolute_drift / initial.abs().max(f64::MIN_POSITIVE),
63    })
64}
65
66#[cfg(test)]
67mod tests {
68    use super::*;
69
70    fn oscillator_energy(state: PhaseState<1>) -> f64 {
71        0.5 * (state.q.data[0].powi(2) + state.p.data[0].powi(2))
72    }
73
74    #[test]
75    fn poisson_bracket_is_antisymmetric() {
76        let fq = Vector::new([1.0, 2.0]);
77        let fp = Vector::new([3.0, 4.0]);
78        let gq = Vector::new([-2.0, 1.0]);
79        let gp = Vector::new([0.5, 3.0]);
80        let fg = canonical_poisson_bracket(fq, fp, gq, gp).unwrap();
81        let gf = canonical_poisson_bracket(gq, gp, fq, fp).unwrap();
82        assert_eq!(fg, -gf);
83    }
84
85    #[test]
86    fn stormer_verlet_is_time_reversible_and_bounds_energy_drift() {
87        let initial = PhaseState {
88            q: Vector::new([1.0]),
89            p: Vector::new([0.0]),
90        };
91        let mut state = initial;
92        for _ in 0..10_000 {
93            stormer_verlet_step(&mut state, 0.01, |q| q, |p| p).unwrap();
94        }
95        let drift = invariant_drift(oscillator_energy(initial), oscillator_energy(state)).unwrap();
96        assert!(drift.absolute_drift < 2e-5);
97
98        for _ in 0..10_000 {
99            stormer_verlet_step(&mut state, -0.01, |q| q, |p| p).unwrap();
100        }
101        assert!(state.q.distance(initial.q).unwrap() < 1e-12);
102        assert!(state.p.distance(initial.p).unwrap() < 1e-12);
103    }
104}
qualia_core_db/solvers/calculus/mechanics.rs

qualia_core_db/solvers/calculus/
mechanics.rs
