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

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
``````

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

``````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]

I think you really need is

``````real.(r)
``````

not

``````abs.(r)
``````
2 Likes

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)
`````` 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(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:
 min(::Complex{Float64}, ::Float64) at .\operators.jl:431
 expand_extrema! at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\axes.jl:335 [inlined]
 expand_extrema!(::Plots.Axis, ::Array{Complex{Float64},1}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\axes.jl:358
 expand_extrema!(::Plots.Subplot{Plots.GRBackend}, ::Dict{Symbol,Any}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\axes.jl:388
 _expand_subplot_extrema(::Plots.Subplot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Symbol) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\pipeline.jl:366
 _process_seriesrecipe(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\pipeline.jl:402
 _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
 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
 plot(::Array{Complex{Float64},1}, ::Array{Complex{Float64},1}) at C:\Users\hermesr\.julia\packages\Plots\qZHsp\src\plot.jl:51
 top-level scope at REPL:1
``````

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

``````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)
`````` where Section 6 talks about the method you linked. TL;DR: just use `exp`.