# Registering Finite Difference Derivatives with Symbolics.jl

**URL:** <https://discourse.julialang.org/t/registering-finite-difference-derivatives-with-symbolics-jl/116636>\
**Category:** Modelling & Simulations\
**Tags:** question, package, modelingtoolkit, symbolics\
**Created:** [July 4, 2024, 6:10pm UTC](https://discourse.julialang.org/t/registering-finite-difference-derivatives-with-symbolics-jl/116636 "2024-07-04T18:10:45Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [July 4, 2024, 6:10pm UTC](https://discourse.julialang.org/t/registering-finite-difference-derivatives-with-symbolics-jl/116636/1 "2024-07-04T18:10:45Z")

</div>

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.

```Julia
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`.

```julia
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.

```julia
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.

```

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [July 4, 2024, 6:19pm UTC](https://discourse.julialang.org/t/registering-finite-difference-derivatives-with-symbolics-jl/116636/2 "2024-07-04T18:19:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [July 6, 2024, 4:31pm UTC](https://discourse.julialang.org/t/registering-finite-difference-derivatives-with-symbolics-jl/116636/3 "2024-07-06T16:31:18Z")

</div>

A [comment](https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/50#issuecomment-1397363491) on an [issue](https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/50) 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`.

```julia
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)

```