|
1 | 1 | # ContinuumArrays.jl
|
2 |
| -A package for representing continuum arrays |
| 2 | +A package for representing quasi arrays with continuous dimensions |
| 3 | + |
| 4 | + |
| 5 | +A quasi array as implemented in [QuasiArrays.jl](https://github.com/JuliaApproximation/QuasiArrays.jl) is a |
| 6 | +generalization of an array that allows non-integer indexing via general axes. This package adds support for |
| 7 | +infinite-dimensional axes, including continuous intervals. Thus it plays the same role as [InfiniteArrays.jl](https://github.com/JuliaArrays/InfiniteArrays.jl) does for standard arrays but now for quasi arrays. |
| 8 | + |
| 9 | +A simple example is the identity function on the interval `0..1`. This can be created using `Inclusion(d)`, |
| 10 | +which returns `x` if `x in d` is true, otherwise throws an error: |
| 11 | +```julia |
| 12 | +julia> using ContinuumArrays, IntervalSets |
| 13 | + |
| 14 | +julia> x = Inclusion(0..1.0) |
| 15 | +Inclusion(0.0..1.0) |
| 16 | + |
| 17 | +julia> size(x) # uncountable (aleph-1) |
| 18 | +(ℵ₁,) |
| 19 | + |
| 20 | +julia> axes(x) # axis is itself |
| 21 | +(Inclusion(0.0..1.0),) |
| 22 | + |
| 23 | +julia> x[0.1] # returns the input |
| 24 | +0.1 |
| 25 | + |
| 26 | +julia> x[1.1] # throws an error |
| 27 | +ERROR: BoundsError: attempt to access Inclusion(0.0..1.0) |
| 28 | + at index [1.1] |
| 29 | +Stacktrace: |
| 30 | + [1] throw_boundserror(::Inclusion{Float64,Interval{:closed,:closed,Float64}}, ::Tuple{Float64}) at ./abstractarray.jl:538 |
| 31 | + [2] checkbounds at /Users/solver/Projects/QuasiArrays.jl/src/abstractquasiarray.jl:287 [inlined] |
| 32 | + [3] getindex(::Inclusion{Float64,Interval{:closed,:closed,Float64}}, ::Float64) at /Users/solver/Projects/QuasiArrays.jl/src/indices.jl:158 |
| 33 | + [4] top-level scope at REPL[14]:1 |
| 34 | +``` |
| 35 | + |
| 36 | +An important usage is representing bases and function approximation, and this package contains |
| 37 | +a basic implementation of linear splines and heaviside functions. For example, we can construct splines |
| 38 | +with evenly spaced nodes via: |
| 39 | +```julia |
| 40 | +julia> L = LinearSpline(0:0.2:1); |
| 41 | + |
| 42 | +julia> size(L) # uncountable (alepha-1) by 11 |
| 43 | +(ℵ₁, 6) |
| 44 | + |
| 45 | +julia> axes(L) # The interval 0.0..1.0 by 1:6. |
| 46 | +(Inclusion(0.0..1.0), Base.OneTo(6)) |
| 47 | + |
| 48 | +julia> L[[0.15,0.25,0.45],1:6] # can index like an array |
| 49 | +3×6 Array{Float64,2}: |
| 50 | + 0.25 0.75 0.0 0.0 0.0 0.0 |
| 51 | + 0.0 0.75 0.25 0.0 0.0 0.0 |
| 52 | + 0.0 0.0 0.75 0.25 0.0 0.0 |
| 53 | +``` |
| 54 | +Functions in this basis are represented by a lazy multiplication by a basis |
| 55 | +and a vector of coefficients: |
| 56 | +```julia |
| 57 | +julia> f = L*[1,2,3,4,5,6] |
| 58 | +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])) |
| 59 | + |
| 60 | +julia> axes(f) |
| 61 | +(Inclusion(0.0..1.0),) |
| 62 | + |
| 63 | +julia> f[0.1] |
| 64 | +1.5 |
| 65 | +``` |
| 66 | + |
| 67 | +Creating a finite element method is possible using standard array terminology. |
| 68 | +We always take the Lebesgue inner product associated with an axes, so in this |
| 69 | +case the mass matrix is just `L'L`. Combined with a derivative operator allows |
| 70 | +us to form the weak Laplacian. |
| 71 | +```julia |
| 72 | +julia> B = L[:,2:end-1]; # drop boundary terms to impose zero Dirichlet |
| 73 | + |
| 74 | +julia> Δ = (D*B)'D*B # weak Laplacian |
| 75 | +4×4 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}: |
| 76 | + 10.0 -5.0 ⋅ ⋅ |
| 77 | + -5.0 10.0 -5.0 ⋅ |
| 78 | + ⋅ -5.0 10.0 -5.0 |
| 79 | + ⋅ ⋅ -5.0 10.0 |
| 80 | + |
| 81 | +julia> B'f # right-hand side |
| 82 | +4-element Array{Float64,1}: |
| 83 | + 0.4 |
| 84 | + 0.6 |
| 85 | + 0.8 |
| 86 | + 1.0 |
| 87 | + |
| 88 | + julia> c = Δ \ B'f # coefficients of Poisson |
| 89 | +4-element Array{Float64,1}: |
| 90 | + 0.24 |
| 91 | + 0.4 |
| 92 | + 0.43999999999999995 |
| 93 | + 0.3199999999999999 |
| 94 | + |
| 95 | +julia> u = B*c; # expand in basis |
| 96 | + |
| 97 | +julia> u[0.1] # evaluate at 0.1 |
| 98 | +0.12 |
| 99 | +``` |
| 100 | + |
| 101 | + |
0 commit comments