1 #ifndef STAN_MATH_PRIM_FUN_MATRIX_EXP_2X2_HPP
2 #define STAN_MATH_PRIM_FUN_MATRIX_EXP_2X2_HPP
3
4 #include <stan/math/prim/fun/Eigen.hpp>
5 #include <stan/math/prim/fun/cosh.hpp>
6 #include <stan/math/prim/fun/exp.hpp>
7 #include <stan/math/prim/fun/sinh.hpp>
8 #include <stan/math/prim/fun/sqrt.hpp>
9 #include <stan/math/prim/fun/square.hpp>
10 #include <cmath>
11
12 namespace stan {
13 namespace math {
14
15 /**
16 * Return the matrix exponential of a 2x2 matrix. Reference for
17 * algorithm: http://mathworld.wolfram.com/MatrixExponential.html
18 * Note: algorithm only works if delta > 0;
19 *
20 * @tparam EigMat type of the matrix
21 * @param[in] A 2x2 matrix to exponentiate.
22 * @return Matrix exponential of A.
23 */
24 template <typename EigMat, require_eigen_t<EigMat>* = nullptr>
25 Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic>
matrix_exp_2x2(const EigMat & A)26 matrix_exp_2x2(const EigMat& A) {
27 using std::cosh;
28 using std::exp;
29 using std::sinh;
30 using std::sqrt;
31
32 using T = value_type_t<EigMat>;
33 T a = A(0, 0), b = A(0, 1), c = A(1, 0), d = A(1, 1), delta;
34 delta = sqrt(square(a - d) + 4 * b * c);
35
36 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> B(2, 2);
37 T half_delta = 0.5 * delta;
38 T cosh_half_delta = cosh(half_delta);
39 T sinh_half_delta = sinh(half_delta);
40 T exp_half_a_plus_d = exp(0.5 * (a + d));
41 T Two_exp_sinh = 2 * exp_half_a_plus_d * sinh_half_delta;
42 T delta_cosh = delta * cosh_half_delta;
43 T ad_sinh_half_delta = (a - d) * sinh_half_delta;
44
45 B(0, 0) = exp_half_a_plus_d * (delta_cosh + ad_sinh_half_delta);
46 B(0, 1) = b * Two_exp_sinh;
47 B(1, 0) = c * Two_exp_sinh;
48 B(1, 1) = exp_half_a_plus_d * (delta_cosh - ad_sinh_half_delta);
49
50 return B / delta;
51 }
52
53 } // namespace math
54 } // namespace stan
55
56 #endif
57