67 lines
2.3 KiB
Rust
67 lines
2.3 KiB
Rust
//! The scalar type an ODE is integrated in.
|
|
//!
|
|
//! All of our problem types: state, time and step size, is generic over
|
|
//! [`Scalar`]. Instantiating a problem at `f64` solves it. Instantiating the
|
|
//! same problem at a forward-mode dual number solves it *and* propagates exact
|
|
//! derivatives alongside the state, at machine precision, with no finite
|
|
//! differencing and no extra right-hand-side evaluations beyond the ones the
|
|
//! dual arithmetic performs anyway.
|
|
//!
|
|
//! [`Scalar`] extends the vocabulary [`num_dual`] already defines rather than
|
|
//! restating it: [`DualNum`] carries the arithmetic, the transcendentals and
|
|
//! `From<f64>`, and [`DualStruct`] carries [`re`](DualStruct::re). Both are
|
|
//! implemented for `f64` as well as for every dual type, so one blanket impl
|
|
//! covers all of them, and a crate depending on this one alongside another
|
|
//! built on `num-dual` can satisfy both with a single bound.
|
|
|
|
use nalgebra::RealField;
|
|
use num_dual::{DualNum, DualStruct};
|
|
use num_traits::NumOps;
|
|
|
|
/// A scalar that an ODE's state, time and step size may take.
|
|
///
|
|
/// `re`, which discards derivative information, is inherited from
|
|
/// [`DualStruct`]; this trait adds only [`lift`](Scalar::lift).
|
|
pub trait Scalar:
|
|
RealField + Copy + NumOps<f64, Self> + DualNum<Primitive = f64> + DualStruct<Real = f64>
|
|
{
|
|
/// A real constant as a value of this type, with zero derivative.
|
|
fn lift(x: f64) -> Self {
|
|
Self::from(x)
|
|
}
|
|
}
|
|
|
|
impl<T> Scalar for T where
|
|
T: RealField + Copy + NumOps<f64, T> + DualNum<Primitive = f64> + DualStruct<Real = f64>
|
|
{
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use super::*;
|
|
|
|
#[test]
|
|
fn test_f64_is_a_scalar() {
|
|
assert_eq!(f64::lift(3.5), 3.5);
|
|
assert_eq!(3.5f64.re(), 3.5);
|
|
}
|
|
|
|
#[test]
|
|
fn test_dual_projects_to_its_real_part() {
|
|
use num_dual::Dual64;
|
|
|
|
// A lifted constant has no derivative; a seeded variable does, and both
|
|
// project back to the same real value.
|
|
let c = <Dual64 as Scalar>::lift(3.5);
|
|
assert_eq!(c.re, 3.5);
|
|
assert_eq!(c.eps, 0.0);
|
|
|
|
let x = Dual64::from_re(3.5).derivative();
|
|
assert_eq!(x.re(), 3.5);
|
|
assert_eq!(x.eps, 1.0);
|
|
|
|
// d/dx x^3 = 3x^2 = 36.75 at x = 3.5
|
|
assert_eq!(x.powi(3).eps, 3.0 * 3.5 * 3.5);
|
|
}
|
|
}
|