16
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Have you tried NDSolve? $\endgroup$ – Michael E2 Jan 10 '17 at 14:16
  • $\begingroup$ @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 $\endgroup$ – Gelasio Salazar Jan 10 '17 at 14:19
  • $\begingroup$ It does what you want, doesn't it? $\endgroup$ – Michael E2 Jan 10 '17 at 14:19
  • $\begingroup$ @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 $\endgroup$ – Gelasio Salazar Jan 10 '17 at 14:23
  • $\begingroup$ 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, ...$? $\endgroup$ – Michael E2 Jan 10 '17 at 14:30
13
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

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]

Mathematica graphics

$\endgroup$
15
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

res=Integrate[%, x]

Mathematica graphics

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]

Mathematica graphics

So the behavior is simply that of $\ln(\pi-x)$ as $x$ approaches $\pi$?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Gelasio Salazar Jan 10 '17 at 15:13
6
$\begingroup$

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}} $$

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

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 - .3`100, 
    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[%]

enter image description here

The line on there is -Log[eps]/Sqrt[2] 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..

$\endgroup$
3
$\begingroup$

Here is your function

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

Mathematica graphics

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]))(-2 ArcSinh[Sqrt[2/3] Cot[δ/2]] + 
  2 ArcSinh[Sqrt[2/3] Cot[ϵ/2]] + 
  Log[-(((Sqrt[2] Cot[δ/2] + 
       2 Sqrt[3 + 2 Cot[δ/2]^2]) (Sqrt[2] Cot[ϵ/2] - 
       2 Sqrt[3 + 2 Cot[ϵ/2]^2]))/((-Sqrt[2] Cot[δ/
         2] + 2 Sqrt[3 + 2 Cot[δ/2]^2]) (Sqrt[2]
         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]

Mathematica graphics

You can see that it grows as ϵ approaches zero.

Take the limit for a particular value of δ

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

(* ∞ *)
$\endgroup$

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.