add readme

Cargo.toml

@@ -1,5 +1,5 @@
 [package]
-name = "diffsol"
+name = "diffeq"
 version = "0.1.0"
 edition = "2021"
 

README.md

@@ -0,0 +1,119 @@
+# Diffeq
+
+Diffeq is a library for solving ordinary differential equations (ODEs) or
+semi-explicit differential algebraic equations (DAEs) in Rust. You can use it
+out-of-the-box with vectors and matrices from the
+[nalgebra](https://nalgebra.org) crate, or you can implement your own types by
+implementing the various vector and matrix traits in `diffeq`.
+
+**Note**: This library is still in the early stages of development and is not
+*ready for production use. The API is likely to change in the future.
+
+## Features
+
+Currently only one solver is implemented, the Backward Differentiation Formula
+(BDF) method. This is a variable step-size implicit method that is suitable for
+stiff ODEs and semi-explicit DAEs and is similar to the BDF method in MATLAB's
+`ode15s` solver or the `bdf` solver in SciPy's `solve_ivp` function.
+
+Users can specify the equations to solve in the following ODE form:
+
+```math
+M \frac{dy}{dt} = f(t, y, p)
+```
+
+where `M` is a mass matrix, `y` is the state vector, `t` is the time, `p` is a
+vector of parameters, and `f` is the right-hand side function. The mass matrix
+`M` is optional (assumed to be the identity matrix if not provided).
+
+Both the RHS function `f` and the mass matrix `M` can be specified as functions
+that take the state vector `y`, the parameter vector `p`, and the time `t` as
+arguments. The action of the jacobian `J` of `f` must also be specified as a
+function that takes the state vector `y`, the parameter vector `p`, the time `t`
+and a vector `v` and returns the product `Jv`.
+
+## Installation
+
+Add the following to your `Cargo.toml` file:
+
+```toml
+[dependencies]
+diffeq = "0.1"
+```
+
+or add it on the command line:
+
+```sh
+cargo add diffeq
+```
+
+## Usage
+
+This example solves the Robertson chemical kinetics problem, which is a stiff ODE with the equations:
+
+```math
+\begin{align*}
+\frac{dy_1}{dt} &= -0.04 y_1 + 10^4 y_2 y_3 \\
+\frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3 \times 10^7 y_2^2 \\
+\frac{dy_3}{dt} &= 3 \times 10^7 y_2^2
+\end{align*}
+```
+
+with the initial conditions `y_1(0) = 1`, `y_2(0) = 0`, and `y_3(0) = 0`. We set
+the tolerances to `1.0e-4` for the relative tolerance and `1.0e-8`, `1.0e-6`,
+and `1.0e-6` for the absolute tolerances of `y_1`, `y_2`, and `y_3`,
+respectively. We set the problem up with the following code:
+
+```rust
+type T = f64;
+type V = nalgebra::DVector<T>;
+let p = V::from_vec(vec![0.04.into(), 1.0e4.into(), 3.0e7.into()]);
+let mut problem = OdeSolverProblem::new_ode(
+    // The rhs function `f` 
+    | x: &V, p: &V, _t: T, y: &mut V | {
+        y[0] = -p[0] * x[0] + p[1] * x[1] * x[2];
+        y[1] = p[0] * x[0] - p[1] * x[1] * x[2] - p[2] * x[1] * x[1];
+        y[2] = p[2] * x[1] * x[1];
+    },
+    // The jacobian function `Jv`
+    | x: &V, p: &V, _t: T, v: &V, y: &mut V | {
+        y[0] = -p[0] * v[0] + p[1] * v[1] * x[2] + p[1] * x[1] * v[2];
+        y[1] = p[0] * v[0] - p[1] * v[1] * x[2] - p[1] * x[1] * v[2]  - 2.0 * p[2] * x[1] * v[1];
+        y[2] = 2.0 * p[2] * x[1] * v[1];
+    },
+    // The initial condition
+    | _p: &V, _t: T | {
+        V::from_vec(vec![1.0, 0.0, 0.0])
+    },
+    p,
+);
+problem.rtol = 1.0e-4;
+problem.atol = Rc::new(V::from_vec(vec![1.0e-8, 1.0e-6, 1.0e-6]));
+
+let mut solver = Bdf::default();
+```
+
+We can then solve the problem up to time `t = 0.4` with the following code:
+
+```rust
+let t = 0.4;
+let _y = solver.solve(&problem, t).unwrap();
+```
+
+Or if we want to explicitly step through time and have access to the solution
+state, we can use the following code:
+
+```rust
+let mut state = OdeSolverState::new(&problem);
+while state.t <= t {
+    solver.step(&mut state).unwrap();
+}
+let _y = solver.interpolate(&state, t);
+```
+
+Note that `step` will advance the state to the next time step as chosen by the
+solver, and `interpolate` will interpolate the solution back to the exact time
+`t` that is requested.
+
+
+

src/lib.rs

@@ -31,3 +31,56 @@ use op::{NonLinearOp, LinearOp, ConstantOp};
 use matrix::{Matrix, MatrixViewMut, MatrixCommon};
 use solver::SolverProblem;
 use linear_solver::{lu::LU, LinearSolver};
+use ode_solver::{OdeSolverProblem, OdeSolverState, bdf::Bdf, OdeSolverMethod};
+
+
+#[cfg(test)]
+#[cfg(test)]
+mod tests {
+    use std::rc::Rc;
+
+    use crate::{OdeSolverProblem, Bdf, OdeSolverState, OdeSolverMethod};
+
+    // WARNING: if this test fails and you make a change to the code, you should update the README.md file as well!!!
+    #[test]
+    fn test_readme() {
+        type T = f64;
+        type V = nalgebra::DVector<T>;
+        let p = V::from_vec(vec![0.04.into(), 1.0e4.into(), 3.0e7.into()]);
+        let mut problem = OdeSolverProblem::new_ode(
+            | x: &V, p: &V, _t: T, y: &mut V | {
+                y[0] = -p[0] * x[0] + p[1] * x[1] * x[2];
+                y[1] = p[0] * x[0] - p[1] * x[1] * x[2] - p[2] * x[1] * x[1];
+                y[2] = p[2] * x[1] * x[1];
+            },
+            | x: &V, p: &V, _t: T, v: &V, y: &mut V | {
+                y[0] = -p[0] * v[0] + p[1] * v[1] * x[2] + p[1] * x[1] * v[2];
+                y[1] = p[0] * v[0] - p[1] * v[1] * x[2] - p[1] * x[1] * v[2]  - 2.0 * p[2] * x[1] * v[1];
+                y[2] = 2.0 * p[2] * x[1] * v[1];
+            },
+            | _p: &V, _t: T | {
+                V::from_vec(vec![1.0, 0.0, 0.0])
+            },
+            p,
+        );
+        problem.rtol = 1.0e-4;
+        problem.atol = Rc::new(V::from_vec(vec![1.0e-8, 1.0e-6, 1.0e-6]));
+
+        let mut solver = Bdf::default();
+
+        let t = 0.4;
+        let _y = solver.solve(&problem, t).unwrap();
+
+        let mut state = OdeSolverState::new(&problem);
+        while state.t <= t {
+            solver.step(&mut state).unwrap();
+        }
+        let _y = solver.interpolate(&state, t);
+    }
+}
+
+
+
+
+
+

src/ode_solver/mod.rs

@@ -14,15 +14,17 @@ pub trait OdeSolverMethod<CRhs: NonLinearOp, CMass: LinearOp<M = CRhs::M, V = CR
     fn set_problem(&mut self, state: &mut OdeSolverState<CRhs::M>, problem: OdeSolverProblem<CRhs, CMass, CInit>);
     fn step(&mut self, state: &mut OdeSolverState<CRhs::M>) -> Result<()>;
     fn interpolate(&self, state: &OdeSolverState<CRhs::M>, t: CRhs::T) -> CRhs::V;
-    fn solve(&mut self, problem: OdeSolverProblem<CRhs, CMass, CInit>, t: CRhs::T) -> Result<CRhs::V> {
+    fn solve(&mut self, problem: &OdeSolverProblem<CRhs, CMass, CInit>, t: CRhs::T) -> Result<CRhs::V> {
+        let problem = problem.clone();
         let mut state = OdeSolverState::new(&problem);
         self.set_problem(&mut state, problem);
         while state.t <= t {
             self.step(&mut state)?;
         }
         Ok(self.interpolate(&state, t))
     }
-    fn make_consistent_and_solve<RS: NonLinearSolver<FilterCallable<CRhs>>>(&mut self, problem: OdeSolverProblem<CRhs, CMass, CInit>, t: CRhs::T, root_solver: &mut RS) -> Result<CRhs::V> {
+    fn make_consistent_and_solve<RS: NonLinearSolver<FilterCallable<CRhs>>>(&mut self, problem: &OdeSolverProblem<CRhs, CMass, CInit>, t: CRhs::T, root_solver: &mut RS) -> Result<CRhs::V> {
+        let problem = problem.clone();
         let mut state = OdeSolverState::new_consistent(&problem, root_solver)?;
         self.set_problem(&mut state, problem);
         while state.t <= t {
@@ -43,7 +45,7 @@ pub struct OdeSolverState<M: Matrix> {
 }
 
 impl <M: Matrix> OdeSolverState<M> {
-    fn new<CRhs, CMass, CInit>(ode_problem: &OdeSolverProblem<CRhs, CMass, CInit>) -> Self
+    pub fn new<CRhs, CMass, CInit>(ode_problem: &OdeSolverProblem<CRhs, CMass, CInit>) -> Self
     where
         CRhs: NonLinearOp<M = M, V = M::V, T = M::T>, 
         CMass: LinearOp<M = M, V = M::V, T = M::T>,