Lecture note · Fig. 9 · control systems
Driven oscillations: how a system answers a push it does not control
Numa Koudsie · for school students and first- and second-year undergraduates
01A push from outside
Push a child on a swing. Push at the wrong moments and the swing barely moves; push in time with it and, with gentle pushes, it soon flies high.
The swing has a rhythm of its own, and your hand is a second rhythm that the swing does not control. Everything interesting about driven systems comes from the meeting of these two rhythms. A washing machine shakes hardest for a few seconds while its drum speeds up through one particular rate of turning. A loudspeaker cone follows the music only because it was designed so that its own rhythm lies far from the notes it plays. London’s Millennium Bridge swayed on its opening day in 2000 because thousands of footsteps fed its sideways motion.
In the first two cases a system with its own natural frequency is driven by a force that depends on time and on nothing else: the driver does not care where the system is. That is the idealisation of this page. The bridge shows its limit: the walkers adjusted their steps to the sway, so the force depended on the very motion it was driving, a feedback loop that needs a larger model. In control engineering this force is the input, the motion is the output, and the system in between is the plant. Our plant is the simplest one with a rhythm of its own (Figure 1).
02The model
Three forces act along the floor. The spring pulls back in proportion to the stretch, \(-kx\) (Hooke’s law). The damper, a piston in oil, resists in proportion to the speed, \(-c\dot x\). The driver pushes with \(F(t)\). Newton’s second law adds them:
\[ m\ddot x + c\dot x + kx = F_0\cos\Omega t . \tag{1} \]It is tidier to describe the plant by two numbers with a direct meaning. The natural frequency \(\omega_0 = \sqrt{k/m}\) is the undamped natural angular frequency; with underdamping, free ringing instead occurs at \(\omega_d=\omega_0\sqrt{1-\zeta^2}\); the damping ratio \(\zeta = c/(2\sqrt{km})\) says how quickly the ringing dies, from \(\zeta = 0\) (forever) to \(\zeta = 1\) (no ringing at all). Dividing (1) by \(m\),
\[ \ddot x + 2\zeta\omega_0\,\dot x + \omega_0^2\,x = \frac{F_0}{m}\cos\Omega t . \tag{2} \]On the page \(m = 1\ \mathrm{kg}\) and the sliders set \(f_0 = \omega_0/2\pi\), \(\zeta\), \(F_0\) and the drive frequency \(f = \Omega/2\pi\). With the defaults, \(f_0 = 0.5\ \mathrm{Hz}\) and \(\zeta = 0.05\), so \(k = 9.87\ \mathrm{N/m}\) and \(c = 0.314\ \mathrm{kg/s}\). The block starts at rest in its rest position, \(x(0) = 0\), \(\dot x(0) = 0\).
Notice what sets this equation apart from the double pendulum of Fig. 2: the time \(t\) appears on its own, inside \(\cos\Omega t\). The future depends not only on where the block is, but also on the clock, because the clock belongs to someone else.
03The exact solution★
Equation (2) is linear, so its solutions add. The general solution is one particular solution that keeps pace with the driver plus any solution of the free equation (right-hand side zero), chosen to satisfy the starting conditions.
3.1The steady state
Guess that, once things settle, the block moves at the driver’s own frequency. Complex numbers make the guess painless: write the force as the real part of \(\tfrac{F_0}{m}e^{i\Omega t}\) and try \(x = \operatorname{Re}\,Z e^{i\Omega t}\). Each time derivative becomes a factor \(i\Omega\), and (2) turns into algebra:
\[ \left(\omega_0^2 - \Omega^2 + 2i\zeta\omega_0\Omega\right) Z = \frac{F_0}{m} \quad\Longrightarrow\quad Z = \frac{F_0/m}{\omega_0^2 - \Omega^2 + 2i\zeta\omega_0\Omega} = X e^{-i\varphi} . \tag{3} \]At the drive frequency, acceleration is \(-\Omega^2\) times displacement. Spring and inertia therefore contribute the opposing real terms \(k\) and \(-m\Omega^2\); the damper contributes \(ic\Omega\), a quarter cycle out of phase. The physical denominator \(k-m\Omega^2+ic\Omega\) keeps these contributions together. Dividing by \(m\) gives the form in (3).
The modulus and argument of \(Z\) give the amplitude and the phase lag of the steady motion \(x_\text{p}(t) = X\cos(\Omega t - \varphi)\):
\[ X = \frac{F_0/m}{\sqrt{\left(\omega_0^2-\Omega^2\right)^2 + \left(2\zeta\omega_0\Omega\right)^2}}, \qquad \varphi = \operatorname{atan2}\!\left(2\zeta\omega_0\Omega,\ \omega_0^2-\Omega^2\right) . \tag{4} \]3.2The transient
The free equation has solutions that ring at the damped frequency \(\omega_d = \omega_0\sqrt{1-\zeta^2}\) and die away with the time constant \(\tau = 1/(\zeta\omega_0)\). Adding one to the steady state and demanding \(x(0) = \dot x(0) = 0\) fixes its two constants:
\[ \begin{aligned} x(t) &= X\cos(\Omega t-\varphi) + e^{-t/\tau}\left(C_1\cos\omega_d t + C_2\sin\omega_d t\right), \\ C_1 &= -X\cos\varphi, \qquad C_2 = \frac{\zeta\omega_0\,C_1 - X\Omega\sin\varphi}{\omega_d} . \end{aligned} \tag{5} \]Read (5) as a story in two acts. At first the block rings at its own frequency, because a sudden push excites its own rhythm. That ringing dies within a few \(\tau\), and afterwards only the driver’s rhythm remains. With the defaults, \(\tau = 6.37\ \mathrm{s}\), \(X = 0.482\ \mathrm{m}\) and \(\varphi = 25.3^\circ\).
04Resonance and phase
Formula (4) contains the whole of Figure 2. Drive slowly, \(\Omega \ll \omega_0\), and the block simply follows the force, stretched by the static amount \(F_0/k\). Drive fast, \(\Omega \gg \omega_0\), and the block hardly moves: its inertia wins. In between, near \(\Omega = \omega_0\), the spring term and the inertia term in (3) cancel, only the small damping term is left to limit the motion, and the amplitude soars. This is resonance.
The height of the peak measures how weakly damped the plant is. At \(\Omega = \omega_0\) exactly, (4) gives
\[ X(\omega_0) = \frac{F_0}{2\zeta k} = Q\,\frac{F_0}{k}, \qquad Q = \frac{1}{2\zeta}, \tag{6} \]so the quality factor \(Q\) is the amplification at resonance. The default \(\zeta = 0.05\) gives \(Q = 10\): a force that would stretch the spring by \(0.101\ \mathrm{m}\) if applied slowly moves the block \(1.01\ \mathrm{m}\) at resonance. The true maximum sits slightly below \(\omega_0\), at \(\Omega_r = \omega_0\sqrt{1-2\zeta^2}\), which is \(0.4987\ \mathrm{Hz}\) for the defaults; for \(\zeta \ge 1/\sqrt2\) there is no peak at all.
The phase tells the swing story in numbers. Below resonance the block moves with the force (\(\varphi \approx 0\)). At resonance it lags by exactly a quarter period, \(\varphi = 90^\circ\): the force is largest when the block is fastest, which is precisely when a push does the most work. Far above resonance the block moves against the force (\(\varphi \to 180^\circ\)).
Set the drive frequency equal to the natural frequency, \(f = f_0 = 0.50\ \mathrm{Hz}\), and press Solve. Wait until the transient has died, then compare Measured, last 10 s with \(Q F_0/k\) from (6), and read Phase lag φ. Now raise the damping ratio to 0.2, repeat, and find the factor by which the amplitude fell.
05The control engineer’s view★
Engineers compress everything above into one function. Take the Laplace transform of (1) with the block starting at rest: each time derivative becomes a factor \(s\), and the output is the input multiplied by the transfer function \(G(s)\),
\[ X(s) = G(s)\,F(s), \qquad G(s) = \frac{1}{ms^2 + cs + k} . \tag{7} \]The zeros of the denominator, the poles \(s = -\zeta\omega_0 \pm i\omega_d\), are the plant’s own rhythm in disguise: their real part is the decay rate \(1/\tau\), their imaginary part the ringing frequency \(\omega_d\). Both lie in the left half of the complex plane, which is what makes the plant stable: left alone, it returns to rest.
Setting \(s = i\Omega\) gives the response to a steady sinusoidal input, and comparing with (3), \(G(i\Omega) = X e^{-i\varphi}/F_0\). Plotting \(|G(i\Omega)|\) and \(\varphi\) against the drive frequency on logarithmic axes gives the Bode plot, the second figure on this page. It is the plant’s identity card: from it an engineer reads how strongly each frequency gets through and how late it arrives. Designing a controller, a device that measures \(x\) and adds a correcting force, means reshaping this card. Fig. 10 does that for a wind turbine.
06The start: beats and decay
Figure 4 is the run the page stores with its default settings. The drive at \(0.45\ \mathrm{Hz}\) is close to the natural \(0.50\ \mathrm{Hz}\), so during the first seconds two nearby rhythms are present at once: the driver’s, and the block’s own ringing from equation (5). Like two guitar strings slightly out of tune, they drift in and out of step once every \(1/|f-f_d| = 20.25\ \mathrm{s}\), where \(f_d = \omega_d/2\pi = 0.499\ \mathrm{Hz}\) is the damped ringing frequency (with \(f_0\) in its place the estimate is 20 s, close enough at this weak damping). The amplitude swells and shrinks with that rhythm; the first swell reaches \(0.596\ \mathrm{m}\) at \(t = 8.0\ \mathrm{s}\). Because the ringing decays with \(\tau = 6.4\ \mathrm{s}\), the second beat is already weak, and after half a minute only the steady oscillation at \(\pm X = \pm 0.482\ \mathrm{m}\) is left.
Lower the damping ratio to 0.01 and press Solve. The beats now last for the whole minute. Count the time between two swells and compare with \(1/|f - f_0|\). Then set \(f = 0.20\ \mathrm{Hz}\), far from resonance, and use equation (5) to explain why the transient hardly shows.
07Where the energy goes
The driver does work at the rate \(P_\text{in} = F\dot x\); the damper turns energy into heat at the rate \(P_\text{loss} = c\dot x^2\). The energy stored in the oscillator, \(E = \tfrac12 m\dot x^2 + \tfrac12 kx^2\), changes by the difference:
\[ E(T) - E(0) = \underbrace{\int_0^T F\dot x\,\mathrm{d}t}_{W_\text{in}} - \underbrace{\int_0^T c\,\dot x^2\,\mathrm{d}t}_{W_\text{loss}} . \tag{8} \]In the steady state the stored energy only sloshes between spring and block, so on average every joule the driver puts in is lost in the damper. Averaging over one period, \(\langle P_\text{in}\rangle = \langle P_\text{loss}\rangle = \tfrac12 c\,\Omega^2 X^2\), which is \(0.292\ \mathrm{W}\) for the defaults. Over the full 60 s run the page finds \(W_\text{in} = 19.03\ \mathrm{J}\), \(W_\text{loss} = 17.93\ \mathrm{J}\) and \(E(T) - E(0) = 1.11\ \mathrm{J}\). Equation (8) then holds to \(1.3\times10^{-12}\) of \(W_\text{in}\): this is the Energy balance readout, a test of the whole computation.
08Solve first, then play
Unlike the earlier figures, this page does not integrate while you watch. Because the driver’s force is known for all time in advance, the whole future is fixed the moment you press Solve. The page uses that.
- Solve. The fourth-order Runge–Kutta method advances \((x, \dot x)\) from \(t = 0\) to \(60\ \mathrm{s}\) in 60 000 steps of \(h = 1\ \mathrm{ms}\), carrying \(W_\text{in}\) and \(W_\text{loss}\) along in the same steps. Measured in headless Chromium, the default solve took 9.5–12.4 ms; the time depends on the machine.
- Store. Every tenth state is written into the page’s dataset, together with \(F\), \(E\), \(P_\text{in}\) and \(P_\text{loss}\): 6 001 rows, one every \(10\ \mathrm{ms}\). The exact solution (5) is stored next to it only to measure the error. With the defaults the largest difference over the run is \(5\times10^{-12}\ \mathrm{m}\), the RK4 vs exact readout.
- Play. Each frame converts the playback time into the nearest stored row, \(n = \operatorname{round}(t_\text{play}/10\,\mathrm{ms})\), and draws that row as it is. Nothing is interpolated or recomputed, so the picture and the data panel always show the same numbers, and the Stored run slider can jump to any instant.
When you move a slider after solving, the page says so and keeps playing the old run until you solve again: the data on screen always belong to the settings named in the status line.
09On screen
The block’s rest position is at \(x_\mathrm{wall}+0.42W_0+b_\mathrm{block}/2\), where \(W_0\) is the width between the wall and the right margin. The largest displacement is drawn as \(0.11W_0\) either side of rest. Measured on 770 px and 356 px canvases, rest lies at 47.8 % and 50.0 % of the width, and the full displacement span is 20.1 % and 19.4 %. The burgundy arrow on the block has a length proportional to \(F/F_0\). The dial at the top right is the driver: its arm turns at \(\Omega\), and its horizontal shadow is proportional to the force. Below, the stored run is drawn over its full 60 s, dark up to the present sample and pale after it, with \(F(t)/k\) dashed in burgundy and the steady amplitude \(\pm X\) dashed in rose. The Bode plot marks the present drive frequency with a filled dot (equation (4)) and the amplitude measured from the last 10 s of the run with a ring around it: when the run has settled, the two coincide.
10What the model leaves out
- Real springs obey Hooke’s law only for small stretches; at large amplitude the restoring force stops being proportional to \(x\), and the resonance peak bends and can jump.
- Real friction is often dry rather than viscous: it does not grow with speed. The linear damper is the version that keeps the equation solvable.
- The driver is ideal: it delivers \(F_0\cos\Omega t\) whatever the block does. A real motor slows down when the load is heavy, and at resonance the load is heaviest.
- The page covers the underdamped range \(\zeta \le 0.9\) only, where the exact solution (5) applies; \(\zeta \ge 1\) needs a different form of the transient.
10.1Further reading
- J. R. Taylor, Classical Mechanics, University Science Books, 2005: chapter 5, oscillations, driven damped oscillators and resonance.
- A. P. French, Vibrations and Waves, W. W. Norton, 1971: chapter 4, forced vibrations and resonance.
- K. J. Åström and R. M. Murray, Feedback Systems: An Introduction for Scientists and Engineers, Princeton University Press, 2008: transfer functions, frequency response and Bode plots.