In my efforts to write a Julia program that solves the time-dependent Schrodinger equation for helium, I’ve gotten as far as solving the time-independent case for the ground state, with some much-appreciated help from @kahliburke. I’m interested in exploring finite difference methods on a non-uniform grid, as this could allow me to overcome the FD methods’ limitation that increasing grid density in one region requires increasing it everywhere, with corresponding costs in terms of operator matrix size and iteration time. I’ve been attempting to locate Julia packages that will help me construct the second-derivative operator on a non-uniform grid and let me access it as a matrix. For example, MethodOfLines.jl seems to support such operator constructions, but only in the context of discretizing and solving the entire PDE. Due to the structure of the problem (detailed here, forgive me for not wanting to type out all that LaTeX again), I need to be able to access the 1D second derivative operator in order to build the full matrix operator to use with KrylovKit.jl 's eigsolve(), so packages that try to handle the full discretization+solution or wrap the operator in a special type won’t work. DiffEqOperators.jl seems promising, but it’s deprecated, so I’m wary of relying on it. CompactBases.jl promises to do what I want and then some, but it’s very unmaintained. FiniteDifferences.jl runs into the special type problem. With all the various differential equations packages in Julia, I suspect there’s something out there that does what I want, but I’m having a hard time finding it.
DiffEqOperators.jl seems promising, but it’s deprecated
Why not use SciMLOperators.jl? That’s what it was deprecated in favor of.
I think you could probably construct a matrix for this using SummationByPartsOperators.jl. I wrote something similar myself a couple years ago, but I’m quite sure their code is more rigorously written and tested than mine. 
There might be an all-in-one function for this, but I am not aware of it, so you’ll probably be interested in fornberg(), which can be used to generate coefficients for an nth order derivative at zero given whatever set of points you provide.
import SummationByPartsOperators
x = [-0.2, -0.1, 0.0, 0.5, 1.0]
#generate coefficients for 2nd derivative @ x=0
SummationByPartsOperators.fornberg(x, 2)
#=
5-element Vector{Float64}:
41.66666666666667
-60.606060606060595
13.999999999999986
5.333333333333334
-0.39393939393939403
=#
Looping through smaller chunks of your grid and offsetting them such that each point is zero, you can put together a derivative matrix for an arbitrary grid.
Have you considered switching to a finite-element method (for which there are many packages)? That will allow you to use a fully unstructured mesh, with higher resolution localized wherever you want. And FEM packages all assemble an explicit sparse matrix for you.
With finite-difference methods, you will typically be restricted to Cartesian products of 1d grids, so the regions of higher resolution are less localized. (And FD methods often try to express matrices as matrix-free linear operators rather than as explicit sparse matrices—this is still compatible with iterative methods and can actually be faster, but it prevents you from using sparse-direct methods).
Sorry for the delayed response—you make a fair point. Is there a particular package you’d recommend from among Julia’s many options? As far as the matrix-free operators go, I’d be open to exploring that down the line; it seems, eg, SciMLOperators.jl will wrap an existing matrix operator.
Looking at the docs, it’s because SciMLOperators.jl won’t construct a derivative matrix operator given a grid or set of points, it’ll just construct a matrix-free operator out of a function or existing matrix.
SciMLOperators.jl won’t construct a derivative matrix operator given a grid or set of points, it’ll just construct a matrix-free operator out of a function or existing matrix.
You can concretize by convert(Abstractmatrix, operator), but you might get better results by using Krylov.jl instead of KrylovKit.jl which accepts arbitrary matrix operators instead of AbstractMatrices. It’s rare in the PDE world to actually want the full matrix.
That could be a useful route—I found SBP Operators previously, but it was hard to follow the usage/applications, and it wasn’t clear that it was actually meant for what I was trying to do. Looking closer, however, it might actually be what I need.
Right, but I still need to construct the matrix first, even if I wrap it as an operator later. This is because the full matrix has some funky structure due to the underlying equation (see the second bullet point in the linked post). That summed V_{F0} term, combined with the fact that the operator spans all instances of the index i, means I need to explicitly construct the matrix operator for the PDE. I can do all the steps of the construction as long as I can get the 1D second-derivative operator on one axis of grid points, but doing that with FD methods for a non-uniform grid is where I’m getting hung up.
We currently use Gridap.jl in my group, but its successor package GalerkinToolkit.jl is coming along nicely. Many people also like Ferrite.jl
Perhaps I’m missing something, but I don’t see how to do what I want with those packages. The nearest I can find is Gridap’s Gridap.Fields.Laplacian() function, which operates on a Gridap Field. At that point, it’s not at all clear what a Field is or how to make one, such that I can get my desired 1D Laplacian operator out.
Regarding your mention of Cartesian products of 1D grids, that’s the whole idea, I think; I make the 1D Laplacian, add the centrifugal barrier term at each grid point ( l(l+1)/r^2), then use the Kronecker product to turn it into a 2D operator.
Switching to finite-element methods (with any software package) requires you learn a little bit of a different approach and terminology. If you start with a differential operator like the Schrödinger operator \hat{H} = -\nabla^2 + V, you translate it to something called a “weak form” a(u,v) = \langle u, \hat{H} v \rangle = \int_{\Omega} (\nabla u \cdot \nabla v + V(x) u^* v) d\Omega (this is very similar to a scalar Helmholtz equation, for which there are some examples in the Gridap manual). Given the weak form a, Gridap (or any FEM package) provides tools to assemble the corresponding Hermitian matrix A that you can use to evolve/solve the Schrödinger equation (in a “Galerkin” formulation).
This is the analogue of a finite-difference matrix approximation of \hat{H}, but notice that the weak form has “integrated by parts” so that it only uses gradient operators and not an explicit Laplacian operator.
There are lots of sources (many, many books etc.) to learn about these concepts. It requires a bit of investment of effort, but IMO it is well-worth it to grasp a new class of approaches that is much more flexible about discretization than finite-difference methods, and has a huge array of flexible software packages (in multiple languages) that is easily applicable to a lot of problems (once you learn the basic concepts).