numra-ode 0.1.4

ODE and DAE solvers for Numra: DoPri5, Tsit5, Verner 6/7/8, Radau5, ESDIRK 3/4/5, BDF, plus forward sensitivity analysis.
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

//!
//! Author: Moussa Leblouba
//! Date: 26 April 2026
//! Modified: 2 May 2026

use numra_ode::events::{EventAction, EventDirection, EventFunction};
use numra_ode::{DoPri5, OdeProblem, Solver, SolverOptions};

struct ZeroCrossing {
    direction: EventDirection,
    action: EventAction,
}

impl EventFunction<f64> for ZeroCrossing {
    fn evaluate(&self, _t: f64, y: &[f64]) -> f64 {
        y[0]
    }

    fn direction(&self) -> EventDirection {
        self.direction
    }

    fn action(&self) -> EventAction {
        self.action
    }
}

#[test]
fn harmonic_zero_crossings_are_recorded_without_terminating() {
    let problem = OdeProblem::new(
        |_t: f64, y: &[f64], dydt: &mut [f64]| {
            dydt[0] = y[1];
            dydt[1] = -y[0];
        },
        0.0,
        10.0,
        vec![1.0, 0.0],
    );
    let options = SolverOptions::default()
        .rtol(1e-9)
        .atol(1e-11)
        .event(Box::new(ZeroCrossing {
            direction: EventDirection::Both,
            action: EventAction::Continue,
        }));

    let result = DoPri5::solve(&problem, 0.0, 10.0, &[1.0, 0.0], &options).unwrap();

    assert!(!result.terminated_by_event);
    assert!(result.events.len() >= 3, "events={:?}", result.events);
    for (event, expected) in result.events.iter().zip([
        std::f64::consts::FRAC_PI_2,
        3.0 * std::f64::consts::FRAC_PI_2,
        5.0 * std::f64::consts::FRAC_PI_2,
    ]) {
        assert!(
            (event.t - expected).abs() < 2e-2,
            "expected crossing near {expected}, got {}",
            event.t
        );
        assert!(event.y[0].abs() < 2e-4, "event state={:?}", event.y);
    }
}

#[test]
fn terminal_event_stops_at_first_analytic_crossing() {
    let problem = OdeProblem::new(
        |_t: f64, y: &[f64], dydt: &mut [f64]| {
            dydt[0] = y[1];
            dydt[1] = -y[0];
        },
        0.0,
        10.0,
        vec![1.0, 0.0],
    );
    let options = SolverOptions::default()
        .rtol(1e-9)
        .atol(1e-11)
        .event(Box::new(ZeroCrossing {
            direction: EventDirection::Falling,
            action: EventAction::Stop,
        }));

    let result = DoPri5::solve(&problem, 0.0, 10.0, &[1.0, 0.0], &options).unwrap();

    assert!(result.terminated_by_event);
    assert_eq!(result.events.len(), 1);
    let event = &result.events[0];
    assert!((event.t - std::f64::consts::FRAC_PI_2).abs() < 2e-2);
    assert!(event.y[0].abs() < 2e-4);
    assert!(
        event.y[1] < 0.0,
        "falling crossing should have negative velocity"
    );
}

#[test]
fn near_initial_crossing_is_located_with_small_step() {
    let problem = OdeProblem::new(
        |_t: f64, _y: &[f64], dydt: &mut [f64]| {
            dydt[0] = 1.0;
        },
        0.0,
        1.0,
        vec![-1e-6],
    );
    let options = SolverOptions::default()
        .h0(1e-4)
        .h_max(1e-4)
        .rtol(1e-9)
        .atol(1e-12)
        .event(Box::new(ZeroCrossing {
            direction: EventDirection::Rising,
            action: EventAction::Stop,
        }));

    let result = DoPri5::solve(&problem, 0.0, 1.0, &[-1e-6], &options).unwrap();

    assert!(result.terminated_by_event);
    let event = &result.events[0];
    assert!(
        (event.t - 1e-6).abs() < 1e-8,
        "expected event at 1e-6, got {}",
        event.t
    );
    assert!(event.y[0].abs() < 1e-8);
}