added skalna06 algorithm for PILS (#115)

* added skalna06 algorithm for PILS

* allow rhs to be non-parametric

* allow apa to be called with single vectors

* renamed @linvars --> @affinevars

* update docstring

* simplified apa structures and updated FEM docs

* updated docstring and FEM example

* retrigger workflow

Author: lucaferranti

docs/literate/applications/FEM_example.jl

@@ -56,10 +56,10 @@
 function unitaryStiffnessMatrix( coordFirstNode, coordSecondNode  )
   diff      = (coordSecondNode - coordFirstNode)
   length   = sqrt( diff' * diff )
-  c        = diff[1] / length ;
-  s        = diff[2] / length ;
-  Qloc2glo = [ c -s 0 0 ; s c 0 0 ; 0 0 c -s ; 0 0 s c ] ;
-  Kloc     = [ 1 0 -1 0 ; 0 0 0 0 ; -1 0 1 0 ; 0 0 0 0 ] ;
+  c        = diff[1] / length
+  s        = diff[2] / length
+  Qloc2glo = [ c -s 0 0 ; s c 0 0 ; 0 0 c -s ; 0 0 s c ]
+  Kloc     = [ 1 0 -1 0 ; 0 0 0 0 ; -1 0 1 0 ; 0 0 0 0 ]
   Kglo     = Qloc2glo * Kloc * transpose(Qloc2glo)
   return     Kglo, length
 end
@@ -79,11 +79,11 @@ E = 2e11 ; # Young modulus
 A = 5e-3 ; # Cross-section area
 # The coordinate matrix is given by
 ## coordinates   x  y
-nodesCMatrix = [ 0. 0. ;
-                 1. 1. ;
-                 2. 0. ;
-                 3. 1. ;
-                 4. 0. ];
+nodesCMatrix = [ 0.0 0.0 ;
+                 1.0 1.0 ;
+                 2.0 0.0 ;
+                 3.0 1.0 ;
+                 4.0 0.0 ];
 # the connectivity matrix is given by
 ## connectivity  start end
 connecMatrix = [ 1     2 ;
@@ -92,54 +92,54 @@ connecMatrix = [ 1     2 ;
                  2     4 ;
                  3     4 ;
                  3     5 ;
-                 4     5 ];
+                 4     5 ]
 # and the fixed degrees of freedom (supports) are defined by the vector
-fixedDofs     = [2 9 10 ];
+fixedDofs = [2 9 10 ]
 #
 # calculations
-numNodes = size( nodesCMatrix )[1]; # compute the number of nodes
-numElems = size( connecMatrix )[1]; # compute the number of elements
-freeDofs = zeros(Int8, 2*numNodes-length(fixedDofs));
-indDof  = 1 ; counter = 0 ;
+numNodes = size( nodesCMatrix )[1]  # compute the number of nodes
+numElems = size( connecMatrix )[1]  # compute the number of elements
+freeDofs = zeros(Int8, 2*numNodes-length(fixedDofs))
+indDof  = 1 ; counter = 0
 while indDof <= (2*numNodes)
   if !(indDof in fixedDofs)
-    global counter = counter + 1 ;
-    freeDofs[ counter ] = indDof ;
+    global counter = counter + 1
+    freeDofs[ counter ] = indDof
   end
-  global indDof = indDof + 1 ;
+  global indDof = indDof + 1
 end
 #
 # assembly
-KG = zeros( 2*numNodes, 2*numNodes );
-FG = zeros( 2*numNodes );
+KG = zeros( 2*numNodes, 2*numNodes )
+FG = zeros( 2*numNodes )
 for elem in 1:numElems
   print(" assembling stiffness matrix of element ", elem , "\n")
   indexFirstNode  = connecMatrix[ elem, 1 ]
   indexSecondNode = connecMatrix[ elem, 2 ]
   dofsElem = [2*indexFirstNode-1 2*indexFirstNode 2*indexSecondNode-1 2*indexSecondNode ]
   KGelem, lengthElem = unitaryStiffnessMatrix( nodesCMatrix[ indexSecondNode, : ], nodesCMatrix[ indexFirstNode, : ] )
-  stiffnessParam = E * A / lengthElem ;
+  stiffnessParam = E * A / lengthElem
   for i in 1:4
     for j in 1:4
       KG[ dofsElem[i], dofsElem[j] ] = KG[ dofsElem[i], dofsElem[j] ] + stiffnessParam * KGelem[i,j]
     end
   end
 end
 FG[4] = -1e4 ;
-KG = KG[ freeDofs, : ] ;
-KG = KG[ :, freeDofs ] ;
-FG = FG[ freeDofs ];
+KG = KG[ freeDofs, : ]
+KG = KG[ :, freeDofs ]
+FG = FG[ freeDofs ]
 
 u = KG \ FG
-UG = zeros( 2*numNodes );
-UG[ freeDofs ] = u ;
+UG = zeros( 2*numNodes )
+UG[ freeDofs ] = u
 #
 # #### Deformed structure
 #
 # The reference (dashed blue line) and deformed (solid red)  configurations of the structure are ploted. Since the displacements are very small, a `scaleFactor` is considered to amplify the deformation and ease the visualization.
 #
 using Plots
-scaleFactor = 2e3 ;
+scaleFactor = 2e3
 plot();
 for elem in 1:numElems
   indexFirstNode  = connecMatrix[ elem, 1 ];
@@ -165,5 +165,59 @@ savefig("deformed.png") # hide
 # ![](deformed.png)
 #
 # ### Problem with interval parameters
-#
-# TO DO
+
+md"""
+Suppose now we have a 10% uncertainty for the stiffness $s_{23}$ associated with the third
+element. To model the problem, we introduce the symbolic variable `s23` using the IntervalLinearAlgebra macro `@linvars`.
+"""
+
+using IntervalLinearAlgebra
+@affinevars s23
+
+# now we can construct the matrix as before
+
+KGp = zeros(AffineExpression{Float64}, 2*numNodes, 2*numNodes );
+for elem in 1:numElems
+  print(" assembling stiffness matrix of element ", elem , "\n")
+  indexFirstNode  = connecMatrix[ elem, 1 ]
+  indexSecondNode = connecMatrix[ elem, 2 ]
+  dofsElem = [2*indexFirstNode-1 2*indexFirstNode 2*indexSecondNode-1 2*indexSecondNode ]
+  KGelem, lengthElem = unitaryStiffnessMatrix( nodesCMatrix[ indexSecondNode, : ], nodesCMatrix[ indexFirstNode, : ] )
+  if elem == 3
+    stiffnessParam = s23
+  else
+    stiffnessParam = E * A / lengthElem
+  end
+  for i in 1:4
+    for j in 1:4
+      KGp[ dofsElem[i], dofsElem[j] ] = KGp[ dofsElem[i], dofsElem[j] ] + stiffnessParam * KGelem[i,j]
+    end
+  end
+end
+KGp = KGp[ freeDofs, : ]
+KGp = KGp[ :, freeDofs ]
+
+# Now we can construct the [`AffineParametricArray`](@ref)
+
+KGp = AffineParametricArray(KGp)
+
+# The range of the stiffness is
+
+srange = E * A / sqrt(2) ± 0.1 * E * A / sqrt(2)
+
+# To solve the system, we could of course just subsitute `srange` into the parametric matrix
+# `KGp` and solve the "normal" interval linear system
+
+usimple = solve(KGp(srange), Interval.(FG))
+
+# This approach, however soffers from the [dependency problem](https://en.wikipedia.org/wiki/Interval_arithmetic#Dependency_problem)
+# and hence the compute displacement will be an overestimation of the true displacement.
+
+# To mitigate this issue, algorithms to solve linear systems with parameters have been developed.
+# In this case we use the algorithm presented in [[SKA06]](@ref)
+
+uparam = solve(KGp, FG, srange)
+
+# We can now compare the naive and parametric solution
+
+hcat(usimple, uparam)/1e-6

docs/src/api/algorithms.md

@@ -27,4 +27,12 @@ LinearKrawczyk
 ```@docs
 LinearOettliPrager
 NonLinearOettliPrager
+```
+
+## Parametric Solvers
+
+### Direct Methods
+
+```@docs
+Skalna06
 ```

src/IntervalLinearAlgebra.jl

@@ -22,7 +22,7 @@ export
     is_H_matrix, is_strongly_regular, is_strictly_diagonally_dominant, is_Z_matrix, is_M_matrix,
     rref,
     eigenbox, Rohn, Hertz, verify_eigen, bound_perron_frobenius_eigenvalue,
-    AffineExpression, @linvars,
+    AffineExpression, @affinevars,
     AffineParametricArray, AffineParametricMatrix, AffineParametricVector
 
 
@@ -35,8 +35,9 @@ include("multiplication.jl")
 include("utils.jl")
 include("classify.jl")
 include("rref.jl")
-include("pils/linexpr.jl")
+include("pils/affine_expressions.jl")
 include("pils/affine_parametric_array.jl")
+include("pils/pils_solvers.jl")
 
 include("eigenvalues/interval_eigenvalues.jl")
 include("eigenvalues/verify_eigs.jl")

src/pils/linexpr.jl src/pils/affine_expressions.jlsrc/pils/linexpr.jl renamed to src/pils/affine_expressions.jl

@@ -8,7 +8,7 @@ Data structure to represent affine expressions, such as ``x+2y+z+4``.
 ### Examples
 
 ```jldoctest
-julia> @linvars x y z
+julia> @affinevars x y z
 3-element Vector{AffineExpression{Int64}}:
  x
  y
@@ -45,32 +45,32 @@ function _get_vars(x...)
 end
 
 """
-    @linvars(x...)
+    @affinevars(x...)
 
 Macro to construct the variables used to represent [`AffineExpression`](@ref).
 
 ### Examples
 
 ```jldoctest
-julia> @linvars x
+julia> @affinevars x
 1-element Vector{AffineExpression{Int64}}:
  x
 
-julia> @linvars x y z
+julia> @affinevars x y z
 3-element Vector{AffineExpression{Int64}}:
  x
  y
  z
 
-julia> @linvars x[1:4]
+julia> @affinevars x[1:4]
 4-element Vector{AffineExpression{Int64}}:
  x1
  x2
  x3
  x4
 ```
 """
-macro linvars(x...)
+macro affinevars(x...)
     vars = _get_vars(x...)
     _vars_dict[:vars] = vars
 

src/pils/affine_parametric_array.jl

@@ -1,99 +1,99 @@
-struct AffineParametricArray{T, N, MT<:AbstractArray{T, N}, VAR} <: AbstractArray{T, N}
+struct AffineParametricArray{T, N, MT<:AbstractArray{T, N}} <: AbstractArray{T, N}
     coeffs::Vector{MT}
-    vars::Vector{VAR}
 end
 
 const AffineParametricMatrix{T, MT} = AffineParametricArray{T, 2, MT} where {T, MT<:AbstractMatrix{T}}
 
 const AffineParametricVector{T, VT} = AffineParametricArray{T, 1, VT} where {T, VT <: AbstractVector{T}}
 
 # apas are callable
-function (apa::AffineParametricArray)(p::AbstractVector)
+function (apa::AffineParametricArray)(p)
     length(p) + 1 == length(apa.coeffs) || throw(ArgumentError("dimension mismatch"))
     return sum(apa.coeffs[i] * p[i] for i in eachindex(p)) + apa.coeffs[end]
 end
 
-function (apa::AffineParametricArray)()
-   return apa.coeffs[end] + sum(apa.coeffs[i] * apa.vars[i] for i in eachindex(apa.vars))
-end
-
 # Array interface
 IndexStyle(::Type{<:AffineParametricArray}) = IndexLinear()
 size(apa::AffineParametricArray) = size(apa.coeffs[1])
 
-getindex(apa::AffineParametricArray, idx...) = getindex(apa(), idx...)
+function getindex(apa::AffineParametricArray, idx...)
+    nvars = length(_vars_dict[:vars])
+    vars = [AffineExpression(Vector(c)) for c in eachcol(Matrix(I, nvars + 1, nvars))]
+    tmp =  sum(getindex(c, idx...) * v for (v, c) in zip(vars, apa.coeffs[1:end-1]))
+    return getindex(apa.coeffs[end], idx...) + tmp
+end
+
+function setindex!(apa::AffineParametricArray, ae::AffineExpression, idx...)
+    @inbounds for i in eachindex(ae.coeffs)
+        setindex!(apa.coeffs[i], ae.coeffs[i], idx...)
+    end
+end
 
-# operations
-function ==(apa1::AffineParametricArray, apa2::AffineParametricArray)
-    return apa1.vars == apa2.vars && apa1.coeffs == apa2.coeffs
+function setindex!(apa::AffineParametricArray, num::Number, idx...)
+    setindex!(apa.coeffs[end], num, idx...)
 end
 
+==(apa1::AffineParametricArray, apa2::AffineParametricArray) = apa1.coeffs == apa2.coeffs
+
 # unary operations
 +(apa::AffineParametricArray) = apa
--(apa::AffineParametricArray) = AffineParametricArray(-apa.coeffs, apa.vars)
+-(apa::AffineParametricArray) = AffineParametricArray(-apa.coeffs)
 
 # addition subtraction
 for op in (:+, :-)
     @eval function $op(apa1::AffineParametricArray, apa2::AffineParametricArray)
-        return AffineParametricArray($op(apa1.coeffs, apa2.coeffs), apa1.vars)
+        return AffineParametricArray($op(apa1.coeffs, apa2.coeffs))
     end
 
-    @eval function $op(apa::AffineParametricArray{T, N, MT, VAR}, B::MS) where {T, N, MT, VAR, S, MS<:AbstractArray{S, N}}
+    @eval function $op(apa::AffineParametricArray{T, N, MT}, B::MS) where {T, N, MT, S, MS<:AbstractArray{S, N}}
         MTS = promote_type(MT, MS)
         coeffs = similar(apa.coeffs, MTS)
         coeffs .= apa.coeffs
         coeffs[end] = $op(coeffs[end], B)
-        return AffineParametricArray(coeffs, apa.vars)
+        return AffineParametricArray(coeffs)
     end
 
-    @eval function $op(B::MS, apa::AffineParametricArray{T, N, MT, VAR}) where {T, N, MT, VAR, S, MS<:AbstractArray{S, N}}
+    @eval function $op(B::MS, apa::AffineParametricArray{T, N, MT}) where {T, N, MT, S, MS<:AbstractArray{S, N}}
         MTS = promote_type(MT, MS)
         coeffs = similar(apa.coeffs, MTS)
         coeffs .= $op(apa.coeffs)
         coeffs[end] = $op(B, apa.coeffs[end])
-        return AffineParametricArray(coeffs, apa.vars)
+        return AffineParametricArray(coeffs)
     end
 end
 
 # multiplication, backslash
 for op in (:*, :\)
     @eval function $op(A::AbstractMatrix, apa::AffineParametricArray)
         coeffs = [$op(A, coeff) for coeff in apa.coeffs]
-        return AffineParametricArray(coeffs, apa.vars)
+        return AffineParametricArray(coeffs)
     end
 end
 
 function *(apa::AffineParametricMatrix, B::AbstractMatrix)
     coeffs = [coeff*B for coeff in apa.coeffs]
-    return AffineParametricArray(coeffs, apa.vars)
+    return AffineParametricArray(coeffs)
 end
 
 function *(apa::AffineParametricMatrix, v::AbstractVector)
     coeffs = [coeff*v for coeff in apa.coeffs]
-    return AffineParametricArray(coeffs, apa.vars)
+    return AffineParametricArray(coeffs)
 end
 
-*(apa::AffineParametricArray, n::Number) = AffineParametricArray(apa.coeffs * n, apa.vars)
-*(n::Number, apa::AffineParametricArray) = AffineParametricArray(apa.coeffs * n, apa.vars)
-/(apa::AffineParametricArray, n::Number) = AffineParametricArray(apa.coeffs / n, apa.vars)
+*(apa::AffineParametricArray, n::Number) = AffineParametricArray(apa.coeffs * n)
+*(n::Number, apa::AffineParametricArray) = AffineParametricArray(apa.coeffs * n)
+/(apa::AffineParametricArray, n::Number) = AffineParametricArray(apa.coeffs / n)
 
 # with AffineExpression
-function AffineParametricArray(A::AbstractArray{<:AffineExpression}, vars::Vector{<:AffineExpression})
-    ncoeffs = length(vars) + 1
+function AffineParametricArray(A::AbstractArray{<:AffineExpression})
+    ncoeffs = length(first(A).coeffs)
     coeffs = [map(a -> a.coeffs[i], A) for i in 1:ncoeffs]
-    return AffineParametricArray(coeffs, vars)
+    return AffineParametricArray(coeffs)
 end
 
-function AffineParametricArray(A::AbstractArray{<:Number}, vars::Vector{<:AffineExpression})
-    ncoeffs = length(vars) + 1
+function AffineParametricArray(A::AbstractArray{<:Number})
+    ncoeffs = length(_vars_dict[:vars]) + 1
     coeffs = [zero(A) for _ in 1:ncoeffs]
     coeffs[end] = A
-    # vars = [AffineExpression(Vector(c)) for c in eachcol(Matrix{T}(I, ncoeffs, ncoeffs-1))]
-    return AffineParametricArray(coeffs, vars)
-end
-
-function setindex!(apa::AffineParametricArray, ae::AffineExpression, idx...)
-    @inbounds for i in eachindex(ae.coeffs)
-        setindex!(apa.coeffs[i], ae.coeffs[i], idx...)
-    end
+    return AffineParametricArray(coeffs)
 end