Electromagnetic Theory (comments to: moloney@nextwork.rose-hulman.edu)

This calculates the approximate surface charge on a long thin wire, where the wire radius is small compared to its overall length. The method is to use a set of hollow cylindrical sections, each of which has its own constant charge density. When the potential at the center of each section on the axis is specified (it will be the same for a conducting wire), one may solve N equations for N unknowns (the charge densities), and then plot the charge density as a function of distance along the wire.

It ultimately determines the charges on each subsection of the wire by solving N equations in N unknowns, by specifying that the potential on the axis of each hollow cylindrically- shaped element be the same V.

restart; with(plots):

Potential due to a charged section of a cylinder when the ends of the cylinder are located at Z1 and Z2

V:=Int(Int(k*rho*2*Pi*r/sqrt(r^2+z^2),r=0..R),z=Z1..Z2);
V:=value(V); # Evaluate the double integral.
Vsub:=subs({k=1,R=1},V); # Set the value of the k constant to 1, and the radius of the cylinder to 1.

L:=5; # The length of a section
Vs:=1; # The potential measured at each section >

Find the potential at different segment distances away from a charged segment.

d:=(dist)->evalf(subs({Z1=-L/2+L*dist,Z2=L/2+L*dist},Vsub));

Make a function for the potential of the section at "dist" sections away from the charged section with rho="charge".

Vsec:=(dist,charge)->subs(rho=charge,d(dist));

Set up 4 equations to find the 4 unknown charges on the 4 sections.

V1:=Vs=Vsec(0,rho1)+Vsec(1,rho2)+Vsec(2,rho3)+Vsec(3,rho4); V2:=Vs=Vsec(1,rho1)+Vsec(0,rho2)+Vsec(1,rho3)+Vsec(2,rho4); V3:=Vs=Vsec(2,rho1)+Vsec(1,rho2)+Vsec(0,rho3)+Vsec(1,rho4); V4:=Vs=Vsec(3,rho1)+Vsec(2,rho2)+Vsec(1,rho3)+Vsec(0,rho4);

fsolve({V.(1..4)},{rho.(1..4)}); # Solve four equations for four unknowns

Now do the same problem with more segments in the cylinder.

NumSeg:=40; # Number of segments in cylinder.

Use a loop to create the equations.

for m from 1 to NumSeg do
V.m:=Vs=convert([seq(Vsec(abs(m-n),rho.n),n=1..NumSeg)],`+`): od:

sol:=fsolve({V.(1..NumSeg)},{rho.(1..NumSeg)}); # Solve the 40 equations for the 40 unknowns.
assign(sol); # Assign the solutions to p1-p40
plot([seq([n,rho.n],n=1..NumSeg)],title=`Charge distribution across the bar`); # plot the results.

Download cylinder.ms [Use 'Save File']

mjm (in packets\Emag) rev 3

(Legendr2.ms)

Objective: To derive the first eight Legendre polynomials in two different ways:
1) Using Maple's built in Legendre function
2) By using the formula given: 1/(2^n*(n!))*(d/dx)^n[(x^2 - 1)^n]

Software: Maple V Release 3.0 Authors: Perry Peters and Greg Williby
Concepts Used: Legendre Polynomials Load the orthopoly package

with(orthopoly); Generate the first eight Legendre polynomials

p1 := P(1,x); p2 := P(2,x); p3 := P(3,x); p4 := P(4,x);
p5 := P(5,x); p6 := P(6,x); p7 := P(7,x); p8 := P(8,x);

By using a simple FOR loop, the formula given at the top

can be evaluated for several values of n.

for n from 1 to 8 do
P.n:= simplify(1/(2^n*(n!))*diff((x^2 - 1)^n, x$n)); od;

Plot the first eight Legendre polynomials

plot({P.(1..8)},x=-1.5..1.5,y=-1.5..1.5,title=`Plot of first eight Legendre Polynomials`,color=black);

Download legendr2.ms [Use 'Save File']

pulse.ms

restart;
with(plots):

The idea here is to create a moving magnetic field by successively energizing coils one after the other. The first coil current ramps up to a max, and then as it ramps down, the second coil's current is ramping up to a max. Then as the second coil current ramps down, the current in the 3rd coil ramps up, and so on.

For a loop of wire, the magnetic field a distance x along its axis is:

B:=mu[0]*i*R^2/(2*(x^2+R^2))^(3/2);

where x is the axial distance from the center of the coil, and R is the radius of the coil. Substituting in values for the constants gives:

Bsub:=subs({mu[0]=1,R=1},B);

Next we generate a triangular pulse to represent the current in an individual coil.

Pulse:=(t)->(1-abs(t))*Heaviside(1-abs(t));
plot(Pulse(t),t=-2..2);

Position five coils at locations x=-1, x=-1/2, x=0, x=1/2, and x=1. The coil at x = -1 is pulsed first, and successive coils are pulsed one second after the start of the pulse in the previous coil.

for n from -2 to 2 do
B.n:=subs({x=x-n/2,i=Pulse(t-n)},Bsub);
od;

Draw the individual coils' contributions to the axial B field in red.

f:=animate({B.(-2..2)},x=-2..2,t=-3..3,color=red,frames=50):

Draw the resultant axial B field in black.

g:=animate(convert({B.(-2..2)},`+`),x=-2..2,t=-3..3,color=black,frames=50):

Animate the individual contributions and resultant fields.

display({f,g});

Download pulse.ms [Use 'Save File']

File: QCAP.MS
Rose-Hulman Institute of Technology
Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney
Software: Maple V Release 4

A constant charge resides on the parallel plate capacitor shown at right. The dielectric material is free to move frictionlessly between the two plates. As the dielectric moves, the capacitance changes, the voltage(V=Q/C) changes, and the potential energy of the system changes as well. This potential energy change with capacitor position results in a force exerted on the dielectric. The nature of this force, and the resulting motion are studied in this exercise. The worksheet QCAP.MS does four things.

a) Find the potential energy of the system as a function of x, w, l, d, Q, k1 and k2 (dielectric constants), and the permittivity of free space (epsilon0).
b) Find the force on the dielectric as a function of x for w = 6 cm, l = 10 cm, d = 1 cm, Q=1 microcoulomb, k1 = 2, and k2 = 1.
c) Plot the resulting motion of the dielectric if x0 = -9 cm, and m = 0.100 kg.
d) Plot the resulting motion of the dielectric if x0 = -10.1 cm, and v0 = 1 cm/s.

Potential energy of the capacitor

> U:=1/(2*C)*Q^2;

Model the capacitor as two capacitors in parallel

> C:=C1+C2;

Find the capacitance of each capacitor

> C1:=k1*epsilon[0]*A1/d;
> C2:=k2*epsilon[0]*A2/d;

Express the area of the two capacitors in terms of l, w, and x.

> A:=l*w;
> A2:=min(w*abs(x),l*w);
> A1:=A-A2;

Potential energy of the system.

> U:=simplify(U);

Part B

> Fx:=-diff(U,x);

Part C Enter in values for the variables.

> w:= 0.06; l:= 0.1; d:= 0.01; Q:= 1*10^(-6); k1:= 2; k2:= 1;
> epsilon[0]:=8.8541878*10^(-12);
> plot(U,x=-0.3..0.3,title =`Potential Energy vs Position`);
> m:=0.1;
> plot(Fx,x=-0.2..0.2,title=`Force vs x`);

put the differential equation in rational form.

> DE:=convert(subs(x=x(t),Fx)=m*diff(x(t),t$2),rational);
> sol:=dsolve({DE,x(0)=-0.09,D(x)(0)=0},x(t),numeric);

A graph can be made using odeplot.

> with(plots):
> odeplot(sol,[t,x(t)],0..1,labels=[t,x],title=`Position vs time`);

Look at the form of the solution.

> sol(2);

A plot can also be made in the usual fashion by making a function that just returns the x(t) part of sol.

> xt:=t->rhs(sol(t)[2]);
> xt(2);
> plot({xt},0..1,title=`Position vs time`);

The output looks sinusoidal. Check this by plotting our solution along with a cosine function on the same plot.

> h1:=plot(xt,0..1,color=red):
> h2:=plot(0.09*cos(9.15*t-3.14),t=0..1,color=blue):
> display({h1,h2});



Part D

Let's see what happens if we start the dielectric moving with an initial velocity from a position completely outside the plates.

>sol2:=dsolve({DE,x(0)=-0.11,D(x)(0)=0.01},x(t),numeric);
> odeplot(sol2,[t,x(t)],0..2,labels=[t,x],title=`Position vs time`);

Download qcap.ms [Use 'Save File']

(slot.ms)

slot.ms This is virtually Example 3. p. 127 in Griffiths text.

A 'slot' is formed by two infinite grounded metal plates parallel to one another, and by a third metal plate. The one grounded metal plate is in the xz plane, from the origin to positive infinity in x and plus/minus infinity in z. The second plate is parallel to the first, located at y=Pi. The third plate runs between the first two from y=0 to y=Pi, and from z=- infinity to to +infinity.

We first plot the fourier sine series solution for this case, then the analytic solution.

restart;with(plots):
s:=1/k*exp(-k*x)*sin(k*y);
V:=4/Pi*sum(subs(k=2*k-1,s),k=1..20);
plot3d(V,x=0..4,y=0..Pi,grid=[40,40],title=`Graph of potential using 20 terms in the series`);
V:=2*1/Pi*arctan(sin(y)/sinh(x));
plot3d(V,x=0..4,y=0..Pi,grid=[40,40],title=`Graph of potential using an explicit function`);

Download slot.ms [Use 'Save File']

vcap.ms

PHYSICS RESOURCE PACKETS PROJECT
File: VCAP.MS
Rose-Hulman Institute of Technology
Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney
Software: Maple V Release 4



A constant voltage is applied across the parallel plate capacitor shown at right. The dielectric material is free to move frictionlessly between the two plates. As the dielectric moves, the capacitance changes. The voltage source maintains a constant plate voltage on the capacitor, either supplying or accepting charge as necessary. The potential energy of the system changes with the position of the dielectric, resulting in a force exerted on the dielectric. The nature of this force, and the resulting motion are studied in this exercise. The worksheet VCAP.MS does four things.


a) Find the potential energy of the system as a function of x, w, l, d, V, k1 and k2 (dielectric constants), and the permittivity of free space (epsilon0). Use 0 as a reference for the initial potential energy of the battery. The origin of our coordinate system is shown in the figure at right.
b) Find the force on the dielectric as a function of x for w = 6 cm, l = 10 cm, d = 1 cm, V = 10000 volts, k1 = 2, and k2 = 1.
c) Plot the resulting motion of the dielectric if x0 = -9 cm, and m = 0.100 kg.
d) Plot the resulting motion of the dielectric if x0 = -10.1 cm, and v0 = 1 mm/s.


Part a:

restart;

The total potential energy of the system includes the potential energy of the charged capacitor and the energy stored in the battery.

U:=Cap_U+Battery_U;

The potential energy of the capacitor is given by

Cap_U:=1/2*C*V^2;

The capacitor breaks down into two regions of interest, one with dielectric between the plates and the remainder without dielectric. We model these regions of the capacitor as two capacitors in parallel

C:=C1+C2;

Find the capacitance of each capacitor

C1:=k1*epsilon[0]*A1/d;
C2:=k2*epsilon[0]*A2/d;

Express the area of the two capacitors in terms of l, w, and x.

A:=l*w;
A2:=min(w*abs(x),l*w);
A1:=A-A2;

Expression for the potential energy of the battery

Battery_U:=deltaQ*V;

Charge flow

deltaQ:=Q0-Q(x);
Q0:=simplify(subs(x=0,C*V));
Q(x):=simplify(C*V);
deltaQ:=simplify(deltaQ);

Potential energy of the system.

U:=simplify(U);

Part b:

Fx:=-diff(U,x);

Part c:

Enter in values for the variables.

w:= 0.06; l:= 0.1; d:= 0.01; V:= 10000; k1:= 2; k2:= 1; epsilon[0]:=8.8541878*10^(-12); plot(Fx,x=-0.2..0.2,title=`Force vs x`);
m:=0.1;

Put the differential equation in rational form.

DE:=convert(subs(x=x(t),Fx)=m*diff(x(t),t$2),rational);
sol:=dsolve({DE,x(0)=-0.09,D(x)(0)=0},x(t),numeric);

A graph can be made using odeplot.

with(plots):
odeplot(sol,[t,x(t)],0..30,labels=[t,x],title=`Position vs time`);

Look at the form of the solution.

sol(2);

A plot can also be made in the usual fashion by making a function that just returns the x(t) part of sol.

xt:=t->rhs(sol(t)[2]);
xt(2);
plot({xt},0..30,title=`Position vs time`);



The output looks sinusoidal. Check this by plotting our solution along with a cosing function on the same plot.

> h1:=plot(xt,0..10,color=red):
> h2:=plot(0.09*cos(.603*t-3.14),t=0..10,color=blue):
> display({h1,h2});

Part d:

Let's see what happens if we start the dielectric moving with an initial velocity from a position completely outside the plates.

sol2:=dsolve({DE,x(0)=-0.11,D(x)(0)=0.001},x(t),numeric);
odeplot(sol2,[t,x(t)],0..30,labels=[t,x],title=`Position vs time`);

Download vcap.ms [Use 'Save File']

(wireskin.ms)

PHYSICS RESOURCE PACKETS PROJECT File: wireskin.MS Rose-Hulman Institute of Technology
Author: mjm (moloney@nextwork.rose-hulman.edu)             Software: Maple V Release 3



BACKGROUND.

Skin depth effects in a long conducting wire at high frequency are readily modelled using numerical integration. This deals directly with the physics without referring to the formidable properties of bessel functions whose argument is complex. The numerical integration scheme actually generates the kelvin functions (ber, bei) and their derivatives, so their features can be examined if desired.

A wire carrying oscillating current is subject to the well-known 'skin depth' effect of current decreasing from the outside of the wire toward the center at high frequencies. This classic problem involves the complication of a pure imaginary term in a bessel-type differetialequation. Texts treating the problem1,2,3 discuss it in terms of ber and bei ('kelvin') functions, a rather daunting approach for most undergraduates.

Numerical integration solves this problem without any reference to bessel functions. It provides one more illustration of the way computers are gradually altering the way physics problems can be approached by students. A cylindrical wire of radius a carrying a current along its axis (the z-axis) is assumed to have a current density J with only a z component. J is proportional to E, and Ez depends only on the radial coordinate r, in addition to its sinusoidal variation in time

[1]                   E(r,t) = E(r) e^(iwt) .

E(r) has real and imaginary parts

[2]       E(r) = ER(r) + i EI(r) .

To make clear that both parts of [2] contribute to the overall current density, we write the real part of [1], Re ( E(r,t) ) = ER(r) cos(wt) - EI(r) sin(wt) .The magnitude of current density in the wire will then be

                           | Jz | = s ( ER(r)^2 + EI(r)^2 )^1/2 , (s is the electrical conductivity)

and its phase will be arctan (EI/ER).

It is readily shown[1] that E(r) is subject to the following equation in a good conductor ( msw >> mew^2, m = magnetic permeability, e = dielectric constant) :

                            r E''(r) + E'(r) = i msw r E(r) .

Since msw has dimensions of inverse length squared, we may define L^2 = 1/ msw. A change of variables to x=r/L gives

[3]                       x E''(x) + E'(x) = i x E(x) .

Using the real and imaginary parts of Er, [3] becomes two coupled equations with real coefficients, ready for numerical integration

[4a]                       x ER''(x) + ER'(x) = - x EI(x) , and

[4b]                       x EI''(x) + EI'(x) = + x ER(x) .

A Runge-Kutta integration scheme for integrating [4a] and [4b] is discussed below, but may as easily (or more easily) be accomplished in symbolic packages like Maple or Mathematica. (The author's students carried out this calculation using symbolic algebra, and then readily plotted the magnitude of E(r) vs. x. ) As variables in a Runge-Kutta (RK) integration scheme[4] we may take

y[1] = ER(x)

dy[1]/dx = y[2]

y[2] = ER'(x)

dy[2]/dx = -y[3] -y[2]/x

y[3] = EI(x)

dy[2]/dx = y[4]

y[4] = EI'(x)

dy[2]/dx = +y[1] -y[4]/x

with initial conditions

y[1] = 1,

y[2] = y[3] = y[4] = 0,

dy[1]/dx = dy[2]/dx = dy[3]/dx = 0, and

dy[4]/dx = 1.

(Taking EI = 0 on the axis amounts to adjusting the time so that the phase of E(r) is zero on the axis.) The integration proceeds from x=0 on the axis, with a unit-valued pure real E(r), until the surface of the wire is reached at an x value depending on L and the wire radius R.

Table I was producedin a few seconds on a 25 Mhz 386 PC, using a 4th order RK scheme with a step size of 5x10^-4. It shows how the magnitude grows steadily from the axis to whatever value of x represents the radius of the wire.

Note that the phase of E(r) changes, along with its magnitude. ( While attempting to check the results in Table I, it was discovered that ER and EI correspond exactly to the ber and bei 'kelvin' functions listed in ref. 6, p. 430. The 'skin depth' d is conventionally[1-3] defined as d= (2)^1/2 L = ( 2/ msw)^1/2 . This is less natural in the differential equation than L, but is conveniently interpreted as the distance over which field amplitude falls by a factor of e.).

In a copper at a frequency of 100 kHz, L = 0.147 mm so a wire radius of 0.5 mm equals 3.41 L. (The wire radius then occurs at x=3.41.) At this frequency, Table I shows that the current density at the wire surface is about 2.6 times greater than along the wire axis.

At a much lower frequency of 1 kHz, however, L=1.47 mm. At this frequency, a 0.5 mm wire radius corresponds to x=0.341, and Table I shows that the current density is nearly constant across the cross-section of the wire.

Based on the kelvin function asymptotic forms (ref. 5, p. 381), exp[ x/(2)^1/2 ]/ (2 Pi x)^1/2 ,one can say that when the wire radius is a very large number of 'lengths' L , the current density approaches doubling (increasing by a factor of 2.03) for each additional length L. We see this toward the bottom of Table I, where an increase of 1.05 L nearly doubles the surface electric field.

Table I. Electric field |E(x)| at wire radius R = 0.50 mm

(  |E(r=0)| = 1.0   )

N. B. Conventional skin depth 'd' is sqrt(2) L.
x =3.5 means radius = 3.5 L, or radius = 3.5 skin depths as defined in this paper.
BUT, x=3.5 also means radius = 3.5/sqrt(2) d, or radius = 2.5 d = 2.5 conventional skin depths

[1] W. T. Scott, The Physics of Elecricity and Magnetism, 2nd Ed. (Wiley, New York, 1966), Sec. 10.8.

[2] M. A. Heald and J. B. Marion, Classical Electromagnetic Radiation, 3rd Ed. (Saunders College Publishing, Fort Worth, 1995) Sec 5.6.

[3] W. R. Smythe, Static and Dynamic Electricity, 3rd Ed. (McGraw-Hill, New York, 1968), Sec 10.04.

[4] W. Press et. al, Numerical Recipes In Pascal (Cambridge U. Press, Cambridge, 1989), Sec. 15.1. [5] M. Abramowitz, and I. Stegun, Handbook of Mathematical Functions (U. S. Government Printing Office, 1966).

---

SOLUTION:

restart;with(plots):

Need to plot a DE evaluated numerically. With real and imag parts, let y(x) = r(x) + I s(x),
then equate real and imag parts of the DE. This gives (where Dop stands for the differential operator)

Dop(r(x)) = -b s(x), and
Dop(s(x)) = +b r(x).

r(0)=1, r'(0)=0, s(0)=0 s'(0)=0

NOTE: To integrate for 'conventional' skin depths,
one must have a factor of two on the right-hand side of er and ei.

er:= x*diff(r(x),x$2)+diff(r(x),x)=-x*s(x);
ei:= x*diff(s(x),x$2)+diff(s(x),x)=+x*r(x);
ers:=er,ei;
fcns:= {r(x), s(x)};

sol2:=dsolve( {ers, r(0)=1, D(r)(0)=0, s(0)=0, D(s)(0)=0 },fcns, numeric);

Substitute in values from problem.

odeplot(sol2,[x, sqrt(r(x)^2+s(x)^2)],0..7,numpoints=100,title=`|E| vs wire radius in skin depths `);

Download wireskin.ms [Use 'Save File']

Back to the Electromagnetic Theory Resource Page