# Can't compute definite integral

Consider a scalar field (in polar co-ordinates), $f(r) = l-r$.

Now I want to evaluate the field integral over a circular region of radius $b$, centered at a distance of $x$ from the origin. According the law of cosines, The arc of radius $r$ from the origin subtends an angle of $\theta = 2 \arccos(\frac{r^2+x^2-b^2}{2 r x})$ in the circular region. Assuming $x>b$, the integral over the region should be: $$\int_{x-b}^{x+b} r f(r) 2\cos^{-1}[\frac{r^2+x^2-b^2}{2 r x} ] dr$$ Or, $$\int_{x-b}^{x+b} r (l-r) 2\cos^{-1}[\frac{r^2+x^2-b^2}{2 r x} ] dr$$

Since this is a real valued integral with real limits, we should be able to compute the integral. In Mathematica, I have:

Integrate[
2 ArcCos[(r^2 + x^2 - b^2)/(2 r x)] r (l - r), {r, x - b, x + b},
Assumptions -> l > 0 &&  x > b && b > 0 ]


which never evaluates!

If, however, I try to do the indefinite integral:

expr = Integrate[2 ArcCos[(r^2 + x^2 - b^2)/(2 r x)] r (l - r), r,
Assumptions -> l > 0 &&  x > b && b > 0 && x - b <= r <= x + b ]


I get the answer

1/18 r ((-9 l + 4 r) x Sqrt[-((
b^4 + (r^2 - x^2)^2 - 2 b^2 (r^2 + x^2))/(r^2 x^2))] -
6 r (-3 l + 2 r) ArcCos[(-b^2 + r^2 + x^2)/(
2 r x)]) - (2 I (b - x)^2 (7 b^2 + x^2) Sqrt[(
b^2 - r^2 - 2 b x + x^2)/(b - x)^2] Sqrt[(
b^2 - r^2 + 2 b x + x^2)/(b + x)^2]
EllipticE[
I ArcSinh[r Sqrt[-(1/(b + x)^2)]], (b + x)^2/(b - x)^2] +
b (9 b l Sqrt[-(1/(b + x)^2)]
Sqrt[-b^4 - (r^2 - x^2)^2 + 2 b^2 (r^2 + x^2)]
ArcTan[(b^2 - r^2 + x^2)/
Sqrt[-b^4 - (r^2 - x^2)^2 + 2 b^2 (r^2 + x^2)]] -
4 I (b - x)^2 (3 b - x) Sqrt[(
b^2 - r^2 - 2 b x + x^2)/(b - x)^2] Sqrt[(
b^2 - r^2 + 2 b x + x^2)/(b + x)^2]
EllipticF[
I ArcSinh[r Sqrt[-(1/(b + x)^2)]], (b + x)^2/(b -
x)^2]))/(9 r x Sqrt[-(1/(b + x)^2)]
Sqrt[-((b^4 + (r^2 - x^2)^2 - 2 b^2 (r^2 + x^2))/(r^2 x^2))])


which looks ugly, but since there should be no discontinuities, I should be able to evaluate the definite integral as

(expr /. {r -> x + b}) - (expr /. {r -> x - b})


Which gives another ugly answer. Now, if I try to evaluate for some numerical values:

N[((expr /. {r -> x + b}) - (expr /. {r -> x - b})) /. {x -> 5,
b -> 1, l -> 10}]


I get errors such as

Power::infy: Infinite expression 1/Sqrt[0] encountered.


Questions:

1. Why can't we evaluate the definite integral directly?
2. Since the integrand is real and finite, why does the integral give infinity expressions?
-
FWIW: it is my opinion that Mathematica's handling of elliptic integrals is very far from optimal. More often than not, I've managed to produce better results from manual evaluation than from letting Mathematica lose at my integral. –  Ｊ. Ｍ. Jul 5 '12 at 2:11
Keep in mind that ReplaceAll is not a mathematical operation, but just pure term replacement. To see the difference try Limit[ArcTan[(b^2-r^2+x^2)/Sqrt[-b^4-(r^2-x^2)^2+2b^2(r^2+x^2)]],r->x-b] and Simplify[(b^2-r^2+x^2)/Sqrt[-b^4-(r^2-x^2)^2+2b^2(r^2+x^2)]]/.r->x-b] (where the argument is the term pointed out by @SjoerdC.deVries below). I think (but haven't tried) that by taking the limit you will be able to evaluate the resulting expression. However, the approaches outlined by Simon Woods and b.gatessucks below are most probably better. –  sebhofer Jul 5 '12 at 8:17
You might be able to get a result by using Limit instead of ReplaceAll. But it could be incorrect if there are jump discontinuities in the indefinite integral. This can happen even if integrand is continuous along integration path, reason being integrand may have branch points. To make matters worse, presence or absence of branch cuts may depend on actual values later used for your symbolic parameters. Which quite possibly is why Integrate punted on that definite integral. –  Daniel Lichtblau Jul 5 '12 at 13:52
@J.M., perhaps that means Mathematica can improve it. Do they take suggestions like this? Does someone from Wolfram read this? –  highBandWidth Jul 5 '12 at 17:08
Maybe, but bear in mind that what's easily done by hand is not always easy to do for a computer following a strict algorithm, and vice-versa... –  Ｊ. Ｍ. Jul 5 '12 at 17:12

## 3 Answers

If the circle is centered at {x0, y0} and has a radius b then it can be parametrized with

$$x = x0 + b \cos(\theta')$$

$$y = y0 + b \sin(\theta')$$

$$\theta' \in [0, 2 \pi]$$

Then your integral is

$$int = \int_{0}^{2 \pi} (l- \sqrt{(x0 + b \cos(\theta'))^2 + (y0 + b \sin(\theta'))^2}) d\theta '$$

int = Integrate[(l - Sqrt[(x0 + b Cos[tp])^2 + (y0 + b Sin[tp])^2 // Expand]), {tp, 0, 2 \[Pi]}, Assumptions -> {x0 > 0, y0 > 0, b > 0}]

(* ConditionalExpression[2 l \[Pi]-4 Sqrt[b^2+x0^2+y0^2-2 b Sqrt[x0^2+y0^2]] EllipticE[-((4 b Sqrt[x0^2+y0^2])/(b^2+x0^2+y0^2-2 b Sqrt[x0^2+y0^2]))],Sqrt[x0^2+y0^2]!=b] *)


The simplest check seems ok :

int /. {x0 -> 0, y0 -> 0}

(* ConditionalExpression[-2 Sqrt[b^2] \[Pi] + 2 l \[Pi], b != 0] *)

-
You have integrated along the circle boundary, rather than the circular disk. Perhaps we need a second integral Integrate[ . , {r,0,b}]? –  highBandWidth Jul 5 '12 at 17:06
@highBandWidth Added another answer, if that is what you need I'll delete this one. –  b.gatessucks Jul 5 '12 at 20:03

I don't have an answer for your questions; below is just a way to calculate the integral.

You want to integrate $f(r)=l-r$ over the disk of radius $b$ centered at $\{x_0,y_0\}$. If you imagine a line from the origin through the center of the disk then, since $f(r)$ is symmetric about that line, we can consider the top half-disk only and multiply the corresponding integral by $2$.

The top half-disk is described by $\{x = x0 + r' \cos(\theta'), y = y0 + r' \sin(\theta ')\}$ (thanks to Heike who pointed this out to me) where

$$r' \in \left[0, b\right]$$ $$\theta' \in [\theta_0, \theta_0 + \pi]$$ $$\tan(\theta_0) = \frac{y0}{x0}$$

Your integral is then

$$\int_{Disc(\{ x0, y0\}, b)} d\theta \ dr \ r \ f(r) = 2 \int_{\theta_0}^{\theta_0+\pi} d\theta' \ \int_{0}^{b} dr' r |J| (l - r) \equiv one - two$$

where $$J = \frac{1}{r} (r' + x0 \cos(\theta') + y0 \sin(\theta'))$$ $$r = \sqrt{(x0 + r' \cos(\theta'))^2 + (y0 + r' \sin(\theta '))^2}$$

th0 = Arctan[x0, y0];

one = 2 Integrate[ l (rp + x0 Cos[tp] + y0 Sin[tp]), {tp, th0, th0 + \[Pi]},
{rp, 0, b}] /. {y0 Cos[th0] -> x0 Sin[th0]}

(* b^2 l \[Pi] *)

aux = Integrate[(1 + \[Beta]^2 + 2 \[Beta] Cos[tp - \[Alpha]])^(3/2), tp,
Assumptions -> \[Beta] \[Element] Reals]

auxTwo[\[Alpha]_, \[Beta]_] = Simplify[(aux /. (tp -> \[Alpha] + \[Pi])) - (aux /. (tp -> \[Alpha])), Assumptions -> {\[Alpha] \[Element] Reals, \[Beta] \[Element] Reals}]

two = Simplify[2/3 ((x0^2 + y0^2)^(3/2) (auxTwo[th0, b/Sqrt[x0^2 + y0^2]]) - \[Pi]
(x0^2 + y0^2)^(3/2)), Assumptions -> {b > 0, x0 \[Element] Reals, y0 \[Element] Reals}]

result = one - two /. {y0 Cos[Arctan[x0, y0]] -> x0 Sin[Arctan[x0, y0]]}

(* b^2 l \[Pi] - 2/9 (x0^2 + y0^2)^(3/2) (-3 \[Pi] + 2 (1 + b/Sqrt[x0^2 + y0^2]) (4 (1 + b^2/(x0^2 + y0^2)) EllipticE[(
4 b Sqrt[x0^2 + y0^2])/(b + Sqrt[x0^2 + y0^2])^2] - (-1 + b/
Sqrt[x0^2 + y0^2])^2 EllipticK[(
4 b Sqrt[x0^2 + y0^2])/(b + Sqrt[x0^2 + y0^2])^2])) *)


One simple check is that it depends only on x0^2 + y0^2 as it should by symmetry. Using this we can simplify the result; with x0^2 + y0^2 -> r0^2 :

Simplify[result /. {x0^2 + y0^2 ->  r0^2}, Assumptions -> {r0 > 0}]

(* b^2 l \[Pi]-2/9 r0^3 (-3 \[Pi]+(2 (b+r0) (4 (b^2+r0^2) EllipticE[(4 b r0)/(b+r0)^2]-
(b-r0)^2 EllipticK[(4 b r0)/(b+r0)^2]))/r0^3) *)


Another simple check is to compare to the case when the circle is centered at the origin :

Integrate[r (l - r), {r, 0, b}, {\[Theta], 0, 2 \[Pi]}] - Limit[Limit[result, y0 -> 0], x0 -> 0]

(* 0 *)

-
Your domain of integration doesn't actually describe a circle. Try for example With[{x0 = 5, y0 = 0, b = 1}, ParametricPlot[{r Cos[th], r Sin[th]}, {r, Sqrt[x0^2 + y0^2] - b, Sqrt[x0^2 + y0^2] + b}, {th, 0, ArcSin[b/Sqrt[x0^2 + y0^2]]}]] to see what it looks like. To describe the circle, the bounds for r would have to depend on theta (or vice versa if you swap the order of integration). –  Heike Jul 6 '12 at 7:09
@Heike You are right, I messed up my notation. Will fix it soon. Many thanks for checking. –  b.gatessucks Jul 6 '12 at 8:24

If I do a FullSimplify on the output of the indefinite Integrate I get this:

Please examine the part of this expression marked in red and you'll see why the substitution of r -> x + b and r -> x - b will cause some problems.

-
Thanks for the insight! Now if we go back to the problem however, the integrand is always finite, even if $r=x+b$ or $r=x-b$. So shouldn't the integral also be finite? I think the terms in red actually cancel with the other such terms in the limit. How should this be evaluated? –  highBandWidth Jul 4 '12 at 17:48
@highBandWidth Well Limit[(b^2 - r^2 + x^2)/( Sqrt[b - r - x] Sqrt[(-b - r + x) (b - r + x) (b + r + x)]), r -> x + b, Assumptions -> {x > b && b > 0}] gives me DirectedInfinity[I]. r -> x - b yields -Infinity. –  Sjoerd C. de Vries Jul 4 '12 at 18:55
@SjoerdC.deVries But of course the limit of the ArcTan does exist. So your comment may be misleading... –  sebhofer Jul 5 '12 at 8:08