Potential and E-Field of a Uniform Sphere of Charge using Laplace’s and Poisson’s Equations

(I know,....the problem is trivial with Gauss’s Law and Integration from Infinity)

PROBLEM

 A sphere of radius R has a uniform charge density ρ and total charge Q.   We seek the potential and E-Field inside and outside of the sphere and wish to use Laplace’s and Poisson’s Equation (for the shear joy of it).

(a)   Taking advantage of the symmetry of this problem, choose an appropriate coordinate system and write down the appropriate form of Laplace’s OR Poisson’s Equation for finding the potential outside of the sphere.

<First enter answer here in this text cell>

The simplest coordinate system is spherical coordinates with the origin at the center of the sphere.  By the usual arguments, V will be a function of r only.  For r > R, i.e., OUTSIDE the sphere, ρ = 0, so we are concerned with Laplace’s Equation.

We can view the full Laplacian in spherical coordinates by executing this statement:

(*Execute code below*)
Laplacian[f[r,θ,ϕ],{r,θ,ϕ},"Spherical"]//Expand

NOTATION:  Derivative[0,0,2][f][r,θ,ϕ]  in M signifies the 2nd derivative of f with respect to ϕ.  

Thus, the first three terms are derivatives with respect to the Second and Third Variables (θ and ϕ).  The last two terms are derivatives wrt r (the First Variable).  

For f a function of r only, these first 3 terms would yield zero.  Replacing f[r] with Vout[r], we therefore can write Laplace’s Equation for our problem as:

 (2/r)Vout'[r] + Vout''[r] == 0.      (holds for r > R)

(b)  Using DSolve, find a ‘general’ solution for Vout[r]  (before introducing boundary conditions).

(*Input code below*)
ClearAll["`*"]  (*Leave the ClearAll["`*"] statement*)
solout = DSolve[Vout''[r] + (2/r)Vout'[r] == 0, Vout[r],r]
Vout[r_] = Vout[r]/.solout[[1]];
Print["Vout[r] = ",Vout[r]]

As expected (for a second order DE), we get two constants of integration, C[1] and C[2].

(c)  To determine the arbitrary constants C[1] and C[2], we need to employ the boundary conditions.

First consider BC1:  We take our V = 0 point to correspond to the limit r → ∞.  

Using BC1 Determine one of the two constants

<Input derivation in this text cell>

BC1:   V = 0 as r → ∞:   The first term of Vout[r] will go to zero as r → ∞ even for non-zero C[1].  This means we must set C[2] = 0 to guarantee that V→ 0.

(d)  Defining a second BC2:   You should recognize that the remaining term in Vout[r] ∝ (1/r) which is the potential of a point charge.  From earlier work we have shown that a uniformly charged sphere has the potential of a point charge Q where Q is the total charge contained in the sphere. 

Here is the tap dance:  Using V[r] = k (∫_(τ'))ρ(r')(1/(|r - r'|))dτ' we showed that the potential OUTSIDE the sphere due to a spherical shell of charge q is equal to the potential due to a point charge at the center of the sphere.  The filled sphere (constant ρ), by superposition of shells gives a Vout equal to Vout[r] = ((k Q)/r).  

So:  equating our “general solution” to ((k Q)/r), we can find the remaining C[1].  With lightning speed, do so and execute your final definition of Vout[r] to get it defined.   This whole mini-exercise gives us BC2:

<Input derivation in this text cell>
BC2 is essentially    ((k Q)/r) = - ((C[1] )/r) which yields  C[1] = - k Q where k = (1/(4 π ϵo))

Here we use our solution (above) for Vout[r] (after defining Q) to get Vout[r] in a final form.

(*Assuming you use M for the math, input your code below*)
Q = Volumesphere*ρ;
Volumesphere = ((4 )/3)π (R^3);   (*  obviously, you can do this by hand *)
Q
k= (1/(4 π ϵo));

(* Define your final potential function below. I call it VVout[r_]. *)
(* M wants the assignments for these C's to be done using a Rule  (via the -> operator)  *)
VVout[r_] = Vout[r]/.{C[1] -> -(4/3) k π (R^3) ρ,C[2]-> 0}

(*Execute your potential function again below to be absolutely sure it is what you expect*)
VVout[r](* just to check *)

(e)  Now consider inside the sphere (r < R).  For V, we want to use Poisson’s Equation because ρ ≠ 0.  Exploiting the symmetry of the problem, derive the appropriate DE and solve it with DSolve.

<Describe your plan in this text cell>
Poisson’s Equation in general is:  (∇^2)V =- ((ρ )/(ϵo)).    For our problem, where r < R, we can exploit the fact that ρ (a constant) shows spherical symmetry; Vin should not depend on θ or ϕ.  Vin will depend on r (think using Gauss’s law to find E[r].  Integrating E from infinity to r < R will yield a function Vin[r].)  We can then write the rth derivative terms of  (∇^2)Vin[r] =- ((ρ )/(ϵo)) and insert into DSolve.

(*Input code below*)
solin = DSolve[Vin''[r] + (2/r)Vin'[r] == -(ρ/(ϵo)), Vin[r],r]
Vin[r_] = Vin[r]/.solin[[1]]

Again, as expected (for a second order DE), we get two constants of integration, C[1] and C[2].

(f)  To determine the arbitrary constants C[1] and C[2], we need to employ the boundary conditions.

For BC1:  we require that V NOT blow up at r = 0;   So what is C[1]?

For BC2:  we require Vin[R] = Vout[R]   (the continuity of V at an interface with no charge)  Use this to find C[2].

(*Input code to set C[1] to the correct value in Vin[]*)
Vin[r_] = Vin[r]/.{C[1]-> 0}    (*  To keep Vin from blowing up at r = 0 *)
Print["Vin[r] = ",Vin[r]]

(*Input code neccesary to solve for C[2].  Set C[2] in Vin[] to this value*)
solin2 = Solve[Vin[R] == VVout[R], C[2]]
 solin2[[1]]
VVin[r_]= Vin[r]/.%//Simplify

(g) Now plot Vin and Vout in the same graph (suggest you use Show).  You will have to assign values to R, k, ρ, and ϵo.

(*Input code below*)
Clear[R, k, ρ, ϵo]  (*Leave Clear[]statement*)
Print["Vin[r] = ",VVin[r]]
Print["Vout[r] = ", VVout[r]]
R =1;
k = 1;
ρ = 1;
ϵo = 1;
Vplotin = Plot[VVin[r],{r,0,R}];
Vplotout = Plot[VVout[r],{r,R, 6R}];
line =Graphics[Line[{{1,0},{R,0.5}}]];
RR = Graphics[Text["R",{R,0.55}]];
Show[{Vplotin,Vplotout,line,RR}, PlotRange-> {{0,6R},{0,.55}},AxesOrigin-> {0,0},Frame-> True, FrameLabel->{"r", "V[r]"}]

Note:  Although we forced it to happen, it’s important to note that V is continuous at r = R (the vertical line).

(h) Remembering that for a V that is dependent on r only, the only nonzero component of the resulting E field is the radial component.  

Use M’s Grad operator (in spherical coordinates) to generate the E field both inside and outside the sphere.  Plot the magnitude of EEradial versus r.

Remember that for a V dependent on r only, the resulting E field will only have a radial component (derivatives wrt θ and ϕ will be zero and insure that only the r component remains).

(* Input code below *)
Clear[R, k, ρ, ϵo] (* Leave the Clear[] statement *)
EEradialin[r_] = -Grad[VVin[r],{r, θ, ϕ},"Spherical"]
EEradialout[r_] = -Grad[VVout[r],{r, θ, ϕ},"Spherical"]
R =1;
k = 1;
ρ = 1;
ϵo = 1;
EEplotin = Plot[EEradialin[r][[1]],{r,0,R},PlotRange-> {{0,6R},{0,.4}}];
EEplotout = Plot[EEradialout[r][[1]],{r,R, 6R},PlotRange-> {{0,6R},{0,.4}}];
line =Graphics[Line[{{1,0},{R,0.35}}]];
RR = Graphics[Text["R",{R,0.38}]];
Show[{ EEplotin,EEplotout,line, RR}, PlotRange-> {{0,6R},{0,.4}},AxesOrigin-> {0,0},Frame-> True, FrameLabel->{"r", "EE[r]"}]

As expected, (OverVector[E]) = EE[r] (OverHat[r]), both inside and outside the sphere.

Notice that because there is no surface charge at r = R, (σ[r=R] = 0), EE[r] is continuous at r = R:  EE[r] = (Subscript[EE, ⊥])[r], where   (Subscript[EE, ⊥])[r] is the normal component of the field at the interface.

Finally, if you compare these results of using Laplace’s and Poisson’s equation for this simple charge distribution, we get the same potentials and E-fields that we would get by employing Gauss’s Law (for (OverVector[E])) and the integral:  V[r] = – (∫_∞^r) (OverVector[E])[(OverVector[r])]·d(OverVector[r]), where we take V → 0 as r → ∞.