metaphysics 0.1.0

Generic extendable first-order differential equation solver
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

use glam::{Vec2, Vec3};
use std::ops::{Add, Deref, DerefMut, Mul};

pub trait Parameter: Sized + Copy + Default {
    type Derivative: Sized
        + Copy
        + Add<Output = Self::Derivative>
        + Mul<f32, Output = Self::Derivative>
        + Default;
    fn step(self, deriv: Self::Derivative, dt: f32) -> Self;
}

/// Independent variable.
#[derive(Clone, Copy, Default, Debug)]
pub struct Var<P: Parameter> {
    value: (P, P),
    deriv: (P::Derivative, P::Derivative),
}

impl<P: Parameter> Var<P> {
    pub fn new(value: P) -> Self {
        Var {
            value: (value, Default::default()),
            deriv: Default::default(),
        }
    }
    pub fn into_value(self) -> P {
        self.value.0
    }
    pub fn value(&self) -> &P {
        &self.value.0
    }
    pub fn value_mut(&mut self) -> &mut P {
        &mut self.value.0
    }

    pub fn add_deriv(&mut self, d: P::Derivative) {
        self.deriv.0 = self.deriv.0 + d;
    }
    pub fn deriv(&self) -> &P::Derivative {
        &self.deriv.0
    }
    pub fn deriv_mut(&mut self) -> &mut P::Derivative {
        &mut self.deriv.0
    }
}

impl<P: Parameter> Deref for Var<P> {
    type Target = P;
    fn deref(&self) -> &Self::Target {
        self.value()
    }
}
impl<P: Parameter> DerefMut for Var<P> {
    fn deref_mut(&mut self) -> &mut Self::Target {
        self.value_mut()
    }
}

pub trait Visitor {
    fn apply<P: Parameter>(&mut self, v: &mut Var<P>);
}

pub trait System {
    fn compute_derivs(&mut self, dt: f32);
    fn visit_vars<V: Visitor>(&mut self, visitor: &mut V);
}

impl Parameter for f32 {
    type Derivative = f32;
    fn step(self, deriv: f32, dt: f32) -> Self {
        self + deriv * dt
    }
}
impl Parameter for Vec2 {
    type Derivative = Vec2;
    fn step(self, deriv: Vec2, dt: f32) -> Self {
        self + deriv * dt
    }
}
impl Parameter for Vec3 {
    type Derivative = Vec3;
    fn step(self, deriv: Vec3, dt: f32) -> Self {
        self + deriv * dt
    }
}

/// Temporal differential equation solver.
///
/// Currently uses the Runge–Kutta method (RK4).
pub struct Solver;

struct Rk4Step {
    stage: u32,
    dt: f32,
}

impl Rk4Step {
    fn dt(&self) -> f32 {
        match self.stage {
            0 => self.dt / 2.0,
            1 => self.dt / 2.0,
            2 => self.dt,
            3 => self.dt,
            _ => unreachable!(),
        }
    }
}

impl Visitor for Rk4Step {
    fn apply<P: Parameter>(&mut self, p: &mut Var<P>) {
        let x = &mut p.value;
        let dx = &mut p.deriv;
        let dt = self.dt();

        match self.stage {
            0 => {
                (x.1, dx.1) = (x.0, dx.0);
                x.0 = x.1.step(dx.0, dt);
            }
            1 => {
                dx.1 = dx.1 + dx.0 * 2.0;
                x.0 = x.1.step(dx.0, dt);
            }
            2 => {
                dx.1 = dx.1 + dx.0 * 2.0;
                x.0 = x.1.step(dx.0, dt);
            }
            3 => {
                dx.1 = (dx.1 + dx.0) * (1.0 / 6.0);
                x.0 = x.1.step(dx.1, dt);
            }
            _ => unreachable!(),
        };

        dx.0 = Default::default();
    }
}

impl Solver {
    pub fn solve_step<S: System>(&self, system: &mut S, dt: f32) {
        for stage in 0..4 {
            let mut step = Rk4Step { stage, dt };
            system.compute_derivs(step.dt());
            system.visit_vars(&mut step);
        }
    }
}