Symbolic matrix exponential

General Usage
,
1

Hello All! I am fairly new to Julia.

I want to find the matrix exponential of a symbolic matrix. How to do this?
When I use the Symbolics and LinearAlgebra packages together, it doesn’t work:

using Symbolics
using LinearAlgebra

function run()
@variables a b c t
M = [ -a b 0
     1 c -a
     0 1 0]
Mt = M*t
expMt = exp(Mt)
print(expMt[1,1])
end

@time run()

Any help will be very appreciated!

2

Open an issue on Symbolics.jl. This isn’t implemented, and I’m not convinced there is a general symbolic solution that can be given on this operation. It’s worth a dev discussion.

1 Like
3

Okay. I did that. I hope this gets implemented. Thanks.

4

I’m not convinced it’s possible in general :sweat_smile: .

I just took a look and saw SymPy can only do it if the Jordan normal form exists for A. That’s probably a fine enough restriction.

1 Like
5

Wolfram Alpha is able ot spit out a result.
https://www.wolframalpha.com/input?i2d=true&i=matrix+exponential+%7B%7B-a%2Cb%2C0%7D%2C%7B1%2Cc%2C-a%7D%2C%7B0%2C1%2C0%7D%7D

image
image1166×107 5.03 KB

I can’t say if it’s correct. It looks hideous - but I think that’s to be expected.

6

The Jordan form exists for all matrices. The only question is if it can be computed symbolically; this might require giving a name to solutions of polynomials of degree 5+ and work with them. Can Symbolics.jl do it?

1 Like
7

That makes sense.

Do you think Wolfram Alpha is spitting out the right results then? I was confused as to what “omega” meant in my screenshot of Wolfram Alpha - and in light of your comment, I suspect it’s a solution to some >5 degree polynomial.

8

From the image it look like ω are the solutions (plural) to a 3rd degree polynomial. The expressions in ω are being summed over all solutions to this polynomial. Didn’t check it, but this is probably the characteristic polynomial of the matrix M. For 3x3 matrix there is a closed form of the solution. For >=5 matrices this would be more problematic and will need exactly the sort of notation Wolfram Alpha uses.

9

This case would be a degree 3 polynomial. But Symbolics hasn’t specialized the <5 cases just because this breaks down at 5x5 so no one has had the motivation yet to write a scheme that works for and 4x4 and smaller matrix only.

10

Approximating the symbolic expression for exp(Mt)[1,1] might be enough for the application. In that case:

using TaylorSeries, LinearAlgebra, Symbolics

@variables a b c t m i
M = [ -a b 0
      1 c -a
      0 1 0 ]
Mt = M*t

k = Taylor1(Rational{Int}, 4)
ek = convert(Taylor1{Rational{Int}},exp(k))

With these:

julia> Mt11 = simplify(substitute(substitute(expand(i*ek(m)), Dict(i => I)), Dict(m => Mt))[1,1])
1 + (1//24)*((t^4)*((b + a^2)^2) + (c*(t^2) - a*(t^2))*(b*c*(t^2) - a*b*(t^2)) - a*b*(t^4)) + (1//6)*(b*t*(c*(t^2) - a*(t^2)) - a*(b + a^2)*(t^3)) + (1//2)*(b + a^2)*(t^2) - a*t

julia> simplify(substitute(Mt11, Dict(a => 1, b => 1//2, c => 1//3, t => 1//5)))
446951//540000

The number of terms in the Taylor expansion is controlled by definition of the k Taylor1 variable.

A formal bound on the approximation quality can be calculated using spectral norm of M (in the usual Taylor series method).