# Is there any faster alternative for matrix exponential?

Lets define matrix $$M$$ as

M = {{ Cos[
t ]^2, - (2 Cos[t ] +
I  Sin[t ]) ((E^(-I t ) - E^(I 2 t )) (s + 2))}, { -8 I  Sin[
2 t ], (2  Cos[t ] -
I  Sin[t ]) ((E^(I t ) - E^(-I 2 t )) (-s + 2))}};


I'd like to find the matrix exponential of $$B$$, which can be done by MatrixExp. It works for the above matrix, and the result is long and ugly. But the matrix of my problem is much more complicated. It is a 6 by 6 matrix, and each element includes the summations and subtractions of trigonometric and exponential terms, similar to (but more complicated than) the above matrix. Even worse, after finding the matrix exponential, I need to take an integral over each element, say, from $$t_1$$ to $$t_2$$, where $$t_2>t_1$$. If I use MatrixExp it takes a very long time, and if I want to do the integration it becomes much longer. Especially, if I write Integrate[MatrixExp[M],{t,t1,t2}], I'm not sure whether it yields the results.

Do you know a faster alternative method that is suitable for more complicated functions?

Revise:

To clarify my problem well, I have a system of linear differential equations with "variable" coefficients: $$\frac{d}{dt}c=M(t) c+\lambda(t)+n(t)$$ In other words, the elements of $$M$$ are time dependent.

I already know how to find the solution if the coefficients are constants, but for the variable case the only method I know is this: Here $$\mathcal{T}$$ is called time-ordering operator and guarantees that earlier time acts first. As you see, this requires both integration and matrix exponential. Also, one should note that in the second term $$\tau$$ is the integration variable. So, if we can assign a number to $$t$$, we cannot do that for $$\tau$$.

• Unfortunately MatrixExp is already quite optimized. There are some implementations lying around the site, but ultimately you’re just doing different approximations of a matrix exponential that MatrixExp would do already. However, I will say two things: 1. Use numerical values, this speeds up evaluation by an incredible amount. And 2. Consider if you can apply a Magnus Expansion, and go that route. Much fewer instances of using MatrixExp and larger time steps. – CA Trevillian Dec 26 '19 at 8:14
• Do you actually need the matrix exponential, or do you just need the integration? If the latter, you can hold off on calculating the matrix exponential until you've assigned numbers within Integrate[ ] and calculate then. Many orders of magnitude faster. – MikeY Dec 26 '19 at 16:06
• Sometimes the heuristics for method selection used by MatrixExp do not give the "best" possibility. You could try doing it by hand using e.g. a Jordan decomposition, to see if that gives a better result. Also, would integrating first help? Or termwise integration of the Taylor expansion for the matrix exponential? – Daniel Lichtblau Dec 26 '19 at 16:43
• Saeid, are you trying to use MatrixExp as an integrator? – CA Trevillian Dec 26 '19 at 17:36
• @CATrevillian I have revised my question. See what I'm looking for. – Saeid Dec 26 '19 at 20:27

Just solve the differential equation without worrying about matrix exponentiation:

NDSolve[{c'[t] == M[t].c[t] + λ[t] + n[t], c == {1, 0}}, c, {t, 0, 1}]

• Is it possible to uncover the methods used within a solution done in this manner? Can you clarify this point? – CA Trevillian Dec 27 '19 at 5:48
• @CATrevillian there's a lot of info on methods on the NDSolve docu page. You can even define your own integration method and use it as a plugin. – Roman Dec 27 '19 at 11:17
• I thought I understood NDSolve, and now I know I know enough to not know enough! Thanks for the link to the plugin system. I’m overall curious how this compares to a Magnus Expansion, so I’ll have to do some investigating. – CA Trevillian Dec 27 '19 at 17:38

Not a full answer, since not sure exactly what you need, but laying it out here...

Assuming your s eventually takes a value, let s=1. Then you can get the numerical integration for the terms of your matrix rapidly by holding off on the matrix exponential until you've set numbers. Define your matrix B as a function:

bb[t_?NumericQ] := {{Cos[t]^2,
-(2 Cos[t] + I Sin[t]) ((E^(-I t) - E^(I 2 t)) (s + 2))},
{-8 I Sin[2 t],
(2 Cos[t] - I Sin[t]) ((E^(I t) - E^(-I 2 t)) (-s + 2))}}


Then NIntegrate. I am using Hold to delay picking off the element of the matrix until you have a matrix, otherwise you get a message that you can ignore anyway.

s = 1;
NIntegrate[Hold@(MatrixExp@bb[t])[[2, 1]], {t, 0, 2}]//Timing
(* 0., 864.408 + 1079.91 I *)


This is actually faster than getting the matrix exponential first in symbolic form, and then integrating.

meb = MatrixExp[B];
s = 1;
NIntegrate[meb[[2, 1]], {t, 0, 1}] // Timing
(* {0.03125, -780.294 + 999.49 I}  *)

• Can you clarify what this is actually doing? Awesome! – CA Trevillian Dec 26 '19 at 17:10