Symbolic integration of potential over a disc : branch cut problem?
up vote
11
down vote
favorite
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
asked Nov 28 at 7:06
chris
12.2k440108
add a comment |
up vote
11
down vote
favorite
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
asked Nov 28 at 7:06
chris
12.2k440108
2
Indeed,Integrateappears to have a problem with the repulsive part;Integrate[ 1/((x - x0)^2 + (y - y0)^2), {x0, y0} [Element] Disk[{0, 0}, 1]]returns0which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.
– Henrik Schumacher
Nov 28 at 8:16
add a comment |
up vote
11
down vote
favorite
up vote
11
down vote
favorite
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
asked Nov 28 at 7:06
chris
12.2k440108
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
asked Nov 28 at 7:06
chris
12.2k440108
asked Nov 28 at 7:06
chris
12.2k440108
edited Nov 28 at 8:19
asked Nov 28 at 7:06
chris
12.2k440108
asked Nov 28 at 7:06
chris
12.2k440108
asked Nov 28 at 7:06
chris
12.2k440108
12.2k440108
2
Indeed,Integrateappears to have a problem with the repulsive part;Integrate[ 1/((x - x0)^2 + (y - y0)^2), {x0, y0} [Element] Disk[{0, 0}, 1]]returns0which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.
– Henrik Schumacher
Nov 28 at 8:16
add a comment |
2
Indeed,Integrateappears to have a problem with the repulsive part;Integrate[ 1/((x - x0)^2 + (y - y0)^2), {x0, y0} [Element] Disk[{0, 0}, 1]]returns0which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.
– Henrik Schumacher
Nov 28 at 8:16
2
2
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.– Henrik Schumacher
Nov 28 at 8:16
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.– Henrik Schumacher
Nov 28 at 8:16
add a comment |
1 Answer
1
active
oldest
votes
up vote
8
down vote
accepted
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 at 8:53
chris
12.2k440108
answered Nov 28 at 8:12
David
1035
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}]]
– chris
Nov 28 at 8:31
@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.
– David
Nov 28 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',
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
up vote
8
down vote
accepted
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 at 8:53
chris
12.2k440108
answered Nov 28 at 8:12
David
1035
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}]]
– chris
Nov 28 at 8:31
@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.
– David
Nov 28 at 10:33
add a comment |
up vote
8
down vote
accepted
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 at 8:53
chris
12.2k440108
answered Nov 28 at 8:12
David
1035
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}]]
– chris
Nov 28 at 8:31
@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.
– David
Nov 28 at 10:33
add a comment |
up vote
8
down vote
accepted
up vote
8
down vote
accepted
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 at 8:53
chris
12.2k440108
answered Nov 28 at 8:12
David
1035
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 at 8:53
chris
12.2k440108
answered Nov 28 at 8:12
David
1035
edited Nov 28 at 8:53
chris
12.2k440108
edited Nov 28 at 8:53
chris
12.2k440108
edited Nov 28 at 8:53
chris
12.2k440108
12.2k440108
answered Nov 28 at 8:12
David
1035
answered Nov 28 at 8:12
David
1035
answered Nov 28 at 8:12
David
1035
1035
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}]]
– chris
Nov 28 at 8:31
@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.
– David
Nov 28 at 10:33
add a comment |
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}]]
– chris
Nov 28 at 8:31
@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.
– David
Nov 28 at 10:33
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}]]– chris
Nov 28 at 8:31
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}]]– chris
Nov 28 at 8:31
@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.
– David
Nov 28 at 10:33
@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.
– David
Nov 28 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.
Some of your past answers have not been well-received, and you're in danger of being blocked from answering.
Please pay close attention to the following guidance:
- 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.
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
Indeed,
Integrateappears to have a problem with the repulsive part;Integrate[ 1/((x - x0)^2 + (y - y0)^2), {x0, y0} [Element] Disk[{0, 0}, 1]]returns0which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support.– Henrik Schumacher
Nov 28 at 8:16