12
$\begingroup$

Context

I am trying to explore the geometry of a crystal made of irregular bubbles.

enter image description here

See animation here.

very vaguely in the spirit of this post (it is in fact motivated by cosmology and galaxy formation).

So I give myself an interaction potential (which is both attractive and repulsive at large and small distances resp.)

pot[r_] = 1/r^2 + r^2

looking like this

Plot[pot[r], {r, 0.1, 5}]

enter image description here

and I integrate it over a Disk

int= Integrate[ pot[Sqrt[(x - x0)^2 + (y - y0)^2]], {x0, y0} ∈ 
  Disk[{0, 0}, 1]]

(* π (x^2+y^2+1/2)  *)

which incidentally looks suspicious, because it is lacking a repulsion near the disc.

But if I take a specific value for {x,y}

rxy = Thread[{x, y} -> {2, 3}]

and carry out the integration numerically

NIntegrate[
 pot[Sqrt[(x - x0)^2 + (y - y0)^2]] /. rxy, {x0, y0} ∈ 
  Disk[{0, 0}, 1], PrecisionGoal -> 6]

(* 42.663 *)

I get a different answer from

  int /. rxy 

(* 42.4115 *)

Indeed if I do the replacement First

Integrate[pot[Sqrt[(x - x0)^2 + (y - y0)^2]] /. rxy, {x0, y0} ∈ Disk[{0, 0}, 1]] 

(* π (27/2+log(13/12)) *)

N[%]

(* 42.663 *)

So mathematica seems to be doing the general integration wrong.

Questions

Is this a bug? Any workaround?

Check

Indeed I can check by integrating numerically radially away from the edge of the disk that the potential generated by the disc is repulsive at close distance:

dat = ParallelTable[
  NIntegrate[
   pot[Sqrt[(x - x0)^2 + (y - y0)^2]] /. {x -> r Cos[t], 
      y -> r Sin[t]} /. t -> Pi/4, {x0, y0}∈ 
    Disk[{0, 0}, 1], PrecisionGoal -> 8],
  {r, 1.01, 2, 0.025}];
  dat // ListLinePlot

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Indeed, Integrate appears to have a problem with the repulsive part; Integrate[ 1/((x - x0)^2 + (y - y0)^2), {x0, y0} \[Element] Disk[{0, 0}, 1]] returns 0 which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support. $\endgroup$ – Henrik Schumacher Nov 28 '18 at 8:16
9
$\begingroup$

It's worth noting that the integrals will evaluate separately!

totalPot[x_]=Integrate[(x-x0)^2+(y0)^2,{x0,y0}∈ Disk[{0,0},1]]+
Integrate[1/((x-x0)^2+(y0)^2),{x0,y0}∈ Disk[{0,0},1],Assumptions->{x>1}];
N[totalPot[Sqrt[2^2 + 3^2]]]
(* 42.663 *)

The exact form of the potential being

$$\frac{1}{2} \pi \left(2 r^2-2 \log \left(r^2-1\right)+4 \log (r)+1\right)$$

Where I made sure to use the manifest rotational symmetry to put y=0, and also added an assumption that x is greater than 1 to avoid any issues with divergences in the 1/r^2 case.

Since it's of physical interest, to put units back in, if I take the potential to be an energy density $k_1 r^2+k_2/r^2$ and the disk is of radius $R$, I find:

$$E(r)=k_1 \frac{\pi}{2}(R^4+2 R^2 r^2)-k_2 \pi \log(1-\frac{R^2}{r^2})$$

As noted by Henrik in the comments, this looks like a bug & should be reported to wolfram support.

$\endgroup$
  • $\begingroup$ thanks! Would you know how to do the integral over an elliptic disc? int= Integrate[ pot[Sqrt[(x - x0)^2 + (y - y0)^2]], {x0, y0} ∈ Disk[{0, 0}, {1,2}]] $\endgroup$ – chris Nov 28 '18 at 8:31
  • $\begingroup$ @chris hmm I don't. In 3D for the 1/r potential and an ellipsoid I know that the result is pretty complicated with no nice answer. But maybe it's easier here. $\endgroup$ – David Nov 28 '18 at 10:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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