Electricity and Magnetism (comments to: moloney@nextwork.rose-hulman.edu)
MAPLE
text and input lines
PHYSICS RESOURCE PACKETS PROJECT File: CIRCLOOP.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
A circular loop of wire of radius R is centered at the origin in the x-y plane and carries a current I counterclockwise from the x-axis to the y-axis. Perform an integration using the Biot-Savart law to determine the magnetic field components (Bx,By,Bz) at a point (0,0,H) due to half of the loop --- that half of the loop where x is greater than zero.
Solution
restart; with(linalg):
P:=vector([0,0,H]); # Vector to point in 3d space.
s:=vector([R*cos(theta),R*sin(theta),0]); # Vector to point on wire loop.
ds:=map(diff,s,theta); # Direction around wire loop.
r:=evalm(P-s); # Vector from loop to 3d point.
dB:=mu[0]*i*crossprod(ds,r)/(4*Pi*sqrt(dotprod(r,r))^3); # Biot-Savart
law
dBvect:=evalm(dB); # Distribute constant over the three components of the
vector.
solsymb:=map(int,dBvect,theta=-Pi/2..Pi/2); # Integrate around half loop
to find B.
Download circloop.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: CSPHERE.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
| The electric field between 2 concentric spheres is E=Q/(4pi eo r^2). where Q is the charge on the inner sphere of radius A. Plot the capacitance of these spheres (let the radius of the smaller sphere be 0.25 m) as a function of the radius of the outer sphere. In the limit as this radius goes to infinity this might be called the capacitance of a single sphere. (This is like clapping with one hand behind your back.) Can you do the same limiting calculation for flat plates or for cylinders? [No, only works for concentric spheres] |
![]() |
Solution
restart; with(plots):
Efield:=q/(4*Pi*epsilon[0]*r^2); # Equation for electric field between
the two spheres
The capacitance of these two spheres is C=q/dV, where dV is the change in potential of the two spheres. Find dV by integrating E.
dV:=int(Efield,r=R1..R2);
Cap:=simplify(q/dV);
Plot the capacitance as a function of the radius of the outer sphere. In order to do this, Eo and R1 must be defined.
epsilon[0]:=8.854*10^(-12);
R1:=.25;
plot(Cap*10^9,R2=0.25..3,C=0..1/2,numpoints=200,title=`Capacitance (nF)
as a function of R2`);
Download csphere.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: DIELCAP.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
The electric field inside concentric conducting spheres of radius R1 and R2 (R2>R1) is reduced by a factor of 2 when a particular dielectric material fills the entire space between R1 and R2. This means that the potential difference between the spheres is also reduced by a factor of 2, and the capacitance is increased by a factor of 2. (The charge Q on the spheres is kept the same.) When the dielectric doesn't fill the entire space between the spheres, the capacitance is increased by less than a factor of 2.
a) How thick would a dielectric spherical shell of thickness d, extending from R1 to R1+d have to be in order to increase the capacitance by a factor 1.5?
b) Would this thickness be the same if it extended from R2-d to R2?
Part A: Determine kappa
Determine k by using the relationship the E=Eo/kappa.
restart; eq1:=E=Eo/kappa;
eq2:=subs(E=Eo/2,eq1); # Since the field is reduced by a factor
of 2, the this equation is true:
kappa=solve(eq2,kappa); # Now solve for kappa.
Part B: Determine the dielectric's thickness
Efield:=q/(4*Pi*epsilon[0]*r^2); # Electric field between the spheres:
V is the difference in potential between the two conducting spheres. V can be expressed in terms of a constant, Q, R1 and R2 by the following formula:
V:=int(Efield,r=R1..R2);
C:=simplify(q/V); # The capacitance without the dielectric
When a dielectric is added, the potential difference becomes:
Delta*'V'=Delta*'Vo'/kappa;
Find the potential between the two spheres again, but this time divide the potential by kappa for r=R1..R1+d, and then add on the regular potential from r=R1+d..R2.
Vdie:=int(Efield,r=R1..R1+d)/kappa+int(Efield,r=R1+d..R2);
Cdie:=simplify(q/Vdie); # The capacitance with the dielectric
Now solve for the thickness of the dielectric (d) that causes the capacitance to increase by a factor of 1.5.
sol:=d=simplify(solve(Cdie=3/2*C,d));
Substitute in the value of kappa obtained in Part A.
sub:=subs(kappa=2,sol);
Perform a sanity check. If R1 was 1 and R2 was 2, then d should be between 0 and 1.
subs({R1=1,R2=2},sub);
Part C
Now find out what the result would be if the dielectric extended from R2-d to R2.
Use the regular potential for r=R1..R2-d, and then divide by kappa for r=R2-d..R2.
Vdie2:=int(Efield,r=R1..R2-d)+int(Efield,r=R2-d..R2)/kappa;
Cdie2:=simplify(q/Vdie2); # The capacitance with the dielectric on the outer shell
Now solve for the thickness of the dielectric (d) that causes the capacitance to increase by a factor of 1.5.
sol2:=d=simplify(solve(Cdie2=3/2*C,d));
Substitute in the value of kappa obtained in Part A.
sub2:=subs(kappa=2,sol2);
Perform a sanity check. If R1 was 1 and R2 was 2, then d should be between 0 and 1.
subs({R1=1,R2=2},sub2);
(The thickness of the dielectric is seen to be dependent on where it is placed between the two spheres. )
Download dielcap.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: EBFIELDS.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
The file 'ebfields' (ebfields.ms or ebfields.ma) contains a general solution of the motion of a particle in constant electric and magnetic fields. For any (Ex, Ey, Ez), (Bx, By, Bz), one can plot the motion of a particle of mass M, charge q and velocity (vx, vy, vz). From this general case, it specializes to the motion of a particle in 'crossed' E and B fields such that the particle travels in a straight line. The particular choices are q =1.61 x 10^-19 C , M = 1.67 x 10^-27 kg , Bz = 0.50 T , and vx = 8.4 x10^6 m/s.
a) Before opening one of these files, figure out what values are needed for (Ex,Ey,Ez) so the particle will travel in a straight line.
b) Now imagine reversing the particle's velocity, and predict what sort of motion it will undergo. (Straight line, motion in a plane, shape of motion etc.) Now make the changes and compare the results to your predictions.
Going Further: c) Predict what will happen in the crossed fields when a particle is released from rest. Try it out and compare results to your predictions.
d) Given the height H bet_ween the top and bottom limits of the y-motion, what can you say about the velocity of the particle at the top and bottom y-limits ofthe motion? [v=0 at bottom where it is released. At top the kinetic energy is mvx^2 which is equal to the work done on the particle by Ey, W = qEy H.]
Execute all the cells in this notebook until you get to the problem input section.
restart;
with(linalg):
with(plots):
Tell maple that all of the components are real to ease the calculations.
assume(Ex,real); assume(Ey,real); assume(Ez,real); assume(Bx,real); assume(By,real); assume(Bz,real); assume(vx0,real); assume(vy0,real); assume(vz0,real); assume(t,real); assume(q,real); assume(m,real);
Define vectors for general solution
Efield:=vector([Ex,Ey,Ez]);
Bfield:=vector([Bx,By,Bz]);
v:=vector([vx(t),vy(t),vz(t)]);
v0:=vector([vx0,vy0,vz0]);
Define Equations
Biot-Savart law
F:=evalm(q*Efield+q*crossprod(v,Bfield))=evalm(m*map(diff,v,t));
Initial acceleration of the particle
a0:=evalm((q*Efield+q*crossprod(v0,Bfield))/m);
Split up force vector into its three components.
Fx:=lhs(F)[1]=rhs(F)[1];
Fy:=lhs(F)[2]=rhs(F)[2];
Fz:=lhs(F)[3]=rhs(F)[3];
Solve the three equations for the vx(t), vy(t), and vz(t), and tell maple the initial conditions.
Also simplify the output. This will take a moment.
sol:=simplify(dsolve({Fx,Fy,Fz, vx(0)=v0[1], vy(0)=v0[2], vz(0)=v0[3], D(vx)(0)=a0[1], D(vy)(0)=a0[2], D(vz)(0)=a0[3]}, [vx(t),vy(t),vz(t)]));
Now collect like terms.
This will put the solution in the form: A*t + B*sin(theta) + C*cos(theta) + D
sol2:=map(collect,sol,{sin(q*sqrt(%1)*t/m),cos(q*sqrt(%1)*t/m),t});
General solution to the problem:
************** INPUT SECTION FOR PROBLEM****************
Now substitute in values the problem. Values may be changed as desired. Values in this section can be changed over and over again without executing the top part of the notebook again. Vectors are in x,y,z component form. Execute all cells between this and the next input section.
EfieldVect:=[0,0,0];
BfieldVect:=[0,0,2/10];
VelocityVect:=[0,-1,0];
charge:=1.6*10^(-19);
mass:=1.67*10^(-27);
************************************************************
Substitute in above values into the general solution.
sub:=subs({Ex=EfieldVect[1], Ey=EfieldVect[2], Ez=EfieldVect[3], Bx=BfieldVect[1], By=BfieldVect[2], Bz=BfieldVect[3], vx0=VelocityVect[1], vy0=VelocityVect[2], vz0=VelocityVect[3], q=charge, m=mass} ,sol2);
Change to floating point format.
velocity:=map(evalf,sub);
Integrate the xyz velocity components to find the position.
position:=subs(velocity,[int(vx(t),t),int(vy(t),t),int(vz(t),t)]);
Evaluate the integrals
positionxyz:=map(value,position);
Location on graph where particle is at t=0.
positionxyz0:=map(evalf,subs(t=0,positionxyz));
************* INPUT SECTION FOR GRAPHICAL OUTPUT *************
Values in this section may be repeatedly changed without executing the upper cells again. The three scalar variables control how large the E, B, and velocity vectors will be drawn on the graph. If too large a value is used, the vector will dominate the plot and the path of the particle will be too small to see. If too small a value is used, the vector will be too small to tell its direction.
scalarE:=1/5;
scalarB:=1/5;
scalarV:=1/10;
The path of the particle is plotted from t=0 to the time entered for tmax. If the output on the graph looks like "star" patterns, the value of tmax is probably too large. If the output looks like a straight or nearly straight line, tmax is probably too small. Execute the remaining cells in this notebook to get the solution.
*******************************************************************
tmax:=2*10^(-7);
Generates the curve for the path of the particle and the field & velocity vectors.
Path:=spacecurve({positionxyz, evalm(EfieldVect*scalarE*t), evalm(BfieldVect*scalarB*t), evalm(positionxyz0+VelocityVect*scalarV*t)}, t=0..tmax, numpoints=100):
Locate the position for the end of each vector.
Eend:=evalm(EfieldVect*scalarE*tmax);
Bend:=evalm(BfieldVect*scalarB*tmax);
Vend:=evalm(positionxyz0+VelocityVect*scalarV*tmax);
Draws letters identifying each vector.
Labels:=textplot3d({[Eend[1],Eend[2],Eend[3],`E`], [Bend[1],Bend[2],Bend[3],`B`], [Vend[1],Vend[2],Vend[3],`v`]}, color=red):
Display everything on the same graph.
display({Path,Labels},title=`Motion of particle through E and B fields`);
Download ebfields.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: EDIPOLE.MWS Rose-Hulman Institute of Technology Author: mjm 2/26/96 http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
Along the axis of a finite dipole, the electric field points one way in the region between the two charges, and it points the other way outside the charges. Which field produces the greater effect? The one in between is stronger, but it doesn't last very long.
To answer this question, we set up a line of dipoles, and integrate along the line of the dipoles, although we do not run right over the top of any individual charge.
Set up a line of 4 dipoles in the y-direction and integrate from x=-a to +a and find out the average electric field Will it be + or - ?
A line of + charges at (0,1/8), (0,3/8), (0, -1/8), (0, -3/8) and a line of - charges at x= +1/4. This gives 4 dipoles pointing in -x direction, f rom y=-3/8 to y=+ 3/8 m. Distance is 1/4 in x and y directions.
Predict the sign of Ex before beginning. [It should be mainly negative outside +/- charges and positive between the charges.] Graph Ex from -2 to 2 and see Ex is negative (mostly) outside dipoles,and large positive between them. Who wins? Integrate Ex from -2 to 2 and see. {Easy way to understand result is that V is positive for x<0 since all + chargesare closer, and V is negative for x> 1/4 because all negative charges are closer !}
restart; with(plots): k:=9e9: q:=1e-9:
Ex:=-k*q*(a-x)/((a-x)^2 + b^2)^(3/2); # Ex at point (x,0) from a
unit charge located at (a,b)
Edipole:=subs(a=0,Ex)-subs(a=1/4,Ex); # Ex due to a dipole with a spacing of 1/4
TotalDipoleEx:=0: for k from -3 to 3 by 2 do
TotalDipoleEx:=TotalDipoleEx+subs(b=k/8,Edipole): od:
plot(TotalDipoleEx,x=-1/2..1/2,title=`E_x of a finite dipole`);
evalf(int(TotalDipoleEx,x=-20..20));
Integral of Ex is positive; stronger positive part in middle wins; Dipoles like these would be created (or aligned)by an electric field pointing in the -x direction, so the field of the dipoles opposes the field that created it. This means in capacitors when dipoles are present in the dielectric, the electric field is smaller in magnitude than without the dielectric. Smaller electric field means smaller potential difference, which (for the same charges on the plates) a larger capacitance.
Yes, but what if we make the dipole distances tiny? (substitute a=1/100 instead of 1/4)
DipoleEx:=0:
for k from -3 to 3 by 2 do
DipoleEx:=DipoleEx+subs(a=0,b=k/8,Ex)-subs(a=1/100,b=k/8,Ex):
od:
plot(DipoleEx,x=-1/2..1/2);
evalf(int(DipoleEx,x=-20..20));
Nope, smaller, but still positive. Dipole field always opposes applied field.
Final note: - integral of Ex dx = V(2) - V(1). Since {integral Ex dx} is positive, V(1) > V(2). The potential at the start is higher than the potential at the end. This is because all + charges are closer to the start point than the - charges.
Download edipole.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: ENVELOPE.MWS Rose-Hulman Institute of Technology Author: mjm http://www.rose-hulman.edu/~moloney Software: Maple V Release 4


In AM (amplitude-modulated) radio, the amplitude of a high-frequency signal is modulated by the lower-frequency audio information. AM radios use 'envelope detectors' involving RC circuits in order to remove the high-frequency (the 'carrier wave') and leave the audio signal.
This problem gives you a simple modulated signal. It also gives you the 'rectified' signal we would see after passing through a diode (which lets current flow in one direction only.).
You are to pass the signal through a series RC circuit, first through a resistor R and then a capacitor C to ground. The 'envelope' signal is to be the voltage across the capacitor.
Letting v be the incoming signal, the equation to be solved is v -IR -q/C = 0 .
You must replace the current I by a function of q, the charge on the capacitor. This (call it eq1) will be the generic equation to be solved for q the charge on the capacitor as a function of time. By substituting values of R and C into eq1, you get eq2. This new equation is to be solved numerically for q (the charge on the capacitor) as a function of time.
a) Get the rectified signal first, by using the Heaviside function. The input signal is
v := sin(96 Pi t) *sin(12000 Pi t).
b) Set up the differential equation, using q(t) for the charge on the capacitor C . C is connected to ground and to the resistor R . R has the input voltage on one side and C on the other. [The current will be dq/dt]
c) Substitute values into the DE to make a new DE which will be solved numerically. Start off by setting RC equal to 1/10 the period of the 'carrier frequency'.
d) Solve the new DE numerically. Plot (use odeplot) and see what this gives you for an envelope function Compare this to the plot of the modulated function.
restart; with(plots):
th:=6000*2*Pi*t; # carrier frequency is 6000 hz
z:=sin(th*8/1000)*sin(th); # amplitude-modulated (AM) signal (modulated
at 48 hz)
v:=z*Heaviside(z); # rectified signal after going through a diode
plot(v,t=0..0.015,numpoints=400, title=`plot of modulated signal`);
From here on down, the students should be doing the work.
eq1:=v-R*diff(q(t),t)-q(t)/C=0;
eq2:=subs(R=45,C=1/100000,eq1);
RC is 0.45 ms. carrier pd= 0.17 ms. modulation pd .= 20 ms.
sol:=dsolve({eq2,q(0)=0}, q(t),numeric);
odeplot(sol,[t,q(t)*100000/1],0..15/1000,numpoints=400,title=`signal across capacitor`);
Note that we don't get all of the signal back, nor is the shape a perfect fit to the envelope.
Download envelope.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: FILTRC.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
| R-C circuits are often used as simple frequency filters in electronic circuits. To see how this works, apply an input voltage of v(t) = A cos (wt) to a series resistor-capacitor (the input voltage is applied at the resistor R and the output voltage is taken at the point where the resistor R and capacitor C join together. One side of the capacitor is connected to the resistor and other side is grounded.). V is the voltage across the capacitor, and q is the charge magnitude on either plate, so q = VC, from the definition of capacitance. The time derivative of this relation is : dq/dt = current = I = C dV/dt. a) Show that the loop equation for voltage can be written v(t) = RC dV/dt + V. b) Use dsolve to obtain a solution to the second of these equations for V(t). c) Assume R= 10 k ohms, C = 0.12 microfarads, A = 0.350 V, and w=6328 rad/s. Plot the voltage across the capacitor vs. time. If the voltage amplitude across the capacitor is close to A, then the circuit is 'passing' most of the input voltage through to the output at this frequency. d) Plot the ratio of the amplitude of the voltage across the capacitor to the input voltage amplitude A from frequencies from 10 rad/sec to 2000 rad/sec. e) Based on this plot, would you say this filter is a 'high-pass' or a 'low-pass' filter? | ![]() |
restart; with(plots):
assume(C>0); assume(R>0); # Tell Maple that C and R are positive.
Part A
i:=C*diff(Vout(t),t);
# Current through the capacitor at the output
eq1:=Vr(t)=i*R; #
Voltage acoss resistor
eq2:=subs(Vr(t)=Vin(t)-Vout(t),eq1); # Write
Vr(t) in terms of Vin(t) and Vout(t)
DE:=subs(Vin(t)=A*cos(omega*t),eq2); #
Substitute the input voltages given in the problem for Vin(t)
Part B
sol:=dsolve({DE,Vout(0)=0},Vout(t)); #
Solve the differential equations for Vout(t) using initial conditions.
Vout_t:=rhs(sol); #
Set Vout_t equal to the right side of the equation.
Part C
Vout_tsub:=subs({R=10000,C=0.12*10^(-6),A=0.35,omega=6238},Vout_t); # Substitute in values for the constants. plot(Vout_tsub,t=0..0.005,title=`Vout_tsub vs time`);
Part D
Vout_tsub2:=subs({R=10000,C=12/100*10^(-6),A=35/100},Vout_t);
Vss:=simplify(subs(exp(-2500/3*t)=0,Vout_tsub2));
# Find the steady state response by eliminating the exponential terms.
Max:=maximize(Vss,t,0..1);
Aout:=Max[2]; # The amplitude of Vss is the
term in Max with the radical in the denominator.
plot(Aout/0.35,omega=0..2000,title=`Gain of the circuit vs omega`);
Part F
The circuit is acting like a low pass filter. Now solve the problem in the S domain.
assume(omega>0); # Tell Maple that omega is
positive.
VoutS:=(35/100)/(C*I*omega)/(R+1/(C*I*omega)); # Use voltage division to
find Vout
AoutS:=abs(VoutS); # Find the amplitude of the output.
Find the amplitude of the output using the values for R and C given in the problem.
AoutSsub:=simplify(subs({R=10000,C=12/100*10^(-6)},AoutS));
Aout-AoutSsub; # Compare this answer to the solution from the differential
equation.
Solutions are the same.
Download filtrc.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: HELIX.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
| Use the basic solution in 'ebfields' to find the conditions which will give a particle trajectory which is a helix. | ![]() |
Download ebfields.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: HOLTPLOT.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4


Make a 3D plot of the magnitude of the magnetic field between a pair of helmholtz coils. (See problem 'loopaxis' for the setup of helmholtz coils.) [In parts a) and b) below, try to use the contours of the 3D plot to help you answer the questions.] a) How far from the center of the coils do we have to go laterally before the field magnitude changes by 10%. Express your answer in terms of the coil radius ( 0.25 radii, 0.32 coil radii, or whatever). b) How far from the center do we have to go along the axis before the magnitude of the field changes by 10%?
restart; with(linalg): with(student): with(plots):
P:=vector([x,y,z]); # Vector to point in 3d space.
s:=vector([R*cos(theta),R*sin(theta),0]); # Vector to point on wire loop.
ds:=map(diff,s,theta); # Direction around wire loop.
r:=evalm(P-s); # Vector from loop to 3d point.
dB:=mu[0]*i*crossprod(ds,r)/(4*Pi*sqrt(dotprod(r,r))^3); # Biot-Savart
law
dBvect:=evalm(dB); # Distribute constant over the three components of the
vector. dBvectsub:=subs({R=1,mu[0]=4*Pi*10^(-7),i=1},op(dBvect)); # Substitute
in the values of the known constants.
Evaluate the integrals using Simpson's method. ":" supresses the output.
sol:=map(simpson,dBvectsub,theta=0..2*Pi,20):
Bxyz:=map(value,sol): # Evaluate the Sigma notation in sol.
Bxyz2:=subs(z=z-1,op(Bxyz)): # Position a second coil 1m from the first
coil.
Bxyztot:=evalm(Bxyz+Bxyz2): # Add the fields of the two coils together.
Bmag:=(Bxyztot[1]^2+Bxyztot[2]^2+Bxyztot[3]^2)^(1/2): # Find the magnitude
of the combined field.
Bmagxz:=subs(y=0,Bmag): # Take a cross section
of the B field in the xz plane.
Plot the B field in the xz plane.
plot3d(Bmagxz,x=-1.5..1.5,z=-0.5..1.5,view=0..2*10^(-6),title=`|B| from helmholtz coils`,grid=[40,40]);
Animate the same function as in the 3d plot to "walk" through the shape above. Notice how little the field in the middle changes.
animate(Bmagxz,x=-1.5..1.5,z=0..1,title=`|B| from helmholtz coils`,color=red);
Find the field at the center of the two coils.
field_at_center:=evalf(subs({z=0.5,x=0,y=0},Bmag));
field90:=0.9*field_at_center; # 90% of the field
field110:=1.1*field_at_center; # 110% of the field
Plot the area in the xz plane that has a magnetic field within 10% of the center.
plot3d(Bmagxz,x=-1.5..1.5, z=-0.5..1.5,view=field90..field110,grid=[40,40],title=`B field within 10% of that at the center`);
Bmagz:=subs({y=0,x=0},Bmag): # Magnitude of field
along z axis.
Bmagx:=subs({y=0,z=1/2},Bmag): # Magnitude of field between the two coils.
plot({Bmagx,field90,field110},x=-1.5..1.5,B=0..1.2*10^(-6),title=`|B| along
x when z=0.5`);
Part A
Find lateral distance for field to change by 10% Since R=1, answer is in radii.
fsolve(Bmagx=field90,x=0.6..0.8);
plot({Bmagz,field90,field110},z=-.5..1.5,B=0..1.2*10^(-6),title=`|B| along
z axis`);
Part B
Find distance along z axis for field to change by 10%
sol:=fsolve(Bmagz=field90,z=0.9..1.2);
Subtract 1/2 from this answer to measure the distance from the center. Again, the answer is in radii
dist:=sol-0.5;
Download holtplot.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: HOWCLOSE.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
A rectangular current loop measures 1.50 m by 0.90 m, and carries a current of 10.5 A. At the mid-point of the long (1.50 m) leg inside the loop, how far would we have to be from the wire so that the magnetic field there would be within 2% of the field we would calculate assuming wire to be infinitely long?
restart;
B:=mu[0]*i*(cos(theta)+cos(phi))/(4*Pi*d); # Formula for magnetic
field of a finite length wire.
B field at distance y from the top horizontal 1.5m segment:
Bth:=subs({d=y,cos(theta)=(3/4)/sqrt((3/4)^2+y^2),cos(phi)=(3/4)/sqrt((3/4)^2+y^2)},B);
Blh:=simplify(subs(y=9/10-y,Bth)); # B field for the lower horizontal
1.5m segment:
B field for the left vertical 0.9m segment:
Blv:=(subs({d=3/4,cos(theta)=y/sqrt(y^2+(3/4)^2),cos(phi)=(9/10-y)/sqrt((9/10-y)^2+(3/4)^2)},B));
B field for the right vertical 0.9m segment:
Brv:=Blv;
Bloop:=Bth+Blh+Blv+Brv; # Total B field
For a long wire, the magnetic field at a distance y from the wire is:
Blong:=mu[0]*i/(2*Pi*y);
Find y where the magnetic field from the loop is within 2% of the magnetic field from a "long" wire at the same distance y.
ratio:=expand(Bloop/Blong)=102/100;
Solve the above equation for the distance. (about 1.1 cm)
fsolve(ratio,y);
Download howclose.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: INDUCTOR.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
| A "real" inductor has inductance, L, and resistance, R. Suppose
this component is connected in parallel with an ideal capacitor, C as shown.
(a) Find the magnitude of the effective impedance of the combination as a function of the frequency. (b) Find the phase relationship between the voltage across the combination and the current through it. |
![]() |
| (c) Suppose the combination is connected in series with a pure resistance,
R as shown. A voltage given by Vin=V0 cos(wt) is applied to the entire
circuit. Find Vout(t), the voltage across the parallel combination.
(d) Let V0=5, L=3 mH,C=10 uF,RL=3 ohms, R=1000 ohms. Plot the magnitude of the frequency response of the circuit for f = 0 Hz to 2000 Hz. |
![]() |
Solution:
restart;
assume(omega>0,C>0,L>0,RL>0,V0>0,f>0,t>0);
Part A
capacitor impediance
ZC:=1/(I*omega*C);
inductor impediance
ZL:=I*omega*L;
resistor impediance
ZR:=RL;
impediance of the combination of elements
Zeff:=1/(1/(ZL+ZR)+1/ZC);
Zeff:=simplify(Zeff);
magnitue of the impediance
magZeff:=simplify(abs(Zeff));
Part B
Convert the complex impediance into polar form.
The first term is the magnitude, and the second term is the angle in radians.
sol:=convert(Zeff,polar);
Evaluate the phase angle of the impediance.
evalc(convert(sol,list)[2]);
Simplify the result. Since Z=V/I, "angle" = (angle of V phasor) - (angle of I phasor)
angle:=simplify(");
Part C
Input voltage in phasor form
Vin:=V0;
Use voltage division to find the output voltage
Vout:=Vin*Zeff/(Zeff+R);
Vout:=simplify(Vout);
Find the phase angle of the output voltage.
evalc(convert(sol,list)[2]);
Simplify the result.
angleV:=simplify(");
Find Vout as a function of time.
'Vout(t)'=magVout*cos(omega*t+angleV);
Part D
Find the magnitude of the frequency responce.
magVout:=evalc(abs(Vout));
Vout1:=simplify(subs({V0=5,omega=2*Pi*f,L=3*10^(-3),C=10*10^(-6),RL=3,R=1000},magVout));
plot(Vout1,f=0..2000,title=`Frequency Response of Circuit (Vout vs f)`);
Download inductor.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: INKDROP.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
An ink drop of radius 0.028 mm (density = 1000 kg/m^3) is travelling at a speed of 90 m/s. It passes through deflecting plates which are 8 mm long and separated by 1.5 mm. We assume the print head slews horizontally across the page. Since the drops must be able to cover the whole field of a single character in the vertical direction, they must be able to be deflected at least 0.5 mm up or down. The voltage difference across the deflecting plates is 2200 V.
a) Find the electric field E between the deflecting plates in V/m.
b) Determine the charge Q which must be on a drop in order to deflect it by 0.50 mm by the time it has reached the edge of the deflecting plates. [Sanity-check your answer by noting that the smallest charge which can be applied is 1.6 x 10^-19C ].
c) The width of the deflecting plate is 6 mm. Determine the total charge on the deflecting plate in order create the potential difference of 2200 V and the electric field E. [Assume the charge on each deflecting plate is uniformly distributed over the plate. Note: the total charge on the plate shoud be much larger than the charge on the drop going through the plates if the drop is not to change the field.]
restart;
Part A: The E field
Efield:=V/d; # The E field is related to the voltage difference and distance between plates Efield2:=subs({V=2200,d=0.0015},Efield); # Now substitute in the given values for V and d.
Part B: Particle's Charge
Calculate the particle's mass. Density equals mass over volume, so mass equals density times volume.
mass:=density*volume;
Now substitue in the given values, assuming the inkdrop to be a sphere.
m:=evalf(subs({density=1000,volume=4/3*Pi*(28*10^(-6))^3},mass));
Now calculate the particle's acceleration. The following kinematic formula can be used:
KE:=x=a*t^2/2+t*v[0]+x[0];
xo and yo are both 0 because the ink drop has no initial motion in the vertical direction. x is the distance the particle is to be deflected (0.0005 m), and t is the time it takes for the drop to travel through the plates (0.008m/(90m/s)). Substituting in these values gives:
KEsub:=subs({v[0]=0,x[0]=0,x=0.0005,t=0.008/90},KE);
acc:=solve(KEsub,a); # Solve the above equation for the acceleration of the ink drop.
Setting the two force equations (F=q*E and F=m*a) equal to each other gives:
force:=charge*Efield2=m*acc;
q:=solve(force,charge); # Now solve for the charge of the ink drop.
Part C: Total charge on the plate
By using Gauss's Law, the total charge on the plate can be calculated.
GaussLaw:=E*A=charge/epsilon[0];
Now substitute in known values.
GaussLaw2:=subs({E=Efield2,A=0.008*0.006,epsilon[0]=8.854*10^(-12)},GaussLaw);
charge=solve(GaussLaw2,charge); # Solve for the total charge
Note that the charge on the plate is larger than the charge on the particle by a factor of about 100.
Download inkdrop.mws [Use 'File Save']
PHYSICS RESOURCE PACKETS PROJECT File: LINEINT.MWS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
In a region where the elecrtric field is (i 3.4 x10^3 + j 2.13) N/C, determine the potential difference between the points P1=(1.0,1.5) m and P2=(0.5, 2.5) m:
a) by taking the dot product of the electric field vector, and the vector from P1 to P2.
b) by integrating along the vector from P1 to P2.
c) by first integrating in the x direction and then in the y direction.
restart; with(linalg):
Part A
Efield:=vector([3400,2.13]); # Enter in the electric
field as a vector.
P1:=vector([1,1.5]); > P2:=vector([0.5,2.5]); # Enter in the
positons of the points as vectors.
s:=evalm(P2-P1); # Vector from P1 to P2.
mag_s:=sqrt(dotprod(s,s)); # Find the magnitude of s.
Since the electric field is a constant, the potential difference is simply -E•s
U:=-dotprod(Efield,s);
Part B
The solution can also be found by finding ds and then integrating.
ds:=evalm(s/mag_s); # Unit vector in the direction
of s.
U:=-int(dotprod(Efield,ds),s=0..mag_s);
Part C
Another way to solve the problem is to first integrate in the x direction, and then in the y.
Define fields in the x and y directions.
Efieldx:=Efield[1]; Efieldy:=Efield[2];
Ux:=-int(Efieldx,x=1..0.5); # Integrate in the x direction from x=1.0 to
0.5
Uy:=-int(Efieldy,y=1.5..2.5); # Integrate in the y direction from y=1.5
to 2.5
The potential difference between P1 and P2 is the sum of the change in potentials along the two paths.
U:=Ux+Uy;
Download lineint.mws [Use 'File Save']


PHYSICS RESOURCE PACKETS PROJECT File: LOOPS.MWS Rose-Hulman Institute of Technology Author: mjm http://www.rose-hulman.edu/~moloney Software: Maple V Release 4
Axial magnetic fields from one, two,and five loops
a) Use the starter material and evaluate the formula at the center of the loop.
b) Make a function f1, where you substitute I=1 amp and mu = 4 pi x 10^(-7). Put the coil center at z=2 give it a radius of 1 m , and plot f1 from z=0 to z=4.
c) At the coil center the slope of Bz is zero and again at z= infinity. Locate the value of z where the slope magnitude is a maximum. [Hint: what does this have to do with the second derivative of the function?]
d) Put a second coil at z=3.2 m and plot Bz for the pair of coils.
e) Evaluate the 1st and 3rd derivatives of the B field halfway between the two coils
f) Use a loop to create a 'loose' solenoid of 5 coils. Plot axial magnetic field (should have some wiggles in it.)
g) Estimate the area under the graph of B along the axis by using a rectangle.
restart; with(plots):
f = Bz @(0,0,z) from loop centered @ z=c, radius R. (factor of mu_o I/4pi left off) This is our basic building block.
f:=mu/(4*Pi)*I*2*Pi*R^2/(R^2 + (z-c)^2)^(3/2);
subs(z=0,c=0,f); # check it at center of the loop; want mu I / 2R
f1:=subs(mu=4*Pi*1e-7,I=1,f);Plot this to
see what it looks like (put at z=2 and plot 0 to 4)
plot(subs(R=1,c=2,f1),z=0..4);
twoloops:=subs(R=1, c=0,f1) + subs(R=1,c=12/10,f1); # two loops, spaced 1.2 units apart
d2:=diff(f,z,z); # 2nd derivative of f:
d3:=diff(f,z,z,z); # 3rd derivative of f:
solution gives two answers for where Bz is zero for a single loop, half a radius away from loop center to either side
sol:=solve(d2=0,z);
d1:=diff(twoloops,z); # first derivative:
evalf(subs(z=6/10,d1)); # check it at z =
0.6 (midpoint)
evalf(subs(z=6/10,d3)); # check 3rd deriv
at z = 0.6
plot(twoloops,z=0..1); # plot the axial field of 2 loops, one at zero, one at 1.2.
twoa:=subs(R=1, c=1/4,f1) + subs(R=1,c=3/4,f1);
plot({1e7*twoloops,1e7*twoa},z=0..1); # Example of plotting two 'functions' together Example of how to use a software loop to make a group of coils
g:=0;
for k from 1 to 4 do
g:=g+subs(R=1/8,c=k/6,f1); od:
plot(1e6*g,z=-1..2,title=` B(axis) of 'loose' solenoid`);
This plot is intentionally wiggly; not a very good solenoid. An estimation exercise would be to guesstimate the area under the curve by means of a a rectangle By eyeball, try 6.1 height from 0.0 to 0.85
Here's a rectangle with a height of 6.2 and going from 0.01 to 0.85 . Then we'll plot it with the g function [it needs capitals on Heaviside!! (I found out the hard way)]
rect:= 62/10*(Heaviside(z-1/100)-Heaviside(z-85/100));
plot({rect,1e6*g},z=-1..2,title=` Area estimation via rectangle`);
evalf(int(1e6*g,z=-5..5)); # Lets integrate and see how good this
is (in the neighborhood)
Download loops.mws [Use 'File Save']
Go to the next Electricity and Magetism Maple Page