diff --git a/russell_ode/README.md b/russell_ode/README.md index efae03f0..534c61fb 100644 --- a/russell_ode/README.md +++ b/russell_ode/README.md @@ -127,17 +127,10 @@ use russell_ode::prelude::*; fn main() -> Result<(), StrError> { // system let ndim = 1; - let system = System::new( - ndim, - |f, x, y, _args: &mut NoArgs| { - f[0] = x + y[0]; - Ok(()) - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(ndim, |f, x, y, _args: &mut NoArgs| { + f[0] = x + y[0]; + Ok(()) + }); // solver let params = Params::new(Method::DoPri8); @@ -235,20 +228,21 @@ See the code [simple_system_with_mass.rs](https://github.com/cpmech/russell/tree ```rust use russell_lab::{vec_approx_eq, StrError, Vector}; use russell_ode::prelude::*; -use russell_sparse::CooMatrix; +use russell_sparse::{CooMatrix, Sym}; fn main() -> Result<(), StrError> { // ODE system let ndim = 3; let jac_nnz = 4; - let mut system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = -y[0] + y[1]; - f[1] = y[0] + y[1]; - f[2] = 1.0 / (1.0 + x); - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = -y[0] + y[1]; + f[1] = y[0] + y[1]; + f[2] = 1.0 / (1.0 + x); + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, alpha * (-1.0)).unwrap(); @@ -257,9 +251,6 @@ fn main() -> Result<(), StrError> { jj.put(1, 1, alpha * (1.0)).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // mass matrix @@ -293,6 +284,7 @@ fn main() -> Result<(), StrError> { println!("{}", solver.stats()); Ok(()) } + ``` The output looks like: diff --git a/russell_ode/examples/pde_1d_heat_spectral_collocation.rs b/russell_ode/examples/pde_1d_heat_spectral_collocation.rs index e2e964e3..c003905d 100644 --- a/russell_ode/examples/pde_1d_heat_spectral_collocation.rs +++ b/russell_ode/examples/pde_1d_heat_spectral_collocation.rs @@ -2,7 +2,7 @@ use plotpy::{Curve, Plot, SuperTitleParams}; use russell_lab::algo::{InterpGrid, InterpLagrange, InterpParams}; use russell_lab::math::{GOLDEN_RATIO, PI}; use russell_lab::{mat_vec_mul, StrError, Vector}; -use russell_ode::{no_jacobian, HasJacobian, Method, NoArgs, OdeSolver, Params, System}; +use russell_ode::{Method, NoArgs, OdeSolver, Params, System}; const PATH_KEY: &str = "/tmp/russell_ode/pde_1d_heat_spectral_collocation"; @@ -80,10 +80,6 @@ fn run( dudt[nn] = 0.0; // homogeneous BCs Ok(()) }, - no_jacobian, - HasJacobian::No, - None, - None, ); // ODE solver diff --git a/russell_ode/examples/simple_ode_single_equation.rs b/russell_ode/examples/simple_ode_single_equation.rs index 6dcbcfde..5cd5a9fb 100644 --- a/russell_ode/examples/simple_ode_single_equation.rs +++ b/russell_ode/examples/simple_ode_single_equation.rs @@ -4,17 +4,10 @@ use russell_ode::prelude::*; fn main() -> Result<(), StrError> { // ODE system let ndim = 1; - let system = System::new( - ndim, - |f, x, y, _args: &mut NoArgs| { - f[0] = x + y[0]; - Ok(()) - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(ndim, |f, x, y, _args: &mut NoArgs| { + f[0] = x + y[0]; + Ok(()) + }); // solver let params = Params::new(Method::DoPri8); diff --git a/russell_ode/examples/simple_system_with_mass.rs b/russell_ode/examples/simple_system_with_mass.rs index 817fd89c..fe752203 100644 --- a/russell_ode/examples/simple_system_with_mass.rs +++ b/russell_ode/examples/simple_system_with_mass.rs @@ -1,19 +1,20 @@ use russell_lab::{vec_approx_eq, StrError, Vector}; use russell_ode::prelude::*; -use russell_sparse::CooMatrix; +use russell_sparse::{CooMatrix, Sym}; fn main() -> Result<(), StrError> { // ODE system let ndim = 3; let jac_nnz = 4; - let mut system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = -y[0] + y[1]; - f[1] = y[0] + y[1]; - f[2] = 1.0 / (1.0 + x); - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = -y[0] + y[1]; + f[1] = y[0] + y[1]; + f[2] = 1.0 / (1.0 + x); + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, alpha * (-1.0)).unwrap(); @@ -22,9 +23,6 @@ fn main() -> Result<(), StrError> { jj.put(1, 1, alpha * (1.0)).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // mass matrix diff --git a/russell_ode/src/enums.rs b/russell_ode/src/enums.rs index f7f2cf0d..99215b96 100644 --- a/russell_ode/src/enums.rs +++ b/russell_ode/src/enums.rs @@ -1,14 +1,3 @@ -/// Indicates whether the analytical Jacobian is available or not -pub enum HasJacobian { - /// Indicates that the analytical Jacobian is available - Yes, - - /// Indicates that the analytical Jacobian is not available - /// - /// **Note:** The [crate::no_jacobian] function can be conveniently used with the `No` option. - No, -} - /// Holds information about the numerical method to solve (approximate) ODEs #[derive(Clone, Copy, Debug)] pub struct Information { diff --git a/russell_ode/src/erk_dense_out.rs b/russell_ode/src/erk_dense_out.rs index 66d6366e..fb8682a9 100644 --- a/russell_ode/src/erk_dense_out.rs +++ b/russell_ode/src/erk_dense_out.rs @@ -2,7 +2,6 @@ use crate::StrError; use crate::{Method, System}; use crate::{DORMAND_PRINCE_5_D, DORMAND_PRINCE_8_AD, DORMAND_PRINCE_8_CD, DORMAND_PRINCE_8_D}; use russell_lab::Vector; -use russell_sparse::CooMatrix; /// Handles the dense output of explicit Runge-Kutta methods pub(crate) struct ErkDenseOut { @@ -58,9 +57,9 @@ impl ErkDenseOut { } /// Updates the data and returns the number of function evaluations - pub(crate) fn update<'a, F, J, A>( + pub(crate) fn update<'a, F, A>( &mut self, - system: &System, + system: &System, x: f64, y: &Vector, h: f64, @@ -70,7 +69,6 @@ impl ErkDenseOut { ) -> Result where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { let mut n_function_eval = 0; @@ -240,7 +238,7 @@ impl ErkDenseOut { #[cfg(test)] mod tests { use super::ErkDenseOut; - use crate::{no_jacobian, HasJacobian, Method, System}; + use crate::{Method, System}; use russell_lab::Vector; #[test] @@ -309,21 +307,14 @@ mod tests { count: usize, fail: usize, } - let mut system = System::new( - 1, - |_f: &mut Vector, _x: f64, _y: &Vector, args: &mut Args| { - args.count += 1; - if args.count == args.fail { - Err("STOP") - } else { - Ok(()) - } - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let mut system = System::new(1, |_f: &mut Vector, _x: f64, _y: &Vector, args: &mut Args| { + args.count += 1; + if args.count == args.fail { + Err("STOP") + } else { + Ok(()) + } + }); let mut out = ErkDenseOut::new(Method::DoPri8, system.ndim).unwrap(); let mut args = Args { count: 0, fail: 1 }; let h = 0.1; diff --git a/russell_ode/src/euler_backward.rs b/russell_ode/src/euler_backward.rs index d5dfede2..1615ba2b 100644 --- a/russell_ode/src/euler_backward.rs +++ b/russell_ode/src/euler_backward.rs @@ -1,19 +1,18 @@ use crate::StrError; use crate::{OdeSolverTrait, Params, System, Workspace}; use russell_lab::{vec_copy, vec_rms_scaled, vec_update, Vector}; -use russell_sparse::{numerical_jacobian, CooMatrix, LinSolver, SparseMatrix}; +use russell_sparse::{numerical_jacobian, LinSolver, SparseMatrix}; /// Implements the backward Euler (implicit) solver (implicit, order 1, unconditionally stable) -pub(crate) struct EulerBackward<'a, F, J, A> +pub(crate) struct EulerBackward<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Holds the parameters params: Params, /// ODE system - system: &'a System, + system: &'a System<'a, F, A>, /// Vector holding the function evaluation /// @@ -36,13 +35,12 @@ where solver: LinSolver<'a>, } -impl<'a, F, J, A> EulerBackward<'a, F, J, A> +impl<'a, F, A> EulerBackward<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Allocates a new instance - pub fn new(params: Params, system: &'a System) -> Self { + pub fn new(params: Params, system: &'a System<'a, F, A>) -> Self { let ndim = system.ndim; let jac_nnz = if params.newton.use_numerical_jacobian { ndim * ndim @@ -63,10 +61,9 @@ where } } -impl<'a, F, J, A> OdeSolverTrait for EulerBackward<'a, F, J, A> +impl<'a, F, A> OdeSolverTrait for EulerBackward<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Enables dense output fn enable_dense_output(&mut self) -> Result<(), StrError> { @@ -115,13 +112,13 @@ where // calculate J_new := h J let kk = self.kk.get_coo_mut().unwrap(); - if self.params.newton.use_numerical_jacobian || !self.system.jac_available { + if self.params.newton.use_numerical_jacobian || self.system.jacobian.is_none() { work.stats.n_function += ndim; let w1 = &mut self.k; // workspace let w2 = &mut self.dy; // workspace numerical_jacobian(kk, h, x_new, y_new, w1, w2, args, &self.system.function)?; } else { - (self.system.jacobian)(kk, h, x_new, y_new, args)?; + (self.system.jacobian.as_ref().unwrap())(kk, h, x_new, y_new, args)?; } // add diagonal entries => calculate K = h J_new - I @@ -190,8 +187,9 @@ where #[cfg(test)] mod tests { use super::EulerBackward; - use crate::{HasJacobian, Method, OdeSolverTrait, Params, Samples, System, Workspace}; + use crate::{Method, OdeSolverTrait, Params, Samples, System, Workspace}; use russell_lab::{array_approx_eq, Vector}; + use russell_sparse::Sym; // Mathematica code: // @@ -383,28 +381,22 @@ mod tests { struct Args { count_f: usize, } - let system = System::new( - 1, - |f, _, _, args: &mut Args| { - f[0] = 1.0; - args.count_f += 1; - if args.count_f == 1 { - Err("f: stop (count = 1)") - } else if args.count_f == 4 { - Err("f: stop (count = 4; num-jacobian)") - } else { - Ok(()) - } - }, - |jj, alpha, _x, _y, _args: &mut Args| { - jj.reset(); - jj.put(0, 0, alpha * (0.0)).unwrap(); - Err("jj: stop") - }, - HasJacobian::Yes, - None, - None, - ); + let mut system = System::new(1, |f, _, _, args: &mut Args| { + f[0] = 1.0; + args.count_f += 1; + if args.count_f == 1 { + Err("f: stop (count = 1)") + } else if args.count_f == 4 { + Err("f: stop (count = 4; num-jacobian)") + } else { + Ok(()) + } + }); + system.set_jacobian(None, Sym::No, |jj, alpha, _x, _y, _args: &mut Args| { + jj.reset(); + jj.put(0, 0, alpha * (0.0)).unwrap(); + Err("jj: stop") + }); let params = Params::new(Method::BwEuler); let mut solver = EulerBackward::new(params, &system); let mut work = Workspace::new(Method::BwEuler); diff --git a/russell_ode/src/euler_forward.rs b/russell_ode/src/euler_forward.rs index da0910d0..42b8f3c9 100644 --- a/russell_ode/src/euler_forward.rs +++ b/russell_ode/src/euler_forward.rs @@ -1,19 +1,17 @@ use crate::StrError; use crate::{OdeSolverTrait, Params, System, Workspace}; use russell_lab::{vec_add, vec_copy, Vector}; -use russell_sparse::CooMatrix; /// Implements the forward Euler method (explicit, order 1, conditionally stable) /// /// **Warning:** This method is interesting for didactic purposes only /// and should not be used in production codes. -pub(crate) struct EulerForward<'a, F, J, A> +pub(crate) struct EulerForward<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// ODE system - system: &'a System, + system: &'a System<'a, F, A>, /// Vector holding the function evaluation /// @@ -24,13 +22,12 @@ where w: Vector, } -impl<'a, F, J, A> EulerForward<'a, F, J, A> +impl<'a, F, A> EulerForward<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Allocates a new instance - pub fn new(system: &'a System) -> Self { + pub fn new(system: &'a System<'a, F, A>) -> Self { let ndim = system.ndim; EulerForward { system, @@ -40,10 +37,9 @@ where } } -impl<'a, F, J, A> OdeSolverTrait for EulerForward<'a, F, J, A> +impl<'a, F, A> OdeSolverTrait for EulerForward<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Enables dense output fn enable_dense_output(&mut self) -> Result<(), StrError> { @@ -87,7 +83,7 @@ where #[cfg(test)] mod tests { use super::EulerForward; - use crate::{no_jacobian, HasJacobian, Method, NoArgs, OdeSolverTrait, Params, Samples, System, Workspace}; + use crate::{Method, NoArgs, OdeSolverTrait, Params, Samples, System, Workspace}; use russell_lab::{array_approx_eq, Vector}; #[test] @@ -172,17 +168,10 @@ mod tests { #[test] fn euler_forward_handles_errors() { - let system = System::new( - 1, - |f, _, _, _: &mut NoArgs| { - f[0] = 1.0; - Err("stop") - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(1, |f, _, _, _: &mut NoArgs| { + f[0] = 1.0; + Err("stop") + }); let mut solver = EulerForward::new(&system); let mut work = Workspace::new(Method::FwEuler); let x = 0.0; diff --git a/russell_ode/src/explicit_runge_kutta.rs b/russell_ode/src/explicit_runge_kutta.rs index ca78534a..b0829fdb 100644 --- a/russell_ode/src/explicit_runge_kutta.rs +++ b/russell_ode/src/explicit_runge_kutta.rs @@ -2,7 +2,6 @@ use crate::constants::*; use crate::StrError; use crate::{detect_stiffness, ErkDenseOut, Information, Method, OdeSolverTrait, Params, System, Workspace}; use russell_lab::{format_fortran, vec_copy, vec_update, Matrix, Vector}; -use russell_sparse::CooMatrix; /// Implements several explicit Runge-Kutta methods /// @@ -21,16 +20,15 @@ use russell_sparse::CooMatrix; /// 2. E. Hairer, G. Wanner (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p -pub(crate) struct ExplicitRungeKutta<'a, F, J, A> +pub(crate) struct ExplicitRungeKutta<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Holds the parameters params: Params, /// ODE system - system: &'a System, + system: &'a System<'a, F, A>, /// Information such as implicit, embedded, etc. info: Information, @@ -80,13 +78,12 @@ where dense_out: Option, } -impl<'a, F, J, A> ExplicitRungeKutta<'a, F, J, A> +impl<'a, F, A> ExplicitRungeKutta<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Allocates a new instance - pub fn new(params: Params, system: &'a System) -> Result { + pub fn new(params: Params, system: &'a System<'a, F, A>) -> Result { // Runge-Kutta coefficients #[rustfmt::skip] let (aa, bb, cc) = match params.method { @@ -156,10 +153,9 @@ where } } -impl<'a, F, J, A> OdeSolverTrait for ExplicitRungeKutta<'a, F, J, A> +impl<'a, F, A> OdeSolverTrait for ExplicitRungeKutta<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Enables dense output fn enable_dense_output(&mut self) -> Result<(), StrError> { @@ -380,7 +376,7 @@ where #[cfg(test)] mod tests { use super::ExplicitRungeKutta; - use crate::{no_jacobian, ErkDenseOut, HasJacobian, Method, OdeSolverTrait, Params, Samples, System, Workspace}; + use crate::{ErkDenseOut, Method, OdeSolverTrait, Params, Samples, System, Workspace}; use russell_lab::{approx_eq, array_approx_eq, Vector}; #[test] @@ -389,14 +385,7 @@ mod tests { for method in Method::erk_methods() { println!("\n... {:?} ...", method); let params = Params::new(method); - let system = System::new( - 1, - |_, _, _, _args: &mut Args| Ok(()), - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(1, |_, _, _, _args: &mut Args| Ok(())); let erk = ExplicitRungeKutta::new(params, &system).unwrap(); let nstage = erk.nstage; assert_eq!(erk.aa.dims(), (nstage, nstage)); @@ -698,18 +687,11 @@ mod tests { // // y(x) = tan(x) + x + 1 // ``` - let system = System::new( - 1, - |f, x, y, _: &mut u8| { - let d = y[0] - x - 1.0; - f[0] = d * d + 2.0; - Ok(()) - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(1, |f, x, y, _: &mut u8| { + let d = y[0] - x - 1.0; + f[0] = d * d + 2.0; + Ok(()) + }); // allocate solver let params = Params::new(Method::Fehlberg4); @@ -741,30 +723,23 @@ mod tests { struct Args { count_f: usize, } - let system = System::new( - 1, - |f, _, _, args: &mut Args| { - f[0] = 1.0; - args.count_f += 1; - if args.count_f == 1 { - Err("f: count = 1") - } else if args.count_f == 3 { - Err("f: count = 3") - } else if args.count_f == 4 { - Err("f: count = 4 (dense output)") - } else if args.count_f == 6 { - Err("f: count = 6 (dense output)") - } else if args.count_f == 8 { - Err("f: count = 8 (dense output)") - } else { - Ok(()) - } - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(1, |f, _, _, args: &mut Args| { + f[0] = 1.0; + args.count_f += 1; + if args.count_f == 1 { + Err("f: count = 1") + } else if args.count_f == 3 { + Err("f: count = 3") + } else if args.count_f == 4 { + Err("f: count = 4 (dense output)") + } else if args.count_f == 6 { + Err("f: count = 6 (dense output)") + } else if args.count_f == 8 { + Err("f: count = 8 (dense output)") + } else { + Ok(()) + } + }); assert_eq!( ExplicitRungeKutta::new(Params::new(Method::Radau5), &system).err(), diff --git a/russell_ode/src/lib.rs b/russell_ode/src/lib.rs index 77d32fec..703c9a9d 100644 --- a/russell_ode/src/lib.rs +++ b/russell_ode/src/lib.rs @@ -58,7 +58,7 @@ //! //! where `[J]` is the Jacobian matrix. //! -//! **Note:** The Jacobian function is not required for explicit Runge-Kutta methods (see [Method] and [Information]). Thus, one may simply pass the [no_jacobian] function and set [HasJacobian::No] in the system. +//! **Note:** The Jacobian function is not required for explicit Runge-Kutta methods (see [Method] and [Information]). //! //! The flag [ParamsNewton::use_numerical_jacobian] may be set to true to compute the Jacobian matrix numerically. This option works with or without specifying the analytical Jacobian function. //! @@ -130,22 +130,23 @@ //! Code: //! //! ``` -//! use russell_lab::{vec_max_abs_diff, StrError, Vector}; +//! use russell_lab::{StrError, Vector}; //! use russell_ode::prelude::*; -//! use russell_sparse::CooMatrix; +//! use russell_sparse::{CooMatrix, Sym}; //! //! fn main() -> Result<(), StrError> { //! // DAE system //! let ndim = 3; //! let jac_nnz = 4; // number of non-zero values in the Jacobian -//! let mut system = System::new( -//! ndim, -//! |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { -//! f[0] = -y[0] + y[1]; -//! f[1] = y[0] + y[1]; -//! f[2] = 1.0 / (1.0 + x); -//! Ok(()) -//! }, +//! let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { +//! f[0] = -y[0] + y[1]; +//! f[1] = y[0] + y[1]; +//! f[2] = 1.0 / (1.0 + x); +//! Ok(()) +//! }); +//! system.set_jacobian( +//! Some(jac_nnz), +//! Sym::No, //! move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { //! jj.reset(); //! jj.put(0, 0, alpha * (-1.0)).unwrap(); @@ -154,9 +155,6 @@ //! jj.put(1, 1, alpha * (1.0)).unwrap(); //! Ok(()) //! }, -//! HasJacobian::Yes, -//! Some(jac_nnz), -//! None, //! ); //! //! // mass matrix diff --git a/russell_ode/src/ode_solver.rs b/russell_ode/src/ode_solver.rs index 0a997b33..18ac356c 100644 --- a/russell_ode/src/ode_solver.rs +++ b/russell_ode/src/ode_solver.rs @@ -3,7 +3,6 @@ use crate::{EulerBackward, EulerForward, ExplicitRungeKutta, Radau5}; use crate::{Method, OdeSolverTrait, Params, Stats, System, Workspace}; use crate::{Output, StrError}; use russell_lab::{vec_all_finite, Vector}; -use russell_sparse::CooMatrix; /// Implements a numerical solver for systems of ODEs /// @@ -32,8 +31,7 @@ use russell_sparse::CooMatrix; /// where `[J]` is the scaled Jacobian matrix and `α` is a scaling coefficient. /// /// **Note:** The Jacobian function is not required for explicit Runge-Kutta methods -/// (see [crate::Method] and [crate::Information]). Thus, one may simply pass the [crate::no_jacobian] -/// function and set [crate::HasJacobian::No] in the system. +/// (see [crate::Method] and [crate::Information]). /// /// The flag [crate::ParamsNewton::use_numerical_jacobian] may be set to true to compute the /// Jacobian matrix numerically. This option works with or without specifying the analytical @@ -67,26 +65,21 @@ use russell_sparse::CooMatrix; /// ``` /// use russell_lab::{vec_approx_eq, StrError, Vector}; /// use russell_ode::prelude::*; +/// use russell_sparse::Sym; /// /// fn main() -> Result<(), StrError> { /// // ODE system /// let ndim = 1; /// let jac_nnz = 1; -/// let system = System::new( -/// ndim, -/// |f, x, y, _args: &mut NoArgs| { -/// f[0] = x + y[0]; -/// Ok(()) -/// }, -/// |jj, alpha, _x, _y, _args: &mut NoArgs| { -/// jj.reset(); -/// jj.put(0, 0, alpha * (1.0))?; -/// Ok(()) -/// }, -/// HasJacobian::Yes, -/// Some(jac_nnz), -/// None, -/// ); +/// let mut system = System::new(ndim, |f, x, y, _args: &mut NoArgs| { +/// f[0] = x + y[0]; +/// Ok(()) +/// }); +/// system.set_jacobian(Some(jac_nnz), Sym::No, |jj, alpha, _x, _y, _args: &mut NoArgs| { +/// jj.reset(); +/// jj.put(0, 0, alpha * (1.0))?; +/// Ok(()) +/// }); /// /// // solver /// let params = Params::new(Method::Radau5); @@ -143,12 +136,10 @@ impl<'a, A> OdeSolver<'a, A> { /// The generic arguments here are: /// /// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- generic argument to assist in the `F` and `J` functions. It may be simply [crate::NoArgs] indicating that no arguments are needed. - pub fn new(params: Params, system: &'a System) -> Result + /// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [crate::NoArgs] indicating that no arguments are needed. + pub fn new(params: Params, system: &'a System<'a, F, A>) -> Result where F: 'a + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: 'a + Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, A: 'a, { if system.mass_matrix.is_some() && params.method != Method::Radau5 { @@ -432,8 +423,8 @@ impl<'a, A> OdeSolver<'a, A> { #[cfg(test)] mod tests { use super::OdeSolver; - use crate::{no_jacobian, HasJacobian, NoArgs, OutCount, OutData, Stats}; use crate::{Method, Params, Samples, System}; + use crate::{NoArgs, OutCount, OutData, Stats}; use russell_lab::{approx_eq, array_approx_eq, vec_approx_eq, Vector}; use russell_sparse::Genie; @@ -838,21 +829,14 @@ mod tests { // system let ndim = 1; - let system = System::new( - ndim, - |f: &mut Vector, _x: f64, _y: &Vector, args: &mut Args| { - if args.f_count == args.f_barrier { - return Err("f: artificial error"); - } - f[0] = 1.0; - args.f_count += 1; - Ok(()) - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(ndim, |f: &mut Vector, _x: f64, _y: &Vector, args: &mut Args| { + if args.f_count == args.f_barrier { + return Err("f: artificial error"); + } + f[0] = 1.0; + args.f_count += 1; + Ok(()) + }); // initial values and final x let x0 = 0.0; diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs index b1032b81..b6d83095 100644 --- a/russell_ode/src/radau5.rs +++ b/russell_ode/src/radau5.rs @@ -3,7 +3,7 @@ use crate::{OdeSolverTrait, Params, System, Workspace}; use russell_lab::math::SQRT_6; use russell_lab::{complex_vec_zip, cpx, format_fortran, vec_copy, Complex64, ComplexVector, Vector}; use russell_sparse::{numerical_jacobian, ComplexCscMatrix, CscMatrix}; -use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, CooMatrix, Genie, LinSolver, SparseMatrix}; +use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, Genie, LinSolver, SparseMatrix}; use std::thread; /// Implements the Radau5 method (Radau IIA) (implicit, order 5, embedded) for ODEs and DAEs @@ -24,16 +24,15 @@ use std::thread; /// 2. E. Hairer, G. Wanner (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p -pub(crate) struct Radau5<'a, F, J, A> +pub(crate) struct Radau5<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Holds the parameters params: Params, /// ODE system - system: &'a System, + system: &'a System<'a, F, A>, /// Holds the Jacobian matrix. J = df/dy jj: SparseMatrix, @@ -119,13 +118,12 @@ where dw12: ComplexVector, // packed (dw1, dw2) } -impl<'a, F, J, A> Radau5<'a, F, J, A> +impl<'a, F, A> Radau5<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Allocates a new instance - pub fn new(params: Params, system: &'a System) -> Self { + pub fn new(params: Params, system: &'a System<'a, F, A>) -> Self { let ndim = system.ndim; let mass_nnz = match system.mass_matrix.as_ref() { Some(mass) => mass.get_info().2, @@ -198,7 +196,7 @@ where } else if !self.jacobian_computed { work.stats.sw_jacobian.reset(); work.stats.n_jacobian += 1; - if self.params.newton.use_numerical_jacobian || !self.system.jac_available { + if self.params.newton.use_numerical_jacobian || self.system.jacobian.is_none() { work.stats.n_function += self.system.ndim; let y_mut = &mut self.w0; // workspace (mutable y) let w1 = &mut self.dw0; // workspace @@ -206,7 +204,7 @@ where vec_copy(y_mut, y).unwrap(); numerical_jacobian(jj, 1.0, x, y_mut, w1, w2, args, &self.system.function)?; } else { - (self.system.jacobian)(jj, 1.0, x, y, args)?; + (self.system.jacobian.as_ref().unwrap())(jj, 1.0, x, y, args)?; } self.jacobian_computed = true; work.stats.stop_sw_jacobian(); @@ -329,10 +327,9 @@ where } } -impl<'a, F, J, A> OdeSolverTrait for Radau5<'a, F, J, A> +impl<'a, F, A> OdeSolverTrait for Radau5<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Enables dense output fn enable_dense_output(&mut self) -> Result<(), StrError> { @@ -736,9 +733,9 @@ const TI: [[f64; 3]; 3] = [ #[cfg(test)] mod tests { use super::Radau5; - use crate::{HasJacobian, Method, OdeSolverTrait, Params, Samples, System, Workspace}; + use crate::{Method, OdeSolverTrait, Params, Samples, System, Workspace}; use russell_lab::{format_fortran, format_scientific, Vector}; - use russell_sparse::Genie; + use russell_sparse::{Genie, Sym}; #[cfg(feature = "with_mumps")] use serial_test::serial; @@ -888,28 +885,22 @@ mod tests { struct Args { count_f: usize, } - let system = System::new( - 1, - |f, _, _, args: &mut Args| { - f[0] = 1.0; - args.count_f += 1; - if args.count_f == 1 { - Err("f: stop (count = 1; initialize)") - } else if args.count_f == 4 { - Err("f: stop (count = 4; num-jacobian)") - } else { - Ok(()) - } - }, - |jj, alpha, _x, _y, _args: &mut Args| { - jj.reset(); - jj.put(0, 0, alpha * (0.0)).unwrap(); - Err("jj: stop") - }, - HasJacobian::Yes, - None, - None, - ); + let mut system = System::new(1, |f, _, _, args: &mut Args| { + f[0] = 1.0; + args.count_f += 1; + if args.count_f == 1 { + Err("f: stop (count = 1; initialize)") + } else if args.count_f == 4 { + Err("f: stop (count = 4; num-jacobian)") + } else { + Ok(()) + } + }); + system.set_jacobian(None, Sym::No, |jj, alpha, _x, _y, _args: &mut Args| { + jj.reset(); + jj.put(0, 0, alpha * (0.0)).unwrap(); + Err("jj: stop") + }); let params = Params::new(Method::Radau5); let mut solver = Radau5::new(params, &system); let mut work = Workspace::new(Method::Radau5); diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index 1f93864f..8931f3e9 100644 --- a/russell_ode/src/samples.rs +++ b/russell_ode/src/samples.rs @@ -1,5 +1,5 @@ use crate::StrError; -use crate::{HasJacobian, NoArgs, PdeDiscreteLaplacian2d, Side, System}; +use crate::{NoArgs, PdeDiscreteLaplacian2d, Side, System}; use russell_lab::math::PI; use russell_lab::Vector; use russell_sparse::{CooMatrix, Genie, Sym}; @@ -36,20 +36,13 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J /// * `y_fn_x` -- is a function to compute the analytical solution - pub fn simple_equation_constant() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn simple_equation_constant<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -58,20 +51,18 @@ impl Samples { // system let ndim = 1; let jac_nnz = 1; // CooMatrix requires at least one value (thus the 0.0 must be stored) - let system = System::new( - ndim, - |f: &mut Vector, _x: f64, _y: &Vector, _args: &mut NoArgs| { - f[0] = 1.0; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, _x: f64, _y: &Vector, _args: &mut NoArgs| { + f[0] = 1.0; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, 0.0 * alpha).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values @@ -144,10 +135,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -157,15 +145,11 @@ impl Samples { /// /// * Mathematica, Numerical Solution of Differential-Algebraic Equations: Solving Systems with a Mass Matrix /// - pub fn simple_system_with_mass_matrix( + pub fn simple_system_with_mass_matrix<'a>( symmetric: bool, genie: Genie, ) -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -178,14 +162,15 @@ impl Samples { // system let ndim = 3; let jac_nnz = if triangular { 3 } else { 4 }; - let mut system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = -y[0] + y[1]; - f[1] = y[0] + y[1]; - f[2] = 1.0 / (1.0 + x); - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = -y[0] + y[1]; + f[1] = y[0] + y[1]; + f[2] = 1.0 / (1.0 + x); + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + sym, move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, alpha * (-1.0)).unwrap(); @@ -196,9 +181,6 @@ impl Samples { jj.put(1, 1, alpha * (1.0)).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - if sym == Sym::No { None } else { Some(sym) }, ); // mass matrix @@ -256,10 +238,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -270,12 +249,8 @@ impl Samples { /// * Hairer E, Nørsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I. /// Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series /// in Computational Mathematics, 528p - pub fn brusselator_ode() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn brusselator_ode<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -284,13 +259,14 @@ impl Samples { // system let ndim = 2; let jac_nnz = 4; - let system = System::new( - ndim, - |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = 1.0 - 4.0 * y[0] + y[0] * y[0] * y[1]; - f[1] = 3.0 * y[0] - y[0] * y[0] * y[1]; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = 1.0 - 4.0 * y[0] + y[0] * y[0] * y[1]; + f[1] = 3.0 * y[0] - y[0] * y[0] * y[1]; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, alpha * (-4.0 + 2.0 * y[0] * y[1])).unwrap(); @@ -299,9 +275,6 @@ impl Samples { jj.put(1, 1, alpha * (-y[0] * y[0])).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values @@ -499,10 +472,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is the discrete Laplacian `PdeDiscreteLaplacian2d` + /// * `system` -- the ODE system /// * `t0` -- initial t /// * `yy0` -- initial `Y0 = [U0, V0]ᵀ` /// * `args: PdeDiscreteLaplacian2d` -- the discrete Laplacian @@ -515,15 +485,15 @@ impl Samples { /// 2. Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p - pub fn brusselator_pde( + pub fn brusselator_pde<'a>( alpha: f64, npoint: usize, second_book: bool, ignore_diffusion: bool, ) -> ( System< + 'a, impl Fn(&mut Vector, f64, &Vector, &mut PdeDiscreteLaplacian2d) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut PdeDiscreteLaplacian2d) -> Result<(), StrError>, PdeDiscreteLaplacian2d, >, f64, @@ -547,34 +517,35 @@ impl Samples { }; // system - let system = System::new( - ndim, - move |f, t, yy, fdm: &mut PdeDiscreteLaplacian2d| { - fdm.loop_over_grid_points(|m, x, y| { - let um = yy[m]; - let vm = yy[s + m]; - let um2 = um * um; - f[m] = 1.0 - 4.4 * um + um2 * vm; - f[s + m] = 3.4 * um - um2 * vm; - if !ignore_diffusion { - fdm.loop_over_coef_mat_row(m, |k, amk| { - let uk = yy[k]; - let vk = yy[s + k]; - f[m] += amk * uk; - f[s + m] += amk * vk; - }); - } - if second_book { - if t >= 1.1 { - let dx = x - 0.3; - let dy = y - 0.6; - let inhomogeneity = if dx * dx + dy * dy <= 0.01 { 5.0 } else { 0.0 }; - f[m] += inhomogeneity; - } + let mut system = System::new(ndim, move |f, t, yy, fdm: &mut PdeDiscreteLaplacian2d| { + fdm.loop_over_grid_points(|m, x, y| { + let um = yy[m]; + let vm = yy[s + m]; + let um2 = um * um; + f[m] = 1.0 - 4.4 * um + um2 * vm; + f[s + m] = 3.4 * um - um2 * vm; + if !ignore_diffusion { + fdm.loop_over_coef_mat_row(m, |k, amk| { + let uk = yy[k]; + let vk = yy[s + k]; + f[m] += amk * uk; + f[s + m] += amk * vk; + }); + } + if second_book { + if t >= 1.1 { + let dx = x - 0.3; + let dy = y - 0.6; + let inhomogeneity = if dx * dx + dy * dy <= 0.01 { 5.0 } else { 0.0 }; + f[m] += inhomogeneity; } - }); - Ok(()) - }, + } + }); + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, move |jj, aa, _x, yy, fdm: &mut PdeDiscreteLaplacian2d| { jj.reset(); let mut nnz_count = 0; @@ -598,9 +569,6 @@ impl Samples { assert_eq!(nnz_count, jac_nnz); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // discrete laplacian @@ -667,10 +635,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `x1` -- final x @@ -682,12 +647,8 @@ impl Samples { /// * Hairer E, Nørsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I. /// Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series /// in Computational Mathematics, 528p - pub fn arenstorf() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn arenstorf<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, f64, @@ -701,19 +662,20 @@ impl Samples { // system let ndim = 4; let jac_nnz = 8; - let system = System::new( - ndim, - |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { - let t0 = (y[0] + MU) * (y[0] + MU) + y[1] * y[1]; - let t1 = (y[0] - MD) * (y[0] - MD) + y[1] * y[1]; - let d0 = t0 * f64::sqrt(t0); - let d1 = t1 * f64::sqrt(t1); - f[0] = y[2]; - f[1] = y[3]; - f[2] = y[0] + 2.0 * y[3] - MD * (y[0] + MU) / d0 - MU * (y[0] - MD) / d1; - f[3] = y[1] - 2.0 * y[2] - MD * y[1] / d0 - MU * y[1] / d1; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { + let t0 = (y[0] + MU) * (y[0] + MU) + y[1] * y[1]; + let t1 = (y[0] - MD) * (y[0] - MD) + y[1] * y[1]; + let d0 = t0 * f64::sqrt(t0); + let d1 = t1 * f64::sqrt(t1); + f[0] = y[2]; + f[1] = y[3]; + f[2] = y[0] + 2.0 * y[3] - MD * (y[0] + MU) / d0 - MU * (y[0] - MD) / d1; + f[3] = y[1] - 2.0 * y[2] - MD * y[1] / d0 - MU * y[1] / d1; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, y: &Vector, _args: &mut NoArgs| { let t0 = (y[0] + MU) * (y[0] + MU) + y[1] * y[1]; let t1 = (y[0] - MD) * (y[0] - MD) + y[1] * y[1]; @@ -749,9 +711,6 @@ impl Samples { jj.put(3, 2, -2.0 * alpha).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values and final x @@ -797,10 +756,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -811,12 +767,8 @@ impl Samples { /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p - pub fn hairer_wanner_eq1() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn hairer_wanner_eq1<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -828,20 +780,18 @@ impl Samples { // system let ndim = 1; let jac_nnz = 1; - let system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = L * (y[0] - f64::cos(x)); - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = L * (y[0] - f64::cos(x)); + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, alpha * L).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values @@ -877,10 +827,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -890,12 +837,8 @@ impl Samples { /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p - pub fn robertson() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn robertson<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -903,14 +846,15 @@ impl Samples { // system let ndim = 3; let jac_nnz = 7; - let system = System::new( - ndim, - |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = -0.04 * y[0] + 1.0e4 * y[1] * y[2]; - f[1] = 0.04 * y[0] - 1.0e4 * y[1] * y[2] - 3.0e7 * y[1] * y[1]; - f[2] = 3.0e7 * y[1] * y[1]; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = -0.04 * y[0] + 1.0e4 * y[1] * y[2]; + f[1] = 0.04 * y[0] - 1.0e4 * y[1] * y[2] - 3.0e7 * y[1] * y[1]; + f[2] = 3.0e7 * y[1] * y[1]; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, -0.04 * alpha).unwrap(); @@ -922,9 +866,6 @@ impl Samples { jj.put(2, 1, 6.0e7 * y[1] * alpha).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values @@ -962,10 +903,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `x1` -- final x @@ -976,15 +914,11 @@ impl Samples { /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p - pub fn van_der_pol( + pub fn van_der_pol<'a>( epsilon: f64, stationary: bool, ) -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, f64, @@ -1008,13 +942,14 @@ impl Samples { // system let ndim = 2; let jac_nnz = 3; - let system = System::new( - ndim, - move |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = y[1]; - f[1] = ((1.0 - y[0] * y[0]) * y[1] - y[0]) / eps; - Ok(()) - }, + let mut system = System::new(ndim, move |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = y[1]; + f[1] = ((1.0 - y[0] * y[0]) * y[1] - y[0]) / eps; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, move |jj: &mut CooMatrix, alpha: f64, _x: f64, y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 1, 1.0 * alpha).unwrap(); @@ -1022,9 +957,6 @@ impl Samples { jj.put(1, 1, alpha * (1.0 - y[0] * y[0]) / eps).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // results @@ -1097,10 +1029,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -1110,12 +1039,8 @@ impl Samples { /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p - pub fn amplifier1t() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn amplifier1t<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -1134,18 +1059,19 @@ impl Samples { // system let ndim = 5; let jac_nnz = 9; - let mut system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - let ue = A * f64::sin(OM * x); - let g12 = BETA * (f64::exp((y[1] - y[2]) / UF) - 1.0); - f[0] = (y[0] - ue) / R; - f[1] = (2.0 * y[1] - UB) / S + GAMMA * g12; - f[2] = y[2] / S - g12; - f[3] = (y[3] - UB) / S + ALPHA * g12; - f[4] = y[4] / S; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + let ue = A * f64::sin(OM * x); + let g12 = BETA * (f64::exp((y[1] - y[2]) / UF) - 1.0); + f[0] = (y[0] - ue) / R; + f[1] = (2.0 * y[1] - UB) / S + GAMMA * g12; + f[2] = y[2] / S - g12; + f[3] = (y[3] - UB) / S + ALPHA * g12; + f[4] = y[4] / S; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, aa: f64, _x: f64, y: &Vector, _args: &mut NoArgs| { let h12 = BETA * f64::exp((y[1] - y[2]) / UF) / UF; jj.reset(); @@ -1160,9 +1086,6 @@ impl Samples { jj.put(4, 4, aa * (1.0 / S)).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // mass matrix @@ -1203,10 +1126,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -1216,12 +1136,8 @@ impl Samples { /// /// * Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H, /// Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley - pub fn kreyszig_eq6_page902() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn kreyszig_eq6_page902<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -1230,20 +1146,18 @@ impl Samples { // system let ndim = 1; let jac_nnz = 1; - let system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = x + y[0]; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = x + y[0]; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, 1.0 * alpha).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values @@ -1284,10 +1198,7 @@ impl Samples { /// /// # Output /// - /// * `system: System` with: - /// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- is `NoArgs` + /// * `system` -- the ODE system /// * `x0` -- initial x /// * `y0` -- initial y /// * `args` -- is a placeholder variable with the arguments to F and J @@ -1297,12 +1208,8 @@ impl Samples { /// /// * Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H, /// Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley - pub fn kreyszig_ex4_page920() -> ( - System< - impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, - NoArgs, - >, + pub fn kreyszig_ex4_page920<'a>() -> ( + System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>, f64, Vector, NoArgs, @@ -1311,13 +1218,14 @@ impl Samples { // system let ndim = 2; let jac_nnz = 3; - let system = System::new( - ndim, - |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { - f[0] = y[1]; - f[1] = -10.0 * y[0] - 11.0 * y[1] + 10.0 * x + 11.0; - Ok(()) - }, + let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { + f[0] = y[1]; + f[1] = -10.0 * y[0] - 11.0 * y[1] + 10.0 * x + 11.0; + Ok(()) + }); + system.set_jacobian( + Some(jac_nnz), + Sym::No, |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| { jj.reset(); jj.put(0, 1, 1.0 * alpha).unwrap(); @@ -1325,9 +1233,6 @@ impl Samples { jj.put(1, 1, -11.0 * alpha).unwrap(); Ok(()) }, - HasJacobian::Yes, - Some(jac_nnz), - None, ); // initial values @@ -1368,7 +1273,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1393,7 +1299,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1412,7 +1319,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1433,7 +1341,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, jac_alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, jac_alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, jac_alpha, &mut args, system.function).unwrap(); @@ -1454,7 +1363,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, jac_alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, jac_alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, jac_alpha, &mut args, system.function).unwrap(); @@ -1478,7 +1388,8 @@ mod tests { // compute the analytical Jacobian matrix let jac_alpha = 2.0; let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, jac_alpha, t, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, jac_alpha, t, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, t, &y0, jac_alpha, &mut args, &system.function).unwrap(); @@ -1498,7 +1409,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1523,7 +1435,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1542,7 +1455,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1561,7 +1475,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1580,7 +1495,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1599,7 +1515,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1631,7 +1548,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); @@ -1656,7 +1574,8 @@ mod tests { // compute the analytical Jacobian matrix let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); - (system.jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); + let jacobian = system.jacobian.as_ref().unwrap(); + (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs index f881b6a0..4df2161a 100644 --- a/russell_ode/src/system.rs +++ b/russell_ode/src/system.rs @@ -1,4 +1,4 @@ -use crate::{HasJacobian, StrError}; +use crate::StrError; use russell_lab::Vector; use russell_sparse::{CooMatrix, Sym}; use std::marker::PhantomData; @@ -40,12 +40,11 @@ pub type NoArgs = u8; /// The generic arguments here are: /// /// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` -/// * `J` -- function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` -/// * `A` -- generic argument to assist in the `F` and `J` functions. It may be simply [NoArgs] indicating that no arguments are needed. +/// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [NoArgs] indicating that no arguments are needed. /// /// # Important /// -/// The implementation requires the `alpha` parameter in the Jacobian function `J` +/// The implementation requires the `alpha` parameter in the Jacobian function /// to scale the Jacobian matrix. For example: /// /// ```text @@ -64,10 +63,9 @@ pub type NoArgs = u8; /// 2. E. Hairer, G. Wanner (2002) Solving Ordinary Differential Equations II. /// Stiff and Differential-Algebraic Problems. Second Revised Edition. /// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p -pub struct System +pub struct System<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// System dimension pub(crate) ndim: usize, @@ -76,10 +74,7 @@ where pub(crate) function: F, /// Jacobian function - pub(crate) jacobian: J, - - /// Indicates whether the analytical Jacobian is available or not - pub(crate) jac_available: bool, + pub(crate) jacobian: Option Result<(), StrError> + 'a>>, /// Number of non-zeros in the Jacobian matrix pub(crate) jac_nnz: usize, @@ -94,10 +89,9 @@ where phantom: PhantomData A>, } -impl<'a, F, J, A> System +impl<'a, F, A> System<'a, F, A> where F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>, - J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>, { /// Allocates a new instance /// @@ -105,96 +99,115 @@ where /// /// * `ndim` -- dimension of the ODE system (number of equations) /// * `function` -- implements the function: `dy/dx = f(x, y)` - /// * `jacobian` -- implements the Jacobian: `J = df/dy` - /// * `has_jacobian` -- indicates that the analytical Jacobian is available (input by `jacobian`) - /// * `jac_nnz` -- the number of non-zeros in the Jacobian; use None to indicate a dense matrix (i.e., nnz = ndim * ndim) - /// * `jac_sym` -- specifies the symmetric flag for the Jacobian and mass matrices + /// + /// **Note:** Even if the (analytical) Jacobian function is not configured, + /// a numerical Jacobian matrix may be computed (see [crate::Params] and [crate::ParamsNewton]). /// /// # Generics /// /// The generic arguments here are: /// /// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)` - /// * `J` -- function to compute the Jacobian: `(jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A)` - /// * `A` -- generic argument to assist in the `F` and `J` functions. It may be simply [NoArgs] indicating that no arguments are needed. + /// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [NoArgs] indicating that no arguments are needed. /// /// # Examples /// - /// ## One equation (ndim = 1) without Jacobian + /// ## One equation (ndim = 1) without Jacobian callback function /// /// ```rust - /// # use russell_ode::prelude::*; - /// let system = System::new( - /// 1, - /// |f, x, y, _args: &mut NoArgs| { + /// use russell_ode::prelude::*; + /// use russell_ode::StrError; + /// + /// fn main() -> Result<(), StrError> { + /// let system = System::new(1, |f, x, y, _args: &mut NoArgs| { /// f[0] = x + y[0]; /// Ok(()) - /// }, - /// no_jacobian, - /// HasJacobian::No, - /// None, - /// None, - /// ); + /// }); + /// Ok(()) + /// } /// ``` /// /// ## Two equation system (ndim = 2) with Jacobian /// /// ```rust - /// # use russell_ode::prelude::*; - /// let ndim = 2; - /// let jac_nnz = 4; - /// let system = System::new( - /// ndim, - /// |f, x, y, _args: &mut NoArgs| { + /// use russell_ode::prelude::*; + /// use russell_ode::StrError; + /// use russell_sparse::Sym; + /// + /// fn main() -> Result<(), StrError> { + /// let ndim = 2; + /// let mut system = System::new(ndim, |f, x, y, _args: &mut NoArgs| { /// f[0] = x + 2.0 * y[0] + 3.0 * y[1]; /// f[1] = x - 4.0 * y[0] - 5.0 * y[1]; /// Ok(()) - /// }, - /// |jj, alpha, _x, _y, _args: &mut NoArgs| { + /// }); + /// + /// let jac_nnz = 4; + /// system.set_jacobian(Some(jac_nnz), Sym::No, |jj, alpha, _x, _y, _args: &mut NoArgs| { /// jj.reset(); /// jj.put(0, 0, alpha * (2.0)).unwrap(); /// jj.put(0, 1, alpha * (3.0)).unwrap(); /// jj.put(1, 0, alpha * (-4.0)).unwrap(); /// jj.put(1, 1, alpha * (-5.0)).unwrap(); /// Ok(()) - /// }, - /// HasJacobian::Yes, - /// Some(jac_nnz), - /// None, - /// ); + /// }); + /// + /// Ok(()) + /// } /// ``` - pub fn new( - ndim: usize, - function: F, - jacobian: J, - has_ana_jacobian: HasJacobian, - jac_nnz: Option, - jac_sym: Option, - ) -> Self { - let jac_available = match has_ana_jacobian { - HasJacobian::Yes => true, - HasJacobian::No => false, - }; + pub fn new(ndim: usize, function: F) -> Self { System { ndim, function, - jacobian, - jac_available, - jac_nnz: if let Some(nnz) = jac_nnz { nnz } else { ndim * ndim }, - jac_sym: if let Some(sym) = jac_sym { sym } else { Sym::No }, + jacobian: None, + jac_nnz: ndim * ndim, + jac_sym: Sym::No, mass_matrix: None, phantom: PhantomData, } } + /// Sets a function to calculate the Jacobian matrix (analytical Jacobian) + /// + /// Use `|jj, alpha, x, y, args|` or `|jj: &mut CooMatrix, alpha: f64, x: f64, y: &Vector, args: &mut A|` + /// + /// # Input + /// + /// * `jac_nnz` -- the number of non-zeros in the Jacobian; use None to indicate a dense matrix (i.e., nnz = ndim * ndim) + /// * `jac_sym` -- specifies the symmetric flag for the Jacobian and mass matrices + /// * `callback` -- the function to calculate the Jacobian matrix + pub fn set_jacobian( + &mut self, + jac_nnz: Option, + jac_sym: Sym, + callback: impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError> + 'a, + ) { + self.jac_nnz = if let Some(nnz) = jac_nnz { + nnz + } else { + self.ndim * self.ndim + }; + self.jac_sym = jac_sym; + self.jacobian = Some(Box::new(callback)); + } + /// Initializes and enables the mass matrix /// - /// **Note:** Later, call [System::mass_put] to "put" elements in the mass matrix. + /// **Important:** The Jacobian callback function must be set first. The symmetry + /// of the mass matrix must be the same as the symmetry of the Jacobian matrix. + /// + /// **Note:** Even if the (analytical) Jacobian function is not configured, + /// a numerical Jacobian matrix may be computed (see [crate::Params] and [crate::ParamsNewton]). /// /// # Input /// /// * `max_nnz` -- Max number of non-zero values + /// + /// Use [System::mass_put] to "put" elements into the mass matrix. pub fn init_mass_matrix(&mut self, max_nnz: usize) -> Result<(), StrError> { + if self.jacobian.is_none() { + return Err("the Jacobian function must be enabled first"); + } self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, self.jac_sym).unwrap()); Ok(()) } @@ -226,24 +239,17 @@ where } } -/// Implements a placeholder function for when the analytical Jacobian is unavailable -/// -/// **Note:** Use this function with the [crate::HasJacobian::No] option. -pub fn no_jacobian(_jj: &mut CooMatrix, _alpha: f64, _x: f64, _y: &Vector, _args: &mut A) -> Result<(), StrError> { - Err("analytical Jacobian is not available") -} - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #[cfg(test)] mod tests { - use super::{no_jacobian, System}; - use crate::{HasJacobian, NoArgs}; + use super::System; + use crate::NoArgs; use russell_lab::Vector; use russell_sparse::{CooMatrix, Sym}; #[test] - fn ode_system_most_none_works() { + fn ode_system_works() { struct Args { n_function_eval: usize, more_data_goes_here: bool, @@ -252,20 +258,13 @@ mod tests { n_function_eval: 0, more_data_goes_here: false, }; - let system = System::new( - 2, - |f, x, y, args: &mut Args| { - args.n_function_eval += 1; - f[0] = -x * y[1]; - f[1] = x * y[0]; - args.more_data_goes_here = true; - Ok(()) - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let system = System::new(2, |f, x, y, args: &mut Args| { + args.n_function_eval += 1; + f[0] = -x * y[1]; + f[1] = x * y[0]; + args.more_data_goes_here = true; + Ok(()) + }); assert_eq!(system.get_ndim(), 2); assert_eq!(system.get_jac_nnz(), 4); // call system function @@ -273,13 +272,8 @@ mod tests { let y = Vector::new(2); let mut k = Vector::new(2); (system.function)(&mut k, x, &y, &mut args).unwrap(); - // call jacobian function - let mut jj = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); - let alpha = 1.0; - assert_eq!( - (system.jacobian)(&mut jj, alpha, x, &y, &mut args), - Err("analytical Jacobian is not available") - ); + // check that jacobian function is none + assert!(system.jacobian.is_none()); // check println!("n_function_eval = {}", args.n_function_eval); assert_eq!(args.n_function_eval, 1); @@ -287,7 +281,7 @@ mod tests { } #[test] - fn ode_system_some_none_works() { + fn ode_system_set_jacobian_works() { struct Args { n_function_eval: usize, n_jacobian_eval: usize, @@ -300,27 +294,21 @@ mod tests { more_data_goes_here_fn: false, more_data_goes_here_jj: false, }; - let mut system = System::new( - 2, - |f, x, y, args: &mut Args| { - args.n_function_eval += 1; - f[0] = -x * y[1]; - f[1] = x * y[0]; - args.more_data_goes_here_fn = true; - Ok(()) - }, - |jj, alpha, x, _y, args: &mut Args| { - args.n_jacobian_eval += 1; - jj.reset(); - jj.put(0, 1, alpha * (-x)).unwrap(); - jj.put(1, 0, alpha * (x)).unwrap(); - args.more_data_goes_here_jj = true; - Ok(()) - }, - HasJacobian::Yes, - Some(2), - None, - ); + let mut system = System::new(2, |f, x, y, args: &mut Args| { + args.n_function_eval += 1; + f[0] = -x * y[1]; + f[1] = x * y[0]; + args.more_data_goes_here_fn = true; + Ok(()) + }); + system.set_jacobian(Some(2), Sym::No, |jj, alpha, x, _y, args: &mut Args| { + args.n_jacobian_eval += 1; + jj.reset(); + jj.put(0, 1, alpha * (-x)).unwrap(); + jj.put(1, 0, alpha * (x)).unwrap(); + args.more_data_goes_here_jj = true; + Ok(()) + }); // analytical_solution: // y[0] = f64::cos(x * x / 2.0) - 2.0 * f64::sin(x * x / 2.0); // y[1] = 2.0 * f64::cos(x * x / 2.0) + f64::sin(x * x / 2.0); @@ -335,7 +323,7 @@ mod tests { // call jacobian function let mut jj = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); let alpha = 1.0; - (system.jacobian)(&mut jj, alpha, x, &y, &mut args).unwrap(); + (system.jacobian.as_ref().unwrap())(&mut jj, alpha, x, &y, &mut args).unwrap(); // check println!("n_function_eval = {}", args.n_function_eval); println!("n_jacobian_eval = {}", args.n_jacobian_eval); @@ -347,17 +335,10 @@ mod tests { #[test] fn ode_system_handles_errors() { - let mut system = System::new( - 1, - |f, _, _, _: &mut NoArgs| { - f[0] = 1.0; - Ok(()) - }, - no_jacobian, - HasJacobian::No, - None, - None, - ); + let mut system = System::new(1, |f, _, _, _: &mut NoArgs| { + f[0] = 1.0; + Ok(()) + }); let mut f = Vector::new(1); let x = 0.0; let y = Vector::new(1); @@ -367,5 +348,9 @@ mod tests { system.mass_put(0, 0, 1.0).err(), Some("mass matrix has not been initialized/enabled") ); + assert_eq!( + system.init_mass_matrix(1).err(), + Some("the Jacobian function must be enabled first") + ); } } diff --git a/russell_ode/zscripts/run-examples.bash b/russell_ode/zscripts/run-examples.bash index 2bba7c66..225fd3eb 100755 --- a/russell_ode/zscripts/run-examples.bash +++ b/russell_ode/zscripts/run-examples.bash @@ -1,7 +1,14 @@ #!/bin/bash +INTEL_MKL=${1:-""} # 0 or 1 to use intel_mkl + +FEAT="" +if [ "${INTEL_MKL}" = "1" ]; then + FEAT="--features intel_mkl,local_suitesparse" +fi + for example in examples/*.rs; do filename="$(basename "$example")" filekey="${filename%%.*}" - cargo run --example $filekey + cargo run --release $FEAT --example $filekey done