# How far should I trust Mathematica's Integrator?/ Why do I get two inconsistent answers?

I am trying to solve the following integral: $$\int_0^{2\pi}\frac{\cos(2\theta)}{\sqrt{d^2+2\lambda^2-d^2\cos{2\theta}} (d^2+2\lambda^2+2\lambda^2\cos{2\theta})}\text{d}\theta$$. To do that I can find its primitive and evaluate the difference at the extremes:

$Assumptions = {θ ∈ Reals, d ∈ Reals, λ ∈ Reals}; T = Cos[2θ]/(Sqrt[d^2 + 2 λ^2 -d^2 Cos[2θ]] (d^2 + 2 λ^2 + 2 λ^2 Cos[2θ])); pr = Integrate[T, θ]; Int = pr/. θ-> 2π - pr/. θ-> 0  This gives as a result$0$which means for any$d$,$\lambda$. However if I try a numerical integration for two given values of$d$,$\lambda$, say$d=\lambda=1$, which I should be able to compute as NIntegrate[T /. d-> 1/.λ-> 1, {θ, 0, 2π}, AccuracyGoal -> Infinity]  I get quite a different result ($-0.477437$in this case). Can anyone tell me which result should I trust or what I am doing wrong? Interestingly, when the same procedure above is used to compute other integrals the results are in agreement. • Have you plotted the primitive? I suspect you'll find that it has discontinuities. – mmeent Sep 15 '17 at 21:20 • @semola: No, the numerical result is most likely correct, the analytic one is wrong. For example, consider a constant function - one possible "primitive" is the sawtooth function with period 1 (which also has discontinuities), and anytime the integration limits differ by an integer, you get 0 as a result for the integral, which is obviously wrong. I hope the example is clear. – Lukas Lang Sep 15 '17 at 21:31 • The "symbolic" result cannot be correct. This is because the Newton-Leibniz method of evaluating an antiderivative at the endpoints will not apply if there are path singularities. – Daniel Lichtblau Sep 15 '17 at 21:31 • Plotting a surface over values of$d$and$\lambda$might suggest how complex the solution (if one exists) might be: surface = Flatten[Table[{d, \[Lambda], NIntegrate[T, {\[Theta], 0, 2 \[Pi]}]}, {d, 1/10, 5, 1/10}, {\[Lambda], 1/10, 5, 1/10}], 1]; ListPlot3D[surface, PlotRange -> {Automatic, Automatic, {-1, 1}}]. – JimB Sep 15 '17 at 22:06 • Please see this blog post about this very issue. – march Sep 15 '17 at 22:55 ## 1 Answer In the comments, the OP requested a closed-form solution. This approach appears to provide one. It uses the Rubi integration package for Mathematica (available at rulebasedintegration.org). Note that the Rubi integration command is "Int": T = Cos[2 θ]/(Sqrt[d^2 + 2 λ^2 - d^2 Cos[2 θ]](d^2 + 2 λ^2 + 2 λ^2 Cos[2 θ])); pr = Int[T, θ] $\frac{\sqrt{\frac{d^2 (-\cos (2 \theta ))+d^2+2 \lambda ^2}{\lambda ^2}} F\left(\theta \left|-\frac{d^2}{\lambda ^2}\right.\right)}{2 \sqrt{2} \lambda ^2 \sqrt{d^2 (-\cos (2 \theta ))+d^2+2 \lambda ^2}}-\frac{\left(d^2+2 \lambda ^2\right) \sqrt{\frac{d^2 (-\cos (2 \theta ))+d^2+2 \lambda ^2}{\lambda ^2}} \Pi \left(\frac{4 \lambda ^2}{d^2+4 \lambda ^2};\theta \left|-\frac{d^2}{\lambda ^2}\right.\right)}{2 \sqrt{2} \lambda ^2 \left(d^2+4 \lambda ^2\right) \sqrt{d^2 (-\cos (2 \theta ))+d^2+2 \lambda ^2}}$[$F$is EllipticF,$\Pi$is EllipticPi and$K$, seen in the next result, is EllipticK] Rubi only calculates indefinite integrals. To obtain the definite integral, we substitute the limits of integration: SymbInt = pr /. θ -> 2 π - pr /. θ -> 0 $\frac{\sqrt{2} K\left(-\frac{d^2}{\lambda ^2}\right)}{\left(\lambda ^2\right)^{3/2}}-\frac{\sqrt{2} \left(d^2+2 \lambda ^2\right) \Pi \left(\frac{4 \lambda ^2}{d^2+4 \lambda ^2}|-\frac{d^2}{\lambda ^2}\right)}{\left(\lambda ^2\right)^{3/2} \left(d^2+4 \lambda ^2\right)}\$

Using JimB's nice code from the comment above for the numeric result from NIntegrate, and Plot3D for the symbolic result from Rubi (SymbInt), we can plot and compare the two for a range of values of d and λ. Here is the code used to plot the surfaces resulting from NIntegrate and Rubi, respectively:

surfaceN = Flatten[Table[{d, λ, NIntegrate[T, {θ, 0, 2 π}]}, {d, 1/10, 5, 1/10}, {λ, 1/10, 5, 1/10}], 1];
ListPlot3D[surfaceN, PlotRange -> {{1/10, 5}, {1/10, 5}}]
Plot3D[SymbInt, {d, 1/10, 5}, {λ, 1/10, 5}, PlotRange ->{{1/10, 5}, {1/10, 5}}]


The graphical results appear to be identical over this set of values: One can also compare the results numerically:

TN = Table[NIntegrate[T, {θ, 0, 2 π}], {d,1/10, 5, 1/10}, {λ, 1/10, 5, 1/10}];
TS = Table[SymbInt, {d, 1/10, 5, 1/10}, {λ,1/10, 5, 1/10}] // N;
diff = TS - TN;
{Max@diff, Min@diff}


{3.0914*10^-7, -5.49857*10^-7}

Giving NIntegrate a higher WorkingPrecision brings its results closer to those of the symbolic form:

TN1 = Table[NIntegrate[T, {θ, 0, 2 π}, MinRecursion -> 2, WorkingPrecision -> 20], {d, 1/10, 5, 1/10}, {λ, 1/10, 5, 1/10}];
diff1 = TS - TN1;
{Max@diff1, Min@diff1}


{5.53158*10^-12, -7.72094*10^-12}

The OP will probably want to check the above approach for more general validity.

• I have changed the URL to our new site for Rubi as Albert's site is going offline soon. – halirutan Aug 13 '18 at 17:32