CompactBases.jl
CompactBases.jl is a package for various bases for function approximation with compact support in the framework of ContinuumArrays
Example usage
julia> a,b = 0.0,1.0 # Extents
(0.0, 1.0)
julia> N = 13 # Number of nodes
13
julia> k = 5 # Order of B-splines
5
julia> x = range(a, stop=b, length=1000)
0.0:0.001001001001001001:1.0
julia> Δx = (b-a)/(N+1) # Grid spacing
0.07142857142857142
julia> # Standard, uniform finite-differences
fd = FiniteDifferences(N, Δx)
Finite differences basis {Float64} on 0.0..1.0 with 13 points spaced by Δx = 0.07142857142857142
julia> # Staggered, uniform finite-differences
sfd_uni = StaggeredFiniteDifferences(N, Δx)
Staggered finite differences basis {Float64} on 0.0..0.9642857142857144 with 13 points spaced by ρ = 0.07142857142857142
julia> # Staggered, non-uniform finite-differences
sfd_nonuni = StaggeredFiniteDifferences(0.01, 0.5, 0.1, b)
Staggered finite differences basis {Float64} on 0.0..1.1112312795594823 with 39 points
julia> # Finite-element boundaries
tf = range(a, stop=b, length=N+2)
0.0:0.07142857142857142:1.0
julia> using LazyArrays, FillArrays
julia> # We can vary the polynomial order in each element
forder = Vcat(7, Fill(4,length(tf)-2))
14-element ApplyArray{Int64,1,typeof(vcat),Tuple{Int64,Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}}:
7
4
4
4
4
4
4
4
4
4
4
4
4
4
julia> # By indexing the second dimension, we can implement Dirichlet0
# boundary conditions.
fem = FEDVR(tf, forder)[:,2:end-1]
FEDVR{Float64} basis with 14 elements on 0.0..1.0, restricted to elements 1:14, basis functions 2..45 ⊂ 1..46
julia> tb = ExpKnotSet(k, -2.0, log10(b), N+1)
23-element ExpKnotSet{5,5,5,Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Array{Float64,1}}:
0.0
0.0
0.0
0.0
0.0
0.01
0.014251026703029978
0.020309176209047358
0.028942661247167503
0.04124626382901352
0.05878016072274912
0.0837677640068292
0.11937766417144363
0.17012542798525887
0.24244620170823283
0.3455107294592219
0.4923882631706739
0.7017038286703828
1.0
1.0
1.0
1.0
1.0
julia> splines = BSpline(tb)[:,2:end-1]
BSpline{Float64} basis with ExpKnotSet(Float64) of on order k = 5 (quartic) on 0,0.01..1.0 (14 intervals), restricted to basis functions 2..17 ⊂ 1..18
Reference
CompactBases.AbstractKnotSet
CompactBases.ArbitraryKnotSet
CompactBases.ArbitraryKnotSet
CompactBases.BSpline
CompactBases.BSpline
CompactBases.BSpline
CompactBases.Density
CompactBases.DiagonalOperator
CompactBases.DiagonalOperator
CompactBases.ExpKnotSet
CompactBases.ExpKnotSet
CompactBases.FunctionProduct
CompactBases.FunctionProduct
CompactBases.FunctionProduct
CompactBases.LinearKnotSet
CompactBases.LinearKnotSet
CompactBases.LinearOperator
CompactBases.LinearOperator
CompactBases.ShiftAndInvert
CompactBases.ShiftAndInvert
CompactBases.ShiftAndInvert
CompactBases.ShiftAndInvert
CompactBases.StaggeredFiniteDifferences
CompactBases.UnitVector
Base.copyto!
Base.copyto!
Base.copyto!
Base.first
Base.first
Base.getindex
Base.getindex
Base.last
Base.last
Base.length
Base.length
CompactBases.assert_multiplicities
CompactBases.assert_multiplicities
CompactBases.centers
CompactBases.change_interval!
CompactBases.deBoor
CompactBases.find_interval
CompactBases.find_interval
CompactBases.findelement
CompactBases.lagrangeder!
CompactBases.leftmultiplicity
CompactBases.leftmultiplicity
CompactBases.lgwt
CompactBases.lgwt
CompactBases.local_step
CompactBases.metric
CompactBases.metric_shape
CompactBases.nonempty_intervals
CompactBases.nonempty_intervals
CompactBases.num_quadrature_points
CompactBases.num_quadrature_points
CompactBases.numfunctions
CompactBases.numfunctions
CompactBases.numintervals
CompactBases.numintervals
CompactBases.order
CompactBases.order
CompactBases.rightmultiplicity
CompactBases.rightmultiplicity
CompactBases.within_interval
CompactBases.within_interval
CompactBases.within_support
CompactBases.within_support
LinearAlgebra.mul!
CompactBases.@materialize
CompactBases.LinearOperator
— TypeLinearOperator(A, B⁻¹, temp)
A helper object used to apply the action of a linear operator, whose matrix representation is given by A
, in the basis whose metric inverse is B⁻¹
. For orthogonal bases (such as e.g. FiniteDifferences
and FEDVR
), only multiplication by A
is necessary, but non-orthonal bases (such as BSpline
) necessitates the application of the metric inverse as well.
CompactBases.LinearOperator
— MethodLinearOperator(A, R)
Construct the linear operator whose matrix representation in the basis R
is A
, automatically deducing if application of the metric inverse is also necessary.
CompactBases.ShiftAndInvert
— TypeShiftAndInvert(A⁻¹, B, temp)
Represents a shifted-and-inverted operator, suitable for Krylov iterations when targetting an interior point of the eigenspectrum. A⁻¹
is a factorization of the shifted operator and B
is the metric matrix.
CompactBases.ShiftAndInvert
— TypeShiftAndInvert(A, R[, σ=0])
Construct the shifted-and-inverted operator corresponding to the eigenproblem $(\mat{A}-σ\mat{B})\vec{x} = \lambda\mat{B}\vec{x}$, where $\mat{B}$ is the suitable operator metric, depending on the AbstractQuasiMatrix
R
employed.
CompactBases.ShiftAndInvert
— TypeShiftAndInvert(A, B[, σ=0])
Construct the shifted-and-inverted operator corresponding to the eigenproblem $(\mat{A}-σ\mat{B})\vec{x} = \lambda\mat{B}\vec{x}$.
CompactBases.ShiftAndInvert
— TypeShiftAndInvert(A, ::UniformScaling[, σ=0])
Construct the shifted-and-inverted operator corresponding to the eigenproblem $(\mat{A}-σ\mat{I})\vec{x} = \lambda\vec{x}$.
CompactBases.StaggeredFiniteDifferences
— TypeStaggeredFiniteDifferences
Staggered finite differences with variationally derived three-points stencils for the case where there is Dirichlet0 boundary condition at r = 0. Supports non-uniform grids, c.f.
- Krause, J. L., & Schafer, K. J. (1999). Control of THz Emission from Stark Wave Packets. The Journal of Physical Chemistry A, 103(49), 10118–10125. http://dx.doi.org/10.1021/jp992144
CompactBases.centers
— Methodcenters(B)
Return the locations of the mass centers of all basis functions of B
; for orthogonal bases such as finite-differences and FE-DVR, this is simply locs
, i.e. the location of the quadrature nodes.
CompactBases.metric
— Methodmetric(B)
Returns the metric or mass matrix of the basis B
, equivalent to S=B'B
.
CompactBases.metric_shape
— Methodmetric_shape(B)
Returns the shape of the metric or mass matrix of the basis B
, i.e. Diagonal
, BandedMatrix
, etc.
CompactBases.nonempty_intervals
— Methodnonempty_intervals(t)
Return the indices of all intervals of the knot set t
that are non-empty.
Examples
julia> nonempty_intervals(ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3))
4-element Array{Int64,1}:
1
3
4
5
CompactBases.numfunctions
— Methodnumfunctions(t)
Returns the number of basis functions generated by knot set t
.
Examples
julia> numfunctions(LinearKnotSet(3, 0, 1, 2))
4
julia> numfunctions(LinearKnotSet(5, 0, 1, 2))
6
CompactBases.numintervals
— Methodnuminterval(t)
Returns the number of intervals generated by the knot set t
.
Examples
julia> numintervals(LinearKnotSet(3, 0, 1, 2))
2
CompactBases.order
— Methodorder(t)
Returns the order k
of the knot set t
.
CompactBases.UnitVector
— TypeUnitVector{T}(N, k)
Helper vector type of length N
where the k
th element is one(T)
and all the others zero(T)
.
Base.first
— Methodfirst(t)
Return the first knot of t
.
Base.getindex
— Methodgetindex(t, i)
Return the i
th knot of the knot set t
, accounting for the multiplicities of the endpoints.
Examples
julia> LinearKnotSet(3, 0, 1, 3)[2]
0.0
julia> LinearKnotSet(3, 0, 1, 3, 1, 1)[2]
0.3333333333333333
Base.last
— Methodlast(t)
Return the last knot of t
.
Base.length
— Methodlength(t)
Return the number of knots of t
.
Examples
julia> length(LinearKnotSet(3, 0, 1, 3))
8
julia> length(LinearKnotSet(3, 0, 1, 3, 1, 1))
4
CompactBases.assert_multiplicities
— Methodassert_multiplicities(k,ml,mr,t)
Assert that the multiplicities at the endpoints, ml
and mr
, respectively, are consistent with the order k
. Also check that the amount of knots in the knot set t
are enough to support the requested order k
.
CompactBases.change_interval!
— Methodchange_interval!(xs, ws, x, w[, a=0, b=1, γ=1])
Transform the Gaußian quadrature roots x
and weights w
on the elementary interval [-1,1]
to the interval [γ*a,γ*b]
and store the result in xs
and ws
, respectively. γ
is an optional root of unity, used to complex-rotate the roots (but not the weights).
CompactBases.deBoor
— MethoddeBoor(t, c, x[, i[, m=0]])
Evaluate the spline given by the knot set t
and the set of control points c
at x
using de Boor's algorithm. i
is the index of the knot interval containing x
. If m≠0
, calculate the m
th derivative at x
instead.
CompactBases.find_interval
— Methodfind_interval(t, x[, i=ml])
Find the interval in the knot set t
that includes x
, starting from interval i
(which by default is the first non-zero interval of the knot set). The search complexity is linear, but by storing the result and using it as starting point for the next call to find_interval
, the knot set need only be traversed once.
CompactBases.findelement
— Functionfindelement(B::FEDVR, k[, i=1, m=k])
Find the finite-element of B
that contains the k
th basis function, optionally starting the search from element i
and basis function m
of that element.
CompactBases.lagrangeder!
— Methodlagrangeder!(xⁱ, m, L′)
Calculate the derivative of the Lagrange interpolating polynomial Lⁱₘ(x), given the roots xⁱ
, at the roots, and storing the result in L′
.
∂ₓ Lⁱₘ(xⁱₘ,) = (xⁱₘ-xⁱₘ,)⁻¹ ∏(k≠m,m′) (xⁱₘ,-xⁱₖ)/(xⁱₘ-xⁱₖ), m≠m′, [δ(m,n) - δ(m,1)]/2wⁱₘ, m=m′
Eq. (20) Rescigno2000
CompactBases.leftmultiplicity
— Methodleftmultiplicity(t)
Return the multiplicity of the left endpoint.
CompactBases.lgwt
— Methodlgwt(t, N) -> (x,w)
Generate the N
Gauß–Legendre quadrature roots x
and associated weights w
, with respect to the B-spline basis generated by the knot set t
.
Examples
julia> CompactBases.lgwt(LinearKnotSet(2, 0, 1, 3), 2)
([0.0704416, 0.262892, 0.403775, 0.596225, 0.737108, 0.929558], [0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667])
julia> CompactBases.lgwt(ExpKnotSet(2, -4, 2, 7), 2)
([2.11325e-5, 7.88675e-5, 0.000290192, 0.000809808, 0.00290192, 0.00809808, 0.0290192, 0.0809808, 0.290192, 0.809808, 2.90192, 8.09808, 29.0192, 80.9808], [5.0e-5, 5.0e-5, 0.00045, 0.00045, 0.0045, 0.0045, 0.045, 0.045, 0.45, 0.45, 4.5, 4.5, 45.0, 45.0])
CompactBases.local_step
— Methodlocal_step(B, i)
The step size around grid point i
. For uniform grids, this is equivalent to step
.
CompactBases.num_quadrature_points
— Methodnum_quadrature_points(k, k′)
The number of quadrature points needed to exactly compute the matrix elements of an operator of polynomial order k′
with respect to a basis of order k
.
CompactBases.rightmultiplicity
— Methodrightmultiplicity(t)
Return the multiplicity of the right endpoint.
CompactBases.within_interval
— Methodwithin_interval(x, interval)
Return the indices of the elements of x
that lie within the given interval
.
CompactBases.within_support
— Methodwithin_support(x, t, j)
Return the indices of the elements of x
that lie withing the compact support of the j
th basis function (enumerated 1..n
), given the knot set t
. For each index of x
that is covered, the index k
of the interval within which x[i]
falls is also returned.
CompactBases.@materialize
— Macro@materialize function op(args...)
This macro simplifies the setup of a few functions necessary for the materialization of Applied
objects:
ApplyStyle
, used to ensure dispatch of the applied object to the routines belowcopyto!(dest::DestType, applied_obj::Applied{...,op})
performs the actual materialization ofapplied_obj
into the destination object which has been generated bysimilar
which usually returns a suitable matrixmaterialize
which makes use of the above functions
Example
@materialize function *(Ac::MyAdjointBasis,
O::MyOperator,
B::MyBasis)
MyApplyStyle # An instance of this type will be returned by ApplyStyle
T -> begin # generates similar
A = parent(Ac)
parent(A) == parent(B) ||
throw(ArgumentError("Incompatible bases"))
# There may be different matrices best representing different
# situations:
if ...
Diagonal(Vector{T}(undef, size(B,1)))
else
Tridiagonal(Vector{T}(undef, size(B,1)-1),
Vector{T}(undef, size(B,1)),
Vector{T}(undef, size(B,1)-1))
end
end
dest::Diagonal{T} -> begin # generate copyto!(dest::Diagonal{T}, ...) where T
dest.diag .= 1
end
dest::Tridiagonal{T} -> begin # generate copyto!(dest::Tridiagonal{T}, ...) where T
dest.dl .= -2
dest.ev .= 1
dest.du .= 3
end
end