// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2009 Jitse Niesen // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either // version 3 of the License, or (at your option) any later version. // // Alternatively, you can redistribute it and/or // modify it under the terms of the GNU General Public License as // published by the Free Software Foundation; either version 2 of // the License, or (at your option) any later version. // // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the // GNU General Public License for more details. // // You should have received a copy of the GNU Lesser General Public // License and a copy of the GNU General Public License along with // Eigen. If not, see . #ifndef EIGEN_MATRIX_FUNCTIONS #define EIGEN_MATRIX_FUNCTIONS #include #include #include #include #include #include namespace Eigen { /** \ingroup Unsupported_modules * \defgroup MatrixFunctions_Module Matrix functions module * \brief This module aims to provide various methods for the computation of * matrix functions. * * To use this module, add * \code * #include * \endcode * at the start of your source file. * * This module defines the following MatrixBase methods. * - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine * - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine * - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential * - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions * - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine * - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine * * These methods are the main entry points to this module. * * %Matrix functions are defined as follows. Suppose that \f$ f \f$ * is an entire function (that is, a function on the complex plane * that is everywhere complex differentiable). Then its Taylor * series * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f] * converges to \f$ f(x) \f$. In this case, we can define the matrix * function by the same series: * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f] * */ #include "src/MatrixFunctions/MatrixExponential.h" #include "src/MatrixFunctions/MatrixFunction.h" /** \page matrixbaseextra MatrixBase methods defined in the MatrixFunctions module \ingroup MatrixFunctions_Module The remainder of the page documents the following MatrixBase methods which are defined in the MatrixFunctions module. \section matrixbase_cos MatrixBase::cos() Compute the matrix cosine. \code const MatrixFunctionReturnValue MatrixBase::cos() const \endcode \param[in] M a square matrix. \returns expression representing \f$ \cos(M) \f$. This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos(). \sa \ref matrixbase_sin "sin()" for an example. \section matrixbase_cosh MatrixBase::cosh() Compute the matrix hyberbolic cosine. \code const MatrixFunctionReturnValue MatrixBase::cosh() const \endcode \param[in] M a square matrix. \returns expression representing \f$ \cosh(M) \f$ This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh(). \sa \ref matrixbase_sinh "sinh()" for an example. \section matrixbase_exp MatrixBase::exp() Compute the matrix exponential. \code const MatrixExponentialReturnValue MatrixBase::exp() const \endcode \param[in] M matrix whose exponential is to be computed. \returns expression representing the matrix exponential of \p M. The matrix exponential of \f$ M \f$ is defined by \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] The matrix exponential can be used to solve linear ordinary differential equations: the solution of \f$ y' = My \f$ with the initial condition \f$ y(0) = y_0 \f$ is given by \f$ y(t) = \exp(M) y_0 \f$. The cost of the computation is approximately \f$ 20 n^3 \f$ for matrices of size \f$ n \f$. The number 20 depends weakly on the norm of the matrix. The matrix exponential is computed using the scaling-and-squaring method combined with Padé approximation. The matrix is first rescaled, then the exponential of the reduced matrix is computed approximant, and then the rescaling is undone by repeated squaring. The degree of the Padé approximant is chosen such that the approximation error is less than the round-off error. However, errors may accumulate during the squaring phase. Details of the algorithm can be found in: Nicholas J. Higham, "The scaling and squaring method for the matrix exponential revisited," SIAM J. %Matrix Anal. Applic., 26:1179–1193, 2005. Example: The following program checks that \f[ \exp \left[ \begin{array}{ccc} 0 & \frac14\pi & 0 \\ -\frac14\pi & 0 & 0 \\ 0 & 0 & 0 \end{array} \right] = \left[ \begin{array}{ccc} \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 0 & 0 & 1 \end{array} \right]. \f] This corresponds to a rotation of \f$ \frac14\pi \f$ radians around the z-axis. \include MatrixExponential.cpp Output: \verbinclude MatrixExponential.out \note \p M has to be a matrix of \c float, \c double, \c complex or \c complex . \section matrixbase_matrixfunction MatrixBase::matrixFunction() Compute a matrix function. \code const MatrixFunctionReturnValue MatrixBase::matrixFunction(typename ei_stem_function::Scalar>::type f) const \endcode \param[in] M argument of matrix function, should be a square matrix. \param[in] f an entire function; \c f(x,n) should compute the n-th derivative of f at x. \returns expression representing \p f applied to \p M. Suppose that \p M is a matrix whose entries have type \c Scalar. Then, the second argument, \p f, should be a function with prototype \code ComplexScalar f(ComplexScalar, int) \endcode where \c ComplexScalar = \c std::complex if \c Scalar is real (e.g., \c float or \c double) and \c ComplexScalar = \c Scalar if \c Scalar is complex. The return value of \c f(x,n) should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. This routine uses the algorithm described in: Philip Davies and Nicholas J. Higham, "A Schur-Parlett algorithm for computing matrix functions", SIAM J. %Matrix Anal. Applic., 25:464–485, 2003. The actual work is done by the MatrixFunction class. Example: The following program checks that \f[ \exp \left[ \begin{array}{ccc} 0 & \frac14\pi & 0 \\ -\frac14\pi & 0 & 0 \\ 0 & 0 & 0 \end{array} \right] = \left[ \begin{array}{ccc} \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ 0 & 0 & 1 \end{array} \right]. \f] This corresponds to a rotation of \f$ \frac14\pi \f$ radians around the z-axis. This is the same example as used in the documentation of \ref matrixbase_exp "exp()". \include MatrixFunction.cpp Output: \verbinclude MatrixFunction.out Note that the function \c expfn is defined for complex numbers \c x, even though the matrix \c A is over the reals. Instead of \c expfn, we could also have used StdStemFunctions::exp: \code A.matrixFunction(StdStemFunctions >::exp, &B); \endcode \section matrixbase_sin MatrixBase::sin() Compute the matrix sine. \code const MatrixFunctionReturnValue MatrixBase::sin() const \endcode \param[in] M a square matrix. \returns expression representing \f$ \sin(M) \f$. This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin(). Example: \include MatrixSine.cpp Output: \verbinclude MatrixSine.out \section matrixbase_sinh const MatrixBase::sinh() Compute the matrix hyperbolic sine. \code MatrixFunctionReturnValue MatrixBase::sinh() const \endcode \param[in] M a square matrix. \returns expression representing \f$ \sinh(M) \f$ This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh(). Example: \include MatrixSinh.cpp Output: \verbinclude MatrixSinh.out */ } #endif // EIGEN_MATRIX_FUNCTIONS