Skip to content

A package for representing quasi arrays with continuous indices

License

Notifications You must be signed in to change notification settings

JuliaApproximation/ContinuumArrays.jl

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

ContinuumArrays.jl

A package for representing quasi arrays with continuous dimensions

Build Status codecov docs

A quasi array as implemented in QuasiArrays.jl is a generalization of an array that allows non-integer indexing via general axes. This package adds support for infinite-dimensional axes, including continuous intervals. Thus it plays the same role as InfiniteArrays.jl does for standard arrays but now for quasi arrays.

A simple example is the identity function on the interval 0..1. This can be created using Inclusion(d), which returns x if x in d is true, otherwise throws an error:

julia> using ContinuumArrays, IntervalSets

julia> x = Inclusion(0..1.0)
Inclusion(0.0..1.0)

julia> size(x) # uncountable (aleph-1)
(ℵ₁,)

julia> axes(x) # axis is itself
(Inclusion(0.0..1.0),)

julia> x[0.1] # returns the input
0.1

julia> x[1.1] # throws an error
ERROR: BoundsError: attempt to access Inclusion(0.0..1.0)
  at index [1.1]
Stacktrace:
 [1] throw_boundserror(::Inclusion{Float64,Interval{:closed,:closed,Float64}}, ::Tuple{Float64}) at ./abstractarray.jl:538
 [2] checkbounds at /Users/solver/Projects/QuasiArrays.jl/src/abstractquasiarray.jl:287 [inlined]
 [3] getindex(::Inclusion{Float64,Interval{:closed,:closed,Float64}}, ::Float64) at /Users/solver/Projects/QuasiArrays.jl/src/indices.jl:158
 [4] top-level scope at REPL[14]:1

An important usage is representing bases and function approximation, and this package contains a basic implementation of linear splines and heaviside functions. For example, we can construct splines with evenly spaced nodes via:

julia> L = LinearSpline(0:0.2:1);

julia> size(L) # uncountable (alepha-1) by 11
(ℵ₁, 6)

julia> axes(L) # The interval 0.0..1.0 by 1:6. 
(Inclusion(0.0..1.0), Base.OneTo(6))

julia> L[[0.15,0.25,0.45],1:6] # can index like an array
3×6 Array{Float64,2}:
 0.25  0.75  0.0   0.0   0.0  0.0
 0.0   0.75  0.25  0.0   0.0  0.0
 0.0   0.0   0.75  0.25  0.0  0.0

Functions in this basis are represented by a lazy multiplication by a basis and a vector of coefficients:

julia> f = L*[1,2,3,4,5,6]
QuasiArrays.ApplyQuasiArray{Float64,1,typeof(*),Tuple{Spline{1,Float64},Array{Int64,1}}}(*, (Spline{1,Float64}([0.0, 0.2, 0.4, 0.6, 0.8, 1.0]), [1, 2, 3, 4, 5, 6]))

julia> axes(f)
(Inclusion(0.0..1.0),)

julia> f[0.1]
1.5

Creating a finite element method is possible using standard array terminology. We always take the Lebesgue inner product associated with an axes, so in this case the mass matrix is just L'L. Combined with a differentiation operator allows us to form the weak Laplacian.

julia> B = L[:,2:end-1]; # drop boundary terms to impose zero Dirichlet

julia> D = Derivative(L); # Differentiation operator

julia> Δ = (D*B)'D*B # weak Laplacian
4×4 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 10.0  -5.0          
 -5.0  10.0  -5.0     
      -5.0  10.0  -5.0
           -5.0  10.0

julia> B'f # right-hand side
4-element Array{Float64,1}:
 0.4
 0.6
 0.8
 1.0

 julia> c = Δ \ B'f # coefficients of Poisson
4-element Array{Float64,1}:
 0.24               
 0.4                
 0.43999999999999995
 0.3199999999999999 

julia> u = B*c; # expand in basis

julia> u[0.1] # evaluate at 0.1
0.12