matrix_exp_2x2.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_2X2_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_2X2_HPP
3 
6 
7 namespace stan {
8 namespace math {
9 
10 /**
11  * Return the matrix exponential of a 2x2 matrix. Reference for
12  * algorithm: http://mathworld.wolfram.com/MatrixExponential.html
13  * Note: algorithm only works if delta > 0;
14  *
15  * @tparam T type of scalar of the elements of input matrix.
16  * @param[in] A 2x2 matrix to exponentiate.
17  * @return Matrix exponential of A.
18  */
19 template <typename Mtype>
20 Mtype matrix_exp_2x2(const Mtype& A) {
21  using T = typename Mtype::Scalar;
22  T a = A(0, 0), b = A(0, 1), c = A(1, 0), d = A(1, 1), delta;
23  delta = sqrt(square(a - d) + 4 * b * c);
24 
25  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> B(2, 2);
26  T half_delta = 0.5 * delta;
27  T cosh_half_delta = cosh(half_delta);
28  T sinh_half_delta = sinh(half_delta);
29  T exp_half_a_plus_d = exp(0.5 * (a + d));
30  T Two_exp_sinh = 2 * exp_half_a_plus_d * sinh_half_delta;
31  T delta_cosh = delta * cosh_half_delta;
32  T ad_sinh_half_delta = (a - d) * sinh_half_delta;
33 
34  B(0, 0) = exp_half_a_plus_d * (delta_cosh + ad_sinh_half_delta);
35  B(0, 1) = b * Two_exp_sinh;
36  B(1, 0) = c * Two_exp_sinh;
37  B(1, 1) = exp_half_a_plus_d * (delta_cosh - ad_sinh_half_delta);
38 
39  return B / delta;
40 }
41 } // namespace math
42 } // namespace stan
43 #endif
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:12
double delta
Definition: runWimpSim.h:98
fvar< T > cosh(const fvar< T > &x)
Definition: cosh.hpp:11
Mtype matrix_exp_2x2(const Mtype &A)
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
const double a
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Float_t d
Definition: plot.C:236
static const double A
Definition: Units.h:82
fvar< T > sinh(const fvar< T > &x)
Definition: sinh.hpp:10
const hit & b
Definition: hits.cxx:21
double T
Definition: Xdiff_gwt.C:5