Symbolic integration of potential over a disc : branch cut problem?
$begingroup$
Context
I am trying to explore the geometry of a crystal made of irregular bubbles.
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}]
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
calculus-and-analysis numerical-integration bugs symbolic
$endgroup$
add a comment |
$begingroup$
Context
I am trying to explore the geometry of a crystal made of irregular bubbles.
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}]
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
calculus-and-analysis numerical-integration bugs symbolic
$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]]
returns0
which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.
$endgroup$
– Henrik Schumacher
Nov 28 '18 at 8:16
add a comment |
$begingroup$
Context
I am trying to explore the geometry of a crystal made of irregular bubbles.
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}]
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
calculus-and-analysis numerical-integration bugs symbolic
$endgroup$
Context
I am trying to explore the geometry of a crystal made of irregular bubbles.
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}]
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
calculus-and-analysis numerical-integration bugs symbolic
calculus-and-analysis numerical-integration bugs symbolic
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]]
returns0
which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.
$endgroup$
– Henrik Schumacher
Nov 28 '18 at 8:16
add a comment |
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]]
returns0
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
add a comment |
1 Answer
1
active
oldest
votes
$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.
$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
add a comment |
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
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
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
$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.
$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
add a comment |
$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.
$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
add a comment |
$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.
$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.
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
add a comment |
$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
add a comment |
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.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
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
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
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
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]]
returns0
which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.$endgroup$
– Henrik Schumacher
Nov 28 '18 at 8:16