Add initial Bessel rule support

Author: AshtonSBradley

Project.toml

@@ -16,6 +16,7 @@ HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
@@ -27,6 +28,7 @@ FresnelIntegrals = "0.2.0"
 HypergeometricFunctions = "0.3.28"
 Nemo = "0.51, 0.52, 0.53, 0.54, 0.55"
 PolyLog = "2.6.0"
+SpecialFunctions = "2"
 SymbolicUtils = "4.27"
 Symbolics = "7"
 julia = "1.10"

src/methods/rule_based/general.jl

@@ -1,6 +1,6 @@
 # ===== Special functions from a lot of packages
 
-# function from SpecialFunctions.jl used:
+# functions from SpecialFunctions.jl used:
 # SpecialFunctions.expinti(x)
 # SpecialFunctions.expint(nu, z)
 # SpecialFunctions.gamma(x, y)
@@ -10,6 +10,18 @@
 # SpecialFunctions.loggamma
 # SpecialFunctions.erfi
 # SpecialFunctions.erf
+# SpecialFunctions.besselj(nu, z)
+# SpecialFunctions.bessely(nu, z)
+# SpecialFunctions.besseli(nu, z)
+# SpecialFunctions.besselk(nu, z)
+# SpecialFunctions.hankelh1(nu, z)
+# SpecialFunctions.hankelh2(nu, z)
+# SpecialFunctions.airyai(z)
+# SpecialFunctions.airybi(z)
+# SpecialFunctions.airyaiprime(z)
+# SpecialFunctions.airybiprime(z)
+
+using SpecialFunctions
 
 using Elliptic 
 @register_symbolic Elliptic.F(phi, m) # incomplete first kind
@@ -230,4 +242,6 @@ all_rules_paths = [
 "7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.1 (a+b arctanh(c x^n))^p.jl"
 "7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.2 (d x)^m (a+b arctanh(c x^n))^p.jl"
 "7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.3 (d+e x)^m (a+b arctanh(c x^n))^p.jl"
+
+"8 Special functions/8.1 Bessel functions.jl"
 ]

src/methods/rule_based/rules2/8 Special functions/8.1 Bessel functions.jl

@@ -0,0 +1,93 @@
+file_rules = [
+    #(* ::Subsection::Closed:: *)
+    #(* 8.1 Bessel functions *)
+
+    (
+        "8_1_1",
+        :((~x)^(~!nu) * besselj((~!mu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~mu), (~nu), (~x)) &&
+                eq((~mu), (~nu) - 1) &&
+                !eq((~a), 0) ?
+                (~x)^(~nu) * SpecialFunctions.besselj((~nu), (~a) * (~x)) ⨸ (~a) : nothing
+        ),
+    )
+
+    (
+        "8_1_2",
+        :((~x)^(~!nu) * bessely((~!mu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~mu), (~nu), (~x)) &&
+                eq((~mu), (~nu) - 1) &&
+                !eq((~a), 0) ?
+                (~x)^(~nu) * SpecialFunctions.bessely((~nu), (~a) * (~x)) ⨸ (~a) : nothing
+        ),
+    )
+
+    (
+        "8_1_3",
+        :((~x)^(~!nu) * besseli((~!mu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~mu), (~nu), (~x)) &&
+                eq((~mu), (~nu) - 1) &&
+                !eq((~a), 0) ?
+                (~x)^(~nu) * SpecialFunctions.besseli((~nu), (~a) * (~x)) ⨸ (~a) : nothing
+        ),
+    )
+
+    (
+        "8_1_4",
+        :((~x)^(~!nu) * besselk((~!mu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~mu), (~nu), (~x)) &&
+                eq((~mu), (~nu) - 1) &&
+                !eq((~a), 0) ?
+                -(~x)^(~nu) * SpecialFunctions.besselk((~nu), (~a) * (~x)) ⨸ (~a) : nothing
+        ),
+    )
+
+    (
+        "8_1_5",
+        :((~x)^(~!m) * besselj((~!nu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~m), (~nu), (~x)) &&
+                igt(simplify((~m) - (~nu) - 1), 0) &&
+                ext_iseven(simplify((~m) - (~nu) - 1)) &&
+                !eq((~a), 0) ?
+                (~x)^(~m) * SpecialFunctions.besselj((~nu) + 1, (~a) * (~x)) ⨸ (~a) -
+                simplify((~m) - (~nu) - 1) ⨸ (~a) * ∫((~x)^((~m) - 1) * SpecialFunctions.besselj((~nu) + 1, (~a) * (~x)), (~x)) : nothing
+        ),
+    )
+
+    (
+        "8_1_6",
+        :((~x)^(~!m) * bessely((~!nu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~m), (~nu), (~x)) &&
+                igt(simplify((~m) - (~nu) - 1), 0) &&
+                ext_iseven(simplify((~m) - (~nu) - 1)) &&
+                !eq((~a), 0) ?
+                (~x)^(~m) * SpecialFunctions.bessely((~nu) + 1, (~a) * (~x)) ⨸ (~a) -
+                simplify((~m) - (~nu) - 1) ⨸ (~a) * ∫((~x)^((~m) - 1) * SpecialFunctions.bessely((~nu) + 1, (~a) * (~x)), (~x)) : nothing
+        ),
+    )
+
+    (
+        "8_1_7",
+        :((~x)^(~!m) * besseli((~!nu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~m), (~nu), (~x)) &&
+                igt(simplify((~m) - (~nu) - 1), 0) &&
+                ext_iseven(simplify((~m) - (~nu) - 1)) &&
+                !eq((~a), 0) ?
+                (~x)^(~m) * SpecialFunctions.besseli((~nu) + 1, (~a) * (~x)) ⨸ (~a) -
+                simplify((~m) - (~nu) - 1) ⨸ (~a) * ∫((~x)^((~m) - 1) * SpecialFunctions.besseli((~nu) + 1, (~a) * (~x)), (~x)) : nothing
+        ),
+    )
+
+    (
+        "8_1_8",
+        :((~x)^(~!m) * besselk((~!nu), (~!a) * (~x))) => :(
+            !contains_var((~a), (~m), (~nu), (~x)) &&
+                igt(simplify((~m) - (~nu) - 1), 0) &&
+                ext_iseven(simplify((~m) - (~nu) - 1)) &&
+                !eq((~a), 0) ?
+                -(~x)^(~m) * SpecialFunctions.besselk((~nu) + 1, (~a) * (~x)) ⨸ (~a) +
+                simplify((~m) - (~nu) - 1) ⨸ (~a) * ∫((~x)^((~m) - 1) * SpecialFunctions.besselk((~nu) + 1, (~a) * (~x)), (~x)) : nothing
+        ),
+    )
+
+]

src/methods/rule_based/rules2/9 Miscellaneous/0.1 Integrand simplification rules.jl

@@ -1,6 +1,14 @@
 file_rules = [
 ("0_1_0", :(+(~~a)) => :(sum(map(f -> ∫(f,~x), arguments(~a)))) )
 
+("0_1_1",
+:(*(~~a)) => :(
+let
+    terms = distribute_special_function_product(arguments(~a), ~x)
+    terms === nothing ? nothing : sum(map(term -> ∫(term, ~x), terms))
+end
+))
+
 ("0_1_6",
 :((~!u)*((~!a)*(~v) + (~!b)*(~v) + (~!w))^(~!p)) => :(
     !contains_var((~a), (~b), (~x)) &&
@@ -172,4 +180,4 @@ end
     ext_isinteger((~p)) ?
 1⨸(~c)^(~p)*∫((~u)*((~b)⨸2 + (~c)*(~x)^(~n))^(2*(~p)), (~x)) : nothing))
 
-]
+]

src/methods/rule_based/rules_utility_functions.jl

@@ -478,6 +478,59 @@ function dist(exp1, exp2, x)
     end
 end
 
+const SPECIAL_FUNCTION_OPERATIONS = (
+    SpecialFunctions.expinti,
+    SpecialFunctions.expint,
+    SpecialFunctions.gamma,
+    SpecialFunctions.sinint,
+    SpecialFunctions.cosint,
+    SpecialFunctions.loggamma,
+    SpecialFunctions.erfi,
+    SpecialFunctions.erf,
+    SpecialFunctions.besselj,
+    SpecialFunctions.bessely,
+    SpecialFunctions.besseli,
+    SpecialFunctions.besselk,
+    SpecialFunctions.hankelh1,
+    SpecialFunctions.hankelh2,
+    SpecialFunctions.airyai,
+    SpecialFunctions.airybi,
+    SpecialFunctions.airyaiprime,
+    SpecialFunctions.airybiprime,
+    Elliptic.F,
+    Elliptic.E,
+    Elliptic.Pi,
+    hypergeometric2f1,
+    hypergeometricpFq,
+    PolyLog.reli,
+    FresnelIntegrals.fresnelc,
+    FresnelIntegrals.fresnels,
+    sinhintegral,
+    coshintegral,
+)
+
+function special_function_operation(op)
+    return any(special_op -> op === special_op, SPECIAL_FUNCTION_OPERATIONS)
+end
+
+function contains_special_function(u)
+    u = Symbolics.unwrap(u)
+    !iscall(u) && return false
+    special_function_operation(operation(u)) && return true
+    return any(contains_special_function(arg) for arg in arguments(u))
+end
+
+function distribute_special_function_product(factors, x)
+    sum_index = findfirst(factor -> issum(factor) && contains_var(factor, x), factors)
+    sum_index === nothing && return nothing
+    any(contains_special_function, factors) || return nothing
+
+    sum_factor = Symbolics.unwrap(factors[sum_index])
+    return map(arguments(sum_factor)) do term
+        prod(i == sum_index ? term : factor for (i, factor) in pairs(factors))
+    end
+end
+
 # linear(a+3x,x) true
 # linear((x+1)^2 - x^2 - 1,x) true
 function linear(args...)
@@ -900,4 +953,4 @@ function nonfree_addends(u, x)
     return contains_var(u, x) ? 1 : u
 end
 
-# TODO are all this unwrap needed?
+# TODO are all this unwrap needed?

src/methods/rule_based/translator_of_rules.jl

@@ -280,6 +280,16 @@ function translate_integrand(integrand)
         ("ArcTanh", "atanh"),
         ("ArcCoth", "acoth"),
 
+        ("BesselJ", "besselj", 2),
+        ("BesselY", "bessely", 2),
+        ("BesselI", "besseli", 2),
+        ("BesselK", "besselk", 2),
+        ("HankelH1", "hankelh1", 2),
+        ("HankelH2", "hankelh2", 2),
+        ("AiryAiPrime", "airyaiprime", 1),
+        ("AiryBiPrime", "airybiprime", 1),
+        ("AiryAi", "airyai", 1),
+        ("AiryBi", "airybi", 1),
         ("PolyLog", "PolyLog.reli", 2),
         ("Gamma", "SymbolicUtils.gamma"),
     ]
@@ -354,6 +364,16 @@ function result_substitutions(result, vardefs)
         ("Erf", "SymbolicUtils.erf"),
         ("SinIntegral", "SymbolicUtils.sinint"),
         ("CosIntegral", "SymbolicUtils.cosint"),
+        ("BesselJ", "SpecialFunctions.besselj", 2),
+        ("BesselY", "SpecialFunctions.bessely", 2),
+        ("BesselI", "SpecialFunctions.besseli", 2),
+        ("BesselK", "SpecialFunctions.besselk", 2),
+        ("HankelH1", "SpecialFunctions.hankelh1", 2),
+        ("HankelH2", "SpecialFunctions.hankelh2", 2),
+        ("AiryAiPrime", "SpecialFunctions.airyaiprime", 1),
+        ("AiryBiPrime", "SpecialFunctions.airybiprime", 1),
+        ("AiryAi", "SpecialFunctions.airyai", 1),
+        ("AiryBi", "SpecialFunctions.airybi", 1),
         # taken from other julia packages
         ("EllipticE", "elliptic_e", (1,2)),
         ("EllipticF", "elliptic_f", 2),
@@ -655,4 +675,4 @@ try
 catch e
     println("Error during translation: $e")
     exit(1)
-end
+end

test/methods/rule_based/test_rule2.jl

@@ -1,4 +1,6 @@
 using Test
+import SymbolicIntegration
+using SpecialFunctions
 using Symbolics
 using SymbolicIntegration: eq
 
@@ -53,6 +55,19 @@ end
     rs = :(exp(~x)) => :(~x)
     @test eq(SymbolicIntegration.rule2(rs, exp(x+1)), x+1)
     @test eq(SymbolicIntegration.rule2(rs, ℯ^x), x)
+
+    rbj = :(besselj(~nu, ~z)) => :(~nu, ~z)
+    besselj_match = SymbolicUtils.arguments(SymbolicIntegration.rule2(rbj, SpecialFunctions.besselj(0, x)))
+    @test SymbolicUtils.unwrap_const(besselj_match[1]) == 0
+    @test isequal(besselj_match[2], x)
+
+    rbi = :(besseli(~nu, ~z)) => :(~nu, ~z)
+    besseli_match = SymbolicUtils.arguments(SymbolicIntegration.rule2(rbi, SpecialFunctions.besseli(1, x)))
+    @test SymbolicUtils.unwrap_const(besseli_match[1]) == 1
+    @test isequal(besseli_match[2], x)
+
+    rai = :(airyai(~z)) => :(~z)
+    @test eq(SymbolicIntegration.rule2(rai, SpecialFunctions.airyai(x)), x)
 end
 
 @testset "Segment" begin

test/methods/rule_based/test_special_functions.jl

@@ -0,0 +1,63 @@
+using Test
+using SymbolicIntegration
+using SpecialFunctions
+using Symbolics
+
+@testset "RuleBased special-function integrals" begin
+    @variables x a
+    method = RuleBasedMethod()
+
+    @test isequal(
+        integrate(x * SpecialFunctions.besselj(0, a * x), x, method),
+        x * SpecialFunctions.besselj(1, a * x) / a,
+    )
+    @test isequal(
+        integrate(x^2 * SpecialFunctions.besselj(1, a * x), x, method),
+        x^2 * SpecialFunctions.besselj(2, a * x) / a,
+    )
+    @test isequal(
+        integrate(x * SpecialFunctions.bessely(0, a * x), x, method),
+        x * SpecialFunctions.bessely(1, a * x) / a,
+    )
+    @test isequal(
+        integrate(x * SpecialFunctions.besseli(0, a * x), x, method),
+        x * SpecialFunctions.besseli(1, a * x) / a,
+    )
+    @test isequal(
+        integrate(x * SpecialFunctions.besselk(0, a * x), x, method),
+        -x * SpecialFunctions.besselk(1, a * x) / a,
+    )
+
+    @test isequal(
+        integrate(x^3 * SpecialFunctions.besselj(0, a * x), x, method),
+        x^3 * SpecialFunctions.besselj(1, a * x) / a -
+            2 * x^2 * SpecialFunctions.besselj(2, a * x) / a^2,
+    )
+    @test isequal(
+        integrate((x + x^3) * SpecialFunctions.besselj(0, a * x), x, method),
+        x * SpecialFunctions.besselj(1, a * x) / a +
+            x^3 * SpecialFunctions.besselj(1, a * x) / a -
+            2 * x^2 * SpecialFunctions.besselj(2, a * x) / a^2,
+    )
+    @test isequal(
+        integrate(x * (1 + x^2) * SpecialFunctions.besselj(0, a * x), x, method),
+        x * SpecialFunctions.besselj(1, a * x) / a +
+            x^3 * SpecialFunctions.besselj(1, a * x) / a -
+            2 * x^2 * SpecialFunctions.besselj(2, a * x) / a^2,
+    )
+    @test isequal(
+        integrate((1 + x) * SpecialFunctions.airyai(x), x, method),
+        SymbolicIntegration.∫(SpecialFunctions.airyai(x), x) +
+            SymbolicIntegration.∫(x * SpecialFunctions.airyai(x), x),
+    )
+    @test isequal(
+        integrate(x^3 * SpecialFunctions.besseli(0, a * x), x, method),
+        x^3 * SpecialFunctions.besseli(1, a * x) / a -
+            2 * x^2 * SpecialFunctions.besseli(2, a * x) / a^2,
+    )
+    @test isequal(
+        integrate(x^3 * SpecialFunctions.besselk(0, a * x), x, method),
+        -x^3 * SpecialFunctions.besselk(1, a * x) / a -
+            2 * x^2 * SpecialFunctions.besselk(2, a * x) / a^2,
+    )
+end

test/methods/rule_based/translator_of_testset.jl

@@ -16,6 +16,16 @@ function translate_mathematica_to_julia(expr::String)
         ("Erf", "SymbolicUtils.erf"),
         ("SinIntegral", "SymbolicUtils.sinint"),
         ("CosIntegral", "SymbolicUtils.cosint"),
+        ("BesselJ", "SymbolicIntegration.SpecialFunctions.besselj", 2),
+        ("BesselY", "SymbolicIntegration.SpecialFunctions.bessely", 2),
+        ("BesselI", "SymbolicIntegration.SpecialFunctions.besseli", 2),
+        ("BesselK", "SymbolicIntegration.SpecialFunctions.besselk", 2),
+        ("HankelH1", "SymbolicIntegration.SpecialFunctions.hankelh1", 2),
+        ("HankelH2", "SymbolicIntegration.SpecialFunctions.hankelh2", 2),
+        ("AiryAiPrime", "SymbolicIntegration.SpecialFunctions.airyaiprime", 1),
+        ("AiryBiPrime", "SymbolicIntegration.SpecialFunctions.airybiprime", 1),
+        ("AiryAi", "SymbolicIntegration.SpecialFunctions.airyai", 1),
+        ("AiryBi", "SymbolicIntegration.SpecialFunctions.airybi", 1),
         # taken from other julia packages
         ("EllipticE", "SymbolicIntegration.elliptic_e", (1,2)),
         ("EllipticF", "SymbolicIntegration.elliptic_f", 2),
@@ -180,4 +190,3 @@ catch e
     println("Error during translation: $e")
     exit(1)
 end
-