diffurch 0.0.3

Numerical solver for ordinary and delay differential equations
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

#![allow(incomplete_features)]
#![feature(generic_const_exprs)]

use diffurch::*;

#[test]
fn constant() {
    let eq = equation!(|| [1., -1., 0.]);
    let ic = [0., 0., 0.];
    let solution = |t: f64| [t, -t, 0.];

    Solver::new()
        .stepsize(0.25)
        .on_step(event!(|t, [x, y, z]| assert_eq!([x, y, z], solution(t))))
        .run(eq, ic, 0. ..10.);
}

#[test]
fn time() {
    let eq = equation!(|t| [0., 1., t]);
    let ic = [0., 0., 0.];
    let solution = |t: f64| [0., t, t * t / 2.];

    Solver::new()
        .stepsize(0.5)
        .on_step(event!(|t, [x, y, z]| assert_eq!([x, y, z], solution(t))))
        .run(eq, ic, 0. ..10.);
}

#[test]
fn ode_exponent() {
    let eq = equation!(|[x]| [-x]);
    let ic = [1.];
    let solution = |t: f64| (-t).exp();

    Solver::new()
        .stepsize(0.1)
        .on_step(event!(|t, [x]| {
            assert!((x - solution(t)).abs() < 1e-14)
        }))
        .run(eq, ic, 0. ..10.);
}

#[test]
fn ode_harmonic() {
    let eq = equation!(|[x, dx]| [dx, -x]);
    let ic = [0., 1.];
    let solution = |t: f64| t.sin();

    Solver::new()
        .stepsize(0.1)
        .on_step(event!(|t, [x, _dx]| {
            assert!((x - solution(t)).abs() < 1e-13)
        }))
        .run(eq, ic, 0. ..10.);
}

#[test]
fn ode2_lin() {
    let eq = equation!(|t, [x]| [2. * x / t]);
    let ic = |t: f64| [t * t];
    let sol = |t: f64| t * t;

    Solver::new()
        .stepsize(0.1)
        .on_step(event!(|t, [x]| assert!((x - sol(t)).abs() < 1e-11)))
        .run(eq, ic, 1. ..10.);
}

#[test]
fn dde_sin() {
    let k = 1.;
    let tau: f64 = 1.;

    let a = k / (k * tau).tan();
    let b = -k / (k * tau).sin();

    let eq = equation!(|t, [x], [x_]| [a * x + b * x_(t - tau)]).const_neutral_delay(tau);
    let ic = |t: f64| [(k * t).sin()];
    let sol = |t: f64| (k * t).sin();

    Solver::new()
        .stepsize(0.33)
        .on_step(event!(|t, [x]| assert!((x - sol(t)).abs() < 1e-11)))
        .run(eq, ic, ..10.);
}

#[test]
fn ndde_sin() {
    let k: f64 = 1.;
    let tau: f64 = 1.;

    let a = -k * (k * tau).sin();
    let b = (k * tau).cos();

    let eq = equation!(|t, [_x], [x]| [a * x(t - tau) + b * x.d(t - tau)]);
    let ic = (|t: f64| [(k * t).sin()], |t: f64| [k * (k * t).cos()]);
    let sol = |t: f64| (k * t).sin();

    Solver::new()
        .stepsize(0.25)
        .on_step(event!(|t, [x]| assert!((x - sol(t)).abs() < 1e-10)))
        .run(eq, ic, ..10.);
}

#[test]
fn bouncing_ball() {
    let eq = equation!(|[_x, dx]| [dx, -2.]).max_delay(f64::INFINITY);

    let ic = [0., 2.];
    let range = 0. ..10.;

    let sol = |t: f64| {
        let t = t % 2.;
        t * (2. - t)
    };

    let when_hit = Loc::new(state_fn!(|[x, _]| x)).sign().bisection();
    let bounce = event_mut!(|[x, dx]| {
        *x = 0.;
        *dx = dx.abs();
    });

    Solver::new()
        .stepsize(0.05)
        .on(when_hit, bounce)
        .on_step(event!(|t, [x, _]| assert!((x - sol(t)).abs() < 1e-13)))
        .run(eq, ic, range);
}