diffsol 0.18.0

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

use crate::{Matrix, OdeBuilder, OdeEquationsImplicit, OdeSolverProblem};
use num_traits::{FromPrimitive, One, Pow};

// du/dt = -u + v - 1
// 0 = u - v - v^3
fn nonlinear_algebraic_rhs<M: Matrix>(x: &[M::T], _p: &[M::T], _t: M::T, y: &mut [M::T]) {
    y[0] = -x[0] + x[1] - M::T::one();
    y[1] = x[0] - x[1] - x[1] * x[1] * x[1];
}

fn nonlinear_algebraic_jac_mul<M: Matrix>(
    x: &[M::T],
    _p: &[M::T],
    _t: M::T,
    v: &[M::T],
    y: &mut [M::T],
) {
    let three = M::T::from_f64(3.0).unwrap();
    y[0] = -v[0] + v[1];
    y[1] = v[0] - (M::T::one() + three * x[1] * x[1]) * v[1];
}

fn nonlinear_algebraic_mass<M: Matrix>(
    x: &[M::T],
    _p: &[M::T],
    _t: M::T,
    beta: M::T,
    y: &mut [M::T],
) {
    y[0] = x[0] + beta * y[0];
    y[1] *= beta;
}

fn nonlinear_algebraic_init<M: Matrix>(_p: &[M::T], _t: M::T, y: &mut [M::T]) {
    y[0] = M::T::from_f64(2.0).unwrap();
    y[1] = M::T::from_f64(5.0).unwrap();
}

pub fn nonlinear_algebraic_problem<M: Matrix + 'static>(
) -> OdeSolverProblem<impl OdeEquationsImplicit<M = M, V = M::V, T = M::T, C = M::C>> {
    OdeBuilder::<M>::new()
        .rhs_implicit(
            nonlinear_algebraic_rhs::<M>,
            nonlinear_algebraic_jac_mul::<M>,
        )
        .mass(nonlinear_algebraic_mass::<M>)
        .init(nonlinear_algebraic_init::<M>, 2)
        .build()
        .unwrap()
}

// as above, with mass M(t) = diag(t, 0)
fn nonlinear_algebraic_time_dependent_mass<M: Matrix>(
    x: &[M::T],
    _p: &[M::T],
    t: M::T,
    beta: M::T,
    y: &mut [M::T],
) {
    y[0] = t * x[0] + beta * y[0];
    y[1] *= beta;
}

pub fn nonlinear_algebraic_time_dependent_mass_problem<M: Matrix + 'static>(
) -> OdeSolverProblem<impl OdeEquationsImplicit<M = M, V = M::V, T = M::T, C = M::C>> {
    OdeBuilder::<M>::new()
        .rhs_implicit(
            nonlinear_algebraic_rhs::<M>,
            nonlinear_algebraic_jac_mul::<M>,
        )
        .mass(nonlinear_algebraic_time_dependent_mass::<M>)
        .init(nonlinear_algebraic_init::<M>, 2)
        .build()
        .unwrap()
}

// du/dt = -u
// 0 = u - v^10
fn power_algebraic_rhs<M: Matrix>(x: &[M::T], _p: &[M::T], _t: M::T, y: &mut [M::T]) {
    y[0] = -x[0];
    y[1] = x[0] - x[1].pow(10i32);
}

fn power_algebraic_jac_mul<M: Matrix>(
    x: &[M::T],
    _p: &[M::T],
    _t: M::T,
    v: &[M::T],
    y: &mut [M::T],
) {
    let ten = M::T::from_f64(10.0).unwrap();
    y[0] = -v[0];
    y[1] = v[0] - ten * x[1].pow(9i32) * v[1];
}

fn power_algebraic_init<M: Matrix>(_p: &[M::T], _t: M::T, y: &mut [M::T]) {
    y[0] = M::T::one();
    y[1] = M::T::from_f64(0.45).unwrap();
}

pub fn power_algebraic_problem<M: Matrix + 'static>(
) -> OdeSolverProblem<impl OdeEquationsImplicit<M = M, V = M::V, T = M::T, C = M::C>> {
    OdeBuilder::<M>::new()
        .rhs_implicit(power_algebraic_rhs::<M>, power_algebraic_jac_mul::<M>)
        .mass(nonlinear_algebraic_mass::<M>)
        .init(power_algebraic_init::<M>, 2)
        .build()
        .unwrap()
}