Something specific for pFq([a,a],[a+1,a+1],-x) ?

[lrnv]

Hey,

I do not know deeply what are the algorithms implemente din this package, but I needed the derivative of gamma_inc for my gradient decent and implemented it there, using a call to

pFq([a,a],[a+1,a+1],-x)

where x would always be a positive real, and a could theoretically be any complex, but I only need positives reals too. The serie collapses to the simple :

$$2F_2(a,a,a+1,a+1,z) = \sum{n = 0}^{\infty} \left(\frac{a}{a+n}\right)^2 \frac{z^n}{n!},$$

or even :

$$2F_2(a,a,a+1,a+1,z) = a^2 \sum{n = 0}^{\infty} \left(\frac{1}{a+n}\right)^2 \frac{z^n}{n!},$$

and so i thought that it might be interesting to give it a specific implementation, surely the code for pFq is way too general for this one.

What do you think ? Could you point me in the right direction to implement something specific ?

[MikaelSlevinsky]

This may sound backwards, but I would probably implement the partial derivative of the incomplete gamma function $\gamma(a,z) = \frac{z^a}{a}M(a,a+1,-z)$ with respect to $a$ with $M$ already implemented here as the dual part of evaluating the mathematical definition at dual number a.

julia> using HypergeometricFunctions

julia> using HypergeometricFunctions.DualNumbers

julia> a = Dual(2.0, 1.0)
2.0 + 1.0ɛ

julia> z = 0.5
0.5

julia> inc_gam(a,z) = z^a/a*HypergeometricFunctions.M(a, a+1, -z)
inc_gam (generic function with 1 method)

julia> inc_gam(a,z)
0.09020401043104985 - 0.11289739433586389ɛ

julia> (inc_gam(2.0+sqrt(eps()), z) - inc_gam(2-sqrt(eps()), z))/2sqrt(eps()) # central difference
-0.11289739375934005

julia> imag(inc_gam(2.0+im*eps(), z)/eps()) # complex difference
-0.11289739433586389

[cadojo]

This may sound backwards, but I would probably implement the partial derivative of the incomplete gamma function γ(a,z)=zaaM(a,a+1,−z) with respect to a with M already implemented here...

For anyone interested, I just got this to work with the SciML ecosystem (ModelingToolkit and Symbolics, specifically) by registering HypergeometricFunctions.M, and SpecialFunctions.gamma as symbolic functions with floating point types. This let me take the gradient of the Power Law Cutoff Potential using ModelingToolkit.calculate_gradient!

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)

3-element Vector{Num}:
 (-2((G*m*lowergamma((3//2) - (1//2)*α, (x(t)^2 + y(t)^2 + z(t)^2) / (c^2))*α) / (4(SpecialFunctions.gamma((5//2) - (1//2)*α)^2)*(sqrt(x(t)^2 + y(t)^2 + z(t)^2)^2)))*x(t)*SpecialFunctions.gamma((5//2) - (1//2)*α)) / sqrt(x(t)^2 + y(t)^2 + z(t)^2) + (-2G*m*x(t)*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((-1//2)*α))) / ((c^3)*SpecialFunctions.gamma((3//2) - (1//2)*α)) + (3G*m*x(t)*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((1//2) - (1//2)*α))*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))) / ((c^2)*SpecialFunctions.gamma((5//2) - (1//2)*α)*sqrt(x(t)^2 + y(t)^2 + z(t)^2)) + (-G*m*x(t)*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((1//2) - (1//2)*α))*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))*α) / ((c^2)*SpecialFunctions.gamma((5//2) - (1//2)*α)*sqrt(x(t)^2 + y(t)^2 + z(t)^2)) + (-2x(t)*SpecialFunctions.gamma((5//2) - (1//2)*α)*((-3G*m*lowergamma((3//2) - (1//2)*α, (x(t)^2 + y(t)^2 + z(t)^2) / (c^2))) / (4(SpecialFunctions.gamma((5//2) - (1//2)*α)^2)*(sqrt(x(t)^2 + y(t)^2 + z(t)^2)^2)))) / sqrt(x(t)^2 + y(t)^2 + z(t)^2)
 (-2((G*m*lowergamma((3//2) - (1//2)*α, (x(t)^2 + y(t)^2 + z(t)^2) / (c^2))*α) / (4(SpecialFunctions.gamma((5//2) - (1//2)*α)^2)*(sqrt(x(t)^2 + y(t)^2 + z(t)^2)^2)))*y(t)*SpecialFunctions.gamma((5//2) - (1//2)*α)) / sqrt(x(t)^2 + y(t)^2 + z(t)^2) + (3G*m*y(t)*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((1//2) - (1//2)*α))*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))) / ((c^2)*SpecialFunctions.gamma((5//2) - (1//2)*α)*sqrt(x(t)^2 + y(t)^2 + z(t)^2)) + (-G*m*y(t)*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((1//2) - (1//2)*α))*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))*α) / ((c^2)*SpecialFunctions.gamma((5//2) - (1//2)*α)*sqrt(x(t)^2 + y(t)^2 + z(t)^2)) + (-2y(t)*SpecialFunctions.gamma((5//2) - (1//2)*α)*((-3G*m*lowergamma((3//2) - (1//2)*α, (x(t)^2 + y(t)^2 + z(t)^2) / (c^2))) / (4(SpecialFunctions.gamma((5//2) - (1//2)*α)^2)*(sqrt(x(t)^2 + y(t)^2 + z(t)^2)^2)))) / sqrt(x(t)^2 + y(t)^2 + z(t)^2) + (-2G*m*y(t)*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((-1//2)*α))) / ((c^3)*SpecialFunctions.gamma((3//2) - (1//2)*α))
        (-2SpecialFunctions.gamma((5//2) - (1//2)*α)*((-3G*m*lowergamma((3//2) - (1//2)*α, (x(t)^2 + y(t)^2 + z(t)^2) / (c^2))) / (4(SpecialFunctions.gamma((5//2) - (1//2)*α)^2)*(sqrt(x(t)^2 + y(t)^2 + z(t)^2)^2)))*z(t)) / sqrt(x(t)^2 + y(t)^2 + z(t)^2) + (-2((G*m*lowergamma((3//2) - (1//2)*α, (x(t)^2 + y(t)^2 + z(t)^2) / (c^2))*α) / (4(SpecialFunctions.gamma((5//2) - (1//2)*α)^2)*(sqrt(x(t)^2 + y(t)^2 + z(t)^2)^2)))*SpecialFunctions.gamma((5//2) - (1//2)*α)*z(t)) / sqrt(x(t)^2 + y(t)^2 + z(t)^2) + (3G*m*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((1//2) - (1//2)*α))*z(t)*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))) / ((c^2)*SpecialFunctions.gamma((5//2) - (1//2)*α)*sqrt(x(t)^2 + y(t)^2 + z(t)^2)) + (-G*m*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((1//2) - (1//2)*α))*z(t)*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))*α) / ((c^2)*SpecialFunctions.gamma((5//2) - (1//2)*α)*sqrt(x(t)^2 + y(t)^2 + z(t)^2)) + (-2G*m*z(t)*exp(-((x(t)^2 + y(t)^2 + z(t)^2) / (c^2)))*(((x(t)^2 + y(t)^2 + z(t)^2) / (c^2))^((-1//2)*α))) / ((c^3)*SpecialFunctions.gamma((3//2) - (1//2)*α))