# This integral is divergent. How to use NIntegrate to see how it grows?

I am trying to get information on the following integral: $$\int_{\pi-0.3}^{\pi-\epsilon} \frac{1}{(3+\cos{x})\sqrt{(3+\cos{x})^2-4}}$$

The lower limit is somewhat arbitrary; the point is that this integral is known to be divergent if the interval of integration is from $\pi-\delta$ to $\pi$ (I'm choosing $\delta=0.3$ just as an example).

I need to explore the behavior as $\epsilon$ goes to $0$, so I started making the upper limit closer and closer to $\pi$. However, I stop getting any further improvement when $\epsilon$ reaches $10^{-16}$. If I type the following

NIntegrate[2/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])),{x,Pi-0.3, Pi - 0.0000000000000001}, Method -> "LocalAdaptive",WorkingPrecision -> 100]


then I get the same answer as if I type, say, $10^{-20}$ or $10^{-25}$ instead of $10^{-16}$. Changing Working Precision doesn't have any effect either. The answer I get with the above code is $25.047662\ldots$.

However, as I wrote above, this integral is divergent, so if I make the upper limit closer and closer to $\pi$ I should be getting bigger and bigger numbers; this shouldn't get stuck around $25.047662\ldots$.

I do get a warning from Mathematica, telling me that "...failed to converge to prescribed accuracy after 9 recursive bisections in $x$ near $\{x\} = \{3.14159265...\}\ldots$"

so I know I shouldn't be expecting an accurate answer. I tried to change this behaviour by making Working Precision larger and larger, to no effect.

I would greatly appreciate if someone could point me in the right direction. This calculation takes about 30 seconds in my computer; I wouldn't mind if it takes 10, 20, 60 minutes, as long as I can see the result growing as the upper limit approaches $\pi$.

• Have you tried NDSolve? – Michael E2 Jan 10 '17 at 14:16
• @Michael E2 : Yes, I tried NDSolve, but the problem was that with NDSolve I didn't know how to "see" the results. When I asked it to adjust the domain of the plot over and over, by hand, as the default plot didn't show how the dependent variable was growing and growing – Gelasio Salazar Jan 10 '17 at 14:19
• It does what you want, doesn't it? – Michael E2 Jan 10 '17 at 14:19
• @Michael E2 : Yes, NDSolve at least gives evidence of the growth, but (I must apologize here) I'm pretty new to Mathematica and I don't know how to obtain the solution given by NDSolve in tabular form, so I can manipulate independently (say with GnuPlot, etc.). I spent some time trying to see how to get from Mathematica the solution given by NDSolve into tabular form, but I haven't found the answer – Gelasio Salazar Jan 10 '17 at 14:23
• Did you want a table of the steps taken by NDSolve, the value of the integra at regular intervals, or at $x = \pi - 10^{-k}$, for $k = 1, 2, ...$? – Michael E2 Jan 10 '17 at 14:30

We can calculate an indefinite integral $\int_a f(x) \; dx$ with NDSolve[{y'[x] == f[x], y[a] == 0}, y, {x, a, b}] yields the integral from a to x, up to x <= b or until NDSolve runs into a singularity or stiffness.

{sol} = NDSolve[{y'[x] ==
2/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])), y[Pi - 0.3] == 0},
y, {x, Pi - 0.3, Pi}];


NDSolve::ndsz: At x == 3.141592536731091, step size is effectively zero; singularity or stiff system suspected.

Quick visualization of the steps:

ListPlot[y /. sol, PlotRange -> All] To make a table, let's examine some of the data stored in sol.

dom = y["Domain"] /. sol    (* domain of the solution *)
(*  {{2.8415926535897933, 3.141592536731091}}  *)

xmax = dom[[1, 2]];         (* upper limit of the domain *)

Log10[Pi - xmax]            (* just to show the power of 10 *)
(*  -6.93234  *)


This creates a table of the log of the change in x from Pi and the value of the integral:

tab = Table[{N@k, N[y[Pi - 10^-k] /. sol]}, {k, Range[-Ceiling@Log10[Pi - xmax]]}]
(*
{{1., 0.769312}, {2., 2.39654}, {3., 4.02471},
{4., 5.65288}, {5., 7.28106}, {6., 8.90923}}
*)

ListPlot[tab, Frame -> True, FrameLabel -> {HoldForm[Log10[Δx]], HoldForm[y[π - Δx]]}] Repeating the process with a higher working precision extends the domain closer to Pi. (The step size of the integration process is limited by the precision of the numbers used.)

{sol} = NDSolve[{y'[x] ==
2/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])), y[Pi - 0.3] == 0},
y, {x, Pi - 0.3, Pi},
PrecisionGoal -> 16, WorkingPrecision -> 50] I need to explore the behavior as ϵ goes to 0

I think the behavior is that of $\ln(\pi-x)$ as $x$ approaches $\pi$.

We only need to look at one term.

expr = 1/((3 + Cos[x]) (Sqrt[(3 + Cos[x])^2 - 4]));
Apart[expr] It is only the first term which makes the integral improper at $x=\pi$. All others are ok.

Integrating the first term, but using series form to make it more clear to see (we only want to see what happens at $\pi$)

Normal@Series[Sqrt[5 + 6 Cos[x] + Cos[x]^2]/(8 (1 + Cos[x])), {x, Pi, 6}] res=Integrate[%, x] We see now why it blows up at $x=\pi$ due to the $\ln(\pi-x)$ in the end there. When $x=\pi$ it becomes $\ln(0)$.

 Limit[res, x -> Pi] So the behavior is simply that of $\ln(\pi-x)$ as $x$ approaches $\pi$?

• Thanks a lot to both Michael E2 and Nasser. I had a difficult time deciding which one was the "answer"; I think Michael E2's is "the" answer since it addresses the question (regarding numerical issues). On the other hand Nasser's comment will turn out to be very, very helpful. Thanks a lot. – Gelasio Salazar Jan 10 '17 at 15:13

The exact solution can be obtained in such a way:

j = Integrate[ 1/((3 + Cos[x])*Sqrt[(3 + Cos[x])^2 - 4]), {x, Pi/2, Pi - epsilon},
Assumptions -> epsilon > 0 && epsilon < 1/2];
k = Series[j, {epsilon, 0, 2}, Assumptions -> epsilon > 0 && epsilon < 1/2];
FullSimplify[k]


$$O\left(\text{epsilon}^3\right) +\frac{13 \text{epsilon}^2}{192 \sqrt{2}}+\frac{\log \left(\frac{1}{243} (-32) \left(8 \sqrt{10}-37\right)\right)-2 \log (\text{epsilon})}{4 \sqrt{2}}$$

• you dropped a 2 from the original integral.. fixing that this agrees well with the numeric results. – george2079 Jan 10 '17 at 17:45
• @george2079 : I followed the TEX formula, not the code. So did Nasser. These differ. – user64494 Jan 10 '17 at 18:30

just to show we can do this with NIntegrate , taking care about precision:

Table[
{n, NIntegrate[
2/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])), {x, Pi - .3100,
Pi - 10^(-n)},
WorkingPrecision -> 100, MaxRecursion -> 20] }, {n, 10, 40, 5}]


{{10., 15.4219}, {15., 23.5628}, {20., 31.7037}, {25., 39.8445},

{30., 47.9854}, {35., 56.1263}, {40., 64.2671}}

This throws a working precision warning for n=40, you need to increase working precision with increasing n..

 ListPlot[%] The line on there is -Log[eps]/Sqrt from @user64494's result

BTW The "LocalAdaptive" method seems to malfunction for some reason. Using EvaluationMonitor you see it do a huge number of evals at the wrong end of the interval..

f[x_] := 1/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)]))


It is periodic. Below is a figure of one period.

Plot[f[x], {x, -\[Pi], \[Pi]}] We can compute the integral making some assumptions about δ and ϵ to cause them to lie in the desired range.

Integrate[1/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])),
{x, π - δ, π - ϵ}, Assumptions -> π > δ > ϵ > 0]

(* (1/(4 Sqrt))(-2 ArcSinh[Sqrt[2/3] Cot[δ/2]] +
2 ArcSinh[Sqrt[2/3] Cot[ϵ/2]] +
Log[-(((Sqrt Cot[δ/2] +
2 Sqrt[3 + 2 Cot[δ/2]^2]) (Sqrt Cot[ϵ/2] -
2 Sqrt[3 + 2 Cot[ϵ/2]^2]))/((-Sqrt Cot[δ/
2] + 2 Sqrt[3 + 2 Cot[δ/2]^2]) (Sqrt
Cot[ϵ/2] + 2 Sqrt[3 + 2 Cot[ϵ/2]^2])))]) *)


Define a function to represent this integral:

integralF[δ_, ϵ_] := %


It looks like:

Plot3D[integralF[δ, ϵ], {δ,
0, π}, {ϵ, 0, π}, AxesLabel -> Automatic] You can see that it grows as ϵ approaches zero.

Take the limit for a particular value of δ

Limit[integralF[0.3, ϵ], ϵ -> 0]

(* ∞ *)