Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

8. Electrostatics

8.1 Introduction

Off[General::spell ];  
Off[General::spell1];

Mathematica Commands for All Sections

VEPlot

Needs["Graphics`PlotField`"]

VEPlot[potential_, xlim_, ylim_, opts___] := 
        Module[{plot1, plot2}, 
plot1 = 
            PlotGradientField[-potential, xlim, ylim, 
              ScaleFunction -> (1 &), 
              DisplayFunction -> Identity]; 
plot2 = 
            ContourPlot[potential, xlim, ylim, ContourShading -> False, 
              ContourSmoothing -> True, 
              DisplayFunction -> Identity, 
              PlotPoints -> 50]; 
 
          Show[{plot1, plot2}, opts, 
            DisplayFunction -> $DisplayFunction ]];

    
VEPlot[potential_, opts___] := 
      VEPlot[potential, {x, 1.1, -1.1}, {y, 1.1, -1.1}, 
        opts]

Protect[VEPlot];

Example:  Equipotential surface and electric field of two-point charges

dipole=(1/(Sqrt[(x^2) + (((-y - (1/2)))^2)])) - 
            (1/(Sqrt[(x^2) + ((((1/2) - y))^2)]));

VEPlot[dipole, 
        Epilog -> {{Hue[ .3],  Disk[{0, 1/2}, 0.1]}, {Hue[0.95],  
              Disk[{0, -1/2}, 0.1]}}];

8.2 Point Charges, Multipoles, and Image Charges

Overview of Point Charges, Multipoles, and Image Charges

Mathematica Commands for Section 8.1

Monopole

Monopole[q_:1,r0_:{0,0,0},r_:{x,y,z}]:=
q/Sqrt[ Sum[(r0[[i]]-r[[i]])^2,{i,1,Length[r0]}] ];

Example

Monopole[q,{x0,y0},{x,y}]

Monopole[q,{x0,y0}]

dipole= ( Monopole[+1,{0,0,-a/2}] + 
          Monopole[-1,{0,0,+a/2}] )

TrigToY

TrigToY[expression_,terms_:2]:= 
Module[{el,m,result}, 
  result=
    Sum[
       Integrate[ (-1)^(m) SphericalHarmonicY[el,-m,θ,ϕ]
                  * expression  Sin[θ]
               ,{θ,0,π},{ϕ,0,2 π}]
    * Y[el,m,θ,ϕ]
   ,{el,0,terms},{m,-el,el}];
  Return[result]
];

Example

eq1= Cos[θ]^2  //TrigToY

eq1 //. Y -> SphericalHarmonicY // Simplify

TrigToP

TrigToP[expression_,terms_:2]:= 
Module[{el,result}, 
  result=
    Sum[ Integrate[ LegendreP[el,Cos[θ]] expression  Sin[θ]
               ,{θ,0,π}]
        * (2 el+1)/2 P[el,Cos[θ]]
   ,{el,0,terms}];
  Return[result]
];

Example

eq1 = Sin[θ]^2 // TrigToP

eq1 //. P -> LegendreP // Simplify

Protect[Monopole,TrigToY,TrigToP];

Problem 1: Superposition of point charges

Remarks and outline

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"] 
Needs["Graphics`PlotField3D`"]

Solution

Clear["Global`*"];

Part a

quad=( Monopole[2 q,{0,0,   0} ]+
       Monopole[ -q,{0,0, d/2} ]+
       Monopole[ -q,{0,0,-d/2} ])

triangle=( Monopole[+2q, {   0,0,   0} ]+
           Monopole[ -q, {+d/2,0,-d/2} ]+
           Monopole[ -q, {-d/2,0,-d/2} ]) ;

square=( Monopole[+q,{+d/2,0, d/2} ]+
         Monopole[-q,{+d/2,0,-d/2} ]+
         Monopole[+q,{-d/2,0,-d/2} ]+
         Monopole[-q,{-d/2,0, d/2} ]) ;

Plot3D[
        square  /. {y -> 0, q -> 1, 
            d -> 1}
, {x, -1, 1}, {z, -1, 
          1}
, 
        AxesLabel -> {"x", "z", ""}];

Part b

charges =  {{Hue[0.3],  
            Disk[{0, 1/2}, 
              0.1]}, 
                 {Hue[0.95], 
            Disk[{0, 0}, 
              0.2]}, 
                {Hue[
              0.3],  Disk[{0, -1/2}, 0.1]}};

p1 = 
        PlotGradientField[-quad //. {q -> 1, d -> 1, 
              y -> 0}
, {x, -1.1, 
            1.1}, {z, -1.1, 1.1}
, 
          ScaleFunction -> (1 &), 
          Epilog -> charges];

p2 = 
        ContourPlot[
          quad //. {q -> 1, d -> 1, 
              y -> 0}, 
{x, -1.1, 
            1.1}, {z, -1.1, 
            1.1}, 
ContourShading -> 
            False, 
PlotPoints -> 
            50, 
DisplayFunction -> 
            Identity
];

Show[p1, p2];

VEPlot[
        quad //. {q -> 1, d -> 1, 
            y -> 0}, 
{x, -1.1, 
          1.1}, {z, -1.1, 1.1}, 
Epilog -> 
          charges];

charges =  {{Hue[0.3],  
            Disk[{-1/2, -1/2}, 
              0.1]}, 
                 {Hue[0.95], 
            Disk[{0, 0}, 
              0.2]}, 
                {Hue[
              0.3],  Disk[{1/2, -1/2}, 0.1]}};

VEPlot[
        triangle /. {q -> 1, d -> 1, 
            y -> 0}, 
{x, -1.1, 
          1.1}, {z, -1.1, 1.1}, 
Epilog -> 
          charges];

charges =  {{Hue[0.95],  
            Disk[{-1/2, -1/2}, 
              0.1]}, 
                 {Hue[0.95], 
            Disk[{1/2, 1/2}, 
              0.1]}, 
                  {Hue[0.3],  
            Disk[{-1/2,  1/2}, 
              0.1]}, 
                {Hue[
              0.3],  Disk[{ 1/2, - 1/2}, 0.1]}};

VEPlot[
        square /. {q -> 1, d -> 1, y -> 0}, {x, -1.1, 
          1.1}, {z, -1.1, 1.1}, Epilog -> charges];

Part c

x2rRule = 
      Thread[{x, y, z} -> 
          CoordinatesToCartesian[{r, θ, ϕ}, 
            Spherical]]

seriesQuad =  Series[quad  //. x2rRule , {r, ∞, 2}]  // 
          Simplify // Normal

seriesQuad //TrigToP //Simplify

seriesSquare =  Series[
            square  //. x2rRule , {r, ∞, 2}]   // Simplify  // 
        Normal

seriesSquare //TrigToY //Simplify

seriesTriangle =  Series[
            triangle  //. x2rRule , {r, ∞, 2}]   // 
          FullSimplify // Normal

seriesTriangle //TrigToY//Simplify

Part d

p1 = 
        PlotGradientField3D[
          quad //. {q -> 1, d -> 1}, {x, -1, 1}, {y, -1, 
            1}, {z, -1, 
            1}, 
ScaleFunction -> (1 &), 
VectorHeads -> 
            True, 
PlotPoints -> 
            6, 
BoxRatios -> {1, 1, 
              2}, 
ViewPoint -> {1.3, -2.4, 
              1}];

Problem 2: Point charges and grounded plane

Remarks and outline

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"]

Solution

Clear["Global`*"];

Part a

dipole= ( Monopole[+q,{0,0,+d}] + 
          Monopole[-q,{0,0,-d}] )

dipole/.z->0

Part b

x2rRule = 
      Thread[{x, y, z} -> 
          CoordinatesToCartesian[{r, θ, phi}, Spherical]]

dipoleR =  Series[dipole //. x2rRule , {r, ∞, 2}]  // 
          Normal // Simplify

dipoleR//TrigToP

efieldR= -Grad[dipoleR,Spherical[r,θ,ϕ]]

Part c

charge = -1/(4 π) D[dipole, z] //. z -> 0

Integrate[charge,{y,-∞,+∞},{x,-∞,+∞}]//PowerExpand

ContourPlot[
        charge /. {d -> 1, q -> 1}, {y, -1, 1}, {x, -1, 
          1}];

Problem 3: Point charges and grounded sphere

Remarks and outline

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"]

Solution

Clear["Global`*"];

Part a

eq1 = {(q/(a + d)) + (q1/(a + d1)) == 0, 
          (q/(-a + d)) + (q1/(a - d1)) == 0};

imageRule= Solve[eq1,{d1,q1}] //ExpandAll //Flatten

potential = (Monopole[q, {0, 0, d}] + Monopole[q1, {0, 0, d1}]) //. 
        imageRule

Part b

x2rRule=Thread[{x,y,z}-> CoordinatesToCartesian[{r,θ,ϕ},Spherical]]

potentialR= potential //.x2rRule //Simplify //Expand

potentialR //. r -> a // Simplify // PowerExpand

Series[potentialR, {r, ∞, 2}]  // Normal   // 
      Collect[#, { Cos[θ], r,  a }] &

VEPlot[
        potential /. {a -> 0.2, d -> 0.4, q -> 1, 
            y -> 0}, {x, -1, 1}, {z, -1, 1}, 
        Epilog -> { {Hue[0.3], 
              Disk[{0, 0.1 }, 
                0.05]}, 
                 {Thickness[0.01], 
              Circle[{0, 0}, 
                0.2]}  , 
                  {Hue[0.95], Disk[{0, 0.4}, 0.1]}}];

Part c

chargedensity=(-(1/(4 π)) D[potentialR, r]/.{r -> a}//Simplify//PowerExpand//Simplify)

totalcharge = 
Integrate[
          chargedensity  2 π a^2  Sin[θ], {θ, 
            0, π}]  // Simplify

rule = {(1/((Sqrt[(a - d)])^2)) -> 
            (1/(-a + d)), (Sqrt[(((a - d))^2)]) -> -a + 
              d};

totalcharge  //. rule // PowerExpand  // Simplify

Part d

\(\(\(θMin = \(+\((FindMinimum[\(+chargedensity\) /. {a -> 1, 
                    q -> 1, d -> 2}, {θ, 1}] // 
              First)\)\);\)\(
\)
    \)\)
\(\(\(θMax = \(-\((FindMinimum[\(-chargedensity\) /. {a -> 1, 
                    q -> 1, d -> 2}, {θ, 1}] // 
              First)\)\);\)\(
\)
    \)\)
\(\(\(gray[θ_] = \((\((chargedensity - θMin)\)/\((θMax - θMin)\) /. {a -> 1, q -> 1, 
              d -> 2})\);\)\(
\)
    \)\)
{θMin, θMax}

radius=1;
ParametricPlot3D[{radius Sin[θ] Cos[ϕ]
                ,radius Sin[θ] Sin[ϕ]
                ,radius Cos[θ]
                ,GrayLevel[gray[θ]]}
          ,{θ,0,π}
          ,{ϕ,0,2 π}
          ,Lighting->False
          ,Axes->False
          ,Boxed->False];

Problem 4: Line charge and grounded plane

Remarks and outline

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"]

Solution

Clear["Global`*"];

Part a

potential = -q Log[(x^2) + (((y - y0))^2)] + 
          q Log[(x^2) + (((y + y0))^2)];

potential //.{y->0}

VEPlot[
        potential //. {q -> 1, y0 -> 1}, {x, -1.1, 
          1.1}, {y, -2, 2}, 
        Epilog -> { {Thickness[ .03], 
              Line[{{-1.1, 0}, {1.1, 
                    0}}] }, 
 {Hue[ .35], 
              Disk[{0,  -1 },  .1]} , 
  {Hue[ .95], Disk[{0,  1 },  .1]}}];

Laplacian[ potential,Cartesian[x,y,z]]//Together

Part b

charge=-(1/(4 π)) D[potential, y]/.{y -> 0}

ContourPlot[
        charge /. {q -> 1, y0 -> 1}, {x, -2, 2}, {y, -2, 
          2}];

Integrate[charge,{x,-∞,+∞}] //PowerExpand

Part c

eField= -Grad[potential,Cartesian[x,y,z]] //Simplify

PlotVectorField[
        eField[[{1, 2}]] //. {q -> 1, 
            y0 -> 1}, 
{x, -5, 5}, {y, 0, 
          3}, 
PlotPoints -> 
          10, 
ScaleFunction -> (1 &) , 
AspectRatio -> 8/10 , 
        Prolog -> { {Thickness[ .03], 
              Line[{{-5, 0}, {5, 0}}] }, {Hue[ .95], 
              Disk[{0,  1 },  .2]}}];

Problem 5: Multipole expansion of a charge distribution

Clear["Global`*"];

Remarks and outline

MultipoleSH

MultipoleSH[distribution_, rlimit_, n_]:=Module[{q, V},q:=Integrate[(Integrate[
                            distribution (r^2), {r, 0, rlimit}])*(((-1))^m) sin(θ) (Subscript[Y, el]^(-m))(θ, ϕ),{θ, 0, π},{ϕ, 0, 2 π}];
V=Sum[((4 π q ((((rlimit/r)))^el) Y(el, 
                                m, θ, ϕ))/((2 el + 
                                  1) r)),{el, 0, n},{m, -el, el}];
Return[V]];

density=((Q (r - 
                    a) (sin^2)(θ) (sin^2)(ϕ))/(a^4));

MultipoleSH[density,a,2] //Simplify //Apart

MultipoleP

MultipoleP[distribution_, rlimit_, n_]:=Module[{q, pot},q:=Integrate[2 π*(Integrate[distribution r^2, {r, 0, rlimit}])*((Sqrt[2 el + 1])/(2 (Sqrt[π]))) LegendreP[el, Cos[θ]]*Sin[θ],{θ, 0, π}];
pot=Sum[((2 (Sqrt[(1/(2 el + 1))]) (Sqrt[π]) q ((((rlimit/r)))^el) P(el, cos(θ)))/r),{el, 0, n}];
Return[pot]];

density=((Q (r - a (cos^2)(θ)))/(a^4));

eq1= MultipoleP[density,a,2]

eq2= MultipoleSH[density,a,2]

eq1 == eq2 //. {P -> LegendreP, 
          Y -> SphericalHarmonicY} // Simplify

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"]

Solution

Part a

chargeDensity1=((q0 (a - r) (sin(θ) - (π/4)))/(a^4));
chargeDensity2=(((E^(((a - r)/a))) q0 (a - 
                    r) cos(θ))/(a^4));
chargeDensity3=((q0 (((a - r))^2) (sin^2)(θ) sin(
                  2 ϕ))/(a^5));

r2xRule=Thread[{r,θ,ϕ}->
            CoordinatesFromCartesian[{x,y,z},Spherical]]

values = {q0 -> 1, a -> 1};

ContourPlot[chargeDensity1 /. r2xRule //. 
            values /. {y -> 0}, {x, -1, 1}, {z, -1, 1}, 
        PlotPoints -> 40];

ContourPlot[chargeDensity2 /. r2xRule //. 
            values /. {y -> 0}, {x, -1, 1}, {z, -1, 1}, 
        PlotPoints -> 40];

ContourPlot[chargeDensity3 /. r2xRule //. 
            values /. {z -> 0}, {x, -1, 1}, {y, -1, 1}, 
        PlotPoints -> 40];

Part b

multi1= MultipoleP[chargeDensity1,a,2]

multi2= MultipoleP[chargeDensity2,a,2] //Simplify

multi3= MultipoleSH[chargeDensity3,a,2] //Simplify

Part c

eField1 = -Grad[multi1, Spherical[r, θ, ϕ]]

eField1 //. P -> LegendreP // Simplify

eField2 = -Grad[multi2, Spherical[r, θ, ϕ]] //. 
          P -> LegendreP // Simplify

eField3 = -Grad[multi3, 
            Spherical[r, θ, ϕ]];
eField3 /. {Y -> SphericalHarmonicY} // ExpToTrig // 
      Simplify

Part d

pot1 =  multi1 /. r2xRule /. {P -> LegendreP, 
                Y -> SphericalHarmonicY} /. 
            values /. {y -> 0}  // Simplify

VEPlot[pot1,{x,1.1,-1.1},{z,1.1,-1.1}];

pot2 =  multi2 /. r2xRule /. {P -> LegendreP, 
              Y -> SphericalHarmonicY} /. 
          values /. {y -> 0}

VEPlot[pot2,{x,-1.1,1.1},{z,-1.1,1.1}];

pot3 =  multi3 /. r2xRule //. {P -> LegendreP, 
                Y -> SphericalHarmonicY} //. 
            values /. {z -> 
              0} //FullSimplify

VEPlot[pot3];

8.3 Laplace's Equation in Cartesian and Cylindrical Coordinates

Overview of Cartesian and Cylindrical Coordinates

Problem 1: Separation of variables in Cartesian and cylindrical coordinates

Required packages

Needs["Calculus`VectorAnalysis`"]

Solution

Clear["Global`*"];

Part a

potential= Vx[x] Vy[y] Vz[z];

SetCoordinates[Cartesian[x,y,z]];

lapEq= Laplacian[potential]

eq1 = (lapEq/(Vx[x] Vy[y] Vz[z])) // Expand

eq2 = eq1 /. Plus -> List

fieldEq=eq2=={-(Cx^2), -(Cy^2), (Cz^2)}//Thread

MapThread[
          DSolve, {fieldEq, {Vx[x], Vy[y], Vz[z]},  {x, y, z}}, 1] // 
        Flatten//ColumnForm

Part b

potential= R[r] P[ϕ] Z[z];

SetCoordinates[Cylindrical[r, ϕ, z]];

lapEq= Laplacian[potential] //Expand

eq1 = (lapEq/(R[r] P[ϕ] Z[z])) // Expand // 
        Collect[#, r] &

zEq= Select[eq1,(!FreeQ[#,z])&] == +Cz^2

DSolve[zEq, Z[z], z] // Flatten

DSolve[zEq /. Cz -> 0, Z[z], z] // 
      Flatten

eq2 =  eq1 /. {zEq /. Equal -> Rule}

eq3 = (r^2) eq2 // Expand

ϕEq = 
      Select[eq3, (! FreeQ[#, ϕ]) &] == -C(ϕ^2)

DSolve[ϕEq,  P[ϕ] , ϕ] // Flatten

DSolve[ϕEq /. C(ϕ^2) -> 0,  P[ϕ] , ϕ] // 
      Flatten

eq4 = eq3/.( ϕEq /. 
            Equal -> Rule)

DSolve[eq4==0, R[r] ,r]//Flatten

( DSolve[eq4 == 0 /. Cz -> 0, R[r], r] // Flatten ) /. 
          C[1] -> A + B /. C[2] -> I (A - B)  // 
      FullSimplify

Problem 2: Potential in a rectangular groove

Solution

Clear["Global`*"];

Part a

A[V_,m_]:= 2/a Integrate[V Sin[π m x/a],{x,0,a}];

pot[V_, n_]:=Sum[A[V, m] (E^(-((m π y)/a))) sin(((m π x)/a)),{m, 1, n}];

Part b

pot[Sin[x],3]

Plot3D[
        pot[Sin[x], 30] /. a -> 1 // Evaluate, 
{x, 0, 
          1}, {y, 0, 0.5}, 
        AxesLabel -> { x , y , V }];

pot15 = pot[Sin[x], 15] /. a -> 1;

VEPlot[pot15, {x, 0, 1}, {y, 0, 0.5} ];

Part c

eq1=pot[V0,∞]

potV0= eq1//. {log(1 - a_) -> 
                log(((1 - a)/(a + 1))) + log(a + 1)} // FullSimplify

Plot3D[potV0 /. a -> 1 /. V0 -> 1 , {x, 0.0001, 
          0.9999}, {y, 0, 0.5}, 
        AxesLabel -> { x , y , V }];

Problem  3:   Rectangular conduit

Remarks and outline

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"]

Solution

Clear["Global`*"];

term[m_]=A(
              m) cosh(((m π x)/b)) sin(((m π y)/b));
pot[n_] = 
        Sum[term[m] // Evaluate, {m, 1, n}];
pot[n]

check =  D[ pot[n] , x]  // MapAll[Evaluate, #] &

check  /. {x -> 0}

pot[n] /. {x -> a}

eq1 = V0 == term[m] /. {x -> a}

rhs = ((Integrate[eq1[[2]] Sin[m π y/b], {y, 0, b}] // 
              Simplify) //. {Sin[2 m π] -> 0})

lhs = ((Integrate[
                eq1[[1]] Sin[m π y/b], {y, 0, b}]) //. {Cos[
                  m π] -> (-1)^m}) // Simplify

aSol= Solve[rhs==lhs,A[m]][[1]] //Simplify;
A[m_]= A[m] /. aSol

pot[n]

Plot[
        Table[pot[n], {n, 2, 8, 2}] //. {b -> 1, V0 -> 1, 
              a -> 1, x -> a} // Evaluate, 
{y, 
          0, 1}, 
PlotRange -> {0, 2}];

Plot3D[
        pot[4] //. {b -> 1, V0 -> 1, 
            a -> 1}, 
{y, 0, 1}, {x, 0, 
          1}, 
ViewPoint -> {2, 2, 1}];

VEPlot[
        pot[4] //. {b -> 1, V0 -> 1, 
            a -> 1}, 
{x, 0, 1}, {y, 0, 
          1}];

Problem 4: Potential inside a rectangular box with five sides at zero potential

Solution

Clear["Global`*"];

Part a

mzRule={mz -> c (Sqrt[((mx^2)/(a^2)) + ((my^2)/(b^2))])};

A[mx_, my_]   = (((4 V0 csch(mz π))/(a b))Integrate[  sin(((mx π x)/a)) sin(((my π y)/b)),{x, 0, a},{y, 0, b}] //Simplify[#, Element[{my,  mx}, Integers]] & )/.mzRule

pot[nx_, ny_]=Sum[A[mx, my] sin(((mx π x)/a)) sin(((my π y)/b)) sinh(((mz π z)/c)),{mx, 1, nx},{my, 1, ny}]/.mzRule

pot[1, 1]

Part b

values = {a -> 1, b -> 1, c -> 1, V0 -> 1, 
          x -> 1/2};

Plot3D[
        pot[20, 20] //. values // Evaluate, {y, 0, 1}, {z, 0, 
          1}];

VEPlot[pot[10, 10] //. values //. {x -> 1/2}  // 
          Evaluate, {y, 0, 1}, {z, 0, 1}];

Problem 5: Conducting cylinder with a potential on the surface

Required packages

Needs["Graphics`ParametricPlot3D`"]

Solution

Clear["Global`*"];

Part a

A[n_] =((a^(-n))/π)(Integrate[
                V0  Sin[n ϕ], {ϕ, 0,  π }]  - 
              Integrate[
                V0  Sin[n ϕ], {ϕ, π, 2 π }] )//Simplify[#, Element[n, Integers]] &

potIn[nf_] := 
        Sum[A[n] Sin[n ϕ] (r^n), {n, 1, nf}];

potIn[6]

eq1= potIn[Infinity]

Vin=eq1//. {log(1 - aa_) -> 
                log(((1 - aa)/(aa + 1))) + log(aa + 1)}//FullSimplify

B[ n_] =((a^n)/π)(Integrate[
                V0  Sin[n ϕ], {ϕ, 
                  0,  π }] 
                     - 
              Integrate[
                V0  Sin[n ϕ], {ϕ, π, 2 π }] )//Simplify[#, Element[n, Integers]] &

potOut[nf_] := 
        Sum[B[n] Sin[n ϕ] (r^(-n)), {n, 1, nf}];

potOut[3]

Vout=(potOut[Infinity] // Simplify)//. {log(1 - aa_) -> 
                log(((1 - aa)/(aa + 1))) + log(aa + 1)}// Simplify

Part b

Off[Plot3D::plnc, Plot3D::gval]

pt1 = 
        Plot3D[ Vin /. a -> 1 /. V0 -> 1  // 
            Evaluate, 
{ϕ, 0, 2  π}, {r, 0, 
            1 }, 
PlotPoints -> 
            40, 
DisplayFunction -> 
            Identity];

pt2 = 
        Plot3D[ Vout /. a -> 1 /. V0 -> 1  // 
            Evaluate, 
{ϕ, 0, 2  π}, {r, 1 , 
            2 }, 
PlotPoints -> 
            40, 
DisplayFunction -> 
            Identity];

Off[Graphics3D::nlist3]

Show[{pt1, pt2}, 
        DisplayFunction -> $DisplayFunction];

Off[ParametricPlot3D::pplr]

pt3 = 
        CylindricalPlot3D[Vin /. a -> 1 /. V0 -> 1, {r, 0 , 
            1}, {ϕ, 
            0.01 ,   0.99 π}  , 
PlotPoints -> 20, DisplayFunction -> Identity];

pt4 = 
        CylindricalPlot3D[Vout /. a -> 1 /. V0 -> 1, {r, 1, 
            1.75}, {ϕ, 0.01, 
            0.99 π}, 
PlotPoints -> 20, 
          DisplayFunction -> Identity];

Show[{pt3 , pt4},  
        DisplayFunction ->  $DisplayFunction ];

8.4 Laplace's Equation in Spherical Coordinates

Overview of Spherical Coordinates

Problem 1: A charged ring

Required packages

Needs["Calculus`VectorAnalysis`"]

Solution

Clear["Global`*"];

Part a

potZ=(q/(Sqrt[(a^2) + (r^2)]));

potZin[n_] := Series[potZ, {r, 0, n}] // PowerExpand // 
          Normal;
potZout[n_] := Series[potZ, {r, ∞, n}] // PowerExpand // 
          Normal;

potZin[9]

potZout[9]

Part b

pot[s_]:=Sum[B(n) ((((r/a)))^n) + ((((a/r)))^(n + 1)) A(
                    n)LegendreP[n, Cos[θ]],{n, 0, s}]

Vout[s_]=Sum[((((a/r)))^(n + 1)) A(n)LegendreP[n, Cos[θ]],{n, 0, s}];

coeffOut[s_] := 
      Solve[0 == 
            CoefficientList[(Vout[s] - potZout[s]) /. θ -> 0, 
              1/r] , A  /@  Range[0, s]] // Flatten

coeffOut[3]

potOut[s_]:=Vout[s]/. coeffOut[s]

potOut[6]

Vin[s_]:=Sum[((((r/a)))^n) B(n)LegendreP[n, Cos[θ]],{n, 0, s}]

coeffIn[s_] := 
      Solve[0 == 
            CoefficientList[(Vin[s] - potZin[s]) /. θ -> 0,  
              r] , B /@ Range[0, s]] // Flatten

potIn[s_] := Vin[s] /. coeffIn[s]

potIn[6]

SetCoordinates[Spherical[r, θ, ϕ]];

electout[s_]:= -Grad[potOut[s]];

electin[s_]:= -Grad[potIn[s]];

electout[3]

electin[3]

Part c

r2xRule = {r, θ, ϕ} -> 
          CoordinatesFromCartesian[{x, y, z}, Spherical] // Thread

values = {a -> 1, q -> 1};

potOutXYZ[x_, y_, z_ ]  = potOut[12] //. r2xRule //. 
          values ;
potInXYZ[x_, y_, z_ ]  = potIn[12] //. r2xRule //. 
          values ;

ff[x_, y_, z_]:=If[N[(Sqrt[(x^2) + (y^2) + (z^2)])]≥1,potOutXYZ[x, y, z],potInXYZ[x, y, z]]

pt1 = 
        Plot3D[ff[x, 0, z], {x, -2, 2}, {z, -2, 2}, 
          PlotPoints -> 30 , 
          DisplayFunction -> Identity];

pt2 = 
        ContourPlot[ff[x, 0, z], {x, -2, 2}, {z, -2, 2}, 
          PlotPoints -> 30,   ContourSmoothing -> True, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{pt1, pt2}]];

VEPlot[
        ff[x, 0, z] // Evaluate, 
{x, -2.1, 
          2.1}, {z, -1.6, 1.6}, 
        Epilog -> {{Hue[ .3],  Disk[{1, 0},  .2]},   {Hue[ .3],  
              Disk[{-1, 0 },  .2]}}];

pt3 = 
        Plot3D[ff[x, y, 0], {x, -2, 2}, {y, -2, 2}, 
          PlotPoints -> 30, 
          DisplayFunction -> Identity];

pt4 = 
        ContourPlot[
          ff[x, y, 0], 
{x, -2, 2}, {y, -2, 
            2}, 
PlotPoints -> 
            30, 
ContourSmoothing -> True, 
          Epilog -> {Hue[0.07], Thickness[ .03], 
              Circle[{0, 0}, 1]}, 
DisplayFunction -> 
            Identity];

Show[GraphicsArray[{pt3, pt4}]];

VEPlot[
        ff[x, y, 0], 
{x, -2.1, 2.1}, {y, -2.1, 
          2.1}, Epilog -> {Hue[0.05], Thickness[ .008], 
            Circle[{0, 0}, 1]}];

Problem 2: Grounded sphere in an electric field

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"]

Solution

Clear["Global`*"];

Part a

Asol = { A[1] -> -E0 a,  
          A[i_] -> 0}  ;

pot[s_]:=Sum[(A(
                        i) ((((r/a)))^i) + ((((a/r)))^(i + 1)) B(
                        i) )LegendreP[i, 
                Cos[θ]],{i, 0, s}]/.Asol

Beq[
        s_]:=LogicalExpand[(pot[s] /. r -> a) == 
          0 + O[Cos[θ]]^(s + 1)]

Beq[4]

Bsol[s_] := Solve[Beq[s]] // Flatten ;

Off[Solve::svars]

Solve[Beq[6]]//Flatten

Bsol = {B[1] -> a E0, B[i_] -> 0};

potential =(Sum[(A(
                            i) ((((r/a)))^i) + ((((a/r)))^(i + 1)) B(
                            i)) LegendreP[i, Cos[θ]],{i, 0, 3}])/.Asol/.Bsol

eField = -Grad[potential, Spherical[r, θ, ϕ]] // 
        ExpandAll

Part b

r2xRule = {r, θ, ϕ} -> 
            CoordinatesFromCartesian[{x, y, z}, Spherical] // 
          Thread;

values = {E0 -> 1, a -> 1};

p1 = 
        VEPlot[potential //. r2xRule //. 
              values //. {y -> 0}, 
{x, -2.1, 
            2.1}, {z, -2.1, 
            2.1}, 
DisplayFunction -> 
            Identity];

disk = Graphics[Disk[{0, 0}, 1]];

Show[p1, disk, 
        DisplayFunction -> $DisplayFunction ];

Problem 3: Sphere with an axially symmetric charge distribution

Required packages

Needs["Calculus`VectorAnalysis`"]
 Needs["Graphics`ParametricPlot3D`"]

Solution

Clear["Global`*"];

Part a

potInA[s_]:=Sum[((((r/a)))^i) A(i) LegendreP[i, Cos[θ]],{i, 0, s}]

potOutA[s_]:=Sum[A[i] ((((a/r)))^(i + 1))LegendreP[i, Cos[θ]],{i, 0, s}]

chargeDensity[s_]:=(1/(4 π))(D[potInA[s], r] - D[potOutA[s], r])//.{r -> a}

chargeDensity[2] // Apart

aEquation[s_, Qdensity_] := 
      LogicalExpand[
        chargeDensity[s] == 
          Qdensity + O[Cos[θ]]^(s + 1)]

aRule[s_, Qdensity_] := 
      Solve[aEquation[s, Qdensity], Array[A, s + 1, 0]] // Flatten

aRule[4,((Q (cos^2)(θ))/(a^2))]

potInA[4]/.aRule[4,((Q (cos^2)(θ))/(a^2))]

potOutA[4]/.aRule[4,((Q (cos^2)(θ))/(a^2))]

Part b

surfaceQ=((Q ((((cos^2)(θ))/2) + 
                    cos(θ)))/(a^2))//Expand

potOutQ= potOutA[8] //. aRule[8,surfaceQ]

potInQ= potInA[8]//.aRule[8,surfaceQ]

values = {Q -> 1, a -> 1};

p1 = 
        CylindricalPlot3D[
          potInQ //. values, {r, 0, 1}, {θ, 0, 2 π}, 
          DisplayFunction -> Identity];

p2 = 
        CylindricalPlot3D[
          potOutQ //. values, {r, 1, 3}, {θ, 0, 2  π}, 
          DisplayFunction -> Identity];

Show[p1, 
        p2, 
AxesLabel -> {a, b, 
            c}, 
ViewPoint -> {-1.5, -2, 
            2}, 
BoxRatios -> {1, 1, 1}, 
        DisplayFunction -> $DisplayFunction];

r2xRule = {r, θ, ϕ} -> 
            CoordinatesFromCartesian[{x, y, z}, Spherical] // 
          Thread;

eq1[y_, 
          z_] = potInQ /. r2xRule /. {Q -> 1, a -> 1, 
              x -> 0} // Simplify;

\(\(\( (*\ 
      for\ 1 < r\ *) \)\(
\)\(eq2[y_, 
          z_] = \(potOutQ /. r2xRule\) /. {Q -> 1, a -> 1, 
              x -> 0}\  // Simplify;\)\)\)

pot[y_, z_]:=If[(y^2) + (z^2) < 1,eq1[y, z],eq2[y, z]]

circle = {Thickness[0.02], Circle[{0, 0}, 1]}  // 
          Graphics;

Show[{VEPlot[pot[y, z], {y, -2, 2}, {z, -2, 2}, 
            DisplayFunction -> 
              Identity], 
circle}, 
DisplayFunction -> $DisplayFunction];

Problem 4: Sphere with a given axially symmetric potential

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`ParametricPlot3D`"]

Solution

Clear["Global`*"];

Part a

potIn[V_, s_]:=Sum[((((r/r0)))^n) B(V, n) P[n, Cos[θ]],{n, 0, s}]

potOut[V_, s_]:=Sum[((((r0/r)))^(n + 1)) B(V, n) P[n, Cos[θ]],{n, 0, s}]

B[V_, n_]:=(1/2) (2 n + 1) Integrate[
          V LegendreP[n, Cos[θ]] Sin[θ], {θ, 
            0, π}]

Part b

pOut = potOut[V0 Cos[θ]^2, 4]

pIn = potIn[V0 Cos[θ]^2, 4]

values = {r0 -> 1, V0 -> 1, 
          P -> LegendreP};

p1 = 
        CylindricalPlot3D[pIn /. values, {r, 0, 1}, {θ, 0, 2  π}, 
          DisplayFunction -> Identity];

p2 = 
        CylindricalPlot3D[pOut /. values, {r, 1, 3}, {θ, 0, 2  π}, 
          DisplayFunction -> Identity];

Show[p1, p2, BoxRatios -> {1, 1, 1}, 
        ViewPoint -> {-1.5, -2, 2}, 
        DisplayFunction -> $DisplayFunction];

Part c

surfacePot=((cos(θ))/(1.1 - 
                cos(θ)));

{pIn2, 
          pOut2} = ({potIn[surfacePot, 3], potOut[surfacePot, 3]} /. 
                P -> LegendreP //. {r0 -> 1} // 
            Simplify);

p1 = 
        CylindricalPlot3D[pIn2, {r, 0, 1}, {θ, 0, 2  π}, 
          PlotRange -> All, 
          DisplayFunction -> Identity];

p2 = 
        CylindricalPlot3D[pOut2, {r, 1, 3}, {θ, 0, 2  π}, 
          PlotRange -> All, 
          DisplayFunction -> Identity];

Show[p1, 
        p2, 
BoxRatios -> {1, 1, 
            1}, 
ViewPoint -> {-1.5, -2, 2}, 
        DisplayFunction -> $DisplayFunction];

Problem 5: Sphere with upper hemisphere V0 and lower hemisphere -V0

Required packages

Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`ParametricPlot3D`"]

Solution

Clear["Global`*"];

Part a

potIn[s_]:=Sum[((((r/r0)))^n) B(n) P[n, Cos[θ]],{n, 0, s}]

potOut[s_]:=Sum[((((r0/r)))^(n + 1)) B(n) P[n, Cos[θ]],{n, 0, s}]

B[n_]:=(1/2) (2 n + 1)(+Integrate[
                  V0 LegendreP[n, Cos[θ]] Sin[θ], {θ, 
                    0, π/2}] + 
              Integrate[-V0 LegendreP[n, Cos[θ]] 
                  Sin[θ], {θ, π/2, π}]);

potOut[5]

Part b

values = {r0 -> 1, V0 -> 1, x -> 0};

r2xRule = {r, θ, ϕ} -> 
            CoordinatesFromCartesian[{x, y, z}, Spherical] // 
          Thread;

allRules =  {P -> LegendreP ,  r2xRule, values} // 
          Flatten ;

eq1 =  potIn[11] //. allRules;

eq2 = potOut[11] //. allRules;

potTot=If[(y^2) + (z^2) ≤ 1,eq1,eq2];

circle= Graphics[{Thickness[0.02],Circle[{0,0},1]}];

Show[{
ContourPlot[
            potTot // Evaluate, 
{y, -2, 
              2}, {z, -2, 2}, 
PlotPoints -> 
              20, 
DisplayFunction -> 
              Identity], 
circle}, 
DisplayFunction -> $DisplayFunction];

8.5 Exercises

Exercise 1 Parallel plates and strips

Exercise 2 Parallel plates and a point charge

Exercise 3 Cylinder divided into segments

Exercise 4 Elliptical cylinder in a uniform electric field

Exercise 5 Cylindrical box and n charges

Exercise 6 Distorted sphere with a potential

Exercise 7 Concentric spheres with different potentials on the two hemispheres

Exercise 8 Ring charge

Exercise 9 Laplace's equations in other coordinates