0
$\begingroup$

I have the following (long) function of g1 and ωthat I wish to integrate with respect to ω

f = (200 (7692800000000 g1^4 (1000000 + ω (100 + ω) (400 + ω (100 + ω))) + 601(44100 - 399 ω^2 + 4 ω^4) (400000001 + 40000 ω (200 + ω)) (396054100 + ω (200 + ω) (79601 + 4 ω (200 + ω))) + 800 g1^2 (165388780876943610000 + ω (9974810249137020000 + ω (338251703143216200 + ω (17987201082240200 + ω (541408115948001 + 4 ω (1956801870100 + ω (13896010001 + 40000 ω (1101 + ω))))))))))/(2560000000000000000 g1^4 (1000000 + ω (100 + ω) (400 + ω (100 + ω))) + (1 + 40000 ω^2) (44100 - 399 ω^2 + 4 ω^4) (400000001 + 40000 ω (200 + ω)) (396054100 + ω (200 + ω) (79601 + 4 ω (200 + ω))) - 3200000000 g1^2 (-83117559000000 + ω (100 + ω) (-17690047067700 + ω (100 + ω) (-5302720399 + 4 ω (100 + ω) (3989599 + 40000 ω (100 + ω))))))

I first construct a table and prepend g1 while integrating with respect to ω

int = Table[{g1, 1/(2 π)*NIntegrate[f, {ω, -200, 200}]}, {g1, 0, 100}]

and I am returned with a bunch of warnings, saying something like

"NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in ω near {ω} = {-0.00621666}. NIntegrate obtained 97.19546034578347` and 0.8830472862609543` for the integral and error estimates"

Which I ignore since it's just a warning (I don't know if I should trust the results). I then proceed to do a ListLogLogPlot of int as a function of g1

plot = ListLogLogPlot[int, PlotRange -> {All, All}, Joined -> True, ImageSize -> Large, Frame -> True, FrameTicksStyle -> Directive[Black, 12], LabelStyle -> Directive[Black, 12], PlotRangePadding -> None, PlotStyle -> {Directive[Black]}]

And I am returned with (ignoring the plot label)enter image description here

Now I wish to increase my time step, so under my NIntegrate I change my integration step to 0.1

int = Evaluate@Table[{g1, 1/(2 π)*NIntegrate[f, {ω, -200, 200}]}, {g1, 1/10, 100, 1/10}]

However, the computation time took extremely long and I was returned with this enter image description here

Clearly, there are some numerics that Mathematica was mishandling in the integration process. I think this might have something to do with the warnings but I'm unsure of how to remedy it. The reason I want to push it to 0.1 was because LogLogPlot does not show g1 starting from 0 since there is no such thing as the Log of 0. So pushing my starting point to 0.1 will bring me close to 0. In fact, I intend to push my starting point to 0.01 in steps of 0.01 for g1 but it will take tremendously long. How should I go about this?

Thanks

$\endgroup$
3
$\begingroup$

Why not use ParametricNDSolveValue?:

int = ParametricNDSolveValue[
    {i'[ω] == f/(2 π), i[-200]==0},
    i[200],
    {ω,-200,200},
    g1
]

enter image description here

Plot:

LogLogPlot[int[g], {g, 1, 100}]

enter image description here

Plot with smaller starting value for g1:

LogLogPlot[int[g], {g, .01, 100}]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks, I'm studying ParametricNDSolveValue as we speak. However, I think there is a mistake on the second plot. When g1 is close to 0 the curve should start at around 300 but you have it over 1000. Is it because of numerics? $\endgroup$ – kowalski Mar 12 '19 at 19:24
  • $\begingroup$ @kowalski That's because I didn't include the factor of 1/(2π) in the integral. $\endgroup$ – Carl Woll Mar 12 '19 at 19:36
  • $\begingroup$ Thanks. I have another question: why did you set i[-200] = 0 followed by i[200]? Are you setting a boundary condition? Also the spike at around g1=30 goes all the way up, to near infinity, is that a pole? It is not the case for the smaller time step setting as shown in my first plot. The spike peaks at 10. Why should the spike peak differently in a smaller integrating step? $\endgroup$ – kowalski Mar 12 '19 at 19:39
  • 1
    $\begingroup$ @kowalski The i[-200], i[200] bits are how you translate an integral into a differential equation. Your first plot has a smaller peak because you didn't sample g1 close enough to the singularity. $\endgroup$ – Carl Woll Mar 12 '19 at 19:45
  • $\begingroup$ Thanks! I will have to study this ParametricNDSolveValue function in detail. This is useful to me. $\endgroup$ – kowalski Mar 12 '19 at 19:52
2
$\begingroup$

Using a higher-order integration rule will reduce the need for recursion when the function is analytic in a (complex) neighborhood of the interval of integration:

int = Evaluate@Table[
     {g1, 1/(2 π) NIntegrate[f, {ω, -200, 0, 200}, 
        Method -> {"GaussKronrodRule", "Points" -> 21}]},
     {g1, 10.^Subdivide[-2., 2, 100]}]; // AbsoluteTiming
(*  {1.07156, Null}  *)

plot = ListLogLogPlot[int, PlotRange -> {All, All}, Joined -> True, 
  ImageSize -> Large, Frame -> True, 
  FrameTicksStyle -> Directive[Black, 12], 
  LabelStyle -> Directive[Black, 12], PlotRangePadding -> None, 
  PlotStyle -> {Directive[Black]}]

Mathematica graphics

Or:

LogLogPlot[1/(2 π)* NIntegrate[f, {ω, -200, 0, 200}, 
   Method -> {"GaussKronrodRule", "Points" -> 21}],
 {g1, 0.01, 100}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ I don't have an issue with integrating steps of g1 from 0 to 100 but rather, from 0.01 to 100 in steps of 0.01. The latter is computationally intensive and I need a simpler way to do that. $\endgroup$ – kowalski Mar 12 '19 at 20:41
  • $\begingroup$ @kowalski Oops, copied the wrong code. Fixed. $\endgroup$ – Michael E2 Mar 12 '19 at 21:48
  • $\begingroup$ Thanks! That works really well too. I have a quick question, what does Points -> 21 do? Is it like a PlotPoints option? If so, why can't I just put it under the LogLogPlot function? The closest that I can find on the documentation is PlotPoints and just Point which isn't used under Methods $\endgroup$ – kowalski Mar 13 '19 at 18:05
  • 1
    $\begingroup$ @kowalski It's an option for the Gauss-Kronrod rule (and other interpolatory rules) that specifies how many points to use for the underlying Gauss-Legendre method; Kronrod showed how to add points to these that simultaneously approximates both the integral and the error. When the integrand is analytic, the Gauss-Legendre method converges very rapidly as the number of points increases, provided there are at least a couple of points per oscillation.The integrand has very few oscillations.... $\endgroup$ – Michael E2 Mar 13 '19 at 18:17
  • $\begingroup$ The default strategy is a heuristic balance of points (default 5) and subdivision of the interval. The more points the fewer subdivisions and vice versa. Too much of one or the other usually wastes time. Subdivisions can help with oscillatory regions and other odd behavior. One subdivision that certainly helps is ω == 0 for small g1, where there is a narrow spike; that's why I added a 0 between the ±200. $\endgroup$ – Michael E2 Mar 13 '19 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.