Homework  Assignment II, Phy711  Due Oct 7 (based on book readings Chaps 1, Chaps 2.1-2.6,2.7,2.8-2.10)

Note: Submit solutions as a Mathematica notebook (no handwritten work will be accepted !). No late submissions.

(1)  A cylinder, of radius R, has  charge density ρ(rϕz) = (Subscript[ρ, 0])r, where rϕz are a set of cylindrical coordinates.

(a) Find the electric potential, and field everywhere in space.

(b) Plot, using Mathematica, the above fields.

By symmetry we assume the electric field is always along the radial direction. We use Gauss’ Law to find its magnitude. In the region outside the cylinder the left hand side of the flux integral is

∫(OverVector[E])·d(OverHat[a]) = L 2π r  E, whereas the r.h.s is  (1/ϵ) ∫ρdV =  (1/ϵ) 2π L (∫_0^R)(Subscript[ρ, 0])r  rdr =   (2/ϵ) π L (Subscript[ρ, 0]) (R^3)/3  

Note that for a length L of the cylinder the total charge is ∫ρdV = L (Subscript[ρ, 0]) 2π (R^3)/3   so the charge per unit length
along the cylinder is  λ = L (Subscript[ρ, 0]) 2π (R^3)/3/L = (Subscript[ρ, 0]) 2π (R^3)/3
so

E = (λ/(2π ϵ r)) outside the cylinder  (along the radial direction).  

Also ϕ = -∫ (OverVector[E])·d(OverHat[r]) = - (λ/(2π ϵ ))log r + (Subscript[C, 1])  for r>R, where (Subscript[C, 1])is a constant.

For r<R repeating as above

L 2π r  E =  (1/ϵ) 2π L (∫_0^r)(Subscript[ρ, 0])r  rdr    or

E = ((Subscript[ρ, 0])/ϵ) (r^2) /3 .  We obtain  ϕ = -∫ (OverVector[E])·d(OverHat[r]) = -((Subscript[ρ, 0])/ϵ) (r^3) /9  + (Subscript[C, 2])

Now we would like ϕ to be continuous across the bound r=R and so

we require 

- (λ/(2π ϵ ))log R  + (Subscript[C, 1])=  -((Subscript[ρ, 0])/ϵ) (R^3) /3  + (Subscript[C, 2])

suppose we set (Subscript[C, 2])=0 (so that ϕ at the origin vanishes) then 

- (λ/(2π ϵ ))log R + (Subscript[C, 1])=  -((Subscript[ρ, 0])/ϵ) (R^3) /9  ⟹ (Subscript[C, 1])= (λ/(2π ϵ ))log R  - ((Subscript[ρ, 0])/ϵ) (R^3) /9

To summarize  (for E along radial direction)

E = λ/(2 π ϵ r)    r > R

E = (Subscript[ρ, 0])/ϵ   (r^2)/3    r < R

and

ϕ =  - λ/(2 π ϵ )log r  + λ/(2 π ϵ )log R  -  (Subscript[ρ, 0])/ϵ   (R^3)/9             r > R

ϕ =    - (Subscript[ρ, 0])/ϵ   (r^3)/9        r <  R

(* according to above discussion *)
ϵ=1;
R=2;
ρ0=1;
λ=ρ0 2 π R^3/3;
(* we set some constants to an arbitrary value *)
ϕout[x_,y_]=-λ/(2 π ϵ)Log[(Sqrt[x^2+y^2])] +  λ/(2 π ϵ) Log[R]-((ρ0)/ϵ)R^3/9
Eout[x_,y_]=-{D[ϕout[x,y],x],D[ϕout[x,y],y]}
ϕin[x_,y_]=-((ρ0)/ϵ)(((x^2+y^2))^((3/2)))/9
Ein[x_,y_]=-{D[ϕin[x,y],x],D[ϕin[x,y],y]}

ϕ[x_,y_]:=Piecewise[{{ϕout[x,y],x^2+y^2>R^2},{ϕin[x,y],x^2+y^2<=R^2}}]

ContourPlot[ϕ[x,y],{x,-5,5},{y,-5,5},ContourShading->Off,Contours->20]
(* this does not look that interesting, lets plot the potential as the third coordinate in a 3 D plot *)

g1=Plot3D[ϕ[x,y],{x,-5,5},{y,-5,5}]

(* for the E field let's look at the gradient *)

Efield[x_,y_]:=Piecewise[{{Eout[x,y],x^2+y^2>R^2},{Ein[x,y],x^2+y^2<=R^2}}]

g2=VectorPlot[Efield[x,y],{x,-5,5},{y,-5,5}]

(2) Consider now a situation in which the cylinder given above, is adjacent to a conducting wall. The axis of the cylinder is parallel to the wall at a constant distance d>R.

(a) Use the method of images to find the potential in the region that includes the cylinder and is bounded by the surface of the conductor.

If the cylinder is displaced a distance d along the x-axis then the result given in question (1) above is

 (Subscript[ϕ, 1])(x,y) =  - (λ/(2π ϵ ))log (Sqrt[(  (((x - d))^2) + (y^2) ) ]) Log R  - ((Subscript[ρ, 0])/ϵ) (R^2) /4   in the region outside the cylinder

and

-((Subscript[ρ, 0])/ϵ) (  (((x - d))^2) + (y^2) ) /4   in the region inside the cylinder

We can call this potential (Subscript[ϕ, 1])(x,y). Now we put an image cylinder (with opposite charge density (Subscript[ρ, 0]))  at the position x=-d and so

(Subscript[ϕ, 2])(x,y) =  (λ/(2π ϵ ))log (Sqrt[(  (((x + d))^2) + (y^2) ) ]) Log R  + ((Subscript[ρ, 0])/ϵ) (R^2) /4  ; in the region outside the cylinder  (we need not consider the
case inside the image since we are only concerned with the region x>0)

Therefore

ϕ(x,y)= (Subscript[ϕ, 1])(x,y)+(Subscript[ϕ, 2])(x,y)

We are assured this is the correct answer since (∇^2)ϕ = -ρ/ϵ  in the region x>d and ϕ vanishes on the conductor on the x=0 plane.

ϵ=1;
R=2;
ρ0=1;
λ=ρ0 2 π R^3/3;
d=2R;

ϕ1out[x_,y_]=-λ/(2 π ϵ)Log[(Sqrt[(x-d)^2+y^2])] +  λ/(2 π ϵ) Log[R]-((ρ0)/ϵ)R^3/9;
ϕ1in[x_,y_]=-((ρ0)/ϵ)((((x-d)^2+y^2))^((3/2)))/9;
ϕ2out[x_,y_]=λ/(2 π ϵ)Log[(Sqrt[(x+d)^2+y^2])] -  λ/(2 π ϵ) Log[R]+((ρ0)/ϵ)R^3/4;

ϕtotal[x_,y_]:=Piecewise[{{ϕ1out[x,y]+ϕ2out[x,y],(x-d)^2+y^2>R^2},{ϕ1in[x,y]+ϕ2out[x,y],(x-d)^2+y^2<=R^2}}]

Plot3D[ϕtotal[x,y],{x,0,10},{y,-5,5}]

(* let's look at the electric field *)

-Grad[ϕ1out[x,y]+ϕ2out[x,y],{x,y}];
g3=StreamPlot[%,{x,0,10},{y,-5,5},RegionFunction->Function[{x,y},(x-d)^2+y^2 > R^2]];
-Grad[ϕ1in[x,y]+ϕ2out[x,y],{x,y}];
g4=StreamPlot[%,{x,0,10},{y,-5,5},RegionFunction->Function[{x,y},(x-d)^2+y^2 < R^2],VectorPoints->10,StreamStyle->Red];
Show[g3,g4]

(b) Find the Dirichlet Green's function for this problem.

Since this is a problem involving cylindrical symmetry we require the analog of the point charge system for three dimensions. That
analog is an infinite wire of line charge density λ a distance x=d (from the conductor) along the x-axis and stretched along the z-axis
(coming out of the plane of the page).
We know from Phy 181 that the potential of such a wire, located at x',y', is

ϕ(x,y)= -(λ/(2π ϵ))Log( (Sqrt[ (((x - x'))^2)+(((y - y'))^2)]) ) thus we expect the 2D Greens function to have a similar structure (in analogy to the 3D Coulomb
potential giving us the 3D Green's function) 

Let's re-derive this using the methods discussed in this course. First we note that an infinitely long wire, located at position
x'',y'' and extending into the z' direction can be represented by a charge density

ρ(x',y',z')=λ δ(x'-x'')δ(y'-y'') and using  Eq. (1.42) in the text we have

ϕ(x,y,z) = (1/(4π ϵ))∫ G(xyz,x'y'z') ρ(x',y',z') dx'dy'dz'
where

G = 1/ (Sqrt[(((x-x'))^2)+(((y-y'))^2)+(((z-z'))^2)]) , now we note that ρ does not depend on z' so we can re-write the above expression as

 (1/(4π ϵ))∫ G(xy,x'y' ) ρ(x',y') dx'dy'   (i.e a 2D integral where we defined a 2D green's function
 
 
G(xy,x'y' ) = (∫_(-∞)^∞) dz'/(Sqrt[(((x-x'))^2)+(((y-y'))^2)+(((z-z'))^2)]) =   (∫_(-∞)^∞)du /(Sqrt[(((x-x'))^2)+(((y-y'))^2)+(u^2)])

Unfortunately this integral does not converge, so we use a method called dimensional regularization to tame this integral. Consider the following 
integral
 (∫_(-∞)^∞)du ((((x-x'))^2)+(((y-y'))^2)+(u^2))^(-1/2-ϵ)    for ϵ=0 we have the integral above, Lets take ϵ a finite number and let ϵ ⟶0 after performing the integral
 
 we then get, (looking in a table of integrals)
 
  (∫_(-∞)^∞)du ((a^2)+(u^2))^(-1/2-δ)   = (((Sqrt[π]) (a^(-2 δ)) (Gamma[δ]))/(Gamma[δ+(1/2)]))
  
  taking the limit as  δ ⟶0  we get 
  
   (∫_(-∞)^∞)du ((a^2)+(u^2))^(-1/2-δ) ⟶ 1/δ  + Log(4) - Log((a^2))  =  1/δ  + Log(4) - Log((((x-x'))^2)+(((y-y'))^2))  where we replaced (((x-x'))^2)+(((y-y'))^2)
   with a.
   
   Inserting this into  (1/(4π ϵ))∫ G(xy,x'y' ) ρ(x',y') dx'dy'  the expression for  ϕ becomes
   
   (λ/(4π ϵ)) ( 1/δ + Log(4) - Log( (((x-x''))^2)+(((y-y''))^2)) ) =  - (λ/(2π ϵ)) Log( (Sqrt[(((x-x''))^2)+(((y-y''))^2)]))  +   (λ/(4π ϵ)) ( 1/δ + Log(4))
   
   The first term agrees with the answer we obtained using Phys 181 methods, however the second, constant, term approaches ∞ as
   δ  ->  0. What to do?  Well it's just a constant and the potential is arbitrary up to a constant and we therefore ignore it even though it approaches ∞.  Therefore we define, the regularized 2D Greens function,
   
   G(xy,x'y' )= - Log((((x-x'))^2)+(((y-y'))^2)) 

  putting the image wire at the reflected point -x',y' we have

(Subscript[G, D])(xy,x'y') = -Log( (((x-x'))^2)+(((y-y'))^2))  + Log( (((x+x'))^2)+(((y-y'))^2)) 

Notice this obeys the Dirichlet condition (Subscript[G, D])(x=0 y,x'y')=0

(c) Use the result for (b) to find an expression for the induced charge density on the conductor.


In the text it is pointed out that Eq. (2.18) for the normal derivative of the Dirichlet Greens function
is the same as the surface charge density derived in Eq. (2.5). Is this also true for this case ?

We know that

According to the above discussion

ϕ(x,y)= (1/(4π ϵ))∫(Subscript[G, D])(x,y,x',y')ρ(x',y')

Now we know that  σ = ϵ (OverVector[E]) · (OverHat[n]) = - (∂_x) (1/(4π ))∫(Subscript[G, D])(x,y,x',y')ρ(x',y')= ((-1)/(4π ))∫(∂_x)(Subscript[G, D])(x,y,x',y')ρ(x',y') evaluated at x=0 (the surface)

Now ρ(x,y) is the continuous charge distribution but since we are interested in ϕ(x,y) in the region r > R our result
above shows that we can replace ρ(x',y') with that of a wire located a distance d from the conductor. i.e.  ρ(x',y')=λ δ(x'-d)δ(y')
and so we can greatly simplify our integral.  We get

 σ =  ((-1)/(4π ))∫(∂_x)(Subscript[G, D])(x,y,x',y')λ δ(x'-d)δ(y')=  ((-λ)/(4π )) (∂_x)(Subscript[G, D])(x,y,d,0) =  ((-λ)/(4π )) (∂_x) ( -Log( (((x-d))^2)+(y^2))  + Log( (((x+d))^2)+(y^2)) )
 =   ((-λ)/(4π ))( ((2 (d-x))/((((d-x))^2)+(y^2)))+((2 (d+x))/((((d+x))^2)+(y^2))) ) evaluated at x=0 ⟹   ((-λ)/(   π ))(( d)/((d^2)+(y^2)))


(d) Find the force in which the conductor attracts the cylinder.  According to Jackson an element of area dA on the surface has
a force dF = (σ^2)/2ϵ  dA  , (Can you verify this result?  Where did the extra factor of 2 come from? ) directed along the positive x axis.
Thus

dF= (((λ/ π ))^2) ((((( d)/((d^2)+(y^2)))))^2) dy dL 1/(2ϵ)  Integrating for all y and a strip of length L along the z-axis we get

F=  (λ^2)L/ (4π ϵ d)   or the force/length ⟹  (λ^2)/ (4π ϵ d)     Which is the same as the force of two lines of charge separated by a distance of 2d. Surprise!

(3) Consider the following periodic function

f(x) == S(Cos(x))  where S(z) is the so-called FresnelS function 
defined as 

S(z) = (∫_0^z)sin(π(t^2)/2)dt.

Mathematica includes this function which is called by the command  

FresnelS[z].

(a) What is the period of f(x) ?

(b) Plot f(x) in an interval of  one period.

Plot[FresnelS[Cos[z]],{z,0,8Pi}]

(c) Write a short numerical Mathematica script that represents f(x) in terms of a Fourier expansion
within the interval given in (b). Plot both f(x) and the Fourier representation of f(x) on the same graph.
How many terms did you use in the expansion in order to obtain a reasonable representation of f(x)?

a=2Pi;
f1[x_,m_]:=Sqrt[2/a] Sin[(2 Pi m)/a x];
f2[x_,m_]:=Sqrt[2/a] Cos[(2 Pi m)/a x];

f[x_]=FresnelS[Cos[x]];

acoff[n_]:=(Sqrt[2/a]) NIntegrate[f[x]*f2[x,n],{x,-a/2,a/2},AccuracyGoal->10^(-4)];
bcoff[n_]:=(Sqrt[2/a]) NIntegrate[f[x]*f1[x,n],{x,-a/2,a/2},AccuracyGoal->10^(-4)] ;

nmax=5;
Table[Chop[bcoff[n]],{n,0,nmax+1}]
atable=Table[Chop[acoff[n]],{n,0,nmax+1}]

fourier[x_]:= 1/2atable[[1]]+(Sqrt[a/2 ])Sum[atable[[n+1]] f2[x,n],{n,1,nmax}]

Plot[{fourier[x],f[x]},{x,-a/2,a/2}]

(d) In your plot you might have noted the Gibbs phenomenon in which violent oscillations in your representation
persist. Why is this not really a problem? (Hint: consider the simple function g(x)=x+(x^2). Numerically integrate 
∫g(x) f(x) dx , in the region defined by (b). Now compare your result with that obtained using the Fourier representation
of f(x). )

No Gibbs phenomena here.

(4) Repeat problem (3) using the Legendre polynomials instead of the Fourier functions. Which method seems to be more efficient for this particular problem ?   (Wait for next Hw assignment)