differential_equations/
bvp.rs

1//! Boundary value problem builder API.
2
3use crate::{
4    error::Error,
5    methods::bvp::BVPMethod,
6    ode::{ODE, solve_bvp},
7    solout::{
8        CrossingDirection, CrossingSolout, DefaultSolout, DenseSolout, EvenSolout, Event,
9        EventWrappedSolout, HyperplaneCrossingSolout, Solout, TEvalSolout,
10    },
11    solution::Solution,
12    traits::{DefaultState, Real, State},
13};
14
15/// Builder for boundary value problems.
16///
17/// Use [`BVP::ode`] for an ODE system that implements [`ODE`] and [`Boundary`],
18/// or [`BVP::ode_from_fn`] when closures are more convenient.
19#[derive(Clone, Debug)]
20pub struct BVP<EqType, T: Real, Y: State<T>, Method, SoloutType> {
21    equation: EqType,
22    t0: T,
23    tf: T,
24    y_guess: Y,
25    method: Method,
26    solout: SoloutType,
27}
28
29/// Boundary condition residual for a boundary value problem.
30///
31/// The solver searches for an initial state such that this residual is zero.
32pub trait Boundary<T = f64, Y = DefaultState<T>>
33where
34    T: Real,
35    Y: State<T>,
36{
37    /// Compute the boundary residual `res = g(y_a, y_b)`.
38    fn boundary(&self, y_a: &Y, y_b: &Y, res: &mut Y);
39}
40
41impl<EqType, T: Real, Y: State<T>> Boundary<T, Y> for &EqType
42where
43    EqType: Boundary<T, Y> + ?Sized,
44{
45    fn boundary(&self, y_a: &Y, y_b: &Y, res: &mut Y) {
46        (*self).boundary(y_a, y_b, res);
47    }
48}
49
50/// Closure-backed ODE boundary value problem.
51#[derive(Clone, Debug)]
52pub struct OdeBvpFromFn<F, B> {
53    diff: F,
54    boundary: B,
55}
56
57impl<F, B> OdeBvpFromFn<F, B> {
58    /// Create a closure-backed ODE BVP definition.
59    pub fn new(diff: F, boundary: B) -> Self {
60        Self { diff, boundary }
61    }
62}
63
64impl<F, B, T, Y> ODE<T, Y> for OdeBvpFromFn<F, B>
65where
66    T: Real,
67    Y: State<T>,
68    F: Fn(T, &Y, &mut Y),
69{
70    fn diff(&self, t: T, y: &Y, dydt: &mut Y) {
71        (self.diff)(t, y, dydt);
72    }
73}
74
75impl<F, B, T, Y> Boundary<T, Y> for OdeBvpFromFn<F, B>
76where
77    T: Real,
78    Y: State<T>,
79    F: Fn(T, &Y, &mut Y),
80    B: Fn(&Y, &Y, &mut Y),
81{
82    fn boundary(&self, y_a: &Y, y_b: &Y, res: &mut Y) {
83        (self.boundary)(y_a, y_b, res);
84    }
85}
86
87impl<'a, F, T, Y> BVP<&'a F, T, Y, (), DefaultSolout>
88where
89    T: Real,
90    Y: State<T>,
91    F: ODE<T, Y> + Boundary<T, Y> + ?Sized,
92{
93    /// Create a boundary value problem for an ODE.
94    pub fn ode(system: &'a F, t0: T, tf: T, y_guess: Y) -> Self {
95        Self {
96            equation: system,
97            t0,
98            tf,
99            y_guess,
100            method: (),
101            solout: DefaultSolout::new(),
102        }
103    }
104}
105
106impl<F, B, T, Y> BVP<OdeBvpFromFn<F, B>, T, Y, (), DefaultSolout>
107where
108    T: Real,
109    Y: State<T>,
110    F: Fn(T, &Y, &mut Y),
111    B: Fn(&Y, &Y, &mut Y),
112{
113    /// Create a boundary value problem for an ODE from closures.
114    pub fn ode_from_fn(diff: F, boundary: B, t0: T, tf: T, y_guess: Y) -> Self {
115        Self {
116            equation: OdeBvpFromFn::new(diff, boundary),
117            t0,
118            tf,
119            y_guess,
120            method: (),
121            solout: DefaultSolout::new(),
122        }
123    }
124}
125
126impl<EqType, T: Real, Y: State<T>, Method, SoloutType> BVP<EqType, T, Y, Method, SoloutType> {
127    fn with_method<NextMethod>(
128        self,
129        method: NextMethod,
130    ) -> BVP<EqType, T, Y, NextMethod, SoloutType> {
131        BVP {
132            equation: self.equation,
133            t0: self.t0,
134            tf: self.tf,
135            y_guess: self.y_guess,
136            method,
137            solout: self.solout,
138        }
139    }
140
141    fn with_solout<NextSolout>(self, solout: NextSolout) -> BVP<EqType, T, Y, Method, NextSolout> {
142        BVP {
143            equation: self.equation,
144            t0: self.t0,
145            tf: self.tf,
146            y_guess: self.y_guess,
147            method: self.method,
148            solout,
149        }
150    }
151
152    /// Set the numerical method used to solve the BVP.
153    pub fn method<NextMethod>(
154        self,
155        method: NextMethod,
156    ) -> BVP<EqType, T, Y, NextMethod, SoloutType> {
157        self.with_method(method)
158    }
159
160    /// Set a custom solout function for the final converged trajectory.
161    pub fn solout<NextSolout>(self, solout: NextSolout) -> BVP<EqType, T, Y, Method, NextSolout> {
162        self.with_solout(solout)
163    }
164
165    /// Output evenly spaced points between the initial and final time.
166    ///
167    /// This controls the output of the final converged IVP trajectory.
168    pub fn even(self, dt: T) -> BVP<EqType, T, Y, Method, EvenSolout<T>> {
169        let solout = EvenSolout::new(dt, self.t0, self.tf);
170        self.with_solout(solout)
171    }
172
173    /// Output dense interpolation points between accepted solver steps.
174    ///
175    /// This controls the output of the final converged IVP trajectory.
176    pub fn dense(self, n: usize) -> BVP<EqType, T, Y, Method, DenseSolout> {
177        self.with_solout(DenseSolout::new(n))
178    }
179
180    /// Use provided time points for final trajectory output.
181    pub fn t_eval(self, points: impl AsRef<[T]>) -> BVP<EqType, T, Y, Method, TEvalSolout<T>> {
182        let solout = TEvalSolout::new(points, self.t0, self.tf);
183        self.with_solout(solout)
184    }
185
186    /// Wrap current final-trajectory output with event detection.
187    pub fn event<'a, E>(
188        self,
189        event: &'a E,
190    ) -> BVP<EqType, T, Y, Method, EventWrappedSolout<'a, T, Y, SoloutType, E>>
191    where
192        E: Event<T, Y>,
193        SoloutType: Solout<T, Y>,
194    {
195        BVP {
196            equation: self.equation,
197            t0: self.t0,
198            tf: self.tf,
199            y_guess: self.y_guess,
200            method: self.method,
201            solout: EventWrappedSolout::new(self.solout, event, self.t0, self.tf),
202        }
203    }
204
205    /// Output points where a component crosses a threshold on the final trajectory.
206    pub fn crossing(
207        self,
208        component_idx: usize,
209        threshold: T,
210        direction: CrossingDirection,
211    ) -> BVP<EqType, T, Y, Method, CrossingSolout<T>> {
212        let crossing_solout =
213            CrossingSolout::new(component_idx, threshold).with_direction(direction);
214        self.with_solout(crossing_solout)
215    }
216
217    /// Output points where a projected final trajectory crosses a hyperplane.
218    pub fn hyperplane_crossing<Y1: State<T>>(
219        self,
220        point: Y1,
221        normal: Y1,
222        extractor: fn(&Y) -> Y1,
223        direction: CrossingDirection,
224    ) -> BVP<EqType, T, Y, Method, HyperplaneCrossingSolout<T, Y1, Y>> {
225        let solout =
226            HyperplaneCrossingSolout::new(point, normal, extractor).with_direction(direction);
227        self.with_solout(solout)
228    }
229}
230
231impl<EqType, T, Y, Method, SoloutType> BVP<EqType, T, Y, Method, SoloutType>
232where
233    T: Real,
234    Y: State<T>,
235    EqType: ODE<T, Y> + Boundary<T, Y>,
236    Method: BVPMethod<T, Y>,
237    SoloutType: Solout<T, Y>,
238{
239    /// Solve the boundary value problem.
240    pub fn solve(mut self) -> Result<Solution<T, Y>, Error<T, Y>> {
241        solve_bvp(
242            &mut self.method,
243            &self.equation,
244            self.t0,
245            self.tf,
246            &self.y_guess,
247            &mut self.solout,
248        )
249    }
250}
differential_equations/bvp.rs

differential_equations/
bvp.rs
