130 lines
3.8 KiB
Rust
130 lines
3.8 KiB
Rust
#[cfg(test)]
|
|
mod tests {
|
|
//https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py
|
|
#[test]
|
|
fn exp_static() {
|
|
use nalgebra::{Matrix1, Matrix2, Matrix3};
|
|
|
|
{
|
|
let m = Matrix1::new(1.0);
|
|
|
|
let f = m.exp();
|
|
|
|
assert!(relative_eq!(f, Matrix1::new(1_f64.exp()), epsilon = 1.0e-7));
|
|
}
|
|
|
|
{
|
|
let m = Matrix2::new(0.0, 1.0, 0.0, 0.0);
|
|
|
|
assert!(relative_eq!(
|
|
m.exp(),
|
|
Matrix2::new(1.0, 1.0, 0.0, 1.0),
|
|
epsilon = 1.0e-7
|
|
));
|
|
}
|
|
|
|
{
|
|
let a: f64 = 1.0;
|
|
let b: f64 = 2.0;
|
|
let c: f64 = 3.0;
|
|
let d: f64 = 4.0;
|
|
let m = Matrix2::new(a, b, c, d);
|
|
|
|
let delta = ((a - d).powf(2.0) + 4.0 * b * c).sqrt();
|
|
let delta_2 = delta / 2.0;
|
|
let ad_2 = (a + d) / 2.0;
|
|
let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh());
|
|
let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh();
|
|
let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh();
|
|
let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
|
|
|
|
let f = Matrix2::new(m11, m12, m21, m22) / delta;
|
|
assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
|
|
}
|
|
|
|
{
|
|
// https://mathworld.wolfram.com/MatrixExponential.html
|
|
use rand::{
|
|
distributions::{Distribution, Uniform},
|
|
thread_rng,
|
|
};
|
|
let mut rng = thread_rng();
|
|
let dist = Uniform::new(-10.0, 10.0);
|
|
loop {
|
|
let a: f64 = dist.sample(&mut rng);
|
|
let b: f64 = dist.sample(&mut rng);
|
|
let c: f64 = dist.sample(&mut rng);
|
|
let d: f64 = dist.sample(&mut rng);
|
|
let m = Matrix2::new(a, b, c, d);
|
|
|
|
let delta_sq = (a - d).powf(2.0) + 4.0 * b * c;
|
|
if delta_sq < 0.0 {
|
|
continue;
|
|
}
|
|
|
|
let delta = delta_sq.sqrt();
|
|
let delta_2 = delta / 2.0;
|
|
let ad_2 = (a + d) / 2.0;
|
|
let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh());
|
|
let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh();
|
|
let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh();
|
|
let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
|
|
|
|
let f = Matrix2::new(m11, m12, m21, m22) / delta;
|
|
println!("a: {}", m);
|
|
assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
|
|
break;
|
|
}
|
|
}
|
|
|
|
{
|
|
let m = Matrix3::new(1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0);
|
|
|
|
let e1 = 1.0_f64.exp();
|
|
let e2 = 2.0_f64.exp();
|
|
|
|
let f = Matrix3::new(
|
|
e1,
|
|
3.0 * e1,
|
|
15.0 * (e2 - 2.0 * e1),
|
|
0.0,
|
|
e1,
|
|
5.0 * (e2 - e1),
|
|
0.0,
|
|
0.0,
|
|
e2,
|
|
);
|
|
|
|
assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn exp_dynamic() {
|
|
use nalgebra::DMatrix;
|
|
|
|
let m = DMatrix::from_row_slice(3, 3, &[1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0]);
|
|
|
|
let e1 = 1.0_f64.exp();
|
|
let e2 = 2.0_f64.exp();
|
|
|
|
let f = DMatrix::from_row_slice(
|
|
3,
|
|
3,
|
|
&[
|
|
e1,
|
|
3.0 * e1,
|
|
15.0 * (e2 - 2.0 * e1),
|
|
0.0,
|
|
e1,
|
|
5.0 * (e2 - e1),
|
|
0.0,
|
|
0.0,
|
|
e2,
|
|
],
|
|
);
|
|
|
|
assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
|
|
}
|
|
}
|