Skip to content

Commit 09ead45

Browse files
AxectAxect
committed
Ver 0.39.6
New ODESolver: RKF78
2 parents e20e55f + fef6992 commit 09ead45

File tree

9 files changed

+288
-2
lines changed

9 files changed

+288
-2
lines changed

‎Cargo.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
[package]
22
name = "peroxide"
3-
version = "0.39.5"
3+
version = "0.39.6"
44
authors = ["axect <[email protected]>"]
55
edition = "2018"
66
description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"

‎README.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -189,9 +189,10 @@ Peroxide can do many things.
189189
- Runge-Kutta 5th order
190190
- Embedded integrator
191191
- Bogacki-Shampine 3(2)
192-
- Runge-Kutta-Fehlberg 4(5)
192+
- Runge-Kutta-Fehlberg 5(4)
193193
- Dormand-Prince 5(4)
194194
- Tsitouras 5(4)
195+
- Runge-Kutta-Fehlberg 8(7)
195196
- Implicit integrator
196197
- Gauss-Legendre 4th order
197198
- Numerical Integration

‎RELEASES.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,8 @@
1+
# Release 0.39.6 (2025-05-16)
2+
3+
- New ODESolver: `RKF78`
4+
- Implement `RKF78` method for `ODESolver`
5+
16
# Release 0.39.5 (2025-04-21)
27

38
- New feature `rkyv`

‎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/lib.rs‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,7 @@
6363
//! - Runge-Kutta-Fehlberg 4(5)
6464
//! - Dormand-Prince 5(4)
6565
//! - Tsitouras 5(4)
66+
//! - Runge-Kutta-Fehlberg 7(8)
6667
//! - Implicit
6768
//! - Gauss-Legendre 4th order
6869
//! - Communication with Python

‎src/numerical/ode.rs‎

Lines changed: 206 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@
2424
//! - Runge-Kutta-Fehlberg 4/5th order (RKF45)
2525
//! - Dormand-Prince 4/5th order (DP45)
2626
//! - Tsitouras 4/5th order (TSIT45)
27+
//! - Runge-Kutta-Fehlberg 7/8th order (RKF78)
2728
//! - **Implicit**
2829
//! - Gauss-Legendre 4th order (GL4)
2930
//!
@@ -895,6 +896,211 @@ impl ButcherTableau for TSIT45 {
895896
}
896897
}
897898

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

0 commit comments

Comments
 (0)