TransWikia.com

Numerical integration of MeijerG function with variables

Mathematica Asked by S.Kirubakaran on November 25, 2020

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.

One 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.

enter image description here

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]]

enter image description here

(* 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}]

enter image description here

Answered by bbgodfrey on November 25, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP