Message Boards

WOLFRAM COMMUNITY

6763 Views

4 Replies

1 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Issues on the following numerical integration?

Mike Luntz

Posted 5 years ago

Unless I am totally screwed up, the field inside a surface with uniform charge density should be 0 so the potential should be a constant. But when I attempt to numerically compute the potential interior to a surface I do not get a constant. The integration does return some error messages, but the result of the integration follows a regular shape so I am guessing that the error messages do not have a major impact on the result. The calculation I am attempting is << ComputerArithmetic` re = ImplicitRegion[x^2/4 + y^2 + z^2 == 1.0, {x, y, z}]; v[x0_, y0_] := 1/Sqrt[(x - x0)^2 + (y - y0)^2 + z^2]; pot = Table[{xx, yy, If[yy > Sqrt[1 - xx^2/4], NaN, NIntegrate[v[xx, yy], {x, y, z} \[Element] re]]}, {xx, .05, 1.95, .1}, {yy, .05, 1.95, .1}]; The results of this calculation vary smoothly between a max of about 17 to a min of about 14.7. In addition to this confusing result I think I have detected a bug in the way the integration routine mates with the Table function. Computing NIntegrate[v[.85, .95], {x, y, z} \[Element] re] returns a result. But when part of a table, the arguments .85 and .95 are apparently not passed to the function v. Table[NIntegrate[ v[xx, yy], {x, y, z} \[Element] re], {xx, .75, .95, .1}, {yy, .75, .95, .1}] returns NIntegrate::inumri: The integrand 1/Sqrt[(-0.95+x)^2+(-0.85+y)^2+z^2] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.625,0.75},{0.999999999999999999990527764909997233148869688962838439242343509680,0.999999999999999980708641734607903536052376492842767408867549327610}}. NIntegrate::inumr: The integrand 1/Sqrt[(x-xx)^2+(y-yy)^2+z^2] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1},{0,1}}.

Unless I am totally screwed up, the field inside a surface with uniform charge density should be 0 so the potential should be a constant. But when I attempt to numerically compute the potential interior to a surface I do not get a constant. The integration does return some error messages, but the result of the integration follows a regular shape so I am guessing that the error messages do not have a major impact on the result. The calculation I am attempting is

<< ComputerArithmetic`
re = ImplicitRegion[x^2/4 + y^2 + z^2 == 1.0, {x, y, z}];
v[x0_, y0_] := 1/Sqrt[(x - x0)^2 + (y - y0)^2 + z^2];
pot = Table[{xx, yy, 
    If[yy > Sqrt[1 - xx^2/4], NaN, 
     NIntegrate[v[xx, yy], {x, y, z} \[Element] re]]}, {xx, .05, 
    1.95, .1}, {yy, .05, 1.95, .1}];

The results of this calculation vary smoothly between a max of about 17 to a min of about 14.7.

In addition to this confusing result I think I have detected a bug in the way the integration routine mates with the Table function. Computing

NIntegrate[v[.85, .95], {x, y, z} \[Element] re]

returns a result. But when part of a table, the arguments .85 and .95 are apparently not passed to the function v.

Table[NIntegrate[
  v[xx, yy], {x, y, z} \[Element] 
   re], {xx, .75, .95, .1}, {yy, .75, .95, .1}]

returns

NIntegrate::inumri: The integrand 1/Sqrt[(-0.95+x)^2+(-0.85+y)^2+z^2] has evaluated to Overflow,    Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.625,0.75},{0.999999999999999999990527764909997233148869688962838439242343509680,0.999999999999999980708641734607903536052376492842767408867549327610}}.

NIntegrate::inumr: The integrand 1/Sqrt[(x-xx)^2+(y-yy)^2+z^2] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1},{0,1}}.

POSTED BY: Mike Luntz

4 Replies

Sort By:

Hans Dolhaine

Hans Dolhaine, retired

Posted 5 years ago

Hmmm, I dont't see a mistake, but that says nothing. My Mma (Version 7) doesn't know ImpllicitRegion, so I help myself with Boole. Then Clear[phi] phi[x0_, y0_, z0_] := NIntegrate[Boole[x^2/4 + y^2 + z^2 == 1]/ Sqrt[(x - x0)^2 + (y - y0)^2 + (z - z0)^2], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] and Plot[phi[x, .5, .1], {x, -1, 1}, PlotStyle -> {Thick, Red}] Table[phi[x, .5, .1], {x, -1, 1, .2}] And modifying your code accordingly << ComputerArithmetic` re = Boole[x^2/4 + y^2 + z^2 == 1.0]; v[x0_, y0_] := 1/Sqrt[(x - x0)^2 + (y - y0)^2 + z^2]; pot = Table[{xx, yy, If[yy > Sqrt[1 - xx^2/4], NaN, NIntegrate[v[xx, yy] re, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]]}, {xx, .05, 1.95, .1}, {yy, .05, 1.95, .1}]; pot // MatrixForm the potential is always zero.

Hmmm, I dont't see a mistake, but that says nothing.

My Mma (Version 7) doesn't know ImpllicitRegion, so I help myself with Boole. Then

Clear[phi]
phi[x0_, y0_, z0_] := 
 NIntegrate[Boole[x^2/4 + y^2 + z^2 == 1]/  Sqrt[(x - x0)^2 + (y - y0)^2 + (z - z0)^2],
  {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]

and

Plot[phi[x, .5, .1], {x, -1, 1}, PlotStyle -> {Thick, Red}]
Table[phi[x, .5, .1], {x, -1, 1, .2}]

And modifying your code accordingly

<< ComputerArithmetic`

re = Boole[x^2/4 + y^2 + z^2 == 1.0];
v[x0_, y0_] := 1/Sqrt[(x - x0)^2 + (y - y0)^2 + z^2];
pot = Table[{xx, yy, 
    If[yy > Sqrt[1 - xx^2/4], NaN, NIntegrate[v[xx, yy] re,
      {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]]},
   {xx, .05, 1.95, .1}, {yy, .05, 1.95, .1}];
pot // MatrixForm

the potential is always zero.

POSTED BY: Hans Dolhaine

Hans Dolhaine

Hans Dolhaine, retired

Posted 5 years ago

Sorry. The above method is simply wrong. I think you cannot calculate a surface integral in a manner like this. I think this is due to that NIntegrate choses a lot of points in the region of integration (discretize the region) and almost none fulfills the condition in the Boole-statement. So the result is rubbish. Try to calculate the surface of a (unit) sphere NIntegrate[ Boole[x^2 + y^2 + z^2 == 1.0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] The result is wrong. But when calculating the volume a sufficent amount of points fulfilling the Boole-condition is found and you get a correct result NIntegrate[Boole[x^2 + y^2 + z^2 <= 1.0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] N[4 Pi /3] So forget my proposal.

POSTED BY: Hans Dolhaine

Hans Dolhaine

Hans Dolhaine, retired

Posted 5 years ago

It seems that the potential inside your surface is indeed not constant. I faintly remember that a charge distrubition on a surface with varying curvature is not constant, but I am not sure. Look at ( I exchanged x and z ) xx = {Cos[phi] Sin[th], Sin[phi] Sin[th], 2 Cos[th]}; xx[[1]]^2 + xx[[2]]^2 + xx[[3]]^2/4 // Simplify ParametricPlot3D[xx, {phi, 0, 2 Pi}, {th, 0, Pi}] and e1 = D[xx, th]; e2 = D[xx, phi]; df1 = Cross[e1, e2]; df2 = FullSimplify[Sqrt[df1.df1]] dd1 = xx - {x, y, z}; dd2 = Sqrt[dd1.dd1] // FullSimplify Then pot[x_, y_, z_] := NIntegrate[Sqrt[(5 - 3 Cos[2 th]) Sin[th]^2]/Sqrt[2]/ Sqrt[(z - 2 Cos[th])^2 + (x - Cos[phi] Sin[th])^2 + (y - Sin[phi] Sin[th])^2], {phi, 0, 2 Pi}, {th, 0, Pi}] Obviously the potential calculated in this manner is not constant: ppot = Table[pot[x, 0, 0], {x, -2.5, 2.5, .09}]; ListLinePlot[ppot]

It seems that the potential inside your surface is indeed not constant.

I faintly remember that a charge distrubition on a surface with varying curvature is not constant, but I am not sure.

Look at ( I exchanged x and z )

xx = {Cos[phi] Sin[th], Sin[phi] Sin[th], 2 Cos[th]};

xx[[1]]^2 + xx[[2]]^2 + xx[[3]]^2/4 // Simplify

ParametricPlot3D[xx, {phi, 0, 2 Pi}, {th, 0, Pi}]

and

e1 = D[xx, th];
e2 = D[xx, phi];
df1 = Cross[e1, e2];
df2 = FullSimplify[Sqrt[df1.df1]]

dd1 = xx - {x, y, z};
dd2 = Sqrt[dd1.dd1] // FullSimplify

Then

pot[x_, y_, z_] :=
 NIntegrate[Sqrt[(5 - 3 Cos[2 th]) Sin[th]^2]/Sqrt[2]/
  Sqrt[(z - 2 Cos[th])^2 + (x - Cos[phi] Sin[th])^2 + (y - 
     Sin[phi] Sin[th])^2],
  {phi, 0, 2 Pi}, {th, 0, Pi}]

Obviously the potential calculated in this manner is not constant:

ppot = Table[pot[x, 0, 0], {x, -2.5, 2.5, .09}];
ListLinePlot[ppot]

POSTED BY: Hans Dolhaine

Mike Luntz

Posted 5 years ago

Thanks Hans. I finally found some references to this that confirm that the potential is not constant.

POSTED BY: Mike Luntz

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback