scirs2_integrate/ode/utils/
interpolation.rs

1//! Interpolation for dense output.
2//!
3//! This module provides interpolation functions for generating
4//! dense output from ODE solvers, allowing evaluation at arbitrary points
5//! within the integration range without recomputing the entire solution.
6
7use crate::IntegrateFloat;
8use scirs2_core::ndarray::Array1;
9use scirs2_core::numeric::Float;
10use std::fmt::Debug;
11
12/// Continuous output method for interpolation
13#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
14pub enum ContinuousOutputMethod {
15    /// Simple linear interpolation between points
16    Linear,
17    /// Cubic Hermite interpolation using function values and derivatives
18    #[default]
19    CubicHermite,
20    /// Specialized high-order interpolation for specific ODE methods
21    MethodSpecific,
22}
23
24/// Find the index in a sorted array where the given value would be inserted
25///
26/// # Arguments
27///
28/// * `sorted_array` - Sorted array
29/// * `value` - Value to look for
30///
31/// # Returns
32///
33/// The index where value would be inserted to keep the array sorted
34#[allow(dead_code)]
35pub fn find_index<F: Float>(sortedarray: &[F], value: F) -> usize {
36    // Binary search
37    let mut left = 0;
38    let mut right = sortedarray.len();
39
40    while left < right {
41        let mid = (left + right) / 2;
42        if sortedarray[mid] < value {
43            left = mid + 1;
44        } else {
45            right = mid;
46        }
47    }
48
49    left
50}
51
52/// Linear interpolation
53///
54/// # Arguments
55///
56/// * `x` - x-coordinates
57/// * `y` - y-values at each x-coordinate
58/// * `x_new` - Point where to interpolate
59///
60/// # Returns
61///
62/// Linearly interpolated value at x_new
63#[allow(dead_code)]
64pub fn linear_interpolation<F: IntegrateFloat>(x: &[F], y: &[Array1<F>], xnew: F) -> Array1<F> {
65    let i = find_index(x, xnew);
66
67    if i == 0 {
68        return y[0].clone();
69    }
70
71    if i >= x.len() {
72        return y[x.len() - 1].clone();
73    }
74
75    let x0 = x[i - 1];
76    let x1 = x[i];
77    let y0 = &y[i - 1];
78    let y1 = &y[i];
79
80    // Linear interpolation: y = y0 + (x - x0) * (y1 - y0) / (x1 - x0)
81    let t = (xnew - x0) / (x1 - x0);
82
83    let mut result = y0.clone();
84    for (r, (a, b)) in result.iter_mut().zip(y0.iter().zip(y1.iter())) {
85        *r = *a + t * (*b - *a);
86    }
87
88    result
89}
90
91/// Cubic Hermite interpolation
92///
93/// # Arguments
94///
95/// * `x` - x-coordinates
96/// * `y` - y-values at each x-coordinate
97/// * `dy` - Derivatives at each x-coordinate
98/// * `x_new` - Point where to interpolate
99///
100/// # Returns
101///
102/// Cubic interpolated value at x_new
103#[allow(dead_code)]
104pub fn cubic_hermite_interpolation<F: IntegrateFloat>(
105    x: &[F],
106    y: &[Array1<F>],
107    dy: &[Array1<F>],
108    x_new: F,
109) -> Array1<F> {
110    let i = find_index(x, x_new);
111
112    if i == 0 {
113        return y[0].clone();
114    }
115
116    if i >= x.len() {
117        return y[x.len() - 1].clone();
118    }
119
120    let x0 = x[i - 1];
121    let x1 = x[i];
122    let y0 = &y[i - 1];
123    let y1 = &y[i];
124    let dy0 = &dy[i - 1];
125    let dy1 = &dy[i];
126
127    // Normalized time coordinate
128    let h = x1 - x0;
129    let t = (x_new - x0) / h;
130
131    // Hermite basis functions
132    let h00 =
133        F::from_f64(2.0).unwrap() * t.powi(3) - F::from_f64(3.0).unwrap() * t.powi(2) + F::one();
134    let h10 = t.powi(3) - F::from_f64(2.0).unwrap() * t.powi(2) + t;
135    let h01 = F::from_f64(-2.0).unwrap() * t.powi(3) + F::from_f64(3.0).unwrap() * t.powi(2);
136    let h11 = t.powi(3) - t.powi(2);
137
138    // Compute interpolant: p(t) = h00*y0 + h10*h*dy0 + h01*y1 + h11*h*dy1
139    let mut result = Array1::zeros(y0.dim());
140
141    for i in 0..y0.len() {
142        result[i] = h00 * y0[i] + h10 * h * dy0[i] + h01 * y1[i] + h11 * h * dy1[i];
143    }
144
145    result
146}
scirs2_integrate/ode/utils/interpolation.rs

scirs2_integrate/ode/utils/
interpolation.rs
