Interface sketch

[jagot]

I will add to the below list, as I think of more things.

Derivatives

D = Derivative(axes(B,1))

Since ∂' = -∂ (anti-symmetric), I think it's worth to point out that a weak Laplacian is -B'*D'*D*B. Furthermore, I assume that B'*(D*D*B) is a strong (?) Laplacian, i.e. the right basis is differentiated twice and is then projected on the left basis.

Higher derivatives

This is not very common in my applications, but I guess a D^n operator would be nice, instead of D*D*D*D...*D*D. For weak vs strong derivatives, I guess the user would write B'*(D^a)'*D^b*B.

Diagonal operators

This expresses the kind of operator that only locally depends on the coordinate, e.g. a potential, that when acting on a function, yields another function that's multiplied by the potential: V*f -> g(x) = V(x)*f(x). This is different from function interpolation below, although superficially similar in orthogonal bases (e.g. finite-differences); in e.g. B-splines, such operators become banded matrices. At the moment, my syntax for this is V = Matrix(x -> -inv(x), B) for e.g. the Coulomb potential, or V = Matrix(x -> -inv(x), B_1, B_2) for differing left and right bases (see below).

Function interpolation

Right now, to expand a function on a basis, I have defined B \ f::Function, which evaluates the function f on the quadrature points (or equivalently) and then does the least-square approximation with respect to the basis functions on the same points.

Properties of bases ("Internals")

Some of the below properties are not applicable to all bases, particularly not global bases, I guess, but are nevertheless handy to have for investigation and logging purposes (Which order did I use for this calculation?)

• order polynomial order of basis functions, defined as polynomial degree + 1; usually decides convergence rate of spatial derivatives/integrals. E.g. finite-differences with the [-1 2 -1] stencil for the 2nd derivative, can be viewed as 2nd order piecewise polynomials, i.e. tent functions.
• locs & weights locations and weights of quadrature points used in e.g. B-splines and finite-elements.

Variable transformations

This is something we discussed some time ago; in e.g. atomic physics, it's common to work with log(r) instead of r to increase accuracy of wavefunctions. However, one must then work with another set of equations. It would be cool to stick a variable transformation "between" the equations and the basis, such that neither the user writing the equations, nor the basis implementing basis functions and derivatives etc, need to care about the distribution of grid points. You mentioned that auto-differentiation could possibly be used for this.

Multiple bases

Sometimes it's useful to express two different functions in separate bases. If these functions are coupled via an equation system, the off-diagonal terms will need to go from one basis to another. Right now, I'm working on this for B-splines, where I use two bases of order k and k-1 to avoid spurious eigenstates of the Dirac Hamiltonian. I can do this easily, since I use the same knot set and the same set of quadrature points/weights, but for two completely different bases, e.g. B-splines and finite-differences, how would one accomplish e.g. B_1'*D*B_2, without forcing each basis implementation to know about all other possible bases?

[MikaelSlevinsky]

If you’re interested in preserving (skew-)symmetry of operators, this can be done with the ultraspherical spectral discretization. This isn’t guaranteed just because the problem is formally self-adjoint, it actually takes some work to setup. J. L. Aurentz & R. M. Slevinsky (2019), On symmetrizing the ultraspherical spectral method for self-adjoint problems, arxiv:1903.08538

[dlfivefifty]

Since ∂' = -∂ (anti-symmetric)

No! That's only true for functions that vanish at the boundary.

D^n operator would be nice

Yes, though it can be represented as ApplyQuasiMatrix(^, D, n). (Whether this is better than including n in the def of D is debatable.)

Diagonal operators

Already exist, e.g.,

x = Inclusion(0..1)
QuasiDiagonal(x) 

Function interpolation

As mentioned in email, I think the canonical version should be

x = axes(B,1)
B \ f.(x)

B \ ::Function might make sense as shorthand for above, but it's best to leave that for last. Conflating functions with quasi-vectors seems like a recipe for confusion.

[jagot]

Since ∂' = -∂ (anti-symmetric)

No! That's only true for functions that vanish at the boundary.

Sure, that's true, but for the cases when I do have vanishing boundary conditions, it is important that I negate: -B'D'*D*B, otherwise I get the wrong derivative operator matrix. I mainly thought we could mention this in the README/documentation, since I've forgotten about it more than once :)

Diagonal operators

Already exist, e.g.,

x = Inclusion(0..1)
QuasiDiagonal(x) 

And then I get the matrix realization as B'X*B, for X = QuasiDiagonal(x)? How does this tie in with QuasiArrays.jl being decoupled from IntervalSets.jl (since Inclusion is in QuasiArrays.jl)? And if I want to represent a function that is undefined/singular at one of the endpoints, Inclusion may sound counterintuitive?:

x = Inclusion(Interval{:open,:closed}(0,1))
invX = QuasiDiagonal(inv(x))

This looks slightly weird to me, but inv(QuasiDiagonal(x)) does not look much better.

Function interpolation

As mentioned in email, I think the canonical version should be

x = axes(B,1)
B \ f.(x)

But what is f.(x), if x is continuous? I.e. we need to some proxy object that results from broadcasting over an Inclusion.

[dlfivefifty]

Sure, that's true, but for the cases when I do have vanishing boundary conditions, it is important that I negate: -B'D'DB,

Sure, but it's still. not true that D'*D*B == -D^2*B because even if B vanishes, D*B does not.

How does this tie in with QuasiArrays.jl being decoupled from IntervalSets.jl

Inclusion(d) is the identity function for any domain/set d, that is:

x = Inclusion(d)
x[k] # same as `k in d ? k : throw(BoundsError())`

So it doesn't need to know anything about functions/quasi-arrays.

Since x is a quasi-vector, inv(x) makes no sense. You want:

x = Inclusion(Interval{:open,:closed}(0,1))
invX = QuasiDiagonal(inv.(x))

But inv(QuasiDiagonal(x)) should also work.

But what is f.(x), if x is continuous

It will return a lazy BroadcastQuasiVector.

[dlfivefifty]

I'm closing this. It should be separate issues or a list of checkboxes to keep track of where we are.