[eigen] Matrix exponential for evaluating a solution to a linear system of differential equations |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: [eigen] Matrix exponential for evaluating a solution to a linear system of differential equations*From*: Douglas Bates <bates@xxxxxxxxxxxxx>*Date*: Thu, 16 Feb 2012 16:24:05 -0600*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:sender:date:x-google-sender-auth:message-id:subject :from:to:content-type; bh=VR4CKdlufQsenbZbZDv5kBS/UmHQHKpBRavK2Y1a6sc=; b=F8+LmLTXx16lw8NnrDCYACb+opI8TT/Gc2w2dJTGSd4HUjFUjnJgM6sUqxL4IfI9zo WF59wxNrIODa+IWlmJSuu/q/9JF/0FhuW+29id5fmwOynBti5A6C+UumE51AsKd4H6WK Af0wguMmMUK1i2KFmCbyOesH9PFUCnB++wrMI=

As described in "Nineteen Dubious Ways of Computing the Exponential of a Matrix" and in the "...: Twenty-five Years Later" follow-up paper (SIAM Review, vol. 45, no. 1, 2003) there are many different approaches to computing the matrix exponential. I would appreciate recommendations on Eigen-based approaches to determining the solution to linear systems of ordinary differential equations using the matrix exponential. In the cases I will consider, the system matrix A is reasonably small, usually smaller than 5 by 5 and parameterized. (They are from compartment models in pharmacokinetics.) For given values of the parameters we need to evaluate e^{At} x_0 for several different values of t. We also may want to evaluate the sensitivities of the solution with respect to the parameters. The nature of the systems is such that there are no complex eigenvalues for A. I realize there is a matrix exponential function in the unsupported MatrixFunctions but I would prefer to use a decomposition method if possible. So, assuming that A is diagonalizable, we can use the Eigen::EigenSolver class and proceed happily. However, we need to detect if A is close to being non-diagonalizable. My current approach is static double precision(std::sqrt(std::numeric_limits<double>::epsilon())); Eigen::EigenSolver<Eigen::MatrixXd> es(A); if (es.eigenvalues().imag().squaredNorm() != 0) throw std::invalid_argument("Complex eigenvalues"); Eigen::VectorXd vals(es.eigenvalues().real()); if (es.eigenvectors().imag().squaredNorm() != 0) throw std::invalid_argument("Complex eigenvectors"); Eigen::MatrixXd vecs(es.eigenvectors().real()); Eigen::VectorXd alpha(Eigen::FullPivLU<Eigen::MatrixXd>(vecs).solve(x0)); bool nonsingular = (vecs * alpha).isApprox(x0, precision); but I am confident that can be improved. For one thing, even when A is non-diagonalizable the matrix of eigenvectors is regarded as nonsingular according to this test and that makes me uncomfortable. I could calculate the condition number from an Eigen::JacobiSVD but that seems to be a very long way around. Would it be better to check for a non-diagonalizable A using the T matrix from the Schur decomposition? Is there a simple way of getting the condition number of a triangular matrix or an estimate of it in Eigen (similar to dtrcon in Lapack)?

**Follow-Ups**:

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Aliasing and matrix multiplication** - Next by Date:
**Re: [eigen] Matrix exponential for evaluating a solution to a linear system of differential equations** - Previous by thread:
**Re: [eigen] Aliasing and matrix multiplication** - Next by thread:
**Re: [eigen] Matrix exponential for evaluating a solution to a linear system of differential equations**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |