Fancy3.mws
vector stuff
op
electric field plot
envelope detector
driven oscillator
| E&M Packet Problem ABW solution Potential and field plots of three point charges Learning objectives: Translation of standard 1/r potential to off-origin points Physical intuition regarding potential Identify equilibrium points Relate potential plot and electric field vectors MAPLE: plot3d, grad, fieldplot commands | ![]() |
restart; with(linalg): with(plots):
arbitrary units are used throughout in Coulomb's law
Q1:=1.;Q2:=1;Q3:=-1: # define relative magnitude of three charges
Set plot range; use non-even values to avoid singular points in grids
xmax:=0.755; ymax:=0.755: #plot domain
-xmax<x<xmax and-ymax<y<ymax
coul1:= Q1/((x-0.5)^2+y^2)^0.5; # potential due to 1st charge only
coul2:=Q2/((x+.5)^2+y^2)^0.5; # potential due to 2nd charge only
coul3:=Q3/((x)^2+(y-.7)^2)^0.5; # potential due to 3rd charge only
coul123:=coul1+coul2+coul3; # potential due to all three charges
3d countour plot of potential function
plot3d(coul123, x=-xmax..xmax, y=-ymax/2..2*ymax, grid=[49,49],style=patchcontour,contours=10,tickmarks=[1,1,100],view=-20..20);
Get electric field from negative gradient of potential
(Note: "-grad(U)" doesn't seem to work, so using "grad(-U)"
vect:=grad(-coul123, vector([x,y]));
Create new vector field of unit vectors E/|E|
vect2:=evalm(vect/sqrt(vect[1]^2+vect[2]^2));
fieldplot([vect2[1],vect2[2]],x=-xmax..xmax,y=-ymax..ymax+1/2,title=`unit
vectors showing direct'n of E field `);
Compress dynamic range of E-vector lengths by taking ln(1+|E|)
log is taken twice to put all lengths within range of fieldplot command
vect3:=evalm(vect2*ln(ln(sqrt(vect[1]^2+vect[2]^2)+1)+1));
fieldplot([vect3[1],vect3[2]],x=-xmax..xmax,y=-ymax..ymax+1/2,title=`electric
field vectors`);
Plot the phase and amplitude response of a driven, damped oscillator as a function of frequency.
This problem illustrates dealing with complex numbers, and the use of the 'op' command to parse an output expression.
It is interesting that very heavy damping (b=10, k=10, m=1) crushes the response away from zero frequency.
restart;with(plots):
assume(k,real);assume(m,real);assume(A,real);
assume(omega>0);assume(t,real);assume(b,real);assume(phi,real);
eqn:=A*exp(I*(omega*t+phi))-k*x(t)-b*diff(x(t),t) = m*diff(x(t),t,t);
s:=dsolve(eqn,x(t)); # (takes a little while)
x:=rhs(s);
values:={A=1,m=1,k=10,b=1/10};
xs:=subs(values,x);
steady:=op(1,xs); # keep only the oscillating terms for steady
state response
ev:=evalc(subs(t=0,phi=0,steady));
Phase of the response as a function of frequency. The conventional definition is the phase of the driver minus the phase of the driven system. The relative amplitude response will be taken as exp(-I alph).
assume(alph,real);
cosal:=op(1,ev); # cos(alph)
sinal:=-op(2,ev)/I; # sin(alph)
alph:=arctan(sinal/cosal);
plot(alph,omega=0..10,title=`response phase vs. omega`);

The driver is in phase with the driven system at very low frequencies, and pulls ahead until it is Pi/2 ahead at resonance, and goes on to be Pi ahead at very high frequencies.
Amplitude response as a function of frequency
xcc:=subs(I=-I,steady); # get complex conjugate
xmag:=simplify(steady*xcc); # get magnitude squared
xev:=simplify(evalc(xmag));
plot(sqrt(xev),omega=0..10,title=`Response amplitude vs freq.`);

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.
End of Fancy3.mws