Skip to content

Commit 24e19ba

Browse files
AxectAxect
committed
Merge branch 'features/ode' into dev
2 parents 16e4dbd + 21fe2a1 commit 24e19ba

5 files changed

Lines changed: 275 additions & 0 deletions

File tree

‎example_data/lorenz_rkf78.png‎

329 KB
Loading

‎example_data/rkf78_test.png‎

60.9 KB
Loading

‎examples/lorenz_rkf78.rs‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
use peroxide::fuga::*;
2+
3+
#[allow(unused_variables)]
4+
fn main() -> Result<(), Box<dyn Error>> {
5+
let initial_conditions = vec![10f64, 1f64, 1f64];
6+
let rkf45 = RKF78::new(1e-4, 0.9, 1e-6, 1e-2, 100);
7+
let basic_ode_solver = BasicODESolver::new(rkf45);
8+
let (_, y_vec) = basic_ode_solver.solve(&Lorenz, (0f64, 100f64), 1e-2, &initial_conditions)?;
9+
let y_mat = py_matrix(y_vec);
10+
let y0 = y_mat.col(0);
11+
let y2 = y_mat.col(2);
12+
13+
#[cfg(feature = "plot")]
14+
{
15+
let mut plt = Plot2D::new();
16+
plt.set_domain(y0)
17+
.insert_image(y2)
18+
.set_xlabel(r"$y_0$")
19+
.set_ylabel(r"$y_2$")
20+
.set_style(PlotStyle::Nature)
21+
.tight_layout()
22+
.set_dpi(600)
23+
.set_path("example_data/lorenz_rkf78.png")
24+
.savefig()?;
25+
}
26+
27+
Ok(())
28+
}
29+
30+
struct Lorenz;
31+
32+
impl ODEProblem for Lorenz {
33+
fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
34+
dy[0] = 10f64 * (y[1] - y[0]);
35+
dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];
36+
dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1];
37+
Ok(())
38+
}
39+
}

‎examples/ode_test_rkf78.rs‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
use peroxide::fuga::*;
2+
3+
#[allow(unused_variables)]
4+
fn main() -> Result<(), Box<dyn Error>> {
5+
let initial_conditions = vec![1f64];
6+
let rkf = RKF78::new(1e-4, 0.9, 1e-6, 1e-1, 100);
7+
let basic_ode_solver = BasicODESolver::new(rkf);
8+
let (t_vec, y_vec) = basic_ode_solver.solve(&Test, (0f64, 10f64), 0.01, &initial_conditions)?;
9+
let y_vec: Vec<f64> = y_vec.into_iter().flatten().collect();
10+
11+
#[cfg(feature = "plot")]
12+
{
13+
let mut plt = Plot2D::new();
14+
plt.set_domain(t_vec)
15+
.insert_image(y_vec)
16+
.set_xlabel(r"$t$")
17+
.set_ylabel(r"$y$")
18+
.set_style(PlotStyle::Nature)
19+
.tight_layout()
20+
.set_dpi(600)
21+
.set_path("example_data/rkf78_test.png")
22+
.savefig()?;
23+
}
24+
25+
Ok(())
26+
}
27+
28+
struct Test;
29+
30+
impl ODEProblem for Test {
31+
fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
32+
Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp())
33+
}
34+
}

‎src/numerical/ode.rs‎

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -895,6 +895,208 @@ impl ButcherTableau for TSIT45 {
895895
}
896896
}
897897

898+
/// Runge-Kutta-Fehlberg 7/8th order integrator.
899+
///
900+
/// This integrator uses the Runge-Kutta-Fehlberg 7(8) method, an adaptive step size integrator.
901+
/// It evaluates f(x,y) thirteen times per step, using embedded 7th and 8th order
902+
/// Runge-Kutta estimates to estimate the solution and the error.
903+
/// The 7th order solution is propagated, and the difference between the 8th and 7th
904+
/// order solutions is used for error estimation and step size control.
905+
///
906+
/// # Member variables
907+
///
908+
/// - `tol`: The tolerance for the estimated error.
909+
/// - `safety_factor`: The safety factor for the step size adjustment.
910+
/// - `min_step_size`: The minimum step size.
911+
/// - `max_step_size`: The maximum step size.
912+
/// - `max_step_iter`: The maximum number of iterations per step.
913+
#[derive(Debug, Clone, Copy)]
914+
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
915+
#[cfg_attr(feature = "rkyv", derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize))]
916+
pub struct RKF78 {
917+
pub tol: f64,
918+
pub safety_factor: f64,
919+
pub min_step_size: f64,
920+
pub max_step_size: f64,
921+
pub max_step_iter: usize,
922+
}
923+
924+
impl Default for RKF78 {
925+
fn default() -> Self {
926+
Self {
927+
tol: 1e-7, // Higher precision default for a higher-order method
928+
safety_factor: 0.9,
929+
min_step_size: 1e-10, // Smaller min step for higher order
930+
max_step_size: 1e-1,
931+
max_step_iter: 100,
932+
}
933+
}
934+
}
935+
936+
impl RKF78 {
937+
pub fn new(
938+
tol: f64,
939+
safety_factor: f64,
940+
min_step_size: f64,
941+
max_step_size: f64,
942+
max_step_iter: usize,
943+
) -> Self {
944+
Self {
945+
tol,
946+
safety_factor,
947+
min_step_size,
948+
max_step_size,
949+
max_step_iter,
950+
}
951+
}
952+
}
953+
954+
impl ButcherTableau for RKF78 {
955+
const C: &'static [f64] = &[
956+
0.0,
957+
2.0 / 27.0,
958+
1.0 / 9.0,
959+
1.0 / 6.0,
960+
5.0 / 12.0,
961+
1.0 / 2.0,
962+
5.0 / 6.0,
963+
1.0 / 6.0,
964+
2.0 / 3.0,
965+
1.0 / 3.0,
966+
1.0,
967+
0.0, // k12 is evaluated at x[i]
968+
1.0, // k13 is evaluated at x[i]+h
969+
];
970+
971+
const A: &'static [&'static [f64]] = &[
972+
// k1
973+
&[],
974+
// k2
975+
&[2.0 / 27.0],
976+
// k3
977+
&[1.0 / 36.0, 3.0 / 36.0],
978+
// k4
979+
&[1.0 / 24.0, 0.0, 3.0 / 24.0],
980+
// k5
981+
&[20.0 / 48.0, 0.0, -75.0 / 48.0, 75.0 / 48.0],
982+
// k6
983+
&[1.0 / 20.0, 0.0, 0.0, 5.0 / 20.0, 4.0 / 20.0],
984+
// k7
985+
&[-25.0 / 108.0, 0.0, 0.0, 125.0 / 108.0, -260.0 / 108.0, 250.0 / 108.0],
986+
// k8
987+
&[31.0 / 300.0, 0.0, 0.0, 0.0, 61.0 / 225.0, -2.0 / 9.0, 13.0 / 900.0],
988+
// k9
989+
&[2.0, 0.0, 0.0, -53.0 / 6.0, 704.0 / 45.0, -107.0 / 9.0, 67.0 / 90.0, 3.0],
990+
// k10
991+
&[
992+
-91.0 / 108.0,
993+
0.0,
994+
0.0,
995+
23.0 / 108.0,
996+
-976.0 / 135.0,
997+
311.0 / 54.0,
998+
-19.0 / 60.0,
999+
17.0 / 6.0,
1000+
-1.0 / 12.0,
1001+
],
1002+
// k11
1003+
&[
1004+
2383.0 / 4100.0,
1005+
0.0,
1006+
0.0,
1007+
-341.0 / 164.0,
1008+
4496.0 / 1025.0,
1009+
-301.0 / 82.0,
1010+
2133.0 / 4100.0,
1011+
45.0 / 82.0,
1012+
45.0 / 164.0,
1013+
18.0 / 41.0,
1014+
],
1015+
// k12
1016+
&[
1017+
3.0 / 205.0,
1018+
0.0,
1019+
0.0,
1020+
0.0,
1021+
0.0,
1022+
-6.0 / 41.0,
1023+
-3.0 / 205.0,
1024+
-3.0 / 41.0,
1025+
3.0 / 41.0,
1026+
6.0 / 41.0,
1027+
0.0,
1028+
],
1029+
// k13
1030+
&[
1031+
-1777.0 / 4100.0,
1032+
0.0,
1033+
0.0,
1034+
-341.0 / 164.0,
1035+
4496.0 / 1025.0,
1036+
-289.0 / 82.0,
1037+
2193.0 / 4100.0,
1038+
51.0 / 82.0,
1039+
33.0 / 164.0,
1040+
12.0 / 41.0,
1041+
0.0,
1042+
1.0,
1043+
],
1044+
];
1045+
1046+
// Coefficients for the 7th order solution (propagated solution)
1047+
// BU_i = BE_i (8th order) - ErrorCoeff_i
1048+
// ErrorCoeff_i = [-41/840, 0, ..., 0, -41/840 (for k11), 41/840 (for k12), 41/840 (for k13)]
1049+
const BU: &'static [f64] = &[
1050+
41.0 / 420.0, // 41/840 - (-41/840)
1051+
0.0,
1052+
0.0,
1053+
0.0,
1054+
0.0,
1055+
34.0 / 105.0,
1056+
9.0 / 35.0,
1057+
9.0 / 35.0,
1058+
9.0 / 280.0,
1059+
9.0 / 280.0,
1060+
41.0 / 420.0, // 41/840 - (-41/840)
1061+
-41.0 / 840.0, // 0.0 - (41/840)
1062+
-41.0 / 840.0, // 0.0 - (41/840)
1063+
];
1064+
1065+
// Coefficients for the 8th order solution (used for error estimation)
1066+
// These are from the y[i+1] formula in the MATLAB description
1067+
const BE: &'static [f64] = &[
1068+
41.0 / 840.0,
1069+
0.0,
1070+
0.0,
1071+
0.0,
1072+
0.0,
1073+
34.0 / 105.0,
1074+
9.0 / 35.0,
1075+
9.0 / 35.0,
1076+
9.0 / 280.0,
1077+
9.0 / 280.0,
1078+
41.0 / 840.0,
1079+
0.0,
1080+
0.0,
1081+
];
1082+
1083+
fn tol(&self) -> f64 {
1084+
self.tol
1085+
}
1086+
fn safety_factor(&self) -> f64 {
1087+
self.safety_factor
1088+
}
1089+
fn min_step_size(&self) -> f64 {
1090+
self.min_step_size
1091+
}
1092+
fn max_step_size(&self) -> f64 {
1093+
self.max_step_size
1094+
}
1095+
fn max_step_iter(&self) -> usize {
1096+
self.max_step_iter
1097+
}
1098+
}
1099+
8981100
// ┌─────────────────────────────────────────────────────────┐
8991101
// Gauss-Legendre 4th order
9001102
// └─────────────────────────────────────────────────────────┘

0 commit comments

Comments
 (0)