Electric Potential of a Uniformly Charged Disk of Charge—Off Axis

A disk of radius R normal to the z axis centered at the origin (i.e., lying in the x-y plane) holds a uniform charge density σ;  Find and plot Vfar and Vnear — the off-axis solutions for z > 0.

HINTS:  
(i) Treat as a 2D problem.  Any plane through the z-axis will do — take the x-z plane (i.e., we consider y = 0).  This is equivalent to saying that there is no ϕ dependence.  Since any choice of ϕ will yield the same result, we can choose ϕ = 0 (and/or π).  These choices correspond to the x-z plane.

(ii)  Use the solution for the Potential along the symmetry (z) axis as a boundary condition for the solution of Laplace’s Equation in Spherical Coordinates.

(a) Set up the series solution for z > 0 ,  r > R  (i.e., consider only the (1/(r^(l +1)))(Subscript[P, l])(Cos[θ] terms in the expansion).

<Enter your answer in this text cell>

Treat as a Laplace's Equation Boundary Condition Problem in Spherical Coordinates.
 
		 V[r, θ] = (∑_( l=0)^∞)(  (Subscript[A, l]) (r^l)  +  (Subscript[B, l]) (1/(r^(l+1))))(Subscript[P, l])(Cos[θ])  
 
 The potential of a uniformly charged disk with surface charge density σ and of radius R for points on the normal z axis of the disk is:
 
		 V[z] = (σ/(2 (Subscript[ϵ, o]) )) ((Sqrt[((z^2) + (R^2))])  - |z|) 
			
(Griffiths’ Problem 2.5.  The problem is worked at https://www.youtube.com/watch?v=Ta1q5_d-7yM)  In spherical coordinates, a point with coordinate z on the axis has coordinates r = z and θ = 0 or π.  (ϕ is indeterminate, but the symmetry implies that the potential does not depend on ϕ.)  

For z >0:   { r, 0 }  i.e., θ = 0; changing to 2D spherical coordinates (we let z → r), we write:

	    V[r, 0] =  (σ/(2 (Subscript[ϵ, o]) )) ((Sqrt[((r^2) + (R^2))]) - r)       (r  ≥ 0)

This formula is valid only for points on the z axis, but we use it as a Boundary Condition.  Applying this boundary condition to the general solution of Laplace's equation in spherical coordinates (Legendre polynomials) yields an approximate solution for the potential off the z axis.

We consider first the case of r > R.   In this case, the  general solution of Laplace's Equation FOR r > R (We through away the A terms.) is

		V[r, θ] = (∑_( l=0)^∞) (Subscript[B, l]) (1/(r^(l+1)))(Subscript[P, l])(Cos[θ])   (The A's are all zero  — this is a second BC:   V→ 0 as r→ infinity.)

(b) Use the facts that at θ = 0, that all of the Legendre Polynomials = 1 when θ = 0, and the hint above to find the (Subscript[B, l])’s for the desired solution.  Write Vfar[r,θ] as the sum of the first four non-zero terms of the resulting series.

<Reality check>
If we set θ = 0 what are the(Subscript[P, l])’s?   To find out, execute:   LegendreP[n, Cos[0]].

(* Input code below *)
ClearAll["`*"]  (* Leave the ClearAll statement *)
LegendreP[n, Cos[0]] 
(* For all n *)

Since n is arbitrary, the answer “1” implies that every Legendre polynomial equals one when θ = 0.  This is a property of Legendre polynomials.  Below is an example of how this works out for n = 5:  (63 - 70 + 15)/8 = 1.

LegendreP[5, Cos[θ]]  (*   just an example  *)

For θ = 0, all of the Legendre Polynomials = 1, where, of course, Cos[θ] = 1.

<Enter your derivation of the (Subscript[B, l]) values below>
<Create new code and text cells where it helps>

Using our θ =  BC:  

V[r, 0] =(σ/(2 (Subscript[ϵ, o]) )) ((Sqrt[((r^2) + (R^2))]) - r) =      (∑_( l=0)^∞) (Subscript[B, l])  (1/(r^(l+1))) * (1)
                       lhs			                                 	rhs
To compare the lhs with the rhs and determine the B’s, we need to expand (Sqrt[((r^2) + (R^2))]) = r (Sqrt[(1 + ((R^2)/(r^2)))]) in powers of (1/r)  
First take the radical alone:  (Sqrt[(1 + ((R^2)/(r^2)))])   

You can use M’s Series[] command to expand the radical as a power series.  Enclosing the Series[] function in a Normal[] writes each term in the series in “normal” form.  Looks better.  This is an infinite series — I truncated it at (1/(r^8)) .

Normal[Series[(((1 + (α^2)))^(1/2)),{α,0,9}]]/.α-> (R/r)

Multiplying each term of the  series by r give the original lhs:  ((Sqrt[((r^2) + (R^2))]) - r)    (without the (σ/(2 (Subscript[ϵ, o]) )) factor).  This is the expanded form of our BC, i.e., V[r,0] (along z axis):

lhs =( Expand[r*%]-r)   (* examine this carefully to justify *)

Write rhs out to B[7] (or more) This is our approximate (symbolic) solution of Laplace’s Equation

rhs=   (∑_( l=0)^7) B[l]  (1/(r^(l+1)))

Now equate coefficients of (1/(r^n)).  Here are the resulting B’s (leaving out the factor of (σ/(2 (Subscript[ϵ, o])))).    This step was done by hand.

Clear[R]
B[0] = ((R^2)/2);
B[1] = 0;
B[2] =- ((R^4)/8);
B[3] = 0;
B[4] = ((R^6)/16);
B[5] = 0;
B[6] =-((5 (R^8))/128);
B[7] = 0;

Combining everything to obtain the final Vfar.  (Realize that by keeping only a finite number of terms, our result is an approximation.  Keeping only a finite number of terms is “truncating the series.”

Vfar[r_,θ_]  =(σ/(2 ϵo)) (∑_( l=0)^7) B[l]  (1/(r^(l+1)))LegendreP[l,Cos[θ]]

(c) Convert your Vfar[r,θ] to Cartesian coordinates for plotting.  Plot the result, V[x,z], using Plot3D.  You will have to choose “reasonable” values for R, σ, and ϵo.  I chose R = 1, σ = 1, and ϵo = 1.  (Computers and humans seem to do better with input values closer to 1 than (10^(-12)).)

(* Input code below *)
Vfarxz[x_, z_] =Vfar[r,θ]/.{r->(Sqrt[(x^2)+(z^2)]),Cos[θ]->(z/(Sqrt[(x^2)+(z^2)]))}

R = 1; σ =1; ϵo = 1;
r = (Sqrt[(x^2)+(z^2)]);
Plot3D[If[r > 1,Vfarxz[x,z],0],{x,-2,2},{z,0,5},PlotRange->{0,.3},AxesLabel->{x,z},PlotLabel->"V[x,z] for r > R",ColorFunction->Function[{x,y,z},Hue[z]]]  
(*I've set V = 0 for r < R for plotting purposes *)

(d)  Now repeat the above for r < R to find Vnear[r,θ].  Keep the first six nonzero terms.

<Enter your derivation of the (Subscript[A, l]) values below>
<Create new code and text cells where it helps>

For r < R: 

We use the Laplace’s Eq. solution:  V[r,θ] =  (∑_( l=0)^∞) (Subscript[A, l]) (r^l) (Subscript[P, l])[Cos[θ]      (Omit the(1/r) terms to avoid infinity at  r = 0.)

To find the A’s. WE again use our θ=0 BC AND the fact that LegendreP[Cos[0]] = 1:

V[r, 0] =(σ/(2 (Subscript[ϵ, o]) )) ((Sqrt[((r^2) + (R^2))]) - r) =   (∑_( l=0)^∞)  (Subscript[A, l]) (r^l)  (the B’s = 0)

We therefore want to first expand the (Sqrt[sqrt]) in powers of r:

(Sqrt[((r^2) + (R^2))]) = R (Sqrt[(((r^2)/(R^2)) + 1)]);  So letting (β^2) = ((r^2)/(R^2))

Remove["`*"] (* Remove deletes symbol names from M's memory. *)
Clear[R]
Normal[Series[(((1 + (β^2)))^(1/2)),{β,0,8}]]/.β-> (r/R)

Putting it all together (without the (σ/(2 (Subscript[ϵ, o]) )) factor), we find is the expanded form of our BC, i.e., V[r,0] (along z axis):

lhs =Expand[ R % ]- r

Write the rhs (our solution of Laplace’s Equation) out to A[8] (or more).

rhs=   (∑_( l=0)^8) A[l] (r^l)

Now equate the A’s to the lhs expansion (BC). Here are the resulting A’s (leaving out the factor of (σ/(2 (Subscript[ϵ, o]))))

Clear[R,σ,ϵo]

A[0] = R;
A[1] =-1;
A[2] = (1/(2 R));
A[3] = 0;
A[4] =- (1/(8 (R^3))) ;
A[5] = 0;
A[6] =(1/(16 (R^5)));
A[7] = 0;
A[8] = -(5/(128 (R^7)));

Putting it all together gives Vnear   (for r < R).

Vnear[r_,θ_]  =(σ/(2 ϵo)) (∑_( l=0)^8) A[l] (r^l) LegendreP[l,Cos[θ]]

(e) Convert your Vnear[r,θ] to Cartesian coordinates for plotting.   Plot the result, V[x,z].  You will have to give R, σ, and ϵo reasonable values again.  The old ones have been erased.

(* Input code below *)
Vnearxz[x_, z_] =Vnear[r,θ]/.{r->(Sqrt[(x^2)+(z^2)]),Cos[θ]->(z/(Sqrt[(x^2)+(z^2)]))}

R = 1; σ = 1; ϵo = 1;
Plot3D[If[(Sqrt[(x^2)+(z^2)]) ≤    1,Vnearxz[x,z],0],{x,-2,2},{z,0,4},PlotRange->{0,.6},AxesLabel->{x,z},PlotLabel->"V[x,z] for r < R",ColorFunction->Function[{x,y,z},Hue[z]]]  (* I've set V = 0 for r > R *)

(f)  Now combine the two plots (both r < R and r > R).   Reminder: here we are considering z >0 only.

For my sanity, I copy and paste the two equations for Vfarxz[x,z] and Vnearxz[x,z] from above, not forgetting to redefine R, σ, and ϵo in terms of actual numbers (all 1).

(* Input code below *)
Vfarxz[x_,z_] =((σ (-((5 (R^8) (((105 (z^2))/((x^2)+(z^2)))+((231 (z^6))/((((x^2)+(z^2)))^3))-((315 (z^4))/((((x^2)+(z^2)))^2))-5))/(2048 ((((x^2)+(z^2)))^(7/2))))+(((R^6) (-((30 (z^2))/((x^2)+(z^2)))+((35 (z^4))/((((x^2)+(z^2)))^2))+3))/(128 ((((x^2)+(z^2)))^(5/2))))-(((R^4) (((3 (z^2))/((x^2)+(z^2)))-1))/(16 ((((x^2)+(z^2)))^(3/2))))+((R^2)/(2 (Sqrt[(x^2)+(z^2)])))))/(2 ϵo));
Vnearxz[x_,z_] = (1/(2 ϵo))(R-z-((5 ((((x^2)+(z^2)))^4) (35+((6435 (z^8))/((((x^2)+(z^2)))^4))-((12012 (z^6))/((((x^2)+(z^2)))^3))+((6930 (z^4))/((((x^2)+(z^2)))^2))-((1260 (z^2))/((x^2)+(z^2)))))/(16384 (R^7)))-((((((x^2)+(z^2)))^2) (3+((35 (z^4))/((((x^2)+(z^2)))^2))-((30 (z^2))/((x^2)+(z^2)))))/(64 (R^3)))+((((x^2)+(z^2)) (-1+((3 (z^2))/((x^2)+(z^2)))))/(4 R))+((((((x^2)+(z^2)))^3) (-5+((231 (z^6))/((((x^2)+(z^2)))^3))-((315 (z^4))/((((x^2)+(z^2)))^2))+((105 (z^2))/((x^2)+(z^2)))))/(256 (R^5)))) σ  ;
R = 1; σ = 1; ϵo = 1;
vplot = Plot3D[If[(Sqrt[(x^2)+(z^2)]) ≤    R,Vnearxz[x,z],If[(Sqrt[(x^2)+(z^2)]) ≥ R,Vfarxz[x,z]]],{x,-2,2},{z,0,5},PlotRange->{0,.6},AxesLabel->{x,z},PlotLabel->"V[x,z] for both r > R and r < R combined",BoxRatios->{1,1,.5},ColorFunction->Function[{x,y,z},Hue[z]]];
disky = Graphics3D[{Thick, Brown,Line[{{-R/2,0,0 },{+R/2,0,0}}]}];
Show[vplot,disky]

(g)  Interpret this Plot.

<Enter your interpretation in this text cell>

The coordinates:  x and z are the variables in the base plane of this plot; “Vertical Direction” is the magnitude of the SCALAR potential.  

The "Vertical Direction" is not “y”; it is V[x,z].  Where is the disk???  We might ask first:  Where is y??  The problem is being treated as a 2D problem:

Remember:  “HINT (i) Treat as a 2D problem.  Any plane through the z axis will do—take the x,z plane (i.e., we consider y = 0).  This is equivalent to saying there is no ϕ dependence; this plane corresponds to a constant ϕ (in this case ϕ = 0).”

The disk is in the x-y plane.  If we set y = 0 (for constant ϕ =  0) we could represent the cut this plane would make as a line (shown in the above plot).  Look at the figure at the very top of the notebook and imagine the plane of the “paper” (which is y = 0).

Realize that a reflection of this plot into the z < 0 region would provide the full 2D solution for V[r,θ].  Without dwelling on how I did it (actually, I just let z → Abs[z] where it mattered), here is what it looks like.  Remember z is the symmetry axis normal to the disk.

Finally note the slight mismatch at the boundary where r = R.  This is because we are using only a few terms in the two expansions for Vnear and Vfar.  With more terms, this mismatch would be reduced.

Vnearxz[x_,z_] = (1/2) (1-Sqrt[(z^2)]-((5 ((((x^2)+(z^2)))^4) (35+((6435 (z^8))/((((x^2)+(z^2)))^4))-((12012 (z^6))/((((x^2)+(z^2)))^3))+((6930 (z^4))/((((x^2)+(z^2)))^2))-((1260 (z^2))/((x^2)+(z^2)))))/16384)-(1/64) ((((x^2)+(z^2)))^2) (3+((35 (z^4))/((((x^2)+(z^2)))^2))-((30 (z^2))/((x^2)+(z^2))))+(1/4) ((x^2)+(z^2)) (-1+((3 (z^2))/((x^2)+(z^2))))+(1/256) ((((x^2)+(z^2)))^3) (-5+((231 (z^6))/((((x^2)+(z^2)))^3))-((315 (z^4))/((((x^2)+(z^2)))^2))+((105 (z^2))/((x^2)+(z^2)))));
vplot = Plot3D[If[(Sqrt[(x^2)+(z^2)]) ≤    R,Vnearxz[x,z],If[(Sqrt[(x^2)+(z^2)]) ≥ R,Vfarxz[x,z]]],{x,-2,2},{z,-2,2},PlotRange->{0,.6},AxesLabel->{x,z},PlotLabel->"V[x,z] for both r > R and r < R combined",BoxRatios->{1,1,.5},ColorFunction->Function[{x,y,z},Hue[z]]]

(h) Create a ContourPlot plot of the potential for points in the x-z plane.  Calculate and plot the electric field for the same region using VectorPlot.  Finally, produce a Streamline plot of the electric field (again, for the same region in the x-z plane).

(* Input code below *)
conplot =ContourPlot[If[(Sqrt[(x^2)+(z^2)]) ≤    R,Vnearxz[x,z],If[(Sqrt[(x^2)+(z^2)]) ≥ R,Vfarxz[x,z]]],{x,-2,2},{z,-2,2},ColorFunction->"Rainbow",Contours->20,PlotRange-> {0,1}]
EEnear[x_,z_] = -Grad[Vnearxz[x,z],{x,z}];
EEfar[x_,z_] = - Grad[Vfarxz[x,z],{x,z}];
vecplot = VectorPlot[If[(Sqrt[(x^2)+(z^2)]) ≤    R,EEnear[x,z],If[(Sqrt[(x^2)+(z^2)]) ≥ R,EEfar[x,z]]], {x,-2,2},{z,-2,2},VectorPoints->16]
streamplot = StreamPlot[If[(Sqrt[(x^2)+(z^2)]) ≤    R,EEnear[x,z],If[(Sqrt[(x^2)+(z^2)]) ≥ R,EEfar[x,z]]], {x,-2,2},{z,-2,2}]

(i) Show the ContourPlot and StreamPlot together.

(* Input code below *)
Show[conplot,streamplot]

(i) Interpret your ContourPlot + StreamPlot.  Where is the disk?

<Enter your interpretation in this text cell>

First, I kept the solution that included both ±z.  

Where is the disk?  It is in a plane PERPENDICULAR to the plane of the “paper”.  The y axis is INTO the paper.  Again, in the plane of the paper (y = 0) the cut made by the y-z plane with the disk would be  LINE along the x axis from x = –R to x = R.

One sees the contours of V (use your mouse to see the values of V on the contours), falling off in magnitude as you increase r (r = (Sqrt[(x^2) + (z^2)])).  The streamlines leave the disk (We assumed σ > 0 so the charge on the disk is +.) and heads off to points at infinity.  As we expect E and the streamlines are perpendicular to the contours (the contours marking equipotentials).  

Beat the Dead Horse:  This represents V and E in the y = 0 plane.  If you rotated the above plots about the z axis you could generate the full 3D picture.  So for r >>R, the contour SURFACES would look like spheres and the streamlines would be radial consistent with the disk looking like a point charge if you get far enough away.

Another curious observation:  the potential does not diverge as you approach the center of the disk (the origin).  The can be confirmed using M’s Limit function.  With the above choices for R, σ, ϵo, the limiting value of the potential is 0.5, consistent with the plot.  In contrast, the potential diverges as you approach a point charge, where the local charge density is infinite.  The actual amount of charge at the center of a uniformly charged disk is “infinitesimal.”  Note that the potential does diverge as you approach an infinite line charge (but not as you approach a finite line charge), even thought the charge at any one point on the line is also infinitesimal.  Divergences are painful.  It pays to keep track.