Symbolic integration of potential over a disc : branch cut problem?












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










share|improve this question











$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


















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










share|improve this question











$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
















12












12








12


1



$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










share|improve this question











$endgroup$




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







calculus-and-analysis numerical-integration bugs symbolic






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited Nov 28 '18 at 8:19







chris

















asked Nov 28 '18 at 7:06









chrischris

12.4k442113




12.4k442113








  • 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
















  • 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










2




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






$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












1 Answer
1






active

oldest

votes


















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-1right)+4 log (r)+1right)$$



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.






share|improve this answer











$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





StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "387"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});


}
});














draft saved

draft discarded


















StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f186859%2fsymbolic-integration-of-potential-over-a-disc-branch-cut-problem%23new-answer', 'question_page');
}
);

Post as a guest















Required, but never shown

























1 Answer
1






active

oldest

votes








1 Answer
1






active

oldest

votes









active

oldest

votes






active

oldest

votes









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-1right)+4 log (r)+1right)$$



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.






share|improve this answer











$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
















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-1right)+4 log (r)+1right)$$



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.






share|improve this answer











$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














9












9








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-1right)+4 log (r)+1right)$$



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.






share|improve this answer











$endgroup$



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-1right)+4 log (r)+1right)$$



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.







share|improve this answer














share|improve this answer



share|improve this answer








edited Nov 28 '18 at 8:53









chris

12.4k442113




12.4k442113










answered Nov 28 '18 at 8:12









DavidDavid

1135




1135












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


















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
















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




$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


















draft saved

draft discarded




















































Thanks for contributing an answer to Mathematica Stack Exchange!


  • Please be sure to answer the question. Provide details and share your research!

But avoid



  • Asking for help, clarification, or responding to other answers.

  • Making statements based on opinion; back them up with references or personal experience.


Use MathJax to format equations. MathJax reference.


To learn more, see our tips on writing great answers.




draft saved


draft discarded














StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f186859%2fsymbolic-integration-of-potential-over-a-disc-branch-cut-problem%23new-answer', 'question_page');
}
);

Post as a guest















Required, but never shown





















































Required, but never shown














Required, but never shown












Required, but never shown







Required, but never shown

































Required, but never shown














Required, but never shown












Required, but never shown







Required, but never shown







Popular posts from this blog

A CLEAN and SIMPLE way to add appendices to Table of Contents and bookmarks

Calculate evaluation metrics using cross_val_predict sklearn

Insert data from modal to MySQL (multiple modal on website)