5
$\begingroup$

I have an integral which doesn't give a closed definite expression. The command

Integrate[x DiracDelta[r x - y] Exp[1/g^2 {Cos[x2 - x] + Cos[x2] + Cos[y+x2]}], {x, 0, 2 Pi}, {y, 0, 2 Pi}, {x2, -Pi, Pi}]

returns the expression itself. Therefore I evaluate this numerically using NIntegrate in order to understand the behaviour of the integral with $r$ and $g^2$.

In this situation my preferred way is to to obtain plots of the integral with respect to $g^2$ for a fixed $r$, and also with $r$ for a fixed $g^2$. How do I combine Plot and NIntegrate together and get a plot for the same over some values for $r$ and $g^2$? In other words, I would need to iterate NIntegrate, and construct a list of points so as to plot on the graph. How do I achieve that?

$\endgroup$
3
  • $\begingroup$ You may define a function f[r_?NumericQ, g_?NumericQ]:=NIntegrate[...]. $\endgroup$
    – Alx
    Commented Feb 29, 2020 at 0:48
  • $\begingroup$ @Alx Can you make it an answer? How would I generate the plot afterwards? $\endgroup$ Commented Feb 29, 2020 at 1:22
  • $\begingroup$ Plot3D[f[r,g],{r,...},{g,...}], you need to put in ranges for r and g. $\endgroup$
    – Alx
    Commented Feb 29, 2020 at 1:36

2 Answers 2

6
$\begingroup$

Progress can be made by integrating over y only (with the correction that {…} be replaced by (…)).

Integrate[x DiracDelta[r x - y] Exp[1/g^2 (Cos[x2 - x] + Cos[x2] + Cos[y + x2])], 
    {y, 0, 2 Pi}, Assumptions -> r > 0 && 0 < x < 2 Pi]
(* E^((Cos[x - x2] + Cos[x2] + Cos[r x + x2])/g^2) x HeavisideTheta[2 Pi - r x] *)

In contrast, if r < 0, the result is 0, as one would expect.

The remaining integrals take the form,

Integrate[E^((Cos[x - x2] + Cos[x2] + Cos[r x + x2])/g^2)
    x HeavisideTheta[2 Pi - r x], {x, 0, 2 Pi}, {x2, -Pi, Pi}]

Edit: which can be simplified further to

Integrate[E^((Cos[x - x2] + Cos[x2] + Cos[r x + x2])/g^2) x, 
    {x, 0, 2 Pi Min[1, 1/r]}, {x2, -Pi, Pi}

Now, if g^2 is much greater than 1, then the remaining integrations can be performed with the exponential ignored.

Integrate[x, {x, 0, 2 Pi Min[1, 1/r]}, {x2, -Pi, Pi}, Assumptions -> r > 0]
(* (4 Pi^3 Min[1, 1/r]^2 *)

For g^2 not large, the remaining integrals must be performed numerically along the lines suggested by Alx in a comment above, but a 2D plot is slow.

f[r_?NumericQ, g_?NumericQ] := NIntegrate[E^((Cos[x - x2] + Cos[x2] + Cos[r x + x2])/g^2) x, 
    {x, 0, 2 Pi Min[1, 1/r]}, {x2, -Pi, Pi}, 
    Method -> {Automatic, "SymbolicProcessing" -> False}]

A 1D plot for g == 1 is

Plot[f[r, 1], {r, 0, 4}, MaxRecursion -> 2, ImageSize -> Large, 
    AxesLabel -> {r, None}, LabelStyle -> {15, Bold, Black}]

enter image description here

For a small value of g, the curve is somewhat different.

LogPlot[f[r, 1/10], {r, 0, 4}, MaxRecursion -> 3, ImageSize -> Large, 
    AxesLabel -> {r, None}, LabelStyle -> {15, Bold, Black}, PlotRange -> All]

enter image description here

Peaks lie at r = 1 and very close to r = 0.

Addendum

The further simplification identified by Andreas (+1), equivalent to

f[r_?NumericQ, g_?NumericQ] := 2 Pi NIntegrate[BesselI[0, g^-2 
    Sqrt[3 + 2 Cos[x] + 2 Cos[r x] + 2 Cos[(1 + r) x]]] x, {x, 0, 2 Pi Min[1, 1/r]}, 
    Method -> {Automatic, "SymbolicProcessing" -> False}]

reproduces the plots above about an order of magnitude faster. For completeness, we provide a plot of the values at the two peak values as functions of g (rescaled by Exp[-2.94 g^-2] to keep the curves readable.

LogPlot[{f[.02, g] Exp[-2.94 g^-2], f[1, g] Exp[-2.94 g^-2]}, {g, 1/10, 3},
    ImageSize -> Large, AxesLabel -> {g, None}, LabelStyle -> {15, Bold, Black}, 
    PlotRange -> {{0, 3}, Automatic}]

The blue and orange curves correspond to r = 0.02 and r = 1, respectively.

enter image description here

$\endgroup$
3
  • $\begingroup$ For my case g^2 is not large, so I won't be able to use this. Thanks though. $\endgroup$ Commented Feb 29, 2020 at 5:59
  • $\begingroup$ Even so, doing the y integration symbolically speeds the subsequent numerical computations and may provide insight. I just added a plot for g == 1 to illustrate the plotting process. $\endgroup$
    – bbgodfrey
    Commented Feb 29, 2020 at 6:09
  • $\begingroup$ I see your point. Let me see if I can use this. $\endgroup$ Commented Feb 29, 2020 at 7:49
4
$\begingroup$

if the range of the x2 integration may be shifted to {0,2π} then it can be done and plotted like:

g = 2.; Plot[{NIntegrate[
   E^((Cos[x - x2] + Cos[x2] + Cos[r x + x2])/g^2) x HeavisideTheta[2 Pi - r x], {x, 0, 2 Pi}, {x2, 0, 2 Pi}], 
  2 Pi  NIntegrate [x BesselI[0, 1/g^2 \[Sqrt](3 + 2 Cos[x] + 2 Cos[r x] + 
           2 Cos[(1 + r) x])] HeavisideTheta[2 Pi - r x], {x, 0, 
      2 Pi}] + 10}, {r, 0, 4}]

because

Integrate[Exp[a Cos[x]] Exp[c Sin[x]], {x, 0, 2 Pi}]
is (* 2 Pi BesselI[0,Sqrt[a^2+c^2]]*)

Above is valid if you shift the integration range by any s (for instance s=-Pi):

Integrate[Exp[a Cos[x]] Exp[c Sin[x]], {x, s, 2 Pi+s}]
is (* 2 Pi BesselI[0,Sqrt[a^2+c^2]]*)
$\endgroup$
1
  • $\begingroup$ Thanks, let me see if this works out. $\endgroup$ Commented Mar 1, 2020 at 10:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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