ordinary-diffeq 0.2.3

A library for solving differential equations based on the DifferentialEquations.jl julia library.
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
149

#![allow(dead_code)]

pub mod callback;
pub mod controller;
pub mod integrator;
pub mod ode;
pub mod problem;

pub mod prelude {
    pub use super::callback::{stop, Callback};
    pub use super::controller::PIController;
    pub use super::integrator::dormand_prince::DormandPrince45;
    pub use super::ode::ODE;
    pub use super::problem::{Problem, Solution};
}

#[cfg(test)]
mod tests {
    use crate::prelude::*;
    use approx::assert_relative_eq;
    use nalgebra::{Vector1, Vector2, Vector6};
    use std::f64::consts::PI;

    #[test]
    fn test_readme() {
        // Define the system (parameters, derivative, and initial state)
        type Params = (f64, f64); // Gravity and Length of Pendulum
        let params = (9.81, 1.0);

        fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
            let &(g, l) = p;
            let theta = y[0];
            let d_theta = y[1];
            Vector2::new(d_theta, -(g / l) * theta.sin())
        }

        let y0 = Vector2::new(0.0, PI / 2.0);

        // Set up the problem (ODE, Integrator, Controller, and Callbacks)
        let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
        let controller = PIController::default();

        let value_too_high = Callback {
            event: &|t: f64, _y: Vector2<f64>, _p: &Params| 5.0 - t,
            effect: &stop,
        };

        // Solve the problem
        let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high);
        let solution = problem.solve();

        // Can interpolate solutions to whatever you want
        let _interpolated_answer = solution.interpolate(4.4);
    }

    #[test]
    fn test_correctness() {
        // Define the system (parameters, derivative, and initial state)
        type Params = ();
        let params = ();

        fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
            Vector1::new(5.0 * y[0] - 3.0)
        }

        let y0 = Vector1::new(1.0);

        // Set up the problem (ODE, Integrator, Controller, and Callbacks)
        let ode = ODE::new(&derivative, 2.0, 3.0, y0, params);
        let dp45 = DormandPrince45::new();
        let controller = PIController::default();

        // Solve the problem
        let mut problem = Problem::new(ode, dp45, controller);
        let solution = problem.solve();
        for (time, state) in solution.times.iter().zip(solution.states.iter()) {
            let exact = 0.4 * (5.0 * (time - 2.0)).exp() + 0.6;
            assert_relative_eq!(state[0], exact, max_relative = 1e-7);
        }
    }

    #[test]
    fn test_orbit() {
        // Calculate one period
        let a = 6.7781363e6_f64;
        let mu = 3.98600441500000e14;
        let period = 2.0 * PI * (a.powi(3) / mu).sqrt();

        // Set up the system
        type Params = (f64,);
        let params = (mu,);
        fn derivative(_t: f64, state: Vector6<f64>, p: &Params) -> Vector6<f64> {
            let acc = -(p.0 * state.fixed_rows::<3>(0)) / (state.fixed_rows::<3>(0).norm().powi(3));
            Vector6::new(state[3], state[4], state[5], acc[0], acc[1], acc[2])
        }
        let y0 = Vector6::new(
            4.263868426884883e6,
            5.146189057155391e6,
            1.1310208421331816e6,
            -5923.454461876975,
            4496.802639690076,
            1870.3893008991558,
        );

        // Integrate
        let ode = ODE::new(&derivative, 0.0, 10.0 * period, y0, params);
        let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-12);
        let controller = PIController::new(0.37, 0.04, 10.0, 0.2, 1000.0, 0.9, 0.01);
        let mut problem = Problem::new(ode, dp45, controller);
        let solution = problem.solve();

        assert_relative_eq!(
            solution.times[solution.states.len() - 1],
            10.0 * period,
            max_relative = 1e-12
        );
        assert_relative_eq!(
            solution.states[solution.states.len() - 1][0],
            y0[0],
            max_relative = 1e-9
        );
        assert_relative_eq!(
            solution.states[solution.states.len() - 1][1],
            y0[1],
            max_relative = 1e-9
        );
        assert_relative_eq!(
            solution.states[solution.states.len() - 1][2],
            y0[2],
            max_relative = 1e-9
        );
        assert_relative_eq!(
            solution.states[solution.states.len() - 1][3],
            y0[3],
            max_relative = 1e-9
        );
        assert_relative_eq!(
            solution.states[solution.states.len() - 1][4],
            y0[4],
            max_relative = 1e-9
        );
        assert_relative_eq!(
            solution.states[solution.states.len() - 1][5],
            y0[5],
            max_relative = 1e-9
        );
    }
}