diff --git a/docs/src/lib/types.md b/docs/src/lib/types.md index 64348cc..4354e3f 100644 --- a/docs/src/lib/types.md +++ b/docs/src/lib/types.md @@ -28,3 +28,10 @@ IntervalMatrix ```@docs IntervalMatrixPower ``` + +## Affine interval matrix + +```@docs +AffineIntervalMatrix1 +AffineIntervalMatrix +``` diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl index 53dd7f0..0b25626 100644 --- a/src/IntervalMatrices.jl +++ b/src/IntervalMatrices.jl @@ -13,6 +13,7 @@ using Reexport # Type defining an interval matrix # ================================ include("matrix.jl") +include("affine.jl") # ================================= # Operations for interval matrices @@ -61,6 +62,8 @@ export IntervalMatrixPower, base, index +export AffineIntervalMatrix, + AffineIntervalMatrix1 set_multiplication_mode(config[:multiplication]) diff --git a/src/affine.jl b/src/affine.jl new file mode 100644 index 0000000..eb1ec79 --- /dev/null +++ b/src/affine.jl @@ -0,0 +1,123 @@ +""" + AffineIntervalMatrix1{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT} + +Interval matrix representing the matrix + +```math +A₀ + λA₁, +``` +where ``A₀`` and ``A₁`` are real (or complex) matrices, and ``λ`` is an interval. + +### Fields + +- `A0` -- matrix +- `A1` -- matrix +- `λ` -- interval + +### Examples + +The matrix ``I + [1 1; -1 1] * (0 .. 1)`` is: + +```jldoctest +julia> using LinearAlgebra + +julia> P = AffineIntervalMatrix1(Matrix(1.0I, 2, 2), [1 1; -1 1.], 0 .. 1); + +julia> P +2×2 AffineIntervalMatrix1{Float64, Interval{Float64}, Matrix{Float64}, Matrix{Float64}}: + [1, 2] [0, 1] + [-1, 0] [1, 2] +``` +""" +struct AffineIntervalMatrix1{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT} + A0::MT0 + A1::MT1 + λ::IT + + # inner constructor with dimension check + function AffineIntervalMatrix1(A0::MT0, A1::MT1, λ::IT) where {T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} + @assert checksquare(A0) == checksquare(A1) "the size of `A0` and `A1` should match, got $(size(A0)) and $(size(A1)) respectively" + return new{T, IT, MT0, MT1}(A0, A1, λ) + end +end + +IndexStyle(::Type{<:AffineIntervalMatrix1}) = IndexLinear() +size(M::AffineIntervalMatrix1) = size(M.A0) +getindex(M::AffineIntervalMatrix1, i::Int) = getindex(M.A0, i) + M.λ * getindex(M.A1, i) +function setindex!(M::AffineIntervalMatrix1{T}, X::T, inds...) where {T} + setindex!(M.A0, X, inds...) + setindex!(M.A1, zero(T), inds...) +end +copy(M::AffineIntervalMatrix1) = AffineIntervalMatrix1(copy(M.A0), copy(M.A1), M.λ) + +""" + AffineIntervalMatrix{T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}} <: AbstractIntervalMatrix{IT} + +Interval matrix representing the matrix + +```math +A₀ + λ₁A₁ + λ₂A₂ + … + λₖAₖ, +``` +where ``A₀`` and ``A₁, …, Aₖ`` are real (or complex) matrices, and ``λ₁, …, λₖ`` +are intervals. + +### Fields + +- `A0` -- matrix +- `A` -- vector of matrices +- `λ` -- vector of intervals + +### Notes + +This type is the general case of the [`AffineIntervalMatrix1`](@ref), which only +contains one matrix proportional to an interval. + +### Examples + +The affine matrix ``I + [1 1; -1 1] * (0 .. 1) + [0 1; 1 0] * (2 .. 3)`` is: + +```jldoctest +julia> using LinearAlgebra + +julia> A0 = Matrix(1.0I, 2, 2); + +julia> A1 = [1 1; -1 1.]; A2 = [0 1; 1 0]; + +julia> λ1 = 0 .. 1; λ2 = 2 .. 3; + +julia> P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ1]); + +julia> P[1, 1] +[1, 2] + +julia> P[1, 2] +[0, 2] +``` +""" +struct AffineIntervalMatrix{T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}, VIT<:AbstractVector{IT}} <: AbstractIntervalMatrix{IT} + A0::MT0 + A::MTA + λ::VIT + + # inner constructor with dimension check + function AffineIntervalMatrix(A0::MT0, A::MTA, λ::VIT) where {T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}, VIT<:AbstractVector{IT}} + k = length(A) + @assert k == length(λ) "expected `A` and `λ` to have the same length, got $(length(A)) and $(length(λ)) respectively" + n = checksquare(A0) + @inbounds for i in 1:k + @assert n == checksquare(A[i]) "each matrix should have the same size $n × $n" + end + return new{T, IT, MT0, MT, MTA, VIT}(A0, A, λ) + end +end + +IndexStyle(::Type{<:AffineIntervalMatrix}) = IndexLinear() +size(M::AffineIntervalMatrix) = size(M.A0) +getindex(M::AffineIntervalMatrix, i::Int) = getindex(M.A0, i) + sum(M.λ[k] * getindex(M.A[k], i) for k in eachindex(M.λ)) +function setindex!(M::AffineIntervalMatrix{T}, X::T, inds...) where {T} + setindex!(M.A0, X, inds...) + @inbounds for k in 1:length(M.A) + setindex!(M.A[k], zero(T), inds...) + end +end +copy(M::AffineIntervalMatrix) = AffineIntervalMatrix(copy(M.A0), copy(M.A), M.λ) diff --git a/test/affine.jl b/test/affine.jl new file mode 100644 index 0000000..c9f45ff --- /dev/null +++ b/test/affine.jl @@ -0,0 +1,35 @@ +@testset "Affine interval matrix in one interval" begin + A0 = [1 0; 0 1.] + A1 = [0 1; 1 0.] + λ = 0 .. 1 + P = AffineIntervalMatrix1(A0, A1, λ) + + @test size(P) == (2, 2) + @test P[1, 1] == interval(1) + @test P[1, 2] == interval(0, 1) + P[1, 1] = 2.0 + @test P[1, 1] == interval(2) + + # complex interval + λc = λ + im*λ/2 + Pc = AffineIntervalMatrix1(A0, A1, λc) + @test eltype(Pc) == Complex{Interval{Float64}} +end + +@testset "Affine interval matrix in several intervals" begin + A0 = [1 0; 0 1.] + A1 = [0 1; 1 0.]; A2 = copy(A1) + λ1 = 0 .. 1; λ2 = copy(λ1) + P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ2]) + + @test size(P) == (2, 2) + @test P[1, 1] == interval(1) + @test P[1, 2] == interval(0, 2) + P[1, 1] = 5.0 + @test P[1, 1] == interval(5) + + # complex interval + λ1c = λ1 + im*λ1/2; λ2c = λ2 - im*λ2/2 + Pc = AffineIntervalMatrix(A0, [A1, A2], [λ1c, λ2c]) + @test eltype(Pc) == Complex{Interval{Float64}} +end diff --git a/test/runtests.jl b/test/runtests.jl index 78a497f..9910562 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,3 +11,4 @@ include("constructors.jl") include("arithmetic.jl") include("setops.jl") include("exponential.jl") +include("affine.jl")