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!
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.
Okay. I did that. I hope this gets implemented. Thanks.
I’m not convinced it’s possible in general 
.
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.
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
I can’t say if it’s correct. It looks hideous - but I think that’s to be expected.
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?
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.
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.
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.
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).