Homework  Assignment III , Phy711  Due Oct 23  (based on book readings Chap 3.1-3.6 )  

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

(1) Repeat problem (3), of the previous assignment, using the Legendre polynomials instead of the Fourier functions. Which method seems to be more efficient for this particular problem ?

Plot[FresnelS[Cos[z]],{z,-1,1}]

(c) Write a short numerical Mathematica script that represents f(x) in terms of a Legendre expansion
within the interval given in (b). Plot both f(x) and the 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)?

f1[x_, m_] := LegendreP[m, x];
f[x_] = FresnelS[Cos[x]]; 
(* Use Eq. (3.24) in text for expansion coefficients *)
acoff[n_] := (2 n+1)/2 NIntegrate[f[x]*f1[x, n], {x, -1, 1}, 
     AccuracyGoal -> 10^(-4)]; 
(* try fit for different values of nmax, find that for n=5 get excellent (visual ) fit *)
nmax=5;
Table[Chop[acoff[n]],{n,0,nmax}];
atable=Table[Chop[acoff[n]],{n,0,nmax}]
fit[x_]:= Sum[atable[[n]] f1[x,n-1],{n,1,nmax}]
Plot[{fit[x],f[x]},{x,-1,1}]

(2) Legendre polynomials

(a)  Plot, in the interval  -1  < x < 1, the first 5 Legendre Polynomials.

Table[Plot[LegendreP[m,x],{x,-1,1}],{m,0,4}]

(b) Plot the product of the functions (Subscript[P, 2])(x) (Subscript[P, 6])(x) and  (Subscript[P, 6])(x) (Subscript[P, 6])(x),  then integrate those products over the interval given above.  Do your answers agree with Eq. (3.21) in the text?

Plot[LegendreP[2,x]LegendreP[6,x],{x,-1,1}]
Integrate[LegendreP[2,x]LegendreP[6,x],{x,-1,1}]
Plot[LegendreP[6,x]LegendreP[6,x],{x,-1,1}]
Integrate[LegendreP[6,x]LegendreP[6,x],{x,-1,1}]

So the normalization condition Eq. (3.21) holds for these two choices.

(c) Consider the three Polynomials (Subscript[P, 5]), (Subscript[P, 6]), (Subscript[P, 7]) show explicitly that Eqs. (3.29) are satisfied for this set.

(L+1)LegendreP[L+1,x]-(2L+1) x LegendreP[L,x]+ L LegendreP[L-1,x]
% /. L->6
Simplify[%]

D[LegendreP[L+1,x],x]- x D[LegendreP[L,x],x]- (L+1) LegendreP[L,x]
% /. L->6
Simplify[%]

( x^2-1) D[LegendreP[L,x],x]- L x LegendreP[L,x]+ L LegendreP[L-1,x] 
% /. L->6
Simplify[%]

So identities 3.29 are satisfied.

(3) Write a short Mathematica routine that uses, the Rodrigues formula,  Eq. (3.49)  to obtain expressions for the associated Legendre polynomials. Use that routine to generate those
Polynomials for the case l=9. Are these polynomials, for a fixed l, but different m orthogonal?

Rodfor[L_,m_]:=(-1)^m (1-x^2)^(m/2) D[LegendreP[L,x],{x,m}]

Table[Rodfor[9,m],{m,0,9}];
MatrixForm[%]
(* You can check these by using Mathematica which us an explicit expression using the LegendreP[L,m,x] command *)

(* Let's integrate the L=9,m=0 with L=9, m=9 Polynomial across range *)

Integrate[LegendreP[9,9,x]LegendreP[9,0,x],{x,-1,1}]

(* that looks OK *) 
(* how about L=9 m=1 with L=9 m=9 ? *)
Integrate[LegendreP[9,9,x]LegendreP[9,1,x],{x,-1,1}]

(* so in general for a given L the different m polynomials are not orthogonal, how about if L differs? *)

(4) Consider the problem discussed in section 3.3 in the text (a spherical conducting shell having both hemispheres at equal and opposite potential V). Repeat this problem (a correct numerical solution is acceptable) for the case where 

V(θ)  =  + V    0 <= θ <= π/3
V(θ)  =   - V     π/3 <= θ <= π

(a) Consider the special case where V=100 Volts,  a=10 cm.  Plot the potential contours on the xz plane. Include the regions  r > a and r <a.

(b) Give a stream plot of the gradient of the potential ϕ(r).  Find the value of the electric field at the center of the sphere.

(* for this problem all we need to do is use Eq. (3.33) in the text, but first we need to evaluate the coefficients (Subscript[A, n]) (Subscript[B, n]) *)

In the region r<a we get, according to Eq. (3.34) and (3.35)

V(θ)= V ∑ A(l) (Subscript[P, l])(Cosθ)

(here I have factored out the constant V, for the function V(x), so that the integrals are dimensionless.)

A(n) = (2 n+1)/2 (∫_(-1)^1)f(x) (Subscript[P, n])(x)dx

where

f(x) = + 1   if    Cos(π/3)  < x < 1
f(x) =  -1    if     -1 < x < Cos(π/3)

We can then use the routine for problem 1 above.

f[x_]:=Piecewise[{{1,1/2≤ x<1},{-1,-1 < x < 1/2}}]

Plot[f[x],{x,-1,1}]

f1[x_, m_] := LegendreP[m, x];
acoff[n_] := (2 n+1)/2   NIntegrate[f[x]f1[x,n],{x,-1,1},  AccuracyGoal -> 10^(-6)] ; 
(* try fit for different values of nmax, find that for n=5 get excellent (visual ) fit *)
nmax=22;
Table[Chop[acoff[n]],{n,0,nmax}];
atable=Table[Chop[acoff[n]],{n,0,nmax}]
fit[x_]:= Sum[atable[[n]] f1[x,n-1],{n,1,nmax}]
Plot[{fit[x],f[x]},{x,-1,1}]

(* Though these 22 terms give a reasonable fit to V(x) Mathematica's numerical Integration starts to tank
for Legendre Polynomials of higher order, so let's try finding an analytic expression for the integrals *)


FullSimplify[-Integrate[LegendreP[m,x],{x,-1,1/2}]+ Integrate[LegendreP[m,x],{x,1/2,1}]]

So it looks like it is better if we use this analytic expression for the exapnsion coefficients

Clear[acoff]

acoff[m_]:=Which[m==0,-1/2,m>0,-(2 m+1)/2 (((Sqrt[3]) π LegendreP[m,1,(1/2)]+2 Sin[m π])/(m π+(m^2) π))]

fit[x_,max_]:= Sum[acoff[n] f1[x,n],{n,0,max}]

Plot[{fit[x,100],f[x]},{x,-1,1}]

Note the Gibbs oscillations at θ=π/3

Now that we have the fit we can now express our solution (in the region inside R < a) as. Using length scale in meters a->0.10

a=0.1;
Vin[y_,z_]:=Sum[acoff[n](((y/a)^2+(z/a)^2)^(1/n))   LegendreP[n,z/Sqrt[z^2+y^2]],{n,0,52}]

ContourPlot[Vin[y,z],{y,-0.15,0.15},{z,-0.15,0.15},Contours->{1,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.0,-0.1,-0.2,-0.3,-0.4,-0.5,-0.6,-0.7,-0.8,-0.9,-1.0},RegionFunction->Function[{y,z},Sqrt[z^2+y^2]<a],ContourShading->Off]

Plot3D[Vin[y,z],{y,-0.2,0.2},{z,-0.2,0.2},RegionFunction->Function[{y,z},Sqrt[z^2+y^2]<a]]
(* note the units here are in terms of 100 V since we factored that out in th expression above *)

(* take gradient of *)
eterm[y_,z_,n_]=FullSimplify[Grad[((Sqrt[(y/a)^2+(z/a)^2]))^n LegendreP[n,z/Sqrt[z^2+y^2]] ,{y,z}]]

Ein[y_,z_]:=-Sum[acoff[n]eterm[y,z,n],{n,0,20}]

ingraph=StreamPlot[Ein[y,z],{y,-0.15,0.15},{z,-0.15,0.15},RegionFunction->Function[{y,z},Sqrt[z^2+y^2]<a],StreamPoints->20]

We repeat for the outside region then  (Subscript[A, n])=0 and according to Eq. 3.33 we require

V(θ) = ∑ (Subscript[B, L])/(a^(L+1)) (Subscript[P, L])(Cosθ)  with

We get the same expansion as above except that the factor  (((r/a))^L)  in our expansion gets replaced with (((a/r))^(L+1))

Vout[y_,z_]:=Sum[acoff[n](a/(Sqrt[y^2+z^2]))^(n+1)  LegendreP[n,z/Sqrt[z^2+y^2]],{n,0,52}]

ContourPlot[Vout[y,z],{y,-0.5,0.5},{z,-0.5,0.5},Contours->10,RegionFunction->Function[{y,z},Sqrt[z^2+y^2]>a],ContourShading->Off]

Plot3D[Vout[y,z],{y,-0.5,0.5},{z,-0.5,0.5},RegionFunction->Function[{y,z},Sqrt[z^2+y^2]>a],PlotRange->{-1,1}]
(* note the units here are in terms of 100 V since we factored that out in th expression above *)

(* take gradient of *)
etermout[y_,z_,n_]=FullSimplify[Grad[((Sqrt[(y/a)^2+(z/a)^2]))^(-n-1) LegendreP[n,z/Sqrt[z^2+y^2]],{y,z}]]

Eout[y_,z_]:=-Sum[acoff[n]etermout[y,z,n],{n,0,20}]

outgraph=StreamPlot[Ein[y,z],{y,-0.5,0.5},{z,-0.5,0.5},RegionFunction->Function[{y,z},Sqrt[z^2+y^2]>a],StreamPoints->30]

V[y_,z_]:=Piecewise[{{Vin[y,z],(Sqrt[y^2+z^2])<a},{Vout[y,z],(Sqrt[y^2+z^2])>a}}]

Plot3D[V[y,z],{y,-0.2,0.2},{z,-0.2,0.2}]

Show[outgraph,ingraph]

(c) To find the electric field at the origin we can take just take the leading order terms in our Legendre expansion, since z->0,y->0

we get

a[0]  (Subscript[P, 0])(Cos θ) + a[1] (r/a) (Subscript[P, 1])(Cosθ) + ... 

Now we know that  z=r Cos[θ]= r (Subscript[P, 1])(Cosθ)  therefore the leading order in the potential is 

Φ  ≈ a[1]/a z   and so  (OverVector[E]) = -a[1]/a (OverHat[k])   from our expansion above we know that a[1]= V

acoff[1]

and so

 (OverVector[E]) = -9/8 V/a (OverHat[k ])  = -1125 Volts (OverHat[k ])