98 lines
3.1 KiB
Rust
98 lines
3.1 KiB
Rust
|
#[cfg(test)]
|
||
|
|
mod tests {
|
||
|
|
//https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py
|
||
|
|
#[test]
|
||
|
|
fn exp() {
|
||
|
|
use nalgebra::{Matrix1, Matrix2, Matrix3};
|
||
|
|
|
||
|
|
{
|
||
|
|
let m = Matrix1::new(1.0);
|
||
|
|
|
||
|
|
let f = m.exp();
|
||
|
|
|
||
|
|
assert_eq!(f, Matrix1::new(1_f64.exp()));
|
||
|
|
}
|
||
|
|
|
||
|
|
{
|
||
|
|
let m = Matrix2::new(0.0, 1.0, 0.0, 0.0);
|
||
|
|
|
||
|
|
assert_eq!(m.exp(), Matrix2::new(1.0, 1.0, 0.0, 1.0));
|
||
|
|
}
|
||
|
|
|
||
|
|
{
|
||
|
|
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!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
|
||
|
|
}
|
||
|
|
|
||
|
|
{
|
||
|
|
// 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!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
|
||
|
|
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!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|