# How do I speed up this numerical integration with Bessel functions? (Tried changing Working Precision and Integration Strategies)

I need to compute the function P[b] given as follows:

P[ b_?NumericQ] :=
NIntegrate[(-(1/(k1^2 + b^2)^(3/2)) a*
b *((xt[s, c] Sinh[s *k1] + yt[s, c] Cosh[s *k1]) ) +
1/Sqrt[k1^2 + b^2]
b *( (s Cosh[s *k1] xt[s, c] +
s Sinh[s *k1] yt[s, c]))) BesselJ[1, s *b], {s, 0, 0,
smax}, {c, 0, l}, MaxRecursion -> 20,
Method -> {"ClenshawCurtisRule", "SymbolicProcessing" -> 0},
PrecisionGoal -> 5, AccuracyGoal -> Automatic,
WorkingPrecision -> 50];


for large values of b. Assume k1=-1/2 and smax=100 although I may want to change those values.

Here

xt[s_?NumericQ,
c_?NumericQ] := ((-((N[
MeijerG[{{0, 3/2}, {1}}, {{0}, {}}, 4/(s^2 k1^2)], wp])) -
s (((l^2 - c^2)/(k1^2 + c^2)^(1/2)) BesselJ[1, s*c]))/(Sinh[
s*k1] - Tanh[s*k2]*Cosh[s*k1]));

yt[s_?NumericQ, c_?NumericQ] := -xt[s, c]*Tanh[s*k2];


Where l and k2 are real positive numbers as well (can be considered fixed). It has two very oscillating definite integrals ({c, 0, l} and {s, 0, 0, smax} with bessel functions and the inbuilt special MeijerG functions. Due to the oscillatory nature, other methods/ strategies aren't helping especially for large values of b and i've played around with the working precision, precision goal. Is there any other way that I can improve the speed to evaluate this function? (I need it to be fast as my next step is to numerically integrate a product of P with bessel functions!)

• Simplify the integrand prior to integration. expr2 = Assuming[{l > c >= 0, s > 0, k2 > 0, Element[{a, b, k1}, Reals]}, expr // FunctionExpand // FullSimplify] This will get rid of the MeijerG Commented Dec 22, 2022 at 4:58
• Keeping the argument s_?NumeriQ, fullsimplifying xt does not seem to get rid of the MeijerG... am i missing something ?
– Dave
Commented Dec 23, 2022 at 6:38

THIS IS AN EXTENDED COMMENT RATHER THAN AN ANSWER

Clear["Global*"]


Note that restricting a function's arguments with NumericQ is only appropriate if the function uses a numeric technique that would preclude the function from being able to evaluate with symbolic values.

xt[s_, c_] := ((-((MeijerG[{{0, 3/2}, {1}}, {{0}, {}}, 4/(s^2 k1^2)])) -
s (((l^2 - c^2)/(k1^2 + c^2)^(1/2)) BesselJ[1, s*c]))/(Sinh[s*k1] -
Tanh[s*k2]*Cosh[s*k1]));

yt[s_, c_] := -xt[s, c]*Tanh[s*k2];


The integrand for subsequent calculation is

expr = (-(1/(k1^2 + b^2)^(3/2)) a*
b*((xt[s, c] Sinh[s*k1] + yt[s, c] Cosh[s*k1])) +
1/Sqrt[k1^2 + b^2] b*((s Cosh[s*k1] xt[s, c] +
s Sinh[s*k1] yt[s, c]))) BesselJ[1, s*b];


This can be simplified to

expr2 = Assuming[{l > c >= 0, s > 0, k2 > 0, Element[{a, b, k1}, Reals]},
expr // FunctionExpand // FullSimplify]

(* (1/((b^2 + k1^2)^(
3/2) s))b BesselJ[1,
b s] (-a + (b^2 + k1^2) s Coth[(k1 - k2) s]) (((c - l) (c + l) s^2 BesselJ[
1, c s])/Sqrt[c^2 + k1^2] + (1/k1)
2 Sqrt[π] Sign[k1]^2 (k1 s + Cosh[k1 s] Sign[k1] - Sinh[k1 s])) *)


The integrand is then free of the MeijerG`

• Than you, this does indeed simplify the integrand, however, it doesnt seem to increase the speed of the double integration. Perhaps this is due to the oscillatory nature ?
– Dave
Commented Dec 26, 2022 at 18:18
• What values of all of the parameters are you using to evaluate these? Commented Dec 26, 2022 at 19:36
• k2 = 2.5, k1 = -1/2, b = 10, a = -1/2 , l = 0.5 Tan[11 Pi/24]
– Dave
Commented Dec 27, 2022 at 0:34