Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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. enter image description here 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?
share|improve this question
    
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. –  J. M. 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... –  J. M. Jul 5 '12 at 17:12

4 Answers 4

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] *)
share|improve this answer
    
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

My version of Mathematica has no problems calculating analytically the integral

$\int_{x-b}^{b+x} 2 r (l-r) \cos ^{-1}\left(\frac{s^2+x^2-b^2}{2 s x}\right) \, dr$

The result looks nice (IMHO) rather than "ugly".

$\frac{1}{9} b^3 \left(9 L \pi -2 (2+z) (8+z (2+z)) \text{EllipticE}\left[\frac{4 (1+z)}{(2+z)^2}\right]+2 z (2+z)^2 \text{EllipticK}\left[-\frac{4 (1+z)}{z^2}\right]\right)/.\left\{z\to \frac{x}{b},L\to \frac{l}{b}\right\}$

$Version
"8.0 for Microsoft Windows (64-bit) (October 7, 2011)"

The idea is to remove as many symbols from the expression as possile by making it dimensionless, i.e. letting s -> r/b, z -> x/b-1, L-> l/b, the transformed integral is easily and quickly calculated by MMA:

g[z_, L_] = 
 b^3 Integrate[
   2 ArcCos[(s^2 + z (2 + z))/(2 s (1 + z))] s (L - s), {s, z, z + 2}, 
   Assumptions -> {z > 0, L > 0}]

(*  1/9 b^3 (9 L \[Pi] - 2 (2 + z) (8 + z (2 + z)) EllipticE[(4 (1 + z))/(2 + z)^2] + 
   2 z (2 + z)^2 EllipticK[-((4 (1 + z))/z^2)]) *)

Checking with NIntegrate shows that the result is correct.

By the way, the g as a function of z looks very simple. But observe that the two elliptic integrals don't. Obviously, their non-trivial parts cancel in a "peaceful" manner.

Best regards, Wolfgang

share|improve this answer

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

halfDisc

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 *)
share|improve this answer
    
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:

Mathematica graphics

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.

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.