Add alternative matrix exponential function

[clinssen]

Add contributed expMt function for faster matrix exponentiation.

CI failure will be fixed after merging with #93.

Before merging this, check in with @oscarbenjamin about licensing, citing, and functionality.

[oscarbenjamin]

Feel free to use the code under whichever license is most convenient for ode-toolbox. Presumably that is GPL V2?

I don't know if you would usually put someone's name in the copyright headers but if so feel free to put my name there.

Sorry that I haven't gotten round to including this in SymPy itself. One thing that I haven't quite resolved with this code is handling the conversion to sin/cos rather than exp:

In [2]: M = Matrix([[1, 1], [-1, 1]])

In [3]: M.exp()
Out[3]: 
⎡ℯ⋅cos(1)   ℯ⋅sin(1)⎤
⎢                   ⎥
⎣-ℯ⋅sin(1)  ℯ⋅cos(1)⎦

In [6]: expMt(M, t)
Out[6]: 
⎡            ⎛ 2                 a⋅t⎞              ⎛ 2                         a⋅t⎞⎤
⎢     RootSum⎝w  - 2⋅w + 2, a ↦ ℯ   ⎠       RootSum⎝w  - 2⋅w + 2, a ↦ (1 - a)⋅ℯ   ⎠⎥
⎢     ───────────────────────────────       ───────────────────────────────────────⎥
⎢                    2                                         2                   ⎥
⎢                                                                                  ⎥
⎢        ⎛ 2                         a⋅t⎞              ⎛ 2                 a⋅t⎞    ⎥
⎢-RootSum⎝w  - 2⋅w + 2, a ↦ (1 - a)⋅ℯ   ⎠       RootSum⎝w  - 2⋅w + 2, a ↦ ℯ   ⎠    ⎥
⎢─────────────────────────────────────────      ───────────────────────────────    ⎥
⎣                    2                                         2                   ⎦

In [7]: expMt(M, t).doit()
Out[7]: 
⎡    t⋅(1 - ⅈ)    t⋅(1 + ⅈ)        t⋅(1 - ⅈ)      t⋅(1 + ⅈ)⎤
⎢   ℯ            ℯ              ⅈ⋅ℯ            ⅈ⋅ℯ         ⎥
⎢   ────────── + ──────────     ──────────── - ────────────⎥
⎢       2            2               2              2      ⎥
⎢                                                          ⎥
⎢     t⋅(1 - ⅈ)      t⋅(1 + ⅈ)     t⋅(1 - ⅈ)    t⋅(1 + ⅈ)  ⎥
⎢  ⅈ⋅ℯ            ⅈ⋅ℯ             ℯ            ℯ           ⎥
⎢- ──────────── + ────────────    ────────── + ──────────  ⎥
⎣       2              2              2            2       ⎦

There should be something to handle the complex conjugate root pairs but in general symbolically it is difficult to separate real roots from complex conjugate pairs. When working with actual complex numbers it is most likely best just to stick with exponentials.

[clinssen]

Thank you for the very kind reply!

I'm afraid these internal magics are a little over my head, so I have added a link to this discussion thread too just for reference.