Correct first derivatives for non-uniform StaggeredFiniteDifferences

Author: jagot

src/fd_derivatives.jl

@@ -1,28 +1,43 @@
 # * Derivatives
 # ** Three-point stencils
+
+# α = super-/subdiagonal of second derivative
+# β = -1/2 diagonal of second derivative
+# γ = diagonal of first derivative
+# δ = super-/subdiagonal of first derivative
+
 α(::FiniteDifferences{T}, ::Integer) where T = one(T)
 β(::FiniteDifferences{T}, ::Integer) where T = one(T)
 γ(::FiniteDifferences{T}, ::Integer) where T = zero(T)
+δ(::FiniteDifferences{T}, ::Integer) where T = one(T)
 
 α(B::StaggeredFiniteDifferences,    j::Integer)         = B.α[j]
-β(B::StaggeredFiniteDifferences{T}, j::Integer) where T = B.β[j]
+β(B::StaggeredFiniteDifferences,    j::Integer)         = B.β[j]
 γ( ::StaggeredFiniteDifferences{T},  ::Integer) where T = zero(T)
+δ(B::StaggeredFiniteDifferences,    j::Integer)         = B.δ[j]
 
-α(B::StaggeredFiniteDifferences)            = B.α
-β(B::StaggeredFiniteDifferences)            = B.β
+α(B::StaggeredFiniteDifferences) = B.α
+β(B::StaggeredFiniteDifferences) = B.β
 γ(B::StaggeredFiniteDifferences{T}) where T = zeros(T, length(B.r))
+δ(B::StaggeredFiniteDifferences) = B.δ
 
 α( ::ImplicitFiniteDifferences{T},  ::Integer) where T = one(T)
 β(B::ImplicitFiniteDifferences{T}, j::Integer) where T = one(T) + (j == 1 ? B.δβ₁ : zero(T))
 γ(B::ImplicitFiniteDifferences{T}, j::Integer) where T = (j == 1 ? B.λ : zero(T))
+δ( ::ImplicitFiniteDifferences{T},  ::Integer) where T = one(T)
 
 α(B::AbstractFiniteDifferences) = α.(Ref(B), B.j[1:end-1])
 β(B::AbstractFiniteDifferences) = β.(Ref(B), B.j)
 γ(B::AbstractFiniteDifferences) = γ.(Ref(B), B.j)
+δ(B::AbstractFiniteDifferences) = δ.(Ref(B), B.j[1:end-1])
 
-function α(B̃::RestrictedFiniteDifferences)
-    j = indices(B̃,2)
-    α(parent(B̃))[j[1]:(j[end]-1)]
+for f in [:α, :δ]
+    @eval begin
+        function $f(B̃::RestrictedFiniteDifferences)
+            j = indices(B̃,2)
+            $f(parent(B̃))[j[1]:(j[end]-1)]
+        end
+    end
 end
 
 for f in [:β, :γ]
@@ -33,7 +48,7 @@ for f in [:β, :γ]
     end
 end
 
-for f in [:α, :β, :γ]
+for f in [:α, :β, :γ, :δ]
     @eval begin
         function $f(Ã::BasisOrRestricted{<:AbstractFiniteDifferences},
                     B̃::BasisOrRestricted{<:AbstractFiniteDifferences},
@@ -84,10 +99,10 @@ end
 
         # Central difference approximation
         μ = 2step(B)
-        a = α(B)/μ
-        dest.dl .= -a
+        d = δ(B)/μ
+        dest.dl .= -d
         dest.d .= γ(B)/μ
-        dest.du .= a
+        dest.du .= d
     end
     dest::BandedMatrix{T} -> begin
         A = parent(Ac)

src/finite_differences.jl

@@ -242,13 +242,15 @@ r = 0. Supports non-uniform grids, c.f.
   10118–10125. http://dx.doi.org/10.1021/jp992144
 
 """
-struct StaggeredFiniteDifferences{T,V<:AbstractVector{T},A,B} <: AbstractFiniteDifferences{T,Int}
+struct StaggeredFiniteDifferences{T,V<:AbstractVector{T},A,B,D} <: AbstractFiniteDifferences{T,Int}
     r::V
     Z::T
     δβ₁::T # Correction used for bare Coulomb potentials, Eq. (22) Schafer2009
     α::A
     β::B
+    δ::D
 
+    # Uniform case
     function StaggeredFiniteDifferences(r::V, Z=one(T),
                                         δβ₁=schafer_corner_fix(local_step(r,1), Z)) where {T,V<:AbstractRange{T}}
         issorted(r) ||
@@ -263,9 +265,10 @@ struct StaggeredFiniteDifferences{T,V<:AbstractVector{T},A,B} <: AbstractFiniteD
 
         β[1] += δβ₁ * step(r)^2
 
-        new{T,V,typeof(α),typeof(β)}(r, T(Z), T(δβ₁), α, β)
+        new{T,V,typeof(α),typeof(β),typeof(α)}(r, T(Z), T(δβ₁), α, β, α)
     end
 
+    # Arbitrary case
     function StaggeredFiniteDifferences(r::V, Z=one(T),
                                         δβ₁=schafer_corner_fix(local_step(r,1), Z)) where {T,V<:AbstractVector{T}}
         issorted(r) ||
@@ -274,6 +277,8 @@ struct StaggeredFiniteDifferences{T,V<:AbstractVector{T},A,B} <: AbstractFiniteD
         N = length(r)
         α = zeros(T, N-1)
         β = zeros(T, N)
+        γ = zeros(T, N)
+        δ = zeros(T, N-1)
 
         r̃ = vcat(r, 2r[end]-r[end-1])
         a = 2r[1]-r[2]
@@ -283,19 +288,21 @@ struct StaggeredFiniteDifferences{T,V<:AbstractVector{T},A,B} <: AbstractFiniteD
             if j < N
                 d = r̃[j+2]
                 # Eq. (A13) Krause 1999
-                α[j] = 2/((c-b)*√((d-b)*(c-a)))*((b+c)/2)^2/(b*c)
+                δ[j] = 2/(√((d-b)*(c-a)))*((b+c)/2)^2/(b*c)
+                α[j] = δ[j]/(c-b)
             end
 
             # Eq. (A14) Krause 1999
-            β[j] = 1/(c-a)*(1/(c-b)*((c+b)/2b)^2 +
-                            1/(b-a)*((b+a)/2b)^2)
+            fp = ((c+b)/2b)^2
+            fm = ((b+a)/2b)^2
+            β[j] = 1/(c-a)*(1/(c-b)*fp + 1/(b-a)*fm)
 
             a,b,c = b,c,d
         end
 
         β[1] += δβ₁
 
-        new{T,V,typeof(α),typeof(β)}(r, T(Z), T(δβ₁), α, β)
+        new{T,V,typeof(α),typeof(β),typeof(δ)}(r, T(Z), T(δβ₁), α, β, δ)
     end
 
     # Constructor of uniform grid

test/fd/derivatives.jl

@@ -128,12 +128,14 @@ end
         @test isapprox(ph, order, atol=0.03) || ph > order
     end
     @testset "kind = StaggeredFiniteDifferences" begin
-        let (order,B) = (2,StaggeredFiniteDifferences)
-            @info "$B derivative accuracy"
-            dd = b-a
+        dd = b-a
+        @testset "$(label)" for (order,label,fun) = [(2.0, "Uniform", N -> (N,dd/(N+1))),
+                                                     (2.5, "Non-uniform", N -> (dd*(1 .- cos.(π*range(0,stop=1,length=N+2)[2:end-1]))/2,))]
+            @info "$(label) StaggeredFiniteDifferences derivative accuracy"
             hs,ϵg,ϵh,ϵh′,pg,ph,ph′ = compute_derivative_errors(Ns, x -> f(x-dd/2), x -> g(x-dd/2), x -> h(x-dd/2), 1) do N
-                Δx = dd/(N+1)
-                StaggeredFiniteDifferences(N, Δx, 1.0, 0.0),Δx
+                R = StaggeredFiniteDifferences(fun(N)..., 1.0, 0.0)
+                Δx = maximum(diff(R.r))
+                R, Δx
             end
 
             @test isapprox(pg, order, atol=0.03) || pg > order