Registering Finite Difference Derivatives with Symbolics.jl

Specific Domains Modelling & Simulations
, , ,
1

I’m trying to take the gradient of the Power Law Cutoff Potential using Symbolics and ModelingToolkit. My end goal is to define an ODESystem with this gradient. The Power Law Cutoff Potential equation is shown below. Here, \Gamma is the gamma function, and \gamma is the lower incomplete gamma function.

How can I register numerical derivative functions, so when I take the gradient there is only one independent variable?

\Phi = \frac{G \alpha m \gamma\left(1.5 - \frac{\alpha}{2}, \frac{x^{2} + y^{2} + z^{2}}{r_{c}^{2}}\right)}{2 \sqrt{x^{2} + y^{2} + z^{2}} \Gamma\left(2.5 - \frac{\alpha}{2}\right)} - \frac{3 G m \gamma\left(1.5 - \frac{\alpha}{2}, \frac{x^{2} + y^{2} + z^{2}}{r_{c}^{2}}\right)}{2 \sqrt{x^{2} + y^{2} + z^{2}} \Gamma\left(2.5 - \frac{\alpha}{2}\right)} + \frac{G m \gamma\left(1 - \frac{\alpha}{2}, \frac{x^{2} + y^{2} + z^{2}}{r_{c}^{2}}\right)}{r_{c} \Gamma\left(1.5 - \frac{\alpha}{2}\right)}

I can define this scalar potential with the following code.

using SpecialFunctions

"""
Computes the lower incomplete gamma function.
"""
function lowergamma(a, x)
    p, q = SpecialFunctions.gamma_inc(a, x)
    return p * gamma(a)
end

using Symbolics, ModelingToolkit

@register_symbolic lowergamma(x::AbstractFloat, y::AbstractFloat)
@register_symbolic SpecialFunctions.gamma(x::AbstractFloat)

@variables t
u = @variables x(t) y(t) z(t)
p = @parameters G m a α c

value = G * α * m * lowergamma(3 // 2 - α // 2, (x^2 + y^2 + z^2) / c^2) /
        (2 * sqrt(x^2 + y^2 + z^2) * gamma(5 // 2 - α // 2)) -
        3 * G * m * lowergamma(3 // 2 - α // 2, (x^2 + y^2 + z^2) / c^2) /
        (2 * sqrt(x^2 + y^2 + z^2) * gamma(5 // 2 - α // 2)) +
        G * m * lowergamma(1 - α // 2, (x^2 + y^2 + z^2) / c^2) /
        (c * gamma(3 // 2 - α // 2))

I then want to register the partial derivatives of lowergamma with respect to the two arguments a and x.

using FiniteDiff

function Symbolics.derivative(::typeof(lowergamma), args::NTuple{2, Any}, ::Val{1})
    x, y = args
    return FiniteDiff.finite_difference_derivative(x -> lowergamma(x, y), x)
end

function Symbolics.derivative(::typeof(lowergamma), args::NTuple{2, Any}, ::Val{2})
    x, y = args
    return FiniteDiff.finite_difference_derivative(y -> lowergamma(x, y), y)
end

Finally, I try to take the gradient, and get the following error.

dvs = map(Differential(t), u)

eqs = dvs .~ Symbolics.gradient(value, u)
ERROR: MethodError: no method matching default_relstep(::Val{:central}, ::Type{Any})

Closest candidates are:
  default_relstep(::Val{fdtype}, ::Type{T}) where {fdtype, T<:Number}
   @ FiniteDiff ~/.julia/packages/FiniteDiff/KKEkv/src/epsilons.jl:26
  default_relstep(::Type{V}, ::Any) where V
   @ FiniteDiff ~/.julia/packages/FiniteDiff/KKEkv/src/epsilons.jl:25

Stacktrace:
  [1] finite_difference_derivative(f::Function, x::SymbolicUtils.BasicSymbolic{…}, fdtype::Val{…}, returntype::Type, fx::Nothing, epsilon::Nothing)
    @ FiniteDiff ~/.julia/packages/FiniteDiff/KKEkv/src/derivatives.jl:103
  [2] derivative(::typeof(lowergamma), args::Tuple{SymbolicUtils.BasicSymbolic{…}, SymbolicUtils.BasicSymbolic{…}}, ::Val{2})
    @ Main ./REPL[11]:3
  [3] derivative_idx(O::SymbolicUtils.BasicSymbolic{Real}, idx::Int64)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:337
  [4] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::SymbolicUtils.BasicSymbolic{Real})
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:253
  [5] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::SymbolicUtils.BasicSymbolic{Real}) (repeats 2 times)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:245
  [6] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::Nothing)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:245
  [7] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:171
  [8] gradient(O::Num, vars::Vector{Num}; simplify::Bool)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:443
  [9] gradient(O::Num, vars::Vector{Num})
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:442
 [10] top-level scope
    @ REPL[14]:1
Some type information was truncated. Use `show(err)` to see complete types.
2

Note that if I don’t register the derivatives, the gradient will not evaluate the partial derivatives of the lower incomplete gamma function with respect to x, y, and z, and the ODESystem constructor will fail due to the appearance of multiple independent variables.

3

A comment on an issue on the SpecialFunctions project pointed to a partial derivative expression that uses HypergeometricFunctions.M. By registering that function, and SpecialFunctions.gamma as symbolic functions which only accept AbstractFloat inputs, I was able to construct an ODESystem.

using SpecialFunctions
using HypergeometricFunctions
using Symbolics

"""
Computes the lower incomplete gamma function.
"""
function lowergamma(a, x)
    p, q = SpecialFunctions.gamma_inc(a, x)
    return p * gamma(a)
end

@register_symbolic lowergamma(x::AbstractFloat, y::AbstractFloat)
@register_symbolic SpecialFunctions.gamma(x::AbstractFloat)
@register_symbolic HypergeometricFunctions.M(
    a::AbstractFloat, b::AbstractFloat, x::AbstractFloat)

function Symbolics.derivative(::typeof(lowergamma), args::NTuple{2, Any}, ::Val{1})
    a, x = args
    return (x^a / a) * HypergeometricFunctions.M(a, a + 1, -x) # https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/50#issuecomment-1397363491
end

function Symbolics.derivative(::typeof(lowergamma), args::NTuple{2, Any}, ::Val{2})
    a, x = args
    return -x^(a - 1) * exp(-x) # https://en.wikipedia.org/wiki/Incomplete_gamma_function#Derivatives
end

using ModelingToolkit

model = begin
    @variables t
    u = @variables x(t) y(t) z(t)
    p = @parameters G m a α c

    G * α * m * lowergamma(3 // 2 - α // 2, (x^2 + y^2 + z^2) / c^2) /
    (2 * sqrt(x^2 + y^2 + z^2) * SpecialFunctions.gamma(5 // 2 - α // 2)) -
    3 * G * m * lowergamma(3 // 2 - α // 2, (x^2 + y^2 + z^2) / c^2) /
    (2 * sqrt(x^2 + y^2 + z^2) * SpecialFunctions.gamma(5 // 2 - α // 2)) +
    G * m * lowergamma(1 - α // 2, (x^2 + y^2 + z^2) / c^2) /
    (c * SpecialFunctions.gamma(3 // 2 - α // 2))
end

Symbolics.gradient(model, u)