Exponential matrix calculation with the exp() function vs procedure with the use of the eigen() function

New to Julia
,
1

Using the procedure shown in Matrix Exponentiation , I have performed the calculation of the exponential matrix:
for the case of the first matrix, the result of the procedure matches the result returned by the function exp ():

08:22:59->>2×2 Array{Int64,2}:
 1  -1
 2   4
Julia>exp(A)
08:22:59->>2×2 Array{Float64,2}:
 -5.30742  -12.6965
 25.393     32.782
Julia>F=eigen(A)
08:22:59->>Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
values:
2-element Array{Float64,1}:
 2.0
 3.0
vectors:
2×2 Array{Float64,2}:
 -0.707107   0.447214
  0.707107  -0.894427
Julia>v=F.vectors
08:22:59->>2×2 Array{Float64,2}:
 -0.707107   0.447214
  0.707107  -0.894427
Julia>Λ=exp.(F.values)
08:22:59->>2-element Array{Float64,1}:
  7.38905609893065
 20.085536923187668
Julia>D=diagm(Λ)
08:22:59->>2×2 Array{Float64,2}:
 7.38906   0.0
 0.0      20.0855
Julia>v*D*inv(v)
08:22:59->>2×2 Array{Float64,2}:
 -5.30742  -12.6965
 25.393     32.782

But for the second matrix, the results do not match !!
I would like to know why and what I should do to obtain the expected result

Julia>C=[-0.4 -1;1 0.45]
08:27:43->>2×2 Array{Float64,2}:
 -0.4  -1.0
  1.0   0.45
Julia>exp(C)
08:27:43->>2×2 Array{Float64,2}:
 0.254525  -0.890921
 0.890921   1.01181
Julia>F=eigen(C)
08:27:43->>Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}
values:
2-element Array{Complex{Float64},1}:
 0.02499999999999991 - 0.9051933495115836im
 0.02499999999999991 + 0.9051933495115836im
vectors:
2×2 Array{Complex{Float64},2}:
 -0.30052-0.640068im  -0.30052+0.640068im
 0.707107-0.0im       0.707107+0.0im
Julia>v=F.vectors
08:27:43->>2×2 Array{Complex{Float64},2}:
 -0.30052-0.640068im  -0.30052+0.640068im
 0.707107-0.0im       0.707107+0.0im
Julia>Λ=exp.(F.values)
08:27:43->>2-element Array{Complex{Float64},1}:
 0.6331664487902933 - 0.8064560400308953im
 0.6331664487902933 + 0.8064560400308953im
Julia>D=diagm(Λ)
08:27:43->>2×2 Array{Complex{Float64},2}:
 0.633166-0.806456im       0.0+0.0im
      0.0+0.0im       0.633166+0.806456im
Julia>v*D*inv(v)
08:27:43->>2×2 Array{Complex{Float64},2}:
 0.254525+0.0im          -0.890921+5.55112e-17im
 0.890921-5.55112e-17im    1.01181-1.11022e-16im
1 Like
2

It seems both are correct. Only the later one has imaginary parts which is about rounding error.

1 Like
3 
Julia>exp(C)
09:22:30->>2×2 Array{Float64,2}:
 0.254525  -0.890921
 0.890921   1.01181

Julia>r=v*D*inv(v)
09:22:38->>2×2 Array{Complex{Float64},2}:
 0.254525+0.0im          -0.890921+5.55112e-17im
 0.890921-5.55112e-17im    1.01181-1.11022e-16im

Julia>abs.(r)
09:22:38->>2×2 Array{Float64,2}:
 0.254525  0.890921
 0.890921  1.01181

Julia>abs.(r)==exp(C)
09:22:38->>false

The difference is in the sign of the element X [1,2]

4

I think you really need is

real.(r)

not

abs.(r)
3 Likes
5

using r = real. (…) I can plot the result:

C=[-0.4 -1;1 0.45]
x0=[0 0.22]
∆t = .25
F=eigen(C*∆t)
v=F.vectors
Λ=exp.(F.values)
D=diagm(Λ)
r=v*D*inv(v)
T=real.(r)
x = x0
x1 =x0 
for i = 1:100
x = x*T #repeatedly multiply by T
x1=vcat(x1, x) # & store current x1(t) in the array x1
end
p1=plot( x1)
p2=plot(x1[:,1],x1[:,2])
plot(p1, p2,layout = (2, 1), legend = false)

plot_Real

how to plot the results using it using imaginary numbers?:

C=[-0.4 -1;1 0.45]
x0=[0 0.22]
∆t = .25
F=eigen(C*∆t)
v=F.vectors
Λ=exp.(F.values)
D=diagm(Λ)
r=v*D*inv(v)
x = x0
x1 =x0 
for i = 1:100
x = x*r #repeatedly multiply by T
x1=vcat(x1, x) # & store current x1(t) in the array x1
end
plot(x1)
plot(x1[:,1],x1[:,2])

plot_Real Imag

plot(x1[:,1],x1[:,2])
ERROR: MethodError: no method matching isless(::Float64, ::Complex{Float64})
Closest candidates are:
  isless(::Float64, ::Float64) at float.jl:465
  isless(::Missing, ::Any) at missing.jl:87
  isless(::AbstractFloat, ::AbstractFloat) at operators.jl:165
  ...
Stacktrace:
 [1] min(::Complex{Float64}, ::Float64) at .\operators.jl:431
 [2] expand_extrema! at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\axes.jl:335 [inlined]
 [3] expand_extrema!(::Plots.Axis, ::Array{Complex{Float64},1}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\axes.jl:358
 [4] expand_extrema!(::Plots.Subplot{Plots.GRBackend}, ::Dict{Symbol,Any}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\axes.jl:388
 [5] _expand_subplot_extrema(::Plots.Subplot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Symbol) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\pipeline.jl:366
 [6] _process_seriesrecipe(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\pipeline.jl:402
 [7] _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{Array{Complex{Float64},1},Array{Complex{Float64},1}}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\plot.jl:234
 [8] plot(::Array{Complex{Float64},1}, ::Vararg{Array{Complex{Float64},1},N} where N; kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\plot.jl:57
 [9] plot(::Array{Complex{Float64},1}, ::Array{Complex{Float64},1}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\plot.jl:51
 [10] top-level scope at REPL[154]:1
6

Maybe plot(real.(x1), imag.(x1))?

7 
p3=plot(real.(x1))
p4=plot(real.(x1[:,1]),real.(x1[:,2]))
p5=plot(real.(x1), imag.(x1))
plot(p3, p4,p5,layout = (3, 1), legend = false)

plot_Real 2

8

This is well-known to be a difficult problem numerically. Please read

https://www.researchgate.net/publication/228590312_Nineteen_Dubious_Ways_to_Compute_the_Exponential_of_a_Matrix_Twenty-Five_Years_Later

where Section 6 talks about the method you linked. TL;DR: just use exp.

2 Likes