Apply parity boundary conditions

docs/src/FDM.md

@@ -21,6 +21,16 @@ The eigenvalue ``E`` is an approximation of the exact energy and the eigenvector
 ```
 A uniform grid spacing is used, ``r_{i+1} = r_{i} + \Delta r``. See the API reference for the expression of the matrix ``\pmb{H}``.
 
+## Boundary conditions
+
+For central differences, the wavefunction is extended across the origin with
+
+```math
+\psi(-r)=(-1)^l\psi(r).
+```
+
+Thus even ``l`` uses ``\psi'(0)=0``, while odd ``l`` uses ``\psi(0)=0``. The origin value for even ``l`` is eliminated with a second-order one-sided derivative, so the radial grid can continue to exclude ``r=0``.
+
 ## Usage
 
 Run the following code before each use.

src/FDM.jl

@@ -52,11 +52,24 @@ function matrix(o::RestEnergy, method::FiniteDifferenceMethod)
 end
 
 function matrix(o::Kinetic, method::FiniteDifferenceMethod)
-  D  = FiniteDifferenceMatrices.fdmatrix(Int64(length(method.R)), n=1, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
-  D² = FiniteDifferenceMatrices.fdmatrix(Int64(length(method.R)), n=2, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
+  D, D² = _derivative_matrices(method)
   return -o.hbar^2/2/o.m * (D² + SparseArrays.spdiagm(2 ./ method.R) * D - method.l*(method.l+1) * SparseArrays.spdiagm(1 ./ method.R .^ 2))
 end
 
+function _derivative_matrices(method::FiniteDifferenceMethod)
+  n = Int64(length(method.R))
+  D = FiniteDifferenceMatrices.fdmatrix(n, n=1, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
+  D² = FiniteDifferenceMatrices.fdmatrix(n, n=2, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
+
+  if method.direction == :c && iseven(method.l) && 1 < n
+    D[1,1] = -2 / (3method.Δr)
+    D[1,2] = 2 / (3method.Δr)
+    D²[1,1] = -2 / (3method.Δr^2)
+    D²[1,2] = 2 / (3method.Δr^2)
+  end
+  return D, D²
+end
+
 function matrix(o::PotentialTerm, method::FiniteDifferenceMethod)
   return SparseArrays.spdiagm([V(o, r) for r in method.R])
 end

test/FDM.jl

@@ -1,5 +1,19 @@
 @testset "FDM.jl" begin
 
+  @testset "boundary condition" begin
+    even = FiniteDifferenceMethod(Δr=0.1, rₘₐₓ=0.3, l=0)
+    odd = FiniteDifferenceMethod(Δr=0.1, rₘₐₓ=0.3, l=1)
+    D_even, D²_even = TwoBody._derivative_matrices(even)
+    D_odd, D²_odd = TwoBody._derivative_matrices(odd)
+
+    even_wavefunction = 1 .+ even.R .^ 2
+    odd_wavefunction = collect(odd.R)
+    @test (D_even * even_wavefunction)[1] ≈ 2even.R[1]
+    @test (D²_even * even_wavefunction)[1] ≈ 2
+    @test (D_odd * odd_wavefunction)[1] ≈ 1
+    @test (D²_odd * odd_wavefunction)[1] ≈ 0 atol=1e-12
+  end
+
   # Testing Results
 
   H = Hamiltonian(