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.
Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this communityIntegrate[(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.
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.
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.).
$\endgroup$
May 17 at 17:30
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}} *)
{138.766, {128/(45 \[Pi]), 0.905415}}
.
$\endgroup$
May 18 at 8:03