# Matrix exponential MatrixExp[] vs Sum[MatrixPower[]] doesn't match?

I might be an idiot, but I cannot get the manual expansion of $e^{At}$ to match the MatrixExp[A t] result. For example, I have the following:

k = 1.0; c = 3.0;
A = {{0., 1.}, {-k, -c}};
MatrixExp[A t] /. {t -> 20}


which produces:

$$\left( \begin{array}{cc} 0.000563347 & 0.000215179 \\ -0.000215179 & -0.0000821912 \\ \end{array} \right)$$

But if I do the expansion manually:

eAt = Sum[(MatrixPower[A, k] t^k)/(k!), {k, 0, 100}];
eAt /. {t -> 30}


I get the result:

$\left( \begin{array}{cc} -2.59526\times 10^{30} & -6.79448\times 10^{30} \\ 6.79448\times 10^{30} & 1.77882\times 10^{31} \\ \end{array} \right)$

This is clearly incorrect, as the A matrix has eigenvalues of -1, -3 so it must be stable. Please help.

(see http://en.wikipedia.org/wiki/Matrix_exponential for matrix exponential definition)

-
Catastrophic numerical error. i.stack.imgur.com/KykyF.png –  Rahul May 8 '14 at 1:52
Is it because you are using "k" in your definition of "A" and as summation index? –  BlacKow May 8 '14 at 2:22

Well, it turns out you are doing the computation with low numerical precision. And this error propagates. If you use high enough precision (infinite maybe), the results turns out fine. Also since you're using a series approximation, including more terms also helps. Here it is:

Let's define A:

A = {{0, 1}, {-1, -3}};


Then

Sum[MatrixPower[20 A, s]/s!, {s, 0, 200}] // N


Gives:

{{0.000563346576, 0.000215179245}, {-0.000215179245, -0.0000821911578}}


Which is what you got before. Notice I include the k = 20 in the MatrixPower definition.

Or you can use finite precision but make sure to carry out the calculation with high enough precision as shown:

A = N[{{0, 1}, {-1, -3}}, 10];


Then:

Block[{$MinPrecision = 20,$MaxPrecision = 20}, Sum[MatrixPower[20 A, s]/s!, {s, 0, 200}]]


Gives as before (only with higher precision)

{{0.00056334657611228189882, 0.00021517924462901202477},
{-0.00021517924462901202477, -0.000082191157774754229696}}

-
This is on target but glosses over something important. Even with exact arithmetic, at t=30 it is not sufficient to sum the first hundred terms of the power series. One requires upwards of 250 terms to get something in the ballpark of a machine precision result. This can be seen using In[246]:= nn = 250; N[MatrixPower[Rationalize[A], nn]*30^nn/(nn!)] Out[247]= {{-3.14113879593*10^-20, -8.22360813112*10^-20}, \ {8.22360813112*10^-20, 2.15296855974*10^-19}}. –  Daniel Lichtblau Nov 11 '14 at 23:07