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

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\text{'}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\text{'}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\text{'}\text{'},\; 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\text{'}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\text{'}t(t)$ contains the information to construct $I\text{'}t(t)$ via its Fourier series. If the Fourier series is truncated, the resulting sum is a filtered form of $I\text{'}t(t)$ omitting the high-frequency components corresponding to the truncated terms. Below we show a graph of the Fourier series of $I\text{'}t(t)$ truncated to 7 terms, superimposed on a graph of $I\text{'}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