On May 24, 8:37 pm, luca <luca.frammol...@gmail.com> wrote: > Hi, > > i have the following problem: given a 3x3 real matrix, compute exp(A). > > I need a really fast way to do this. I have searched a bit with > google, but it seems to me that > computing the matrix exponential is not so simple, at least if your > matrix does not have a special > structure (for example A=diagonal matrix). > > I have found a simple method that use the diagonalization of A. If A > has 3 distinct eigenvalues, than compute > A=PDP^-1, where P is the matrix of the eigenvectors, D is a diagonal > matrix (whose diagonal elements are > the eigenvalues of A). Than, exp(A) = P exp(D) P^-1. Since P^-1 is > fast enough and exp(D) is simple > to compute, this should be a fast method. > > But, the problem is: i am not sure that the matrix A will always have > 3 distinct eigenvalues...what happens > if this does not happen? Can i use that formula even if 2 (or all > three) eigenvalues are equal? > > Are there any other ways to compute exp(A) in a fast way? > > Thank you, > Luca
If you have time to do the diagonalization (by similarity), then of course that is a way to get the matrix exponential.
If a matrix is not diagonalizable, then it can be written in Jordan canonical form, i.e. diagonal matrix plus some nilpotent blocks corresponding to eigenvalues with algebraic multiplicity greater than their geometric multiplicity.
The application of power series to such diagonal + nilpotent terms gives (esp. in the 3x3 case) just a slight modification to what you were already looking at for the exponential of a purely diagonal matrix.
Let's consider an illustrative case, where a eigenvalue r has geometric multiplicity 1 but algebraic multiplicity 3. Instead of being diagonalizable, such a matrix is then similar to J =
[ r 1 0 ] [ 0 r 1 ] [ 0 0 r ]
Since rI and the nilpotent part N of this commute, exp(J) = exp(r)I * exp(N). The nice thing about the power series exp(N) is that it has only a (small) finite number of terms:
exp(N) = I + N + (1/2)N^2 =
[ 1 1 0.5 ] [ 0 1 1 ] [ 0 0 1 ]
All the other cases are actually simpler than this one.