Solutions for Problem Set 2

Problem 1.23 in Jackson

4 point calculation of the potential inside a cylindrical capacitior with square cross section

Corner sum

In this subsection, we let   ϕ[i]=ϕ[x,y] and c[i] = c[x,y] =  ϕ[x+h,y]+ϕ[x-h,y]+ϕ[x,y+h]+ϕ[x,y-h].
Then solve 4 ϕ[i] = c[i]

Clear["`*"]

c[1]  :=  2 ϕ[2 ] +  V

c[2]  :=  ϕ[1] + ϕ[3] + V

c[3]  :=  ϕ[2] + ϕ[4]  +  V

c[4]  :=  2  ϕ[3]

eqc := Table[c[i] == 4  ϕ[i], {i, 1, 4}]

vars :=  Table[ϕ[i], {i, 1, 4}]

solc =  Solve[eqc, vars] // First

vars  /.  solc  /.  V -> 100 // N

Side sum

In this section, we let s[i] = ϕ[x+h,y+h]+ϕ[x+h,y-h]+ϕ[x-h,y+h]+ϕ[x-h,y-h].
We then set up the average a[i] = (4c[i]+s[i])/20 and solve  ϕ[i] = a[i]

s[1]  :=  2 V

s[2]  :=  2 V

s[3 ] :=  V +  ϕ[3]

s[4] :=  V

Table[a[i]  =  (4 c[i] + s[i])/20, {i, 1, 4}]

eqa =  Table[a[i] == ϕ[i], {i, 1, 4}]

sola =  Solve[eqa, vars] // First

The resulting potential array (given V=100) is

N[vars  /.  sola  /.  V -> 100, 
      4]

Compare with Jackson's "exact" values  {48.87, 47.18, 38.34, 19.81}

Surface charge on Inner Surface

First determine surface charge density on inner surface. We use a 3-point differentiation formula 
f'[1] = (-3f[1]+4f[2]-f[3])/2h  
to determine normal component of the electric field. The points [1], [2], etc, are points on the inner surface hust above the points [1], [2], ... at which the potential was eveluated.

σ[
        1]  =  -(-3 V + 4  ϕ[1]) (Subscript[ϵ, 0])/(2  
                h) /.  sola

σ[
        2] =  -(-3 V + 4  ϕ[2]) (Subscript[ϵ, 0])/(2  
                h) /.  sola

σ[
        3]  =  -(-3 V + 4  ϕ[3]) (Subscript[ϵ, 0])/(2  
                h) /.  sola

ListPlot[
      Table[h σ[i]/(Subscript[ϵ, 0])  /. V -> 100, {i, 1, 3}], 
      PlotStyle ->  {RGBColor[1, 0, 0], PointSize[0.03]}, 
      PlotRange -> {{0, 5}, {0, 100}}]

We use the Simpson's rule to evaluate the integral Q= charge/length on the inner surface.  We integrate over 1/8th of the surface and multiply by 8

qi =  
      8 h (σ[1] + 4 σ[2] + σ[3])/3

Next, we determine charge/length on outer surface

σ'[
        1] = -(4  ϕ[1] - V) (Subscript[ϵ, 0])/(2  h) /.  
        sola

σ'[
        2] = -(4  ϕ[2] - V) (Subscript[ϵ, 0])/(2  h) /.  
        sola

σ'[
        3] = -(4  ϕ[3] - V) (Subscript[ϵ, 0])/(2  h) /.  
        sola

σ'[
        4] = -(4  ϕ[4] - ϕ[3]) (Subscript[ϵ, 0])/(2  
                h) /.  sola

ListPlot[
      Table[h σ'[i]/(Subscript[ϵ, 0])  /. V -> 100, {i, 1, 4}], 
      PlotStyle ->  {RGBColor[0, 0, 1], PointSize[0.03]}, 
      PlotRange -> {{0, 5}, {-100, 0}}]

Digression on semi-open Newton Cotes formula.

We desire to integrate a function f[x] from 0 to 4h, but the value of f[0] is unknown. To remedy the situation, we run a cubic interpolating polynomial through the points f[h], f[2h], f[3h], f[4h] and integrate the interpolation polynomial from 0 to 4h.

p  =  Product[(x - i h), {i, 1, 4}]

pol  =  0

Do[{num = p/(x - i h), den = num /. x -> i h,  
        pol = pol + num f[i]/den}, {i, 1, 4}]

pol

Integrate[pol, {x, 0, 4  h}]

Total charge on outer surface using above quadrature rule & Capacitance

qo = 
      8 h (4/3) (2 σ'[4] - σ'[3] + 
            2  σ'[2])

Note that the prediction of the total charge on each surbace is different. We use an average

qi/((Subscript[ϵ, 0]) V) // N

qo/((Subscript[ϵ, 0]) V) // N

q  =  (qi - qo)/2

The capcitance/length  (in units of (Subscript[ϵ, 0])) is

N[c  =  q/((Subscript[ϵ, 0]) V), 4]

21 point calculation & more!

Initialize the potential array

The potential is defined on a grid ϕ[i,j], i=0..8, j=0..4.   In the code below,  m =  4 k, k = 1,2,3,...
Put m = 4 to solve the previous 4-point problem,  put m = 8 to solve the  21-point problem, put m = 16 to 
halve the grid once again.

As a first step, we fill in boundary values.

Clear["`*"]

m = 32

n = m/2

Do[ϕ[i, 0] = 0, {i, 0, m + 1}]

Do[ϕ[0, j] = 0, {j, 1, n}]

Do[ϕ[i, n] = V, {i, n, m + 1}]

Do[Do[ϕ[i, j] = ϕ[j, i], {i, 1, j - 1}], {j, 2, 
        n}]

Do[ϕ[m + 1, j] = ϕ[m - 1, j], {j, 1, n - 1}]

Table[ϕ[i, j], {i, 0, m + 1}, {j, 0, n}] // TableForm

Corner sums

Do[c[i, 
          j] = ϕ[i - 1, j] + ϕ[i + 1, j] + ϕ[i, 
            j - 1] + ϕ[i, j + 1], {i, 1, m}, {j, 1, n - 1}]

vars  :=  
      Table[Table[ϕ[i, j], {i, j, m}], {j, 1, n - 1}] // 
        Flatten

Side Sums

Now evaluate the diagonal sums

Do[s[i, 
          j] = ϕ[i - 1, j - 1] + ϕ[i + 1, j + 1] + ϕ[i - 1, 
            j + 1] + ϕ[i + 1, j - 1], {i, 1, m}, {j, 1, 
        n - 1}]

Average the two and solve for potential

Form the average array

Do[a[i, j] = (4  c[i, j] + s[i, j])/20, {i, 1, m}, {j, 1, 
        n - 1}]

Set up the equations

eqa  :=  
      Table[Table[a[i, j] == ϕ[i, j], {i, j, m}], {j, 1, n - 1}] // 
        Flatten

Solve the equations exactly

sola =  Solve[eqa, vars] // First

Look at the solutions

N[Table[ϕ[i, j] /.  sola /.  V ->  100 , {i, 0, 
            m + 1}, {j, 0, n}], 4] // TableForm

ListPlot3D[
      Table[ϕ[i, j] /.  sola /.  V ->  100  // N, {j, 0, 
          n}, {i, 0, m + 1}], BoxRatios -> {8, 4, 5}]

Surface Charges

Use the 3-point differetiator to determine charge on inner surface

Do[σ[
          i]  =  -(-3 V + 
                  4  ϕ[i, n - 1] - ϕ[i, n - 2])/(2  h) /.  
          sola, {i, n, m + 1}]

Look at the inner surface charges

ListPlot[Table[σ[i], {i, n, m + 1}] /.  
          V -> 100  /.  h -> 1, 
      PlotStyle -> {PointSize[0.03], RGBColor[1, 0, 0]}, 
      PlotStyle -> {PointSize[0.03], RGBColor[1, 0, 0]}, 
      PlotRange -> {0, 320/m}]

Integrate to find total inner surface charge

If[m == 4, 
      qi = 8  h (σ[2] +  4  σ[3] + σ[4])/3 , 
      qi = 8  h ((σ[n] +  
                    4  σ[n + 1] + σ[n + 2])/
                3
 +  
              Sum[(σ[i] + 4  σ[i + 1] + σ[i + 2])/
                  3, {i, n + 2, m - 2, 2}]) ]

qi/V // N

Find charge on outer surface (here we find the negative of the cgarge density)

Do[σ'[
          i] = (4  ϕ[i, 1] - ϕ[i, 2])/(2  h)  /.  
          sola, {i, 1, m + 1}]

Look at the charge distribution

ListPlot[Table[σ'[i], {i, 1, m + 1}] /.  
          V -> 100  /.  h -> 1, 
      PlotStyle -> {PointSize[0.03], RGBColor[0, 1, 0]}, 
      PlotRange -> {{0, m + 1},  {0, 240/m}}]

Integrate to find the  total outer surface charge

If[m == 4, 
      qo  =  8  
          h ( 4 (2  σ'[1] - σ'[2] + 
                    2  σ'[3])/3 ), 
      qo  = 8  
          h ( 4 (2  σ'[1] - σ'[2] + 
                      2  σ'[3])/3  +  
              Sum[(σ'[i] + 
                      4  σ'[i + 1] + σ'[i + 2])/
                  3, {i, 4, m - 2, 2}])]

qo/V // N

Find the average charge

q = (qi + qo)/2

Find the Capacitance

Cap = q/V // N

Problem 1.24 in Jackson

Three-point solution

Set up and solve equations

Notation: Potential at bottom-left ϕ[1], bottom-center ϕ[2], center ϕ[3]

Clear["`*"]

h  =  1/4

We choose ρ = 4π in order to compare with Jackson's answers to this problem.

ρ =  4 π

c[1] := 2  ϕ[2]

c[2] := 2  ϕ[1] + ϕ[3]

c[3] := 4 ϕ[2]

s[1] := ϕ[3]

s[2] := 2 ϕ[2]

s[3] := 4 ϕ[1]

Table[a[i] = (4 c[i] + s[i])/5, {i, 1, 3}]

Exact solution

Multiply ρ by ((  12)/10) (h^2)  (If we use c[i] instead, then we must multiply by (h^2))

eqs := 
      Table[4  ϕ[i] == a[i] +   6 h^2 ρ /5 , {i, 1, 
          3}]

var := Table[ϕ[i], {i, 1, 3}]

sol = Solve[eqs, var] // First

t = 
      Table[{{0, 0, 0, 0, 0}, {0, ϕ[1], ϕ[2], ϕ[1], 
              0}, {0, ϕ[2], ϕ[3], ϕ[2], 
              0}, {0, ϕ[1], ϕ[2], ϕ[1], 0}, {0, 0, 0, 0, 0}}] /. 
        sol

ListPlot3D[t]

N[t, 3] // TableForm

Jacobi Solution

Clear["`*"]

h = 1/4

ρ =  4 π

For[{k = 1, ϕ[1] = 1, ϕ[2] = 1, ϕ[3] = 1}, 
      k < 20, k++, {s[1] = (8 ϕ[2] + ϕ[3])/5, 
        s[2] = (2  ϕ[2] + 4 (2 ϕ[1] + ϕ[3]))/5, 
        s[3] = (4 ϕ[1] + 16 ϕ[2])/5, 
        Table[ϕ[i] = (s[i] +  6 h^2 ρ/5 )/4, {i, 1, 
            3}]}]

N[Table[{{0, 0, 0, 0, 0}, {0, ϕ[1], ϕ[2], ϕ[1], 
              0}, {0, ϕ[2], ϕ[3], ϕ[2], 
              0}, {0, ϕ[1], ϕ[2], ϕ[1], 0}, {0, 0, 0, 0, 0}}], 
        3] // TableForm

Gauss-Seidel Solution

Clear["`*"]

h = 1/4

ρ =  4 π

σ =  6 h^2 ρ/5

For[{k = 1, ϕ[1] = 1, ϕ[2] = 1, ϕ[3] = 1}, 
      k < 20, k++, {ϕ[
            1] = ((8 ϕ[2] + ϕ[3])/5 + σ)/4, ϕ[
            2] = ((2  ϕ[2] + 4 (2 ϕ[1] + ϕ[3]))/
                  5 + σ)/4, ϕ[
            3] = ((4 ϕ[1] + 16 ϕ[2])/5 + σ)/
            4}]

N[Table[{{0, 0, 0, 0, 0}, {0, ϕ[1], ϕ[2], ϕ[1], 
              0}, {0, ϕ[2], ϕ[3], ϕ[2], 
              0}, {0, ϕ[1], ϕ[2], ϕ[1], 0}, {0, 0, 0, 0, 0}}], 
        3] // TableForm

General Case

Potential for uniform charge density in a  square region  with V=0 on boundary (Exact Solution)

Initialize potential array

Clear["`*"]

n = 4

h = 1/n

Do[{ϕ[i, 0] = 0, ϕ[i, n] = 0}, {i, 0, n}]

Do[{ϕ[0, j] = 0, ϕ[n, j] = 0}, {j, 1, n - 1}]

Table[ϕ[i, j], {i, 0, n}, {j, 0, n}] // TableForm

Evaluate corner sum

Do[c[i, 
          j] = ϕ[i - 1, j] + ϕ[i + 1, j] + ϕ[i, 
            j - 1] + ϕ[i, j + 1], {i, 1, n - 1}, {j, 1, 
        n - 1}]

Table[c[i, j], {i, 1, n - 1}, {j, 1, n - 1}] // TableForm

Set up equations and variables

eqs  =  
      Table[ϕ[i, j] == (c[i, j] + h^2  ρ)/4, {i, 1, n - 1}, {j, 
            1, n - 1}] // Flatten

vars  =  
      Table[ϕ[i, j], {i, 1, n - 1}, {j, 1, n - 1}] // Flatten

Solve and plot

soln  =  Solve[eqs, vars] // First

ListPlot3D[Table[ϕ[i, j], {i, 0, n}, {j, 0, n}] /.  
          soln /.  ρ -> 1 ]

4 π Max[
          Table[ϕ[i, j] /.  soln /. ρ -> 1, {i, 0, 
                n}, {j, 0, n}] // Flatten] // N

Jacobi Solution

Clear["`*"]

n = 6

h = 1/n

ρ = 4  π

For[{Do[ϕ[i, j] = 0, {i, 0, n}, {j, 0, n}], k = 0}, 
      k < 20, k++, {Do[
          c[i, j] = ϕ[i - 1, j] + ϕ[i + 1, j] + ϕ[i, 
                j - 1] + ϕ[i, j + 1], {i, 1, n - 1}, {j, 1, n - 1}], 
        Do[ϕ[i, j] = (c[i, j] + h^2  ρ)/4, {i, 1, n - 1}, {j, 
            1, n - 1}], 
        ListPlot3D[Table[ϕ[i, j], {i, 0, n}, {j, 0, n}], 
          PlotRange -> {{1, n + 1}, {1, n + 1}, {0, 1}}]}]

Max[
        Table[ϕ[i, j] /.  soln /. ρ -> 1, {i, 0, n}, {j, 
              0, n}] // Flatten] // N

Gauss-Seidel Solution

Clear["`*"]

n = 6

h = 1/n

ρ = 4  π

For[{k = 0, Do[ϕ[i, j] = 0, {i, 0, n}, {j, 0, n}]}, 
      k < 20, k++, {Do[ϕ[i, 
                j] = (ϕ[i - 1, j] + ϕ[i + 1, j] + ϕ[i, 
                      j - 1] + ϕ[i, j + 1] + h^2 ρ)/4, {i, 1, 
              n - 1}, {j, 1, n - 1}] 
          ListPlot3D[Table[ϕ[i, j], {i, 0, n}, {j, 0, n}], 
            PlotRange -> {{1, n + 1}, {1, n + 1}, {0, 1}}]}]

Max[
        Table[ϕ[i, j] /.  soln /. ρ -> 1, {i, 0, n}, {j, 
              0, n}] // Flatten] // N