conspire/math/tensor/rank_2/logarithm/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::math::Quantity;
5use crate::units::Dimensionless;
6
7use super::{
8    super::{Rank2, Tensor, TensorArray, TensorError, rank_4::TensorRank4},
9    TensorRank2,
10    eigen::{find_orthonormal_eigenvectors, reconstruct_symmetric, solve_cubic_symmetric},
11};
12use crate::math::assert::Assert;
13
14impl<I> TensorRank2<3, I, I, Dimensionless> {
15    /// Returns the matrix logarithm of the 3x3 symmetric tensor.
16    pub fn logm(&self) -> Result<Self, TensorError> {
17        if self.is_diagonal() {
18            let mut logm = TensorRank2::zero();
19            logm.iter_mut()
20                .enumerate()
21                .zip(self.iter())
22                .for_each(|((i, logm_i), self_i)| logm_i[i] = self_i[i].ln());
23            Ok(logm)
24        } else {
25            let tensor = self - &TensorRank2::identity();
26            let norm = tensor.norm();
27            if norm < 1e-2 {
28                let num_terms = if norm < 1e-4 {
29                    2
30                } else if norm < 1e-3 {
31                    3
32                } else {
33                    5
34                };
35                let mut logm = tensor.clone();
36                let mut power = tensor.clone();
37                (2..=num_terms).for_each(|k| {
38                    power *= &tensor;
39                    logm += &power * (if k % 2 == 0 { -1.0 } else { 1.0 } / k as f64);
40                });
41                Ok(logm)
42            } else if self.is_symmetric() {
43                let mut eigenvalues = solve_cubic_symmetric(self.invariants())?;
44                if eigenvalues.iter().any(|eigenvalue| eigenvalue <= &0.0) {
45                    panic!("Symmetric matrix has a non-positive eigenvalue")
46                }
47                let eigenvectors = find_orthonormal_eigenvectors(&eigenvalues, self);
48                eigenvalues
49                    .iter_mut()
50                    .for_each(|eigenvalue| *eigenvalue = eigenvalue.ln());
51                Ok(reconstruct_symmetric(eigenvalues, eigenvectors))
52            } else {
53                panic!("Matrix logarithm only implemented for symmetric cases")
54            }
55        }
56    }
57    /// Returns the derivative of the matrix logarithm of the 3x3 symmetric tensor.
58    pub fn dlogm(&self) -> Result<TensorRank4<3, I, I, I, I, Dimensionless>, TensorError> {
59        if self.is_diagonal() {
60            let mut dlogm = TensorRank4::zero();
61            dlogm.iter_mut().enumerate().for_each(|(i, dlogm_i)| {
62                dlogm_i.iter_mut().enumerate().for_each(|(j, dlogm_ij)| {
63                    dlogm_ij.iter_mut().enumerate().for_each(|(k, dlogm_ijk)| {
64                        dlogm_ijk
65                            .iter_mut()
66                            .enumerate()
67                            .filter(|(l, _)| i == k && &j == l)
68                            .for_each(|(_, dlogm_ijkl)| {
69                                *dlogm_ijkl = if Assert::default()
70                                    .eq_within_tols(self[i][i], &self[j][j])
71                                    .is_ok()
72                                {
73                                    1.0 / self[j][j]
74                                } else {
75                                    (self[i][i].ln() - self[j][j].ln()) / (self[i][i] - self[j][j])
76                                }
77                            })
78                    })
79                })
80            });
81            Ok(dlogm)
82        } else if self.is_symmetric() {
83            let eigenvalues = solve_cubic_symmetric(self.invariants())?;
84            if eigenvalues.iter().any(|eigenvalue| eigenvalue <= &0.0) {
85                panic!("Symmetric matrix has a non-positive eigenvalue")
86            }
87            let divided_difference: Self = eigenvalues
88                .iter()
89                .map(|eigenvalue_i| {
90                    eigenvalues
91                        .iter()
92                        .map(|eigenvalue_j| {
93                            if Assert::default()
94                                .eq_within_tols(eigenvalue_i, eigenvalue_j)
95                                .is_ok()
96                            {
97                                1.0 / eigenvalue_j
98                            } else {
99                                (eigenvalue_i.ln() - eigenvalue_j.ln())
100                                    / (eigenvalue_i - eigenvalue_j)
101                            }
102                        })
103                        .collect()
104                })
105                .collect();
106            let eigenvectors = find_orthonormal_eigenvectors(&eigenvalues, self).transpose();
107            Ok(eigenvectors.iter().map(|eigenvector_i|
108                eigenvectors.iter().map(|eigenvector_j|
109                    eigenvectors.iter().map(|eigenvector_k|
110                        eigenvectors.iter().map(|eigenvector_l|
111                            eigenvector_i.iter().zip(eigenvector_k.iter().zip(divided_difference.iter())).map(|(eigenvector_ip, (eigenvector_kp, divided_difference_p))|
112                                eigenvector_j.iter().zip(eigenvector_l.iter().zip(divided_difference_p.iter())).map(|(eigenvector_jq, (eigenvector_lq, divided_difference_pq))|
113                                    eigenvector_ip * eigenvector_kp * divided_difference_pq * eigenvector_jq * eigenvector_lq
114                                ).sum::<Quantity>()
115                            ).sum::<Quantity>()
116                        ).collect()
117                    ).collect()
118                ).collect()
119            ).collect())
120        } else {
121            panic!("Matrix logarithm only implemented for symmetric cases")
122        }
123    }
124}
conspire/math/tensor/rank_2/logarithm/mod.rs

conspire/math/tensor/rank_2/logarithm/
mod.rs
