Introduction
I’m happy to announce the first release of CompactBases.jl, a library for easily approximating functions and matrix representations of linear operators in various bases where the basis functions have compact support, which leads to highly structured matrices which are non-zero only close to the diagonal (for some definition of close). At the moment, these bases are offered:
-
various finite-differences (only three-point stencils as of now ⇒ tridiagonal matrices, but more will come; local operators diagonal),
-
finite-elements, discrete-variable representation (FE-DVR, almost-block-diagonal matrices, local operators diagonal),
-
B-splines (non-orthogonal basis functions ⇒ all operators banded)
A while back, I approached @dlfivefifty with the question as to how one would structure a code such that one could easily swap out the underlying basis with minimal effort. It turned out that he had been thinking along those lines as well, and this led to QuasiArrays.jl and ContinuumArrays.jl, which consider the basis as a quasi-matrix, where the first axis is continuous on an interval, and the second axis is discrete, corresponding to the number of basis functions. A companion package for global bases, i.e. orthogonal polynomials, is OrthogonalPolynomialsQuasi.jl.
Simple usage
Install and load as usual:
julia> ]add CompactBases
julia> using CompactBases
We start out by constructing the basis. In principle, this is the only step that differs between different kinds of bases, ideally, the rest should be exactly the same.
For instance, we can easily construct simple finite-differences, FE-DVR, and B-spline bases like so:
julia> fd = FiniteDifferences(9, 0.1)
Finite differences basis {Float64} on 0.0..1.0 with 9 points spaced by Δx = 0.1
julia> fedvr = FEDVR(range(0, stop=1, length=3), 6)
FEDVR{Float64} basis with 2 elements on 0.0..1.0
julia> bsplines = BSpline(LinearKnotSet(4, 0.0, 1.0, 5))
BSpline{Float64} basis with LinearKnotSet(Float64) of order k = 4 (cubic) on 0.0..1.0 (5 intervals)
We can investigate the axes of the different bases and see that they are indeed continuous along the first dimension, and discrete along the second:
julia> axes(fd)
(Inclusion(0.0..1.0), Base.OneTo(9))
julia> axes(fedvr)
(Inclusion(0.0..1.0), Base.OneTo(11))
julia> axes(bsplines)
(Inclusion(0.0..1.0), Base.OneTo(8))
If we wish to construct the matrix representation of the second derivative in the different bases, we do this quite naturally using the language of linear operators:
julia> x = axes(fd, 1)
Inclusion(0.0..1.0)
julia> D = Derivative(x)
Derivative{Float64,IntervalSets.Interval{:closed,:closed,Float64}}(Inclusion(0.0..1.0))
(where D acts along the continuous dimension, which happens to be the same for all bases in this example)
julia> fd' * D' * D * fd
9×9 LinearAlgebra.SymTridiagonal{Float64,Array{Float64,1}}:
-200.0 100.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
100.0 -200.0 100.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 100.0 -200.0 100.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 100.0 -200.0 100.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 100.0 -200.0 100.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 100.0 -200.0 100.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 100.0 -200.0 100.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 100.0 -200.0 100.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 100.0 -200.0
julia> fedvr' * D' * D * fedvr
3×3-blocked 11×11 BlockBandedMatrices.BlockSkylineMatrix{Float64,Array{Float64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Array{Int64,1},Array{Int64,1},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
-1240.0 579.721 -62.6353 19.3726 -10.2715 │ 5.65685 │ ⋅ ⋅ ⋅ ⋅ ⋅
579.721 -385.83 138.991 -29.0091 13.668 │ -7.26304 │ ⋅ ⋅ ⋅ ⋅ ⋅
-62.6353 138.991 -174.17 98.332 -29.0091 │ 13.6985 │ ⋅ ⋅ ⋅ ⋅ ⋅
19.3726 -29.0091 98.332 -174.17 138.991 │ -44.2898 │ ⋅ ⋅ ⋅ ⋅ ⋅
-10.2715 13.668 -29.0091 138.991 -385.83 │ 409.924 │ ⋅ ⋅ ⋅ ⋅ ⋅
──────────────────────────────────────────────────────────┼───────────────┼──────────────────────────────────────────────────────────
5.65685 -7.26304 13.6985 -44.2898 409.924 │ -1240.0 │ 409.924 -44.2898 13.6985 -7.26304 5.65685
──────────────────────────────────────────────────────────┼───────────────┼──────────────────────────────────────────────────────────
⋅ ⋅ ⋅ ⋅ ⋅ │ 409.924 │ -385.83 138.991 -29.0091 13.668 -10.2715
⋅ ⋅ ⋅ ⋅ ⋅ │ -44.2898 │ 138.991 -174.17 98.332 -29.0091 19.3726
⋅ ⋅ ⋅ ⋅ ⋅ │ 13.6985 │ -29.0091 98.332 -174.17 138.991 -62.6353
⋅ ⋅ ⋅ ⋅ ⋅ │ -7.26304 │ 13.668 -29.0091 138.991 -385.83 579.721
⋅ ⋅ ⋅ ⋅ ⋅ │ 5.65685 │ -10.2715 19.3726 -62.6353 579.721 -1240.0
julia> bsplines' * D' * D * bsplines
8×8 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
-9.0 6.375 2.375 0.25 ⋅ ⋅ ⋅ ⋅
6.375 -7.5 -0.1875 1.25 0.0625 ⋅ ⋅ ⋅
2.375 -0.1875 -3.375 0.166667 0.979167 0.0416667 ⋅ ⋅
0.25 1.25 0.166667 -3.33333 0.625 0.979167 0.0625 ⋅
⋅ 0.0625 0.979167 0.625 -3.33333 0.166667 1.25 0.25
⋅ ⋅ 0.0416667 0.979167 0.166667 -3.375 -0.1875 2.375
⋅ ⋅ ⋅ 0.0625 1.25 -0.1875 -7.5 6.375
⋅ ⋅ ⋅ ⋅ 0.25 2.375 6.375 -9.0
Finally, it is very easy to evaluate the basis functions at arbitrary points, which is useful for function reconstruction (or simply plotting):
julia> xx = range(0, stop=1, length=1000)
0.0:0.001001001001001001:1.0
julia> χfd = fd[xx, :];
julia> χfedvr = fedvr[xx, :];
julia> χbsplines = bsplines[xx, :];
Many more features are available, most of which are described in the documentation.
Known limitations
-
Anything but Dirichlet-0 boundary conditions are awkward at the moment, since we have not yet figured out what the interface should look like. Being initially geared towards atoms, i.e. spherical symmetry, there are some hacks for particular situations arising therein.
-
Slow compile time.
-
Documentation in a bit of disarray, since the finite-differences, FE-DVR, and B-splines started out as separate packages, but nothing should be incorrect, per se. If you find anything, please let me know.