pflow_solver/
methods.rs

1//! Runge-Kutta solver methods (Butcher tableaux).
2
3/// An ODE solver method defined by its Butcher tableau.
4#[derive(Debug, Clone)]
5pub struct Solver {
6    pub name: &'static str,
7    pub order: usize,
8    pub c: Vec<f64>,
9    pub a: Vec<Vec<f64>>,
10    pub b: Vec<f64>,
11    pub b_hat: Vec<f64>,
12}
13
14/// Tsitouras 5/4 Runge-Kutta solver.
15///
16/// A 5th order explicit method with embedded 4th order error estimator.
17pub fn tsit5() -> Solver {
18    Solver {
19        name: "Tsit5",
20        order: 5,
21        c: vec![0.0, 0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0],
22        a: vec![
23            vec![],
24            vec![0.161],
25            vec![-0.008480655492356924, 0.335480655492357],
26            vec![2.8971530571054935, -6.359448489975075, 4.362295432869581],
27            vec![
28                5.325864828439257,
29                -11.748883564062828,
30                7.4955393428898365,
31                -0.09249506636175525,
32            ],
33            vec![
34                5.86145544294642,
35                -12.92096931784711,
36                8.159367898576159,
37                -0.071584973281401,
38                -0.028269050394068383,
39            ],
40            vec![
41                0.09646076681806523,
42                0.01,
43                0.4798896504144996,
44                1.379008574103742,
45                -3.290069515436081,
46                2.324710524099774,
47                0.0,
48            ],
49        ],
50        b: vec![
51            0.09646076681806523,
52            0.01,
53            0.4798896504144996,
54            1.379008574103742,
55            -3.290069515436081,
56            2.324710524099774,
57            0.0,
58        ],
59        b_hat: vec![
60            0.001780011052226,
61            0.000816434459657,
62            -0.007880878010262,
63            0.144711007173263,
64            -0.582357165452555,
65            0.458082105929187,
66            1.0 / 66.0,
67        ],
68    }
69}
70
71/// Dormand-Prince 5(4) Runge-Kutta solver (MATLAB's ode45).
72pub fn rk45() -> Solver {
73    Solver {
74        name: "RK45",
75        order: 5,
76        c: vec![0.0, 1.0 / 5.0, 3.0 / 10.0, 4.0 / 5.0, 8.0 / 9.0, 1.0, 1.0],
77        a: vec![
78            vec![],
79            vec![1.0 / 5.0],
80            vec![3.0 / 40.0, 9.0 / 40.0],
81            vec![44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0],
82            vec![
83                19372.0 / 6561.0,
84                -25360.0 / 2187.0,
85                64448.0 / 6561.0,
86                -212.0 / 729.0,
87            ],
88            vec![
89                9017.0 / 3168.0,
90                -355.0 / 33.0,
91                46732.0 / 5247.0,
92                49.0 / 176.0,
93                -5103.0 / 18656.0,
94            ],
95            vec![
96                35.0 / 384.0,
97                0.0,
98                500.0 / 1113.0,
99                125.0 / 192.0,
100                -2187.0 / 6784.0,
101                11.0 / 84.0,
102            ],
103        ],
104        b: vec![
105            35.0 / 384.0,
106            0.0,
107            500.0 / 1113.0,
108            125.0 / 192.0,
109            -2187.0 / 6784.0,
110            11.0 / 84.0,
111            0.0,
112        ],
113        b_hat: vec![
114            35.0 / 384.0 - 5179.0 / 57600.0,
115            0.0,
116            500.0 / 1113.0 - 7571.0 / 16695.0,
117            125.0 / 192.0 - 393.0 / 640.0,
118            -2187.0 / 6784.0 + 92097.0 / 339200.0,
119            11.0 / 84.0 - 187.0 / 2100.0,
120            -1.0 / 40.0,
121        ],
122    }
123}
124
125/// Classic 4th order Runge-Kutta (fixed step).
126pub fn rk4() -> Solver {
127    Solver {
128        name: "RK4",
129        order: 4,
130        c: vec![0.0, 0.5, 0.5, 1.0],
131        a: vec![vec![], vec![0.5], vec![0.0, 0.5], vec![0.0, 0.0, 1.0]],
132        b: vec![1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0],
133        b_hat: vec![0.0, 0.0, 0.0, 0.0],
134    }
135}
136
137/// Forward Euler method (1st order, fixed step).
138pub fn euler() -> Solver {
139    Solver {
140        name: "Euler",
141        order: 1,
142        c: vec![0.0],
143        a: vec![vec![]],
144        b: vec![1.0],
145        b_hat: vec![0.0],
146    }
147}
148
149/// Heun's method / improved Euler (2nd order).
150pub fn heun() -> Solver {
151    Solver {
152        name: "Heun",
153        order: 2,
154        c: vec![0.0, 1.0],
155        a: vec![vec![], vec![1.0]],
156        b: vec![0.5, 0.5],
157        b_hat: vec![0.0, 0.0],
158    }
159}
160
161/// Midpoint method (2nd order).
162pub fn midpoint() -> Solver {
163    Solver {
164        name: "Midpoint",
165        order: 2,
166        c: vec![0.0, 0.5],
167        a: vec![vec![], vec![0.5]],
168        b: vec![0.0, 1.0],
169        b_hat: vec![0.0, 0.0],
170    }
171}
172
173/// Bogacki-Shampine 3(2) method.
174pub fn bs32() -> Solver {
175    Solver {
176        name: "BS32",
177        order: 3,
178        c: vec![0.0, 0.5, 0.75, 1.0],
179        a: vec![
180            vec![],
181            vec![0.5],
182            vec![0.0, 0.75],
183            vec![2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0],
184        ],
185        b: vec![2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0, 0.0],
186        b_hat: vec![
187            2.0 / 9.0 - 7.0 / 24.0,
188            1.0 / 3.0 - 1.0 / 4.0,
189            4.0 / 9.0 - 1.0 / 3.0,
190            -1.0 / 8.0,
191        ],
192    }
193}
pflow_solver/methods.rs

pflow_solver/
methods.rs
