WIP #88 - Symbolic interval matrix power

Project.toml

@@ -7,6 +7,7 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"
 
 [compat]
 julia = "1"

src/IntervalMatrices.jl

@@ -6,6 +6,8 @@ using LinearAlgebra
 import Random
 using Random: AbstractRNG, GLOBAL_RNG, seed!
 
+import SymEngine
+
 using Reexport
 @reexport using IntervalArithmetic
 
@@ -18,6 +20,7 @@ include("matrix.jl")
 # Arithmetic for interval matrices
 # =================================
 include("arithmetic.jl")
+include("power.jl")
 
 # =======================================================
 # Methods to handle the exponential of an interval matrix

src/power.jl

@@ -0,0 +1,27 @@
+import Base: ^
+
+function ^(M::IntervalMatrix{T}, k::Integer) where {T}
+    n = checksquare(M)
+
+    # create symbolic matrix
+    Ms = Matrix{SymEngine.Basic}(undef, n, n)
+    dict = Dict{SymEngine.Basic, Interval{T}}()
+    for i in 1:n, j in 1:n
+        var = SymEngine.symbols("M_$(i)_$(j)")
+        Ms[i, j] = var
+        dict[var] = M[i, j]
+    end
+
+    # compute symbolic matrix power
+    Msᵏ = Ms^k
+
+    # substitute intervals in symbolic result
+    Mᵏ = similar(M)
+    for i in 1:n, j in 1:n
+        println(Msᵏ[i, j])
+        println(dict)
+        substituted = SymEngine.subs(Msᵏ[i, j], dict)
+        Mᵏ[i, j] = SymEngine.lambdify(substituted)
+    end
+    return Mᵏ
+end