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 nonzero only close to the diagonal (for some definition of close). At the moment, these bases are offered:

various finitedifferences (only threepoint stencils as of now ⇒ tridiagonal matrices, but more will come; local operators diagonal),

finiteelements, discretevariable representation (FEDVR, almostblockdiagonal matrices, local operators diagonal),

Bsplines (nonorthogonal 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 quasimatrix, 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 finitedifferences, FEDVR, and Bspline 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×3blocked 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 Dirichlet0 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 finitedifferences, FEDVR, and Bsplines started out as separate packages, but nothing should be incorrect, per se. If you find anything, please let me know.