Numerical integration of MeijerG function with variables

I applied the integration, How to solve the function which is mentioned below. here, 'theta' is the only variable. q13, q1,q2,p,qi1 are all variables will take value from the loop ranges from 0 to 2, mu is 1.5, L=3, gammak1=1, BI=9, C1=4, Bsd=4, M=2

• a[theta_] = MeijerG[{{-(q13 + q1 + q2 + (2mu)+ p - 1), 1 - (q13 + q1 + qi1 + mu + p + (Lmu)), 1 - (q13 + q1 + mu + p)}, {}}, {{0}, {-(q13 + q1 + mu + p)}}, C1 / (BI*gammak1* (Bsd + ((Sin[Pi/M]^2)/(Sin[theta]^2))))]

• Then NIntegrate[a[theta],{theta,0,Pi/2}]

I am not getting the values when I integrate. Even for a day, the simulation is going on.

• You have some variables names together. Also, what are the values for qi1 and M? Jul 23, 2020 at 14:50
• @bbgodfrey - with M -> 1, a[theta] is independent of theta. Jul 23, 2020 at 15:12
• @BobHanlon qi1 = 0 to 2 and M=2 Jul 23, 2020 at 19:30
• @bbgodfrey then is it possible to Change this Meijer G function to hypergeometric function or any other form, so that I need to integrate the theta variable? Jul 23, 2020 at 19:32
• Have thought more about this. Some progress appears possible. Jul 24, 2020 at 5:32

1 Answer

The root cause of the difficulties encountered by the OP is that MeijerG, as implemented in Mathematica, often is not accurate when evaluated at machine precision. For instance, with the parameters,

params = {p -> 1, q1 -> 1, q2 -> 1, q13 -> 1, qi1 -> 1, mu -> 3/2,
L -> 3, gammak1 -> 1, BI -> 9, C1 -> 4, Bsd -> 4, M -> 2};


merely plotting the integrand in the Question over a portion of its domain

mg = MeijerG[{{-(q13 + q1 + q2 + (2 mu) + p - 1),
1 - (q13 + q1 + qi1 + mu + p + (L mu)),
1 - (q13 + q1 + mu + p)}, {}}, {{0}, {-(q13 + q1 + mu + p)}},
C1/(BI*gammak1*(Bsd + ((Sin[Pi/M]^2)/(Sin[theta]^2))))] /. params;
LogPlot[mg, {theta, 1/2, Pi/2}, PlotRange -> {{0, Pi/2}, All},
ImageSize -> Large, AxesLabel -> {theta, "mg"}, LabelStyle -> {15, Bold, Black}]


is very slow, noisy, and highly inaccurate.

Not surprisingly, trying to NIntegrate this fails. In fact, obtaining a smooth curve requires a very large WorkingPrecision, especially near theta = 0.

lp = LogPlot[mg, {theta, 1/100, Pi/2}, WorkingPrecision -> 2000,
PlotRange -> {{0, Pi/2}, All}, ImageSize -> Large, AxesLabel -> {theta, "mg"},
LabelStyle -> {15, Bold, Black}])
Cases[lp, Line[z_] -> z, Infinity][[1, 2, 1]]


(* 0.022446 *)


Even with WorkingPrecision -> 2000, Plot can reach only down to 0.022446, well above the requested lower bound, 0.01. About 70 seconds was required for this computation. (Note that the first value obtained by Plot has been discarded, because it visibly is inaccurate.) Performing the desired integral requires the integrand to be evaluated all the way to theta = 0. To obtain a reasonable estimate of the function there, Fit the very left of the curve above to a biquadratic in theta.

Rest[First[Cases[lp, Line[z_] -> z, Infinity]]];
fi = Interpolation[%];
val = %%[[1, 1]];
Fit[%%%[[;; 10]], {1, theta^2, theta^4}, theta]
s = Piecewise[{{fi[theta], theta >= val}}, %];
(* 17.877 - 25.3885 theta^2 + 159.428 theta^4 *)


mg at theta = 1/100 can be computed from

Block[{$MaxExtraPrecision = 100000}, N[mg /. theta -> 10^-2, 6]] (* 5.79133*10^7 *)  in about 70 seconds and agrees well with Exp[s /. theta -> 1/100]. Plotting Exp[s] has a credible appearance as well. The integral now is trivial to perform. NIntegrate[Exp[s], {theta, 0, Pi/2}] (* 1.62895*10^7 *)  In general, the same procedure should work for other parameters. Addendum: Faster, simpler code Most of the code above can be replaced by fb = Interpolation@Block[{$MaxExtraPrecision = 100000},
Table[{theta, N[mg, 6]}, {theta, 2/100, Ceiling[Pi/2, 1/100], 1/100}]];
Quiet@NIntegrate[fb[theta], {theta, 0, Pi/2}]
(* 1.6289*10^7 *)


as before. This computation requires about 18 seconds, rather than 70 seconds for the previous approach.

Applying the present approach to a different set of parameters, say,

SeedRandom[1812];
params = Join[Thread[{p, q1, q2, q13, qi1} -> Rationalize[RandomReal[2, 5], 0]],
{mu -> 3/2, L -> 3, gammak1 -> 1, BI -> 9, C1 -> 4, Bsd -> 4, M -> 2}];


and defining mg and fb as before, but with the new params, then gives for the integral

(* 2.0744*10^6 *)


in about 30 seconds. For completeness, the correspond plot of the integrand is

Quiet@LogPlot[fb[theta], {theta, 0, Pi/2}, ImageSize -> Large,
AxesLabel -> {theta, "mg"}, LabelStyle -> {15, Bold, Black}]


• I am really satisfied with your answer, Sir. I have one doubt, If I integrate from 2/100 to Pi/2, will I get the approx answer as from 0 to Pi/2? Jul 26, 2020 at 15:41
• @S.Kirubakaran Integrating over {2/100, Pi/2} is not very accurate. However, NIntegrate[fb[theta], {theta, 0, Pi/2}] is quite accurate, because fb, an InterpolationFunction, automatically extrapolates to 0 for NIntegrate with good accuracy, since mg changes slowly there. See the final plot. Jul 26, 2020 at 17:03
• So, it means that I can directly integrate the mg with fb code. Am I right Sir? Jul 26, 2020 at 17:16
• @S.Kirubakaran If mg is smooth at theta = 0, as it appears to be, then you can integrate mg with fb over all of {0, Pi/2}. As a test, you can extend Table to include 1/100, but doing so can be very slow. By the way, notice that I applied Quiet to NIntegrate to turn off warning messages that extrapolation is being used near theta = 0. Jul 26, 2020 at 17:24