Registering Finite Difference Derivatives with Symbolics.jl

Specific Domains Modelling & Simulations
, , ,
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)