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

use nalgebra::SVector;

use super::super::ode::ODE;
use super::Integrator;

/// Integrator Trait
pub trait RosenbrockIntegrator<'a> {
    const GAMMA: f64;
    const A: &'a [f64];
    const B: &'a [f64];
    const C: &'a [f64];
    const C2: &'a [f64];
    const D: &'a [f64];
}

pub struct Rodas4<const D: usize> {
    k: Vec<SVector<f64,D>>,
    a_tol: f64,
    r_tol: f64,
}

impl<const D: usize> Rodas4<D> where Rodas4<D>: Integrator<D> {
    pub fn new(a_tol: f64, r_tol: f64) -> Self {
        Self {
            k: vec![SVector::<f64,D>::zeros(); Self::STAGES],
            a_tol,
            r_tol,
        }
    }
}

impl<'a, const D: usize> RosenbrockIntegrator<'a> for Rodas4<D> {
    const GAMMA: f64 = 0.25;
    const A: &'a [f64] = &[
        1.544000000000000,
        0.9466785280815826,
        0.2557011698983284,
        3.314825187068521,
        2.896124015972201,
        0.9986419139977817,
        1.221224509226641,
        6.019134481288629,
        12.53708332932087,
        -0.6878860361058950,
    ];
    const B: &'a [f64] = &[
        10.12623508344586,
        -7.487995877610167,
        -34.80091861555747,
        -7.992771707568823,
        1.025137723295662,
        -0.6762803392801253,
        6.087714651680015,
        16.43084320892478,
        24.76722511418386,
        -6.594389125716872,
    ];
    const C: &'a [f64] = &[
        -5.668800000000000,
        -2.430093356833875,
        -0.2063599157091915,
        -0.1073529058151375,
        -9.594562251023355,
        -20.47028614809616,
        7.496443313967647,
        -10.24680431464352,
        -33.99990352819905,
        11.70890893206160,
        8.083246795921522,
        -7.981132988064893,
        -31.52159432874371,
        16.31930543123136,
        -6.058818238834054,
    ];
    const C2: &'a [f64] = &[
        0.0,
        0.386,
        0.21,
        0.63,
    ];
    const D: &'a [f64] = &[
        0.2500000000000000,
        -0.1043000000000000,
        0.1035000000000000,
        -0.03620000000000023,
    ];
}

impl<const D: usize> Integrator<D> for Rodas4<D>
where
    Rodas4<D>: RosenbrockIntegrator,
{
    const STAGES: usize = 6;
    const ADAPTIVE: bool = true;

    // TODO: Finish this
    fn step(&self, ode: &ODE<D>, _h: f64) -> (SVector<f64,D>, Option<f64>) {
        let next_y = ode.y.clone();
        let err = SVector::<f64, D>::zeros();
        (next_y, Some(err.norm()))
    }
}