num_dual/datatypes/
dual2.rs

1use crate::{DualNum, DualNumFloat, DualStruct};
2use num_traits::{FloatConst, FromPrimitive, Inv, Num, One, Signed, Zero};
3#[cfg(feature = "serde")]
4use serde::{Deserialize, Serialize};
5use std::fmt;
6use std::iter::{Product, Sum};
7use std::ops::{
8    Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign,
9};
10
11/// A scalar second order dual number for the calculation of second derivatives.
12#[derive(Copy, Clone, Debug)]
13#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
14pub struct Dual2<T> {
15    /// Real part of the second order dual number
16    pub re: T,
17    /// First derivative part of the second order dual number
18    pub v1: T,
19    /// Second derivative part of the second order dual number
20    pub v2: T,
21}
22
23#[cfg(feature = "ndarray")]
24impl<T: DualNum> ndarray::ScalarOperand for Dual2<T> {}
25
26pub type Dual2_32 = Dual2<f32>;
27pub type Dual2_64 = Dual2<f64>;
28
29impl<T> Dual2<T> {
30    /// Create a new second order dual number from its fields.
31    #[inline]
32    pub fn new(re: T, v1: T, v2: T) -> Self {
33        Self { re, v1, v2 }
34    }
35}
36
37impl<T: One> Dual2<T> {
38    /// Set the derivative part to 1.
39    /// ```
40    /// # use num_dual::{Dual2, DualNum};
41    /// let x = Dual2::from_re(5.0).derivative().powi(2);
42    /// assert_eq!(x.re, 25.0);             // x²
43    /// assert_eq!(x.v1, 10.0);    // 2x
44    /// assert_eq!(x.v2, 2.0);     // 2
45    /// ```
46    ///
47    /// Can also be used for higher order derivatives.
48    /// ```
49    /// # use num_dual::{Dual64, Dual2, DualNum};
50    /// let x = Dual2::from_re(Dual64::from_re(5.0).derivative())
51    ///     .derivative()
52    ///     .powi(2);
53    /// assert_eq!(x.re.re, 25.0);      // x²
54    /// assert_eq!(x.re.eps, 10.0);     // 2x
55    /// assert_eq!(x.v1.re, 10.0);      // 2x
56    /// assert_eq!(x.v1.eps, 2.0);      // 2
57    /// assert_eq!(x.v2.re, 2.0);       // 2
58    /// ```
59    #[inline]
60    pub fn derivative(mut self) -> Self {
61        self.v1 = T::one();
62        self
63    }
64}
65
66impl<T: Zero> Dual2<T> {
67    /// Create a new second order dual number from the real part.
68    #[inline]
69    pub fn from_re(re: T) -> Self {
70        Self::new(re, T::zero(), T::zero())
71    }
72}
73
74/* chain rule */
75impl<T: DualNum> Dual2<T> {
76    #[inline]
77    fn chain_rule(&self, f0: T, f1: T, f2: T) -> Self {
78        Self::new(
79            f0,
80            self.v1.clone() * f1.clone(),
81            self.v2.clone() * f1 + self.v1.clone() * self.v1.clone() * f2,
82        )
83    }
84}
85
86/* product rule */
87impl<T: DualNum> Mul<&Dual2<T>> for &Dual2<T> {
88    type Output = Dual2<T>;
89    #[inline]
90    fn mul(self, other: &Dual2<T>) -> Dual2<T> {
91        Dual2::new(
92            self.re.clone() * other.re.clone(),
93            other.v1.clone() * self.re.clone() + self.v1.clone() * other.re.clone(),
94            other.v2.clone() * self.re.clone()
95                + self.v1.clone() * other.v1.clone()
96                + other.v1.clone() * self.v1.clone()
97                + self.v2.clone() * other.re.clone(),
98        )
99    }
100}
101
102/* quotient rule */
103impl<T: DualNum> Div<&Dual2<T>> for &Dual2<T> {
104    type Output = Dual2<T>;
105    #[inline]
106    fn div(self, other: &Dual2<T>) -> Dual2<T> {
107        let inv = other.re.recip();
108        let inv2 = inv.clone() * inv.clone();
109        Dual2::new(
110            self.re.clone() * inv.clone(),
111            (self.v1.clone() * other.re.clone() - other.v1.clone() * self.re.clone())
112                * inv2.clone(),
113            self.v2.clone() * inv.clone()
114                - (other.v2.clone() * self.re.clone()
115                    + self.v1.clone() * other.v1.clone()
116                    + other.v1.clone() * self.v1.clone())
117                    * inv2.clone()
118                + other.v1.clone()
119                    * other.v1.clone()
120                    * ((T::one() + T::one()) * self.re.clone() * inv2 * inv),
121        )
122    }
123}
124
125/* string conversions */
126impl<T: DualNum> fmt::Display for Dual2<T> {
127    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
128        write!(f, "{} + {}ε1 + {}ε1²", self.re, self.v1, self.v2)
129    }
130}
131
132impl_second_derivatives!(Dual2, [v1, v2]);
133impl_dual!(Dual2, [v1, v2]);
134#[cfg(feature = "nalgebra")]
135impl_nalgebra!(Dual2, [v1, v2]);
num_dual/datatypes/dual2.rs

num_dual/datatypes/
dual2.rs
