Ver 0.39.9

- Change an implementation of GL4

Author: Axect

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.39.8"
+version = "0.39.9"
 authors = ["axect <axect@outlook.kr>"]
 edition = "2018"
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"

RELEASES.md

@@ -1,3 +1,7 @@
+# Release 0.39.9 (2025-07-10)
+
+- Change implementation of Gauss-Legendre 4th order method
+
 # Release 0.39.8 (2025-06-23)
 
 - Implement `LogNormal` distribution

example_data/gl4_test.png

@@ -1155,6 +1155,18 @@ impl ButcherTableau for RKF78 {
 // ┌─────────────────────────────────────────────────────────┐
 //  Gauss-Legendre 4th order
 // └─────────────────────────────────────────────────────────┘
+
+// Correct coefficients for 4th-order Gauss-Legendre method
+const SQRT3: f64 = 1.7320508075688772;
+const C1: f64 = 0.5 - SQRT3 / 6.0;
+const C2: f64 = 0.5 + SQRT3 / 6.0;
+const A11: f64 = 0.25;
+const A12: f64 = 0.25 - SQRT3 / 6.0;
+const A21: f64 = 0.25 + SQRT3 / 6.0;
+const A22: f64 = 0.25;
+const B1: f64 = 0.5;
+const B2: f64 = 0.5;
+
 /// Enum for implicit solvers.
 ///
 /// This enum defines the available implicit solvers for the Gauss-Legendre 4th order integrator.
@@ -1198,7 +1210,7 @@ impl Default for GL4 {
     fn default() -> Self {
         GL4 {
             solver: ImplicitSolver::FixedPoint,
-            tol: 1e-6,
+            tol: 1e-8,
             max_step_iter: 100,
         }
     }
@@ -1219,41 +1231,40 @@ impl ODEIntegrator for GL4 {
     #[inline]
     fn step<P: ODEProblem>(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result<f64> {
         let n = y.len();
-        let sqrt3 = 3.0_f64.sqrt();
-        let c = 0.5 * (3.0 - sqrt3) / 6.0;
-        let d = 0.5 * (3.0 + sqrt3) / 6.0;
+        //let sqrt3 = 3.0_f64.sqrt();
+        //let c = 0.5 * (3.0 - sqrt3) / 6.0;
+        //let d = 0.5 * (3.0 + sqrt3) / 6.0;
         let mut k1 = vec![0.0; n];
         let mut k2 = vec![0.0; n];
-        let mut y1 = vec![0.0; n];
-        let mut y2 = vec![0.0; n];
+
+        // Initial guess for k1, k2.
+        problem.rhs(t, y, &mut k1)?;
+        k2.copy_from_slice(&k1);
 
         match self.solver {
             ImplicitSolver::FixedPoint => {
                 // Fixed-point iteration
+                let mut y1 = vec![0.0; n];
+                let mut y2 = vec![0.0; n];
+
                 for _ in 0..self.max_step_iter {
+                    let k1_old = k1.clone();
+                    let k2_old = k2.clone();
+
                     for i in 0..n {
-                        y1[i] = y[i] + dt * (c * k1[i] + d * k2[i] - sqrt3 * (k2[i] - k1[i]) / 2.0);
-                        y2[i] = y[i] + dt * (c * k1[i] + d * k2[i] + sqrt3 * (k2[i] - k1[i]) / 2.0);
+                        y1[i] = y[i] + dt * (A11 * k1[i] + A12 * k2[i]);
+                        y2[i] = y[i] + dt * (A11 * k1[i] + A12 * k2[i]);
                     }
 
-                    problem.rhs(t + c * dt, &y1, &mut k1)?;
-                    problem.rhs(t + d * dt, &y2, &mut k2)?;
+                    // Compute new k1 and k2
+                    problem.rhs(t + C1 * dt, &y1, &mut k1)?;
+                    problem.rhs(t + C2 * dt, &y2, &mut k2)?;
 
+                    // Check for convergence
                     let mut max_diff = 0f64;
                     for i in 0..n {
-                        max_diff = max_diff
-                            .max(
-                                (y1[i]
-                                    - y[i]
-                                    - dt * (c * k1[i] + d * k2[i] - sqrt3 * (k2[i] - k1[i]) / 2.0))
-                                    .abs(),
-                            )
-                            .max(
-                                (y2[i]
-                                    - y[i]
-                                    - dt * (c * k1[i] + d * k2[i] + sqrt3 * (k2[i] - k1[i]) / 2.0))
-                                    .abs(),
-                            );
+                        max_diff = max_diff.max((k1[i] - k1_old[i]).abs());
+                        max_diff = max_diff.max((k2[i] - k2_old[i]).abs());
                     }
 
                     if max_diff < self.tol {
@@ -1264,43 +1275,28 @@ impl ODEIntegrator for GL4 {
             ImplicitSolver::Broyden => {
                 let m = 2 * n;
                 let mut U = vec![0.0; m];
-
-                // Initial Guess via Fixed-Point Iteration
-                {
-                    let mut k1 = vec![0.0; n];
-                    let mut k2 = vec![0.0; n];
-                    let mut y1 = vec![0.0; n];
-                    let mut y2 = vec![0.0; n];
-
-                    y1.copy_from_slice(y);
-                    y2.copy_from_slice(y);
-                    problem.rhs(t + c * dt, &y1, &mut k1)?;
-                    problem.rhs(t + d * dt, &y2, &mut k2)?;
-                    for i in 0..n {
-                        U[i] = k1[i];
-                        U[n + i] = k2[i];
-                    }
-                }
+                U[..n].copy_from_slice(&k1);
+                U[n..].copy_from_slice(&k2);
 
                 // F_vec = F(U)
                 let mut F_vec = vec![0.0; m];
-                compute_F(problem, t, y, dt, c, d, sqrt3, &U, &mut F_vec)?;
+                compute_F(problem, t, y, dt, &U, &mut F_vec)?;
 
                 // Initialize inverse Jacobian matrix
                 let mut J_inv = eye(m);
 
                 // Repeat Broyden's method
                 for _ in 0..self.max_step_iter {
                     // delta = - J_inv * F_vec
-                    let mut delta = (&J_inv * &F_vec).mul_scalar(-1.0);
+                    let delta = (&J_inv * &F_vec).mul_scalar(-1.0);
 
                     // U <- U + delta
                     U.iter_mut()
-                        .zip(delta.iter_mut())
+                        .zip(delta.iter())
                         .for_each(|(u, d)| *u += *d);
 
                     let mut F_new = vec![0.0; m];
-                    compute_F(problem, t, y, dt, c, d, sqrt3, &U, &mut F_new)?;
+                    compute_F(problem, t, y, dt, &U, &mut F_new)?;
 
                     // If infinity norm of F_new is less than tol, break
                     if F_new.norm(Norm::LInf) < self.tol {
@@ -1312,12 +1308,13 @@ impl ODEIntegrator for GL4 {
 
                     // J_inv * delta_F
                     let J_inv_delta_F = &J_inv * &delta_F;
+
                     let denom = delta.dot(&J_inv_delta_F);
                     if denom.abs() < 1e-12 {
                         break;
                     }
 
-                    // Broyden "good" update
+                    // Broyden's "good" update for the inverse Jacobian
                     // J_inv <- J_inv + ((delta - J_inv * delta_F) * delta^T * J_inv) / denom
                     let delta_minus_J_inv_delta_F = delta.sub_vec(&J_inv_delta_F).to_col();
                     let delta_T_J_inv = &delta.to_row() * &J_inv;
@@ -1326,57 +1323,87 @@ impl ODEIntegrator for GL4 {
                     F_vec = F_new;
                 }
 
-                // Extract k1 and k2
-                for i in 0..n {
-                    k1[i] = U[i];
-                    k2[i] = U[n + i];
-                }
+                k1.copy_from_slice(&U[..n]);
+                k2.copy_from_slice(&U[n..]);
             }
         }
 
         for i in 0..n {
-            y[i] += dt * 0.5 * (k1[i] + k2[i]);
+            y[i] += dt * (B1 * k1[i] + B2 * k2[i]);
         }
 
         Ok(dt)
     }
 }
 
-// Helper function to compute the function F(U) for the implicit solver.
-// y1 = y + dt*(c*k1 + d*k2 - sqrt3/2*(k2-k1))
-// y2 = y + dt*(c*k1 + d*k2 + sqrt3/2*(k2-k1))
+//// Helper function to compute the function F(U) for the implicit solver.
+//// y1 = y + dt*(c*k1 + d*k2 - sqrt3/2*(k2-k1))
+//// y2 = y + dt*(c*k1 + d*k2 + sqrt3/2*(k2-k1))
+//#[allow(non_snake_case)]
+//fn compute_F<P: ODEProblem>(
+//    problem: &P,
+//    t: f64,
+//    y: &[f64],
+//    dt: f64,
+//    c: f64,
+//    d: f64,
+//    sqrt3: f64,
+//    U: &[f64],
+//    F: &mut [f64],
+//) -> Result<()> {
+//    let n = y.len();
+//    let mut y1 = vec![0.0; n];
+//    let mut y2 = vec![0.0; n];
+//
+//    for i in 0..n {
+//        let k1 = U[i];
+//        let k2 = U[n + i];
+//        y1[i] = y[i] + dt * (c * k1 + d * k2 - sqrt3 * (k2 - k1) / 2.0);
+//        y2[i] = y[i] + dt * (c * k1 + d * k2 + sqrt3 * (k2 - k1) / 2.0);
+//    }
+//
+//    let mut f1 = vec![0.0; n];
+//    let mut f2 = vec![0.0; n];
+//    problem.rhs(t + c * dt, &y1, &mut f1)?;
+//    problem.rhs(t + d * dt, &y2, &mut f2)?;
+//
+//    // F = [ k1 - f1, k2 - f2 ]
+//    for i in 0..n {
+//        F[i] = U[i] - f1[i];
+//        F[n + i] = U[n + i] - f2[i];
+//    }
+//    Ok(())
+//}
+/// Helper function to compute the residual F(U) = U - f(y + dt*A*U)
 #[allow(non_snake_case)]
 fn compute_F<P: ODEProblem>(
     problem: &P,
     t: f64,
     y: &[f64],
     dt: f64,
-    c: f64,
-    d: f64,
-    sqrt3: f64,
-    U: &[f64],
+    U: &[f64], // U is a concatenated vector [k1, k2]
     F: &mut [f64],
 ) -> Result<()> {
     let n = y.len();
+    let (k1_slice, k2_slice) = U.split_at(n);
+
     let mut y1 = vec![0.0; n];
     let mut y2 = vec![0.0; n];
 
     for i in 0..n {
-        let k1 = U[i];
-        let k2 = U[n + i];
-        y1[i] = y[i] + dt * (c * k1 + d * k2 - sqrt3 * (k2 - k1) / 2.0);
-        y2[i] = y[i] + dt * (c * k1 + d * k2 + sqrt3 * (k2 - k1) / 2.0);
+        y1[i] = y[i] + dt * (A11 * k1_slice[i] + A12 * k2_slice[i]);
+        y2[i] = y[i] + dt * (A21 * k1_slice[i] + A22 * k2_slice[i]);
     }
+    
+    // F is an output parameter, its parts f1 and f2 are stored temporarily
+    let (f1, f2) = F.split_at_mut(n);
+    problem.rhs(t + C1 * dt, &y1, f1)?;
+    problem.rhs(t + C2 * dt, &y2, f2)?;
 
-    let mut f1 = vec![0.0; n];
-    let mut f2 = vec![0.0; n];
-    problem.rhs(t + c * dt, &y1, &mut f1)?;
-    problem.rhs(t + d * dt, &y2, &mut f2)?;
-
-    // F = [ k1 - f1, k2 - f2 ]
+    // Compute final residual F = [k1 - f1, k2 - f2]
     for i in 0..n {
-        F[i] = U[i] - f1[i];
-        F[n + i] = U[n + i] - f2[i];
+        f1[i] = k1_slice[i] - f1[i];
+        f2[i] = k2_slice[i] - f2[i];
     }
     Ok(())
 }