ordinary_diffeq/
lib.rs

1#![allow(dead_code)]
2
3pub mod callback;
4pub mod controller;
5pub mod integrator;
6pub mod ode;
7pub mod problem;
8
9pub mod prelude {
10    pub use super::callback::{stop, Callback};
11    pub use super::controller::PIController;
12    pub use super::integrator::dormand_prince::DormandPrince45;
13    pub use super::ode::ODE;
14    pub use super::problem::{Problem, Solution};
15}
16
17#[cfg(test)]
18mod tests {
19    use crate::prelude::*;
20    use approx::assert_relative_eq;
21    use nalgebra::{Vector1, Vector2, Vector6};
22    use std::f64::consts::PI;
23
24    #[test]
25    fn test_readme() {
26        // Define the system (parameters, derivative, and initial state)
27        type Params = (f64, f64); // Gravity and Length of Pendulum
28        let params = (9.81, 1.0);
29
30        fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
31            let &(g, l) = p;
32            let theta = y[0];
33            let d_theta = y[1];
34            Vector2::new(d_theta, -(g / l) * theta.sin())
35        }
36
37        let y0 = Vector2::new(0.0, PI / 2.0);
38
39        // Set up the problem (ODE, Integrator, Controller, and Callbacks)
40        let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
41        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
42        let controller = PIController::default();
43
44        let value_too_high = Callback {
45            event: &|t: f64, _y: Vector2<f64>, _p: &Params| 5.0 - t,
46            effect: &stop,
47        };
48
49        // Solve the problem
50        let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high);
51        let solution = problem.solve();
52
53        // Can interpolate solutions to whatever you want
54        let _interpolated_answer = solution.interpolate(4.4);
55    }
56
57    #[test]
58    fn test_correctness() {
59        // Define the system (parameters, derivative, and initial state)
60        type Params = ();
61        let params = ();
62
63        fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
64            Vector1::new(5.0 * y[0] - 3.0)
65        }
66
67        let y0 = Vector1::new(1.0);
68
69        // Set up the problem (ODE, Integrator, Controller, and Callbacks)
70        let ode = ODE::new(&derivative, 2.0, 3.0, y0, params);
71        let dp45 = DormandPrince45::new();
72        let controller = PIController::default();
73
74        // Solve the problem
75        let mut problem = Problem::new(ode, dp45, controller);
76        let solution = problem.solve();
77        for (time, state) in solution.times.iter().zip(solution.states.iter()) {
78            let exact = 0.4 * (5.0 * (time - 2.0)).exp() + 0.6;
79            assert_relative_eq!(state[0], exact, max_relative = 1e-7);
80        }
81    }
82
83    #[test]
84    fn test_orbit() {
85        // Calculate one period
86        let a = 6.7781363e6_f64;
87        let mu = 3.98600441500000e14;
88        let period = 2.0 * PI * (a.powi(3) / mu).sqrt();
89
90        // Set up the system
91        type Params = (f64,);
92        let params = (mu,);
93        fn derivative(_t: f64, state: Vector6<f64>, p: &Params) -> Vector6<f64> {
94            let acc = -(p.0 * state.fixed_rows::<3>(0)) / (state.fixed_rows::<3>(0).norm().powi(3));
95            Vector6::new(state[3], state[4], state[5], acc[0], acc[1], acc[2])
96        }
97        let y0 = Vector6::new(
98            4.263868426884883e6,
99            5.146189057155391e6,
100            1.1310208421331816e6,
101            -5923.454461876975,
102            4496.802639690076,
103            1870.3893008991558,
104        );
105
106        // Integrate
107        let ode = ODE::new(&derivative, 0.0, 10.0 * period, y0, params);
108        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-12);
109        let controller = PIController::new(0.37, 0.04, 10.0, 0.2, 1000.0, 0.9, 0.01);
110        let mut problem = Problem::new(ode, dp45, controller);
111        let solution = problem.solve();
112
113        assert_relative_eq!(
114            solution.times[solution.states.len() - 1],
115            10.0 * period,
116            max_relative = 1e-12
117        );
118        assert_relative_eq!(
119            solution.states[solution.states.len() - 1][0],
120            y0[0],
121            max_relative = 1e-9
122        );
123        assert_relative_eq!(
124            solution.states[solution.states.len() - 1][1],
125            y0[1],
126            max_relative = 1e-9
127        );
128        assert_relative_eq!(
129            solution.states[solution.states.len() - 1][2],
130            y0[2],
131            max_relative = 1e-9
132        );
133        assert_relative_eq!(
134            solution.states[solution.states.len() - 1][3],
135            y0[3],
136            max_relative = 1e-9
137        );
138        assert_relative_eq!(
139            solution.states[solution.states.len() - 1][4],
140            y0[4],
141            max_relative = 1e-9
142        );
143        assert_relative_eq!(
144            solution.states[solution.states.len() - 1][5],
145            y0[5],
146            max_relative = 1e-9
147        );
148    }
149}
ordinary_diffeq/lib.rs

ordinary_diffeq/
lib.rs
