Laplace's Equation with Boundary Conditions in One Dimension

To date we have used Gauss's Law and the Method of Images to find the potential and electric field for rather symmetric geometries.  For more complex geometries, V(x,y,z) can often be found by solving Laplace's equation:  (∇^2) V(x,y,z) = 0.  Although the number of solutions to Laplace's equation is infinite, knowing either the potential or the field along the boundaries of the region of interest ("The Boundary Conditions") normally guarantees a unique solution (the uniqueness theorem).    

Solving differential equations can be relatively difficult.  In this module, we will restrict ourselves to simple 1-D systems for which solutions are readily derived.  After learning what we can from these, we will move on to more complicated systems where we "up the ante".

In this context, a one-dimensional problem does not mean that our fields exist in a one-dimensional vector space.  Our electric fields are vectors in 3-space, as they always are.  Rather, the direction and magnitude of the electric fields only depend on one variable.  In Cartesian space, this could be x.  In this case, the direction and magnitude of the electric field would not depend on y or z.  Similarly, the value of the potential would depend on x and not y or z.

Laplace’s Equation in One Dimension—Infinite Parallel Plates

In the infinite parallel plate geometry, the fields and potentials depend on only one Cartesian variable, say x.  The boundaries of the region of interest are planes parallel to the y-z plane, which we will assume intersect the x-axis at points x = x1 and x = x2.   We will assume that V(x1, y, z) = V1 (a specific value of potential) for all values of y and z.  Similarly V(x2,y,z) = V2.   Between x1 and x2, the potential will be some function V(x).  This is the geometry of a parallel plate capacitor with plates that extend to infinity in the y and z directions.  You will remember from your work with Coulomb’s Law and Gauss’s Law that V(x) in this system is proportional to x and the E field is constant in magnitude and direction (± (OverHat[x])  — in the direction of decreasing potential).

In 1-D Cartesian coordinates, Laplace’s equation takes the form

(((d^2)V[x])/(d (x^2)))  =  0

(a)  Although you probably know the solution, use M’s DSolve[V’’[x] == 0, V[x], x] command to solve this equation.  By NOT inserting the BCs, we will see that we can generate a general solution.   

The double prime after the V in the first argument indicates the second derivative.  (The double prime is not the quotation mark, “, it is two hits on the apostrophe key, ’ ’.)  We seek V[x], the second argument of the command, and the derivatives are taken with respect to x, the third argument.  To refer to the solution later, set the variable sol to the result of the DSolve[] command.

(* Input code below *)
ClearAll["Global`*"]  (* Leave the ClearAll command *)

sol=DSolve[V''[x]==0, V[x],x]      (* do not use " - quotation mark; use ''  two hits on the ' key *)

M returns a list of solutions (in the case ONLY ONE!). The entire list is enclosed in brackets, and each solution is enclosed in its own, additional brackets.  The solution {V[x] -> C[1] + x C[2]} is part of this larger list, so we get nested brackets.

(b)  Execute sol[[1]];  this should get you this single bracket element of sol:

(* Input code below *)

sol[[1]]

Technically, M treats the solution sol[[1]] as a "rule" that can be applied in various ways.  We want to convert this rule to a function.  We can use way the "ReplaceAll" command, usually written "/." (without the quotes) to replace V[x] with the rule given by sol[[1]].

(c)  Type in and execute   V[x]/.sol[[1]]

(* Input code below *)

V[x]/.sol[[1]]

This is still not a function that we can plot.

(d)  To get a function, we need to set a real function (in x_) equal to V[x].  This adds an extra step.  (use % for the above result to achieve this definition.)

(* Input code below *)

V[x_] = %   (* We now define a real function, in x. NOTE we are REDEFINING V[x], as all hackers love to do*)

(e)  Just as a check enter and execute the “new” V[x]

(* Input code below *)

V[x]

As expected, this  solution is a linear function of x.  We started with V’’[x], the second derivative.  M had to integrate twice to find the solution.  C[1] and C[2] are the two constants of integration.

(f) To meet our boundary conditions, we require V[x1] = V1 and V[x2] = V2.  Use these two requirements and M's Solve[] command to determine the values C[1] and C[2].   If the uniqueness theorem applies to this problem (It Does!), only one choice of C[1] and C[2] will work.  Like DSolve[], Solve[] returns a nested list of solutions.  We must specify which solution we want, even if there is only one.   SO:  Use Solve and extract the two constants C[1] and C[2] that satisfy our BCs.  Then substitute these values of C[1] and C[2] into V[x].

(* Input code below *)

sol3=Solve[{V[x1] ==V1,V[x2]==V2}, {C[1],C[2]}]  
V[x]/.sol3[[1]]    (* The ReplaceAll command substitutes the values of C[1] and C[2] into VV[x] *)
Print["V[x] = ",V[x_]=%]      (* The function V[x_] is now defined in terms of these solutions. *)

Handling M’s constants of integration is tricky.  Here is another way to do it.  Introducing C1 and C2 as temporary constants gets around M’s touchiness regarding the symbol “C”

Note that the ReplaceAll trick works only once per M session unless you change the constant names.  We can’t change the names of C[1] and C[2].  If you get an error message, quit M and re-execute only one version of of the ReplaceAll command.  Or enter the values manually.  

The following code will work in place of the method described above.  But you must use one method or the other—not both.  
sol3=Solve[{V[x1] ==V1,V[x2]==V2}, {C[1],C[2]}]  
C1= C[1]/.sol3[[1]]
C2 = C[2]/.sol3[[1]]
V[x_]  = C1 + C2 x

(g)  To verify that this choice of C[1] and C[2] satisfies the boundary conditions, evaluate V[x1] and V[x2] explicitly.  [You are trying to show that  V[x1] = V1 and V[x2] = V2.  

The Simplify[] command may be useful here.  For instance:  Simplify[V[x1]].   (This also works:  V[x1]//Simplify )      //  =  PostFix)

(* Input code below *)

Simplify[V[x1]]
Simplify[V[x2]]

We see that the BC's are satisfied.  If you get something like C[1]+x1 C[2], you need to make sure V[x_] is defined properly in terms of x1, x2, V1 and V2.

(h)  Now use DSolve with BCs inserted (notice how it saves a lot of trouble) to find the potential between the plates where again 
(code will look something like this:  sol2=DSolve[{V’’[x]==0,V[x1] ==V1, V[x2] ==V2}, V[x],x]
Just a reminder:  the plates are located at x = x1 and x = x2.   V(x1, y, z) = V1 and  V(x2,y,z) = V2.

(* Input code below *)
ClearAll["`*"]  (* Leave the ClearAll statement.  We are "starting over". *)

sol2=DSolve[{V''[x]==0,V[x1] ==V1, V[x2] ==V2}, V[x],x]
V[x_] = V[x]/.sol2[[1]]  (* this extracts the solution and assigns it to V[x]  *)

Same as V[x] above.  Here we have let M do all the work.

(i) Griffith's ( in Section 3.1.2) claims that all solutions V[x] of the 1-D Laplace's equation are featureless, and the potential at x is the average of the potentials at x+a and x-a;  it has no minimum or maximum except at the boundaries. [You can think of a being small but it’s not necessary].

We therefore expect:    V[xo] = ((V[xo-a] + V[xo+a])/2)

 Use M to see if Wolfram agrees.  (Use the test:  == ; might want to use Simplify.)

(* Input code below *)

lhs = V[xo]//Simplify
rhs =  (( V[xo-a] + V[xo+a])/2)//Simplify    (*  Should = V[xo]  *)

lhs == rhs   (* This does the TEST of equality  *)

So the statement that the potential at x is the average of the potentials at x+a and x–a is true.

This general result is important, because it shows that the solutions of Laplace’s equation have “no local minima or maxima”.  All the minima and maxima potential values lie on the boundaries.  This is also the basis of the “relaxation method” for numerically solving Laplace’s equation in Cartesian coordinates.

(j) For laughs, show that this is true for some consistent (numerical) choices of x1, x2, xo (located between x1 and x2), V1, and V2, and some arbitrary a.   

I chose:

x1 = 1; x2 = 2; 
V1 = 10;  V2 = 100;
xo = 1.5;  a = 0.02;

(* Input code below *)

x1 = 1; x2 = 2; 
V1 = 10;  V2 = 100;  xo = 1.5;  a = 0.02;
V[xo]
(V[xo - a] + V[xo+a])/2
V[xo]== (V[xo - a] + V[xo+a])/2  (*  True/False test for Equality *)

(k) Now Plot the Solution for x1 < x < x2.

(* Input code below *)

x1 = 1; x2 = 2; 
V1 = 10;  V2 = 100; 
plotcartesian = Plot[V[x], {x,x1,x2},PlotLabel->"Potential vs Distance between Parallel Plates",AxesLabel->{"Displacement (x)" ,"Voltage (V[x])"},PlotRange-> {{0,x2 + 0.5},{0,V2 + 5}}];
edgeone= Graphics[{Red,Dashed,Line[{{x1,0},{x1,V2}}]}];
edgetwo = Graphics[{Blue,Dashed,Line[{{x2,0},{x2,V2}}]}];
texty = Graphics[{Text["V1",{x1-0.1,V1}],Text["x1",{x1,V1-4}],Text["¬",{x1-0.01,V1-0.9}],Text["V2",{x2-0.1,V2}],Text["x2",{x2+ 0.02,V2-4}],Text["¬",{x2-0.01,V2-1}]}];

Show[plotcartesian, edgeone, edgetwo,texty]   (* dashed lines show positions of the plates *)

(l) Although somewhat mundane, interpret your results.  (Is the result consistent with your previous work in this geometry? — think parallel plate capacitor.)

<Enter interpretation in this text cell>

Not exciting, but linear, as expected.  The potential values at the ends (10 V and 100 V) are consistent with the values specified in the boundary conditions.  There are no min or max between the endpoints at the boundaries.

(m)  Find the E Field.  In 3D, EE = –Grad[V];  This yields a 3D Vector.

The electric field amplitude should be MINUS the slope of the potential plot times (–1), which is about –(90 V)/(1 m) or –90 V/m.

(* Input code below *)
Clear [xo, x1, x2, V1, V2, a] (* Leave Clear statement *)

V[x]  (*  just a check  *)
-Grad[V[x],{x,y,z}]

(OverVector[EE])[x] is CONSTANT VECTOR and points from the more positive plate towards the less positive plate (in the –(OverHat[x]) direction).  The E Field is UNIFORM between the plates.

 Everything looks consistent with our previous Coulomb’s Law/Gauss’s Law approach.

A general expression for just the x-component of the E field (the x-component being the only non-zero component) is given by:

(* Execute code below *)

V[x]
EE[x_] = -D[V[x],x]  (* note that this is the - of the slope of V[x] *)

So lets rewrite:    (OverVector[EE])[x] =  ((V2-V1)/(x2-x1))(OverHat[x]).  For V2 > V1 and x2>x1,  (OverVector[EE])[x] points in the + x direction as it should.

Not much more to say here.  BUT:  by going through this simple 1D potential problem using Laplace's Equation we get an understanding of the importance of BCs, the nature of the resulting potential, and the simple connection of the E-field to this V.

Laplace’s Equation in One Dimension—Concentric Cylinders

In systems with cylindrical symmetry, the potential is a function of the radial coordinate, s, and independent of ϕ and z.  We will seek the potential between two circular cylinders which serve as boundaries, where the boundaries are equipotentials.  The boundary conditions then take the form V(s1) = V1 and V(s2) = V2.

Ignoring the terms in ϕ and z, Laplace’s equation in cylindrical coordinates is given by:

(1/s) (∂/(∂s)) (s((∂V)/(∂s)))   =    (((∂^2)V)/(∂(s^2))) + (1/s) ((∂V)/(∂s))    =    0

By convention, we multiply each term by s before attempting a solution.  Factors of s are often easier to work with than factors of (1/s), that is,    s*V’’[s] + V’[s] == 0

(n) Use M to find V[s] using :  DSolve[{s*V''[s] + V'[s]  == 0,  V[s1] ==V1,V[s2]==V2 }, V[s], s].   WE'RE putting in the BCs into DSolve from the start.  This gives us the particular solution.  Is your solution consistent with what you found from Gauss’s Law earlier in the semester?  Comment.

(* Input code below *)
ClearAll["Global`*"] (* Leave the ClearAll statement *)
sol=DSolve[{s*V''[s] + V'[s] ==0,V[s1] ==V1,V[s2]==V2}, V[s],s]
V[s_] = V[s]/.sol[[1]]  (* This extracts the desired solution from the "rule" contained in the curly brackets  *)

<Enter comment in this text cell>

Obviously, the solution V[s]  is not linear, but is proportional to Log[s], as we found with Gauss’s Law.

(o) Verify that these expressions satisfy the boundary conditions (i.e., does V[s1] = V1, etc.??)  Comment.

(* Input code below *)

Simplify[V[s1]]
Simplify[V[s2]]

<Enter comment in this text cell>

We see that V[s1] = V1 and V[s2] = V2.  The BC's are satisfied.

(p) Plot V[s] for the co-axial cylinders.  For plotting, I assumed:
		s1 = 1;    s2 = 2; 
		V1 = 10; V2 = 100;

(* Input code below *)

s1 = 1; s2 = 2; V1 = 10; V2 = 100;
VV[s_] = If[s < s1,Null,If [s1≤ s < s2, V[s], If[s > s2,Null]]];  
(* Our solution V (or VV) is valid only in the region s1 ≤ s ≤ s2 *)
plotty= Plot[VV[s], {s, 0, s2},PlotRange->{{0,1.2 s2},{0,1.1V2}},
PlotLabel-> "The Potential vs. Distance between Co-axial Conducting Cylinders",AxesLabel->{"s", "V[s]"}];
texty = Graphics[{Text["V1",{s1-0.1,V1}],Text["s1",{s1,V1-4}],Text["¬",{s1-0.01,V1-0.9}],Text["V2",{s2-0.1,V2}],Text["s2",{s2+ 0.02,V2-4}],Text["¬",{s2-0.01,V2-.5}]}];

Show[plotty,texty]

Note that for the given constants s1, s2, V1, V2, the slope of the potential function V[s] is monotonically decreasing; as expected; no min or max in V[s] occurs except at the endpoints (the boundaries).

(q)  Clear the constants (s1, s2, V1, V2) and find the ONE DIMENSIONAL E field from V[s].
		Use M’s Grad function in Cylindrical Coordinates.

Here is some help:

M writes it as:  Grad[f[r,θ,z],{r,θ,z},”Cylindrical”], where f is a general, undefined, scalar function in cylindrical coordinates.  Evaluate the following:

(* Evaluate code below *)
Clear[s1,s2, V1, V2]
Grad[f[s,ϕ,z],{s,ϕ,z},"Cylindrical"]

M assumes that EACH component of the gradient can have a dependence on all three coordinates.  The result is a 3D Vector which is easily translated into:
						  {(∂/(∂s))f[s,ϕ,z],(((∂/(∂ϕ))f[s,ϕ,z])/s),(∂/(∂z))f[s,ϕ,z]}  
						  Note M’s notation:  [Derivative[1,0,0][f][s,ϕ,z] = (∂/(∂s))f[s,ϕ,z],  etc.]
						  
For our task, I suggest that you simply enter:  Ecylindrical[s_] = –Grad[V[s],{s,ϕ,z},"Cylindrical"] and execute.  Don't overlook the – sign.  The lack of ϕ and z dependence in V does its magic (i.e., the ϕ and z components come out zero.)

(* Input code below *)

Ecylindrical[s_]= -Grad[V[s],{s,ϕ,z},"Cylindrical"]

Repeat:  Since V[s] is a function of s only, only the s component of the Gradient is non zero ((∂/(∂ϕ))V[s] and (∂/(∂z))V[s] both = 0).

Given the symmetry of the problem we are not surprised to see that Ecylindrical is: 

			A function of s only 
AND 
			The only non-zero component is the (OverHat[s]) component.

(r)  Plot the 2D EE field using Vector Field.  (Since the field does not depend on z, consider z = 0).  REMEMBER:  M requires that we transform EEcylindrical into EEcartesian for plotting.

(* Execute the code below to express the field in terms of the Cartesian coordinates x and y *)

EECartesian[s_]= TransformedField["Cylindrical"->"Cartesian",Ecylindrical[s],{s,ϕ,ζ}->{x,y,z}]    (* does not like a common z!!! so we write ζ instead of z in {s,ϕ,ζ} *)
EECartesian2D[x_,y_] = {EECartesian[s][[1]], EECartesian[s][[2]]}   (* selects out the x,y components *)

Now put in some values for the constants and plot

(* Input code below *)

s1 = 1; s2 = 2; 
V1 = 10; V2 = 100;

vecy = VectorPlot[If[(Sqrt[(x^2)+(y^2)]) ≤ s1 || (Sqrt[(x^2)+(y^2)]) ≥ s2,{0,0},EECartesian2D[x,y]],{x,-3,3},{y,-3,3},VectorScale->0.05,VectorPoints->25];
innerelectrode =  Graphics[{Thick,Red,Circle[{0,0},s1]}];
outerelectrode =  Graphics[{Thick,Purple,Circle[{0,0},s2]}];
Show[vecy,innerelectrode,outerelectrode]

Remember this is a 2D PLOT.  It represents a cut through the cylinders normal to the z axis at (any) particular z, say z = 0.  

The E vectors are all radial and only are defined between the electrodes.  By default, M positions the middle of each arrow at the location of the vector.  Therefore the heads and tails will appear to cross the boundaries at points.  The middle of each arrow lies between the boundaries.

For our assumed constants, E points from the outer more positive electrode (cylinder), inward. E increases in intensity (the vectors get longer) as one moves from s2 into s1.  This is consistent with the  SLOPE of V[s] (largest at s = s1).  (Again, don’t forget the – sign in (OverVector[E]) = – Grad[V]  (so the E vectors point towards smaller, decreasing s).  Here is a copy/paste of the V[s] plot above.

Just to nail our understanding, we will generate the contours of V[s] (between the electrodes), the streamline plot of the E field, and their superposition.  Note: to plot the contour map, we have to convert V[s] into it’s Cartesian form (in this case, it is simply substituting s = (Sqrt[(x^2) + (y^2)])  into V[s] (from above, V[s] = ((V1 Log[s]-V2 Log[s]+V2 Log[s1]-V1 Log[s2])/(Log[s1]-Log[s2])))

I use M’s TransformedField function to do the work;  this time it is transforming a scalar function, V[s].

(* Execute code below *)
Clear[s1,s2,V1,V2] (* Leave the Clear statement *)

VV[x_,y_] =TransformedField["Cylindrical" -> "Cartesian",V[s], {s,θ,ζ}-> {x,y,z}]
s1 = 1; s2 = 2; 
V1 = 10; V2 = 100;
contoury = ContourPlot[If[(Sqrt[(x^2)+(y^2)]) ≤ s1 || (Sqrt[(x^2)+(y^2)]) ≥ s2,Null,VV[x,y]],{x,-3,3},{y,-3,3}] ;   
Show[contoury, innerelectrode,outerelectrode]

As expected the contours of V[s] are circles (s = constants); in 3D they would be concentric cylinders.  If you want to see the magnitudes of V corresponding to each contour, click your mouse inside the plot and point to the contour lines.

Now we plot the streamlines:

(* Execute code below *)
streamy = StreamPlot[If[(Sqrt[(x^2)+(y^2)]) ≤ s1 || (Sqrt[(x^2)+(y^2)]) ≥ s2,{0,0},EECartesian2D[x,y]],{x,-3,3},{y,-3,3},StreamPoints-> 100];

Show[streamy, innerelectrode,outerelectrode]

And finely superimpose the streamlines onto the contourplot.

(* Execute code below *)
Show[ contoury,innerelectrode,outerelectrode,streamy]

One more time — because of the symmetry, this is a ONE D Laplace’s Equation problem.  (Remember, this is a SLICE through the CYLINDERS, e.g., at z = 0).

Laplace’s Equation in One Dimension—Spherical Symmetry

In systems with spherical symmetry, the potential is a function of the radial coordinate, r, and independent of θ and ϕ.  We will seek the potential between two concentric spheres which serve as boundaries, where the boundaries are equipotentials.  The boundary conditions then take the form V(r1) = V1 and V(r2) = V2.

Laplace’s equation for spherical coordinates, ignoring the terms in θ and ϕ, is:  (1/(r^2)) (∂/(∂r))((r^2) ((∂V)/(∂r)))= (2/r) ((∂V)/(∂r)) + (((∂^2)V)/(∂(r^2))) = 0

It is conventional to multiply each term (including the terms in θ and ϕ, which are not shown) by (r^2).  This removes the r’s in the denominator of each term (including those not shown).

(s) Use the M DSolve function with the BCs:  V[r1] = V1 and V[r2] = V2 to find V[r] between the spherical electrodes.  Is your solution consistent with what your found earlier in the semester, using Gauss’s Law?  Comment.

(* Input code below *)
ClearAll["Global`*"] (* Leave the ClearAll statement *)
sol = DSolve[{(2/r)*V'[r] + V''[r] == 0, V[r1] == V1, V[r2] == V2}, V[r], r]
V[r_] = Simplify[V[r] /. sol[[1]]]

<Enter comment in this text cell>

Looking carefully, we see that the solution (V[r]) itself is not linear, but is of the form  K1 + (K2/r), consistent with what we found using Gauss’s Law.

(t)  Check to make sure V[r] does indeed satisfy the boundary conditions, V[r1] = V1 and V[r2] = V2.

(* Input code below *)

Simplify[V[r1]]
Simplify[V[r2]]

(u) Plot your results for the potential as a function of r for the concentric spheres.  For plotting, I assumed:
		r1 = 1;   r2 = 2; 
		V1 = 10; V2 = 100;

(* Input code below *)

r1 = 1; r2 = 2; V1 = 10; V2 = 100;
VV[r_] = If[r < r1,Null,If [r1≤ r < r2, V[r], If[r > r2,Null]]];  
(* Our solution V (or VV) is valid only in the region r1 ≤ r ≤ r2 *)
plotty= Plot[VV[r], {r, 0, r2},PlotRange->{{0,1.2 r2},{0,1.1V2}},
PlotLabel-> "The Potential vs. Distance between Concentric Conducting Spheres",AxesLabel->{"r", "V[r]"}];
texty = Graphics[{Text["V1",{r1-0.1,V1}],Text["r1",{r1,V1-4}],Text["¬",{r1-0.01,V1-0.9}],Text["V2",{r2-0.1,V2}],Text["r2",{r2+ 0.02,V2-4}],Text["¬",{r2-0.01,V2-.5}]}];

Show[plotty,texty]

Note that for the given constants r1, r2, V1, V2, the slope of the potential function V[r] is monotonically decreasing; as expected, no min or max in V[r] occurs except at the endpoints (the boundaries).

(v)  Clear the constants (r1, r2, V1, V2) and find the ONE DIMENSIONAL E field from V[r].
		Use M’s Grad function in Spherical Coordinates.

Here is some help:

M writes it as:  Grad[f[r,θ,z],{r,θ,z},”Spherical”], where f is the scalar function being operated on.
Evaluate the following:

(* Execute code below *)
Clear[r1,r2, V1, V2]
Grad[f[r,θ,ϕ],{r,θ,ϕ},"Spherical"]

Here f[r, θ, ϕ] is a general, undefined, scalar function in spherical coordinates; M assumes that EACH component of the gradient can have a dependence on all three coordinates.  The result is a 3D Vector which is easily translated into:

		  {(∂/(∂s))f[r,θ,ϕ],(((∂/(∂θ))f[r,θ,ϕ])/r),((Csc[θ](∂/(∂ϕ))f[r,θ,ϕ])/r)},    which should look familiar.   
		 Note M’s notation:  Derivative[0,0,1][f][r,θ,ϕ] = (∂/(∂ϕ))f[r,θ,ϕ] ] 
		 
SO:   suggest that you simply enter:  Espherical[r_] = –Grad[V[r],{s,θ,ϕ},”Spherical”] and execute.  Don’t overlook the – sign.  The lack of ϕ and z dependence in V does its magic (i.e., the ϕ and z components come out zero.)

(* Input code below *)

Espherical[r_]= -Grad[V[r],{r,θ,ϕ},"Spherical"]

Broken Record:  This operation yields a 3D vector.  Since V[r] is a function of r only, only the r component of the Gradient is non zero ((∂/(∂θ))V[r] and (∂/(∂ϕ))V[r] both = 0).

Given the symmetry of the problem we are not surprised to see that Espherical is: 

			A function of r only 
AND 
			The only non-zero component is the (OverHat[r]) component.

(w)  Plot the 2D EE field using Vector Field (there is no z dependence; consider letting z = 0).  REMEMBER:  M requires that we transform EEspherical into EEcartesian for plotting.

(* Execute code below *)
EECartesian[x_,y_,z_]= TransformedField["Spherical"->"Cartesian",Espherical[r],{r,θ,ϕ}->{x,y,z}]    
EECartesian2D[x_,y_] = {EECartesian[x,y,0][[1]], EECartesian[x,y,0][[2]]}   (* selects out the x,y components *)

(x)  Now put in some values for the constants and make a VectorPlot of the E-Field.

(* Input code below *)

r1 = 1; r2 = 2; 
V1 = 10; V2 = 100;
vecy = VectorPlot[If[(Sqrt[(x^2)+(y^2)]) ≤ r1 || (Sqrt[(x^2)+(y^2)]) ≥ r2,{0,0},EECartesian2D[x,y]],{x,-3,3},{y,-3,3},VectorScale->0.05,VectorPoints->25];
innerelectrode =  Graphics[{Thick,Red,Circle[{0,0},r1]}];
outerelectrode =  Graphics[{Thick,Purple,Circle[{0,0},r2]}];
Show[vecy,innerelectrode,outerelectrode]

Remember this is a 2D VectorPlot.  It represents a cut through (or cross section of) the spheres.  It’s easiest to think of a plane passing through the polar (z) axis which of course passes through the origin.

The E vectors are all radial and only are defined between the electrodes.  For our assumed constants, E points from the outer more positive sphere, inward. E increases in intensity (the vectors get longer) as one moves from r2 into r1.  This is consistent with the  SLOPE of V[r] (largest at r= r1).  (Again, don’t forget the –sign in (OverVector[E]) = – Grad[V]  (so the E vectors point towards smaller, decreasing r.)

Just to nail our understanding, we will generate the contours of V[r] (between the electrodes), the streamline plot of the E field, and their superposition.  Note: to plot the contour map, we have to convert V[r] into it’s Cartesian form.

I use M’s TransformedField function to do the work;  again, it is transforming a scalar function, V[r].

(* Execute code below *)
Clear[r1,r2,V1,V2] (* Leave Clear statement *)

VV[x_,y_] =TransformedField["Spherical" -> "Cartesian",V[r], {r,θ,ϕ}-> {x,y,z}]/.z-> 0
r1 = 1; r2 = 2; 
V1 = 10; V2 = 100;
contoury = ContourPlot[If[(Sqrt[(x^2)+(y^2)]) ≤ r1 || (Sqrt[(x^2)+(y^2)]) ≥ r2,Null,VV[x,y]],{x,-3,3},{y,-3,3}] ;   
Show[contoury, innerelectrode,outerelectrode]

As expected the contours of V[r] are circles (r= constants); in 3D they would be concentric spheres.

Now we plot the streamlines:

(* Execute code below *)
streamy = StreamPlot[If[(Sqrt[(x^2)+(y^2)]) ≤ r1 || (Sqrt[(x^2)+(y^2)]) ≥ r2,{0,0},EECartesian2D[x,y]],{x,-3,3},{y,-3,3},StreamPoints-> 100];

Show[streamy, innerelectrode,outerelectrode]

And finely superimpose the streamlines onto the contour plot.

(* Execute code below *)
Show[ contoury,innerelectrode,outerelectrode,streamy]

One more time — because of the symmetry, this is a ONE D Laplace’s Equation problem (r dependent only).   Remember, the plots above are slices through the SPHERES.

Clearly, and not surprising, the basic behavior of the 2D plots for the cylinders and the spheres look somewhat similar — BUT, importantly, they are represented by different functions.