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

Question

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!)

Thanks in advance!

Answer 1

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