Electric Circuit Dynamics

This example demonstrates MLAB's differential equation-solving facilities, the use of MLAB's Fourier transform operations, and MLAB graphics in the context of the classic problem of analyzing an LRC circuit. A circuit containing a coil with an inductance of Lhenries, a resistor with a resistance of R ohms, and a capacitor with a capacitance of Cfarads in series is traditionally called an LRC circuit. Consider the LRC circuit shown below which also contains a voltage source component which exhibits a voltage drop of -f(t) at time t across the voltage source component. (This picture was constructed using MLAB.)

When the switch is closed at time 0, a current will begin to flow. Let I(t) be the current flowing in the circuit at time t, measured in amperes. Let Q(t) be the charge on the capacitor at time t, measured in coulombs. We have dQ(t)/dt = I(t).

The voltage drop across the resistor at time t is given by Ohm's law as R . I(t). The voltage drop across the capacitor at time t is Q(t)/C. The voltage drop across the coil at time t is L(dI(t)/dt), and b Kirchhoff's first law, the sum of the voltage drops across each of the circuit components is 0. Thus, we have

L{dI(t) / (dt)} + {Q / C} + RI(t) - f(t) = 0, or
L{d2I(t) / dt2} + R{dI(t) / (dt)}+{I(t) / C}-{df(t) / (dt)} = 0.

Take the initial conditions to be I(0) = 0 and dI(0)/dt = 0, and define f(t) = exp(-(t mod 1)). Fix L = 1, C = 1, and R = .5. Now we may solve the differential equation defining the current flow function I(t) and produce a graph of this function as shown below.

*function i''t(t)=(f't(t)-r*i't -i/c)/l
*initial i(0) = 0
*initial i't(0) = 0
*function f(t)=exp(-mod(t,1))
*l=1; c=1; r=.5
*m = integrate(i't,i''t,0:25!200)
*type odestr
   odestr = t   i't     i't't   i       i't
*draw m col (1,4)
*view

The value of the string variable ODESTR tells us what functions are numerically tabulated in the successive columns of M. The graph of the rate-of-change of the current function I't(t) can also be plotted from the data in M.

*delete w
*draw m col 1:2
*view

We can display the phase diagram graph for I(t) by plotting I't(t) vs. I(t) as follows

*delete w
*draw m col (4,2)
*view

We may ``zoom-in'' to see the neighborhood of the limit cycle more clearly as follows.

*WINDOW -.75 TO -.5, -.1 TO .1
*VIEW

We can use the particular MLAB Fourier transform operator realdft to compute the amplitude and phase-shift spectra of the current function I(t) tabulated in M COL (1,4).

*d=realdft(m col(1,4))
*delete w
*draw d col 1:2
*frame 0 to 1, 0 to .5
*top title "Amplitude Spectrum" size .2 inches
*w1=w

*draw d col(1,3)
*frame 0 to 1, .5 to 1
*top title "Phase-Shift Spectrum" size .2 inches
*view

Ignoring the large DC value at frequency zero, we see that the maximum amplitude occurs at about .2Hertz; this is the resonant frequency of the circuit. To see the amplitude spectrum at higher ``resolution'', we should subtract the DC-value from our signal I(t) and compute the amplitude spectrum of this shifted signal whose mean value is now zero.

The Fourier transform of I(t) contains the information to construct I(t) as a periodic function via its Fourier series. If the Fourier series is truncated, the resulting sum is a filtered form of I(t) omitting the high-frequency components corresponding to the truncated terms. Below we show a graph of the Fourier series of I(t) truncated to 7 terms. Note that Gibbs' phenomenon is exhibited, showing non-uniform convergence to the mid-point of the discontinuities occuring at the points between successive periods of~I(t).

*fct s(t)=sum(i,1,n, d(i,2)*cos(2*pi*d(i,1)*t + d(i,3)) )
*n=7
*q=points(s,0:25!120)
*delete w,w1
*draw m col 1:2 lt dotted
*draw q
*view

We can also use the MLAB Fourier transform operator to compute the amplitude and phase-shift spectra of the current rate-of-change function I't(t) tabulated in M COL (1,2). As above, this amplitude spectrum shows the resonant frequency of the circuit to be about .2Hertz. The peak at 1 Hertz corresponds to the frequency of oscillation of the forcing function f(t).

*d=realdft(m col 1:2)
*delete w,w1
*draw d col 1:2
*frame 0 to 1, 0 to .5
*top title "Amplitude Spectrum" size .2 inches
*w1=w

*draw d col(1,3)
*frame 0 to 1, .5 to 1
*top title "Phase-Shift Spectrum" size .2 inches
*view

Just as before, the Fourier transform of I't(t) contains the information to construct I't(t) via its Fourier series. If the Fourier series is truncated, the resulting sum is a filtered form of I't(t) omitting the high-frequency components corresponding to the truncated terms. Below we show a graph of the Fourier series of I't(t) truncated to 7 terms, superimposed on a graph of I't(t) plotted as a dotted line.

*fct s(t)=sum(i,1,n, d(i,2)*cos(2*pi*d(i,1)*t + d(i,3)) )
*n=7
*q=points(s,0:25!120)
*delete w,w1
*draw m col 1:2 lt dotted
*draw q
*view