Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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? (http://en.wikipedia.org/wiki/Matrix_exponential)

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

1 Answer 1

up vote 7 down vote accepted

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}}
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.