Skip to main content

conspire/math/integrate/field/euler/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{Integrable, RECONSTRUCT_FAILED};
5use crate::math::{
6    Derivative, Differentiable, Quantity, TensorVec,
7    integrate::{IntegrationError, Times},
8};
9use std::ops::Mul;
10
11/// Explicit Euler for a single [`Integrable`], one step per interval of `time`.
12///
13/// ```math
14/// \mathbf{x}_{n+1} = \mathrm{reconstruct}\!\left(\mathbf{x}_n,\ h\,\mathbf{f}(t_n, \mathbf{x}_n)\right)
15/// ```
16pub fn integrate_euler<Field, U, T>(
17    mut rate: impl FnMut(Quantity<T>, &Field::Point) -> Result<Derivative<Field::Increment, T>, String>,
18    time: &[Quantity<T>],
19    initial_condition: Field::Point,
20) -> Result<(Times<T>, U), IntegrationError>
21where
22    Field: Integrable,
23    Field::Point: Clone,
24    Field::Increment: Differentiable<T>,
25    for<'a> &'a Derivative<Field::Increment, T>: Mul<Quantity<T>, Output = Field::Increment>,
26    U: TensorVec<Item = Field::Point>,
27{
28    let mut point = initial_condition;
29    let mut points = U::new();
30    let mut times = Times::new();
31    points.push(point.clone());
32    times.push(time[0]);
33    for step in time.windows(2) {
34        let increment = &rate(step[0], &point)? * (step[1] - step[0]);
35        point = Field::reconstruct(&point, &increment)
36            .map_err(|_| IntegrationError::from(RECONSTRUCT_FAILED.to_string()))?;
37        points.push(point.clone());
38        times.push(step[1]);
39    }
40    Ok((times, points))
41}