Skip to content

Conversation

@clinssen
Copy link
Contributor

@clinssen clinssen commented Nov 26, 2025

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
Copy link

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                 at⎞              ⎛ 2                         at⎞⎤
⎢     RootSumw  - 2w + 2, aRootSumw  - 2w + 2, a ↦ (1 - a)⋅   ⎠⎥
⎢     ───────────────────────────────       ───────────────────────────────────────⎥
⎢                    2                                         2                   ⎥
⎢                                                                                  ⎥
⎢        ⎛ 2                         at⎞              ⎛ 2                 at⎞    ⎥
⎢-RootSumw  - 2w + 2, a ↦ (1 - a)⋅RootSumw  - 2w + 2, a   ⎠    ⎥
⎢─────────────────────────────────────────      ───────────────────────────────    ⎥
⎣                    2                                         2In [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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants