diffsol 0.1.9

A library for solving ordinary differential equations (ODEs) in Rust.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

use anyhow::Result;

use crate::{op::Op, solver::SolverProblem};

pub struct NonLinearSolveSolution<V> {
    pub x0: V,
    pub x: V,
}

impl<V> NonLinearSolveSolution<V> {
    pub fn new(x0: V, x: V) -> Self {
        Self { x0, x }
    }
}

/// A solver for the nonlinear problem `F(x) = 0`.
pub trait NonLinearSolver<C: Op> {
    /// Get the problem to be solved.
    fn problem(&self) -> &SolverProblem<C>;

    /// Set the problem to be solved, any previous problem is discarded.
    fn set_problem(&mut self, problem: &SolverProblem<C>);

    /// Reset the approximation of the Jacobian matrix.
    fn reset_jacobian(&mut self, x: &C::V, t: C::T);

    // Solve the problem `F(x, t) = 0` for fixed t, and return the solution `x`.
    fn solve(&mut self, x: &C::V, t: C::T) -> Result<C::V> {
        let mut x = x.clone();
        self.solve_in_place(&mut x, t)?;
        Ok(x)
    }

    /// Solve the problem `F(x) = 0` in place.
    fn solve_in_place(&mut self, x: &mut C::V, t: C::T) -> Result<()>;

    /// Solve the linearised problem `J * x = b`, where `J` was calculated using [Self::reset_jacobian].
    /// The input `b` is provided in `x`, and the solution is returned in `x`.
    fn solve_linearised_in_place(&self, x: &mut C::V) -> Result<()>;

    // Set the maximum number of iterations for the solver.
    fn set_max_iter(&mut self, max_iter: usize);

    // Get the maximum number of iterations for the solver.
    fn max_iter(&self) -> usize;

    // Get the number of iterations taken by the solver on the last call to `solve`.
    fn niter(&self) -> usize;
}

pub mod convergence;
pub mod newton;
pub mod root;

//tests
#[cfg(test)]
pub mod tests {
    use std::rc::Rc;

    use self::newton::NewtonNonlinearSolver;
    use crate::{
        linear_solver::nalgebra::lu::LU,
        matrix::MatrixCommon,
        op::{closure::Closure, NonLinearOp},
        scale, DenseMatrix, Vector,
    };

    use super::*;
    use num_traits::{One, Zero};

    pub fn get_square_problem<M>() -> (
        SolverProblem<impl NonLinearOp<M = M, V = M::V, T = M::T>>,
        Vec<NonLinearSolveSolution<M::V>>,
    )
    where
        M: DenseMatrix + 'static,
    {
        let jac1 = M::from_diagonal(&M::V::from_vec(vec![2.0.into(), 2.0.into()]));
        let jac2 = jac1.clone();
        let p = Rc::new(M::V::zeros(0));
        let op = Closure::new(
            // 0 = J * x * x - 8
            move |x: &<M as MatrixCommon>::V, _p: &<M as MatrixCommon>::V, _t, y| {
                jac1.gemv(M::T::one(), x, M::T::zero(), y); // y = J * x
                y.component_mul_assign(x);
                y.add_scalar_mut(M::T::from(-8.0));
            },
            // J = 2 * J * x * dx
            move |x: &<M as MatrixCommon>::V, _p: &<M as MatrixCommon>::V, _t, v, y| {
                jac2.gemv(M::T::from(2.0), x, M::T::zero(), y); // y = 2 * J * x
                y.component_mul_assign(v);
            },
            2,
            2,
            p,
        );
        let rtol = M::T::from(1e-6);
        let atol = M::V::from_vec(vec![1e-6.into(), 1e-6.into()]);
        let problem = SolverProblem::new(Rc::new(op), Rc::new(atol), rtol);
        let solns = vec![NonLinearSolveSolution::new(
            M::V::from_vec(vec![2.1.into(), 2.1.into()]),
            M::V::from_vec(vec![2.0.into(), 2.0.into()]),
        )];
        (problem, solns)
    }

    pub fn test_nonlinear_solver<C>(
        mut solver: impl NonLinearSolver<C>,
        problem: SolverProblem<C>,
        solns: Vec<NonLinearSolveSolution<C::V>>,
    ) where
        C: NonLinearOp,
    {
        solver.set_problem(&problem);
        let t = C::T::zero();
        for soln in solns {
            let x = solver.solve(&soln.x0, t).unwrap();
            let tol = x.clone() * scale(problem.rtol) + problem.atol.as_ref();
            x.assert_eq(&soln.x, &tol);
        }
    }

    type MCpu = nalgebra::DMatrix<f64>;

    #[test]
    fn test_newton_cpu_square() {
        let lu = LU::default();
        let (prob, soln) = get_square_problem::<MCpu>();
        let s = NewtonNonlinearSolver::new(lu);
        test_nonlinear_solver(s, prob, soln);
    }
}