Can Mathematica calculate this elliptic, triple integral?

Integrate[(4 a b/Pi) (a^2 + b^2 - 2 a b Cos[c])^(1/2), {a, 0, 1}, {b,
0, 1}, {c, 0, Pi}]


I'm using basic plan, it gives me the result like that. The approximate is 0.9054... but I want exact value with closed form.

• I'm afraid it's impossible
– yode
Commented May 17, 2022 at 8:35
• Actually it is possible, see my answer below, and the result is exactly $128/(45\pi)$. Commented May 17, 2022 at 9:14

Doing the Integrals separately for a, b and c gives the answer $$\frac{128}{45 \pi}\sim 0.905415$$, which agrees with the numerical estimate from NIntegrate:

NIntegrate[(4 a b/Pi) (a^2 + b^2 - 2 a b Cos[c])^(1/2), {a, 0, 1}, {b, 0, 1}, {c, 0, Pi}]


0.905415

In a nutshell, the trick is to do the integrals separately as indefinite and then take the limits properly. This is a somewhat straight-forward integral and the trick works, in more general cases, caution is advised!

In more detail: First do the integral over $$a$$:

inta = Integrate[(4 a b/Pi) (a^2 + b^2 - 2 a b Cos[c])^(1/2), a]


(1/(3 [Pi]))b (Sqrt[a^2 + b^2-2 a b Cos[c]] (4 a^2 + b^2 - 2 a b Cos[c] - 3 b^2 Cos[2 c]) + 6 b^3 Cos[c] Log[a - b Cos[c] + Sqrt[a^2 + b^2 - 2 a b Cos[c]]] Sin[c]^2)

and then the one over b:

intb = Integrate[(inta /. a -> 1) - (inta /. a -> 0), b] // PowerExpand


(1/(15 [Pi]))(b^5 (-1 + 3 Cos[2 c]) + Sqrt[1 + b^2 - 2 b Cos[c]] (1 + 8 b^2 + b^4 - 2 (b + b^3) Cos[c] - 3 (1 + b^4) Cos[2 c]) - 6 b^5 Cos[c] Log[b - b Cos[c]] Sin[c]^2 + 6 Cos[c] Log[b - Cos[c] + Sqrt[1 + b^2 - 2 b Cos[c]]] Sin[c]^2 + 6 b^5 Cos[c] Log[1 - b Cos[c] + Sqrt[1 + b^2 - 2 b Cos[c]]] Sin[c]^2)

Then the one over c:

intc = Limit[intb, b -> 1] - (Series[intb, {b, 0, 0}] // FullSimplify // Normal) //  FullSimplify


-(1/(15 [Pi])) 2 (1 - 5 Sqrt[2 - 2 Cos[c]] + 3 (-1 + Sqrt[2 - 2 Cos[c]]) Cos[2 c] + 2 Cos[c] (Sqrt[2 - 2 Cos[c]] + 3 (Log[1 - Cos[c]] - Log[1 + Sqrt[2 - 2 Cos[c]] - Cos[c]]) Sin[c]^2))

Integrate[intc, {c, 0, Pi}]
% // N


128/(45 [Pi])

0.905415

Note: In some cases one has to use Limit or even a series expansion to extract the correct value of the definite integral.

• I'd like to inform that Integrate[(4 a b/Pi) (a^2 + b^2 - 2 a b Cos[c])^(1/2), {a, 0, 1}, {b, 0, 1}, Assumptions -> c > -Infinity] produces a correct huge conditional expression in several hours as NIntegrate[%[[1]], {c, 0, Pi}] confirms with the result 0.905415.Unfortunately, Integrate[%%[[1]], {c, 0, Pi}] returns the input (I think because of Arg[-E^(-I c) (-1 + E^(I c))^2] and similar expressions.). Commented May 17, 2022 at 17:30
• For those who wonder why "caution is advised" when using this method, read this blog post from a Mathematica developer (which I never miss an opportunity to link to.) Commented May 17, 2022 at 18:40
• @MichaelSeifert Great link thanks! I was playing fast and loose with the calculus in the post, hence just the warning ;-) Commented May 17, 2022 at 19:12

The unevaluated integral that Mathematica produced in the OP can be massaged a bit. First, note that due to symmetry, we can instead consider the integral

(16/π) Integrate[a b (a + b) EllipticE[(4 a b)/(a + b)^2], {a, 0, 1}, {b, 0, a}]


and using the conventional definition for EllipticE, we have

(16/π) Integrate[a b (a + b) Sqrt[1 - (4 a b)/(a + b)^2 Sin[t]^2],
{t, 0, π/2}, {a, 0, 1}, {b, 0, a}]


Due again to symmetry, this is equivalent to the integral

(8/π) Integrate[a b (a + b) Sqrt[1 - (4 a b)/(a + b)^2 Sin[t]^2],
{t, -π/2, π/2}, {a, 0, 1}, {b, 0, a}]


and letting w == Sin[t], we have

(8/π) Integrate[a b (a + b) Sqrt[(1 - (4 a b w^2)/(a + b)^2)/(1 - w^2)],
{w, -1, 1}, {a, 0, 1}, {b, 0, a}]


which happens to evaluate pretty quickly to 128/(45 π).

It should be noted that the equivalent integrals

16/π Integrate[a b (a + b) Sqrt[(1 - (4 a b)/(a + b)^2 w^2)/(1 - w^2)],
{w, 0, 1}, {a, 0, 1}, {b, 0, a}]


and

4/π Integrate[a b (a + b) Sqrt[(1 - (4 a b)/(a + b)^2 w^2)/(1 - w^2)],
{w, -1, 1}, {a, 0, 1}, {b, 0, 1}]


do not seem to finish as quickly (I waited ~ 2 minutes for each before giving up.)

If you regard symmetry of integrand (two times integral from 0<b<a) and rearange to do integration over c at last, you get the result in a non splitted integration in 10 seconds.

Original integrations does not work may be because integrate can not handle the case a == b.

f[a_, b_, c_] = (4 a b/Pi) (a^2 + b^2 - 2 a b Cos[c])^(1/2);

f[a, b, c] == f[b, a, c]     (*   True   *)

{(int = 2*Integrate[f[a, b, c], {c, 0, Pi}, {a, 0, 1}, {b, 0, a}]),
int // N} // Timing

(*   {10.141, {128/(45 \[Pi]), 0.905415}}   *)

• +1. Simply and strongly. On my comp {138.766, {128/(45 \[Pi]), 0.905415}}. Commented May 18, 2022 at 8:03