Skip to main content
Logo image

Introduction to Differential Equations

Section 2.6 Forced oscillations and resonance

Note: 2 lectures, §3.6 in [EP], §3.8 in [BD]
Let us return back to the example of a mass on a spring. We examine the case of forced oscillations, which we did not yet handle. That is, we consider the equation
\begin{equation*} mx'' + cx' + kx = F(t) , \end{equation*}
for some nonzero \(F(t)\text{.}\) The setup is again: \(m\) is mass, \(c\) is friction, \(k\) is the spring constant, and \(F(t)\) is an external force acting on the mass.
We are interested in periodic forcing, such as noncentered rotating parts, or perhaps loud sounds, or other sources of periodic force. Once we learn about Fourier series in Chapter 5, we will see that we cover all periodic functions by simply considering \(F(t) = F_0 \cos (\omega t)\) (or sine instead of cosine, the calculations are essentially the same).

Subsection 2.6.1 Undamped forced motion and resonance

First let us consider undamped (\(c=0\)) motion. We have the equation
\begin{equation*} mx'' + kx = F_0 \cos (\omega t) . \end{equation*}
This equation has the complementary solution (solution to the associated homogeneous equation)
\begin{equation*} x_c = C_1 \cos (\omega_0 t) + C_2 \sin (\omega_0 t) , \end{equation*}
where \(\omega_0 = \sqrt{\nicefrac{k}{m}}\) is the natural frequency (angular). It is the frequency at which the system “wants to oscillate” without external interference.
Suppose that \(\omega_0 \not= \omega\text{.}\) We try the solution \(x_p = A \cos (\omega t)\) and solve for \(A\text{.}\) We do not need a sine in our trial solution as after plugging in we only have cosines. If you include a sine, it is fine; you will find that its coefficient is zero (I could not find a second rhyme).
We solve using the method of undetermined coefficients. We find that
\begin{equation*} x_p = \frac{F_0}{m(\omega_0^2 - \omega^2)} \cos (\omega t) . \end{equation*}
We leave it as an exercise to do the algebra required.
The general solution is
\begin{equation*} \mybxbg{ ~~ x = C_1 \cos (\omega_0 t) + C_2 \sin (\omega_0 t) + \frac{F_0}{m(\omega_0^2 - \omega^2)} \cos (\omega t) . ~~ } \end{equation*}
Written another way
\begin{equation*} x = C \cos (\omega_0 t - \gamma) + \frac{F_0}{m(\omega_0^2 - \omega^2)} \cos (\omega t) . \end{equation*}
The solution is a superposition of two cosine waves at different frequencies.

Example 2.6.1.

Take
\begin{equation*} 0.5 x'' + 8 x = 10 \cos (\pi t), \qquad x(0)=0, \qquad x'(0)=0 . \end{equation*}
Let us compute. First we read off the parameters: \(\omega = \pi\text{,}\) \(\omega_0 = \sqrt{\nicefrac{8}{0.5}} = 4\text{,}\) \(F_0 = 10\text{,}\) \(m=0.5\text{.}\) The general solution is
\begin{equation*} x = C_1 \cos (4 t) + C_2 \sin (4 t) + \frac{20}{16 - \pi^2} \cos (\pi t) . \end{equation*}
Solve for \(C_1\) and \(C_2\) using the initial conditions: \(C_1 = \frac{-20}{16 - \pi^2}\) and \(C_2 = 0\text{.}\) Hence
\begin{equation*} x = \frac{20}{16 - \pi^2} \bigl( \cos (\pi t)- \cos (4 t) \bigr) . \end{equation*}
Figure 2.5. Graph of \(\frac{20}{16 - \pi^2} \bigl( \cos (\pi t)- \cos (4 t) \bigr)\text{.}\)
Do notice the “beating” behavior in Figure 2.5. To see it, use the trigonometric identity
\begin{equation*} 2\sin \left( \frac{A-B}{2} \right) \sin \left( \frac{A+B}{2} \right) = \cos B -\cos A \end{equation*}
to get
\begin{equation*} x = \frac{20}{16 - \pi^2} \left( 2 \sin \left(\frac{4-\pi}{2} t \right) \sin \left( \frac{4+\pi}{2} t \right) \right) . \end{equation*}
The function \(x\) is a high frequency wave modulated by a low frequency wave.
Now suppose \(\omega_0 = \omega\text{.}\) Obviously, we cannot try the solution \(A \cos (\omega t)\) and then use the method of undetermined coefficients. We notice that \(\cos (\omega t)\) solves the associated homogeneous equation. Therefore, we try \(x_p = A t \cos (\omega t) + B t \sin (\omega t)\text{.}\) This time we need the sine term, since the second derivative of \(t \cos (\omega t)\) contains sines. We write the equation
\begin{equation*} x'' + \omega^2 x = \frac{F_0}{m} \cos ( \omega t) . \end{equation*}
Plugging \(x_p\) into the left-hand side we get
\begin{equation*} 2 B \omega \cos (\omega t) - 2 A \omega \sin (\omega t) = \frac{F_0}{m} \cos (\omega t) . \end{equation*}
Hence \(A = 0\) and \(B = \frac{F_0}{2m\omega}\text{.}\) Our particular solution is \(\frac{F_0}{2m\omega} \, t \sin (\omega t)\) and our general solution is
\begin{equation*} x = C_1 \cos (\omega t) + C_2 \sin (\omega t) + \frac{F_0}{2m\omega} \, t \sin (\omega t) . \end{equation*}
The important term is the last one (the particular solution we found). This term grows without bound as \(t \to \infty\text{.}\) In fact it oscillates between \(\frac{F_0 t}{2m\omega}\) and \(\frac{- F_0 t}{2m\omega}\text{.}\) The first two terms only oscillate between \(\pm\sqrt{C_1^2 + C_2^2}\text{,}\) which becomes smaller and smaller in proportion to the oscillations of the last term as \(t\) gets larger. In Figure 2.6 we see the graph with \(C_1=C_2=0\text{,}\) \(F_0 = 2\text{,}\) \(m=1\text{,}\) \(\omega = \pi\text{.}\)
Figure 2.6. Graph of \(\frac{1}{\pi} t \sin (\pi t)\text{.}\)
By forcing the system in just the right frequency we produce very wild oscillations. This kind of behavior is called resonance or perhaps pure resonance. Sometimes resonance is desired. For example, remember when as a kid you could start swinging by just moving back and forth on the swing seat in the “correct frequency”? You were trying to achieve resonance. The force of each one of your moves was small, but after a while it produced large swings.
On the other hand resonance can be destructive. In an earthquake some buildings collapse while others may be relatively undamaged. This is due to different buildings having different resonance frequencies. So figuring out the resonance frequency can be very important.
A common (but wrong) example of destructive force of resonance is the Tacoma Narrows bridge failure. It turns out there was a different phenomenon at play 1 .

Subsection 2.6.2 Damped forced motion and practical resonance

In real life things are not as simple as they were above. There is, of course, some damping. Our equation becomes
\begin{equation} mx'' + cx' + kx = F_0 \cos (\omega t) ,\tag{2.11} \end{equation}
for some \(c > 0\text{.}\) We solved the homogeneous problem before. We let
\begin{equation*} p = \frac{c}{2m}, \qquad \omega_0 = \sqrt{\frac{k}{m}} . \end{equation*}
We replace equation (2.11) with
\begin{equation*} x'' + 2px' + \omega_0^2x = \frac{F_0}{m} \cos (\omega t) . \end{equation*}
The roots of the characteristic equation of the associated homogeneous problem are \(r_1,r_2 = -p \pm \sqrt{p^2 - \omega_0^2}\text{.}\) The form of the general solution of the associated homogeneous equation depends on the sign of \(p^2 - \omega_0^2\text{,}\) or equivalently on the sign of \(c^2 - 4km\text{,}\) as before:
\begin{equation*} x_c = \begin{cases} C_1 e^{r_1 t} + C_2 e^{r_2 t} & \text{if } \; c^2 > 4km , \\ C_1 e^{-p t} + C_2 t e^{-p t} & \text{if } \; c^2 = 4km , \\ e^{-p t} \bigl( C_1 \cos (\omega_1 t) + C_2 \sin (\omega_1 t) \bigr) & \text{if } \; c^2 < 4km , \end{cases} \end{equation*}
where \(\omega_1 = \sqrt{\omega_0^2 - p^2}\text{.}\) In any case, we see that \(x_c(t) \to 0\) as \(t \to \infty\text{.}\)
Let us find a particular solution. There can be no conflicts when trying to solve for the undetermined coefficients by trying \(x_p = A \cos (\omega t) + B \sin (\omega t)\text{.}\) Let us plug in and solve for \(A\) and \(B\text{.}\) We get (the tedious details are left to reader)
\begin{equation*} \bigl((\omega_0^2 - \omega^2)B - 2\omega p A\bigr) \sin (\omega t) + \bigl((\omega_0^2 - \omega^2)A + 2\omega p B\bigr) \cos (\omega t) = \frac{F_0}{m} \cos (\omega t) . \end{equation*}
We solve for \(A\) and \(B\text{:}\)
\begin{equation*} \begin{aligned} & A=\frac{(\omega_0^2-\omega^2) F_0} {m{(2\omega p)}^2+m{(\omega_0^2-\omega^2)}^2} , \\ & B=\frac{2 \omega p F_0} {m{(2\omega p)}^2+m{(\omega_0^2-\omega^2)}^2} . \end{aligned} \end{equation*}
We also compute \(C = \sqrt{A^2+B^2}\) to be
\begin{equation*} C = \frac{F_0}{m \sqrt{{(2\omega p)}^2+{(\omega_0^2-\omega^2)}^2}} . \end{equation*}
Thus our particular solution is
\begin{equation*} x_p = \frac{(\omega_0^2-\omega^2) F_0} {m{(2\omega p)}^2+m{(\omega_0^2-\omega^2)}^2} \cos (\omega t) + \frac{2 \omega p F_0} {m{(2\omega p)}^2+m{(\omega_0^2-\omega^2)}^2} \sin (\omega t) . \end{equation*}
Or in the alternative notation we have amplitude \(C\) and phase shift \(\gamma\) where (if \(\omega \not= \omega_0\))
\begin{equation*} \tan \gamma = \frac{B}{A} = \frac{2\omega p}{\omega_0^2-\omega^2} . \end{equation*}
Hence,
\begin{equation*} \mybxbg{~~ x_p = \frac{F_0}{m \sqrt{{(2\omega p)}^2+{(\omega_0^2-\omega^2)}^2}} \cos ( \omega t - \gamma ) . ~~} \end{equation*}
If \(\omega = \omega_0\text{,}\) then \(A=0\text{,}\) \(B = C = \frac{F_0}{2m\omega p}\text{,}\) and \(\gamma = \nicefrac{\pi}{2}\text{.}\)
For reasons we will explain in a moment, we call \(x_c\) the transient solution and denote it by \(x_{tr}\text{.}\) We call the \(x_p\) from above the steady periodic solution and denote it by \(x_{sp}\text{.}\) The general solution is
\begin{equation*} x = x_c + x_p = x_{tr} + x_{sp} . \end{equation*}
The transient solution \(x_c = x_{tr}\) goes to zero as \(t \to \infty\text{,}\) as all the terms involve an exponential with a negative exponent. So for large \(t\text{,}\) the effect of \(x_{tr}\) is negligible and we see essentially only \(x_{sp}\text{.}\) Hence the name transient. Notice that \(x_{sp}\) involves no arbitrary constants, and the initial conditions only affect \(x_{tr}\text{.}\) Thus, the effect of the initial conditions is negligible after some period of time. We might as well focus on the steady periodic solution and ignore the transient solution. See Figure 2.7 for a graph given several different initial conditions.
Figure 2.7. Solutions with different initial conditions for parameters \(k=1\text{,}\) \(m=1\text{,}\) \(F_0 = 1\text{,}\) \(c=0.7\text{,}\) and \(\omega=1.1\text{.}\)
The speed at which \(x_{tr}\) goes to zero depends on \(p\) (and hence \(c\)). The bigger \(p\) is (the bigger \(c\) is), the “faster” \(x_{tr}\) becomes negligible. So the smaller the damping, the longer the “transient region.” This is consistent with the observation that when \(c=0\text{,}\) the initial conditions affect the behavior for all time (i.e. an infinite “transient region”).
Let us describe what we mean by resonance when damping is present. Since there were no conflicts when solving with undetermined coefficient, there is no term that goes to infinity. We look instead at the maximum value of the amplitude of the steady periodic solution. Let \(C\) be the amplitude of \(x_{sp}\text{.}\) If we plot \(C\) as a function of \(\omega\) (with all other parameters fixed), we can find its maximum. We call the \(\omega\) that achieves this maximum the practical resonance frequency. We call the maximal amplitude \(C(\omega)\) the practical resonance amplitude. Thus when damping is present we talk of practical resonance rather than pure resonance. A sample plot for three different values of \(c\) is given in Figure 2.8. As you can see the practical resonance amplitude grows as damping gets smaller, and practical resonance can disappear altogether when damping is large.
Figure 2.8. Graph of \(C(\omega)\) showing practical resonance with parameters \(k=1\text{,}\) \(m=1\text{,}\) \(F_0 = 1\text{.}\) The top line is with \(c=0.4\text{,}\) the middle line with \(c=0.8\text{,}\) and the bottom line with \(c=1.6\text{.}\)
To find the maximum we need to find the derivative \(C'(\omega)\text{.}\) Computation shows
\begin{equation*} C'(\omega) = \frac{- 2\omega( 2p^2+\omega^2-\omega_0^2)F_0} {m {\bigl({(2\omega p)}^2+{(\omega_0^2-\omega^2)}^2\bigr)}^{3/2}} . \end{equation*}
This is zero either when \(\omega = 0\) or when \(2p^2+\omega^2-\omega_0^2 = 0\text{.}\) In other words, \(C'(\omega) = 0\) when
\begin{equation*} \mybxbg{ ~~ \omega = \sqrt{\omega_0^2 - 2p^2} \quad \text{or} \quad \omega = 0 . ~~ } \end{equation*}
If \(\omega_0^2 - 2p^2\) is positive, then \(\sqrt{\omega_0^2 - 2p^2}\) is the practical resonance frequency (that is the point where \(C(\omega)\) is maximal). This follows by the first derivative test for example as then \(C'(\omega) > 0\) for small \(\omega\) in this case. If on the other hand \(\omega_0^2 - 2p^2\) is not positive, then \(C(\omega)\) achieves its maximum at \(\omega=0\text{,}\) and there is no practical resonance since we assume \(\omega > 0\) in our system. In this case the amplitude gets larger as the forcing frequency gets smaller.
If practical resonance occurs, the frequency is smaller than \(\omega_0\text{.}\) As the damping \(c\) (and hence \(p\)) becomes smaller, the practical resonance frequency goes to \(\omega_0\text{.}\) So when damping is very small, \(\omega_0\) is a good estimate of the practical resonance frequency. This behavior agrees with the observation that when \(c=0\text{,}\) then \(\omega_0\) is the resonance frequency.
Another interesting observation to make is that when \(\omega \to \infty\text{,}\) then \(C \to 0\text{.}\) This means that if the forcing frequency gets too high it does not manage to get the mass moving in the mass-spring system. This is quite reasonable intuitively. If we wiggle back and forth really fast while sitting on a swing, we will not get it moving at all, no matter how forceful. Fast vibrations just cancel each other out before the mass has any chance of responding by moving one way or the other.
The behavior is more complicated if the forcing function is not an exact cosine wave, but for example a square wave. A general periodic function will be the sum (superposition) of many cosine waves of different frequencies. The reader is encouraged to come back to this section once we have learned about the Fourier series.

Example 2.6.2.

Geogebra Activity: Use this Geogebra applet 3  to explore the behaviour of a forced damped harmonic oscillator. Try to find the practical resonance for some choice of parameters.

Subsection 2.6.3 Exercises

Exercise 2.6.1.

Derive a formula for \(x_{sp}\) if the equation is \(m x'' + c x' + kx = F_0 \sin (\omega t)\text{.}\) Assume \(c > 0\text{.}\)
Solution.
We know that in such a case, the particular solution should look like \(x_{sp}=A\cos(\omega t)+B\sin(\omega t)\text{.}\) So differentiating this and plugging it into the ODE we get
\begin{equation*} \begin{aligned} \left((\omega_0^2-\omega^2)B-2\omega pA\right)\sin(\omega t)+\left((\omega_0^2-\omega^2)A+2\omega pB\right)\cos(\omega t)=\frac{F_0}{m}\sin(\omega t) \end{aligned} \end{equation*}
Where \(\omega_0\) and \(p\) are defined as before. Solving for \(A\) and \(B\) we get
\begin{equation*} \begin{aligned} &A=\frac{-F_0(2\omega p)}{m(2\omega p)^2+m(\omega_0^2-\omega^2)^2},\quad B=\frac{(\omega_0^2-\omega^2)^2F_0}{m(2\omega p)^2+m(\omega_0^2-\omega^2)^2} \\ & C=\sqrt{A^2+B^2}=\frac{F_0}{m\sqrt{(2\omega p)^2+(\omega_0^2-\omega^2)^2}} \\ & \tan(\gamma)=\frac{B}{A}=-\frac{\omega_0^2-\omega^2}{2\omega p} \end{aligned} \end{equation*}
Therefore
\begin{equation*} \begin{aligned} x_{sp}=\frac{F_0}{m\sqrt{(2\omega p)^2+(\omega_0^2-\omega^2)^2}}\cos(\omega t-\gamma) \end{aligned} \end{equation*}
Note that only the phase shift changes (which makes sense given that the force also only changed by a phase shift).

Exercise 2.6.2.

Derive a formula for \(x_{sp}\) if the equation is \(m x'' + c x' + kx = F_0 \cos (\omega t) + F_1 \cos (3\omega t)\text{.}\) Assume \(c > 0\text{.}\)
Answer.
We find \(x_{sp}\) for each term on the right hand side.
\begin{equation*} \begin{aligned} x_{sp_0}=\frac{F_0}{m\sqrt{(2\omega p)^2+(\omega_0^2-\omega^2)^2}}\cos(\omega t-\gamma),\quad \tan \gamma_0=\frac{2\omega p}{\omega_0^2-\omega^2} \end{aligned} \end{equation*}
And
\begin{equation*} \begin{aligned} x_{sp_1}=\frac{F_0}{m\sqrt{(2(3\omega)p)^2+(\omega_0^2-(3\omega)^2)^2}}\cos(3\omega t-\gamma_1),\quad \tan \gamma_1=\frac{2(3\omega)p}{\omega_0^2-(3\omega)^2} \end{aligned} \end{equation*}
So finally \(x_{sp}=x_{sp_0}+x_{sp_1}\text{.}\)

Exercise 2.6.3.

Derive a formula for \(x_{sp}\) for \(mx''+cx'+kx = F_0 \cos(\omega t) + A\text{,}\) where \(A\) is some constant. Assume \(c > 0\text{.}\)
Answer.
\(x_{sp} = \frac{(\omega_0^2-\omega^2) F_0} {m{(2\omega p)}^2+m{(\omega_0^2-\omega^2)}^2} \cos (\omega t) + \frac{2 \omega p F_0} {m{(2\omega p)}^2+m{(\omega_0^2-\omega^2)}^2} \sin (\omega t) + \frac{A}{k}\text{,}\) where \(p = \frac{c}{2m}\) and \(\omega_0 = \sqrt{\frac{k}{m}}\text{.}\)

Exercise 2.6.4.

Take \(m x'' + c x' + kx = F_0 \cos (\omega t)\text{.}\) Fix \(m > 0\text{,}\) \(k > 0\text{,}\) and \(F_0 > 0\text{.}\) Consider the function \(C(\omega)\text{.}\) For what values of \(c\) (solve in terms of \(m\text{,}\) \(k\text{,}\) and \(F_0\)) will there be no practical resonance (that is, for what values of \(c\) is there no maximum of \(C(\omega)\) for \(\omega > 0\))?
Solution.
There is no practical resonance if \(\omega_0^2-2p^2\leq 0\) so
\begin{equation*} \begin{aligned} (\frac{k}{m})-\frac{c^2}{2m^2}& \leq 0 \\ c \geq \sqrt{2km} \end{aligned} \end{equation*}

Exercise 2.6.5.

Take \(m x'' + c x' + kx = F_0 \cos (\omega t)\text{.}\) Fix \(c > 0\text{,}\) \(k > 0\text{,}\) and \(F_0 > 0\text{.}\) Consider the function \(C(\omega)\text{.}\) For what values of \(m\) (solve in terms of \(c\text{,}\) \(k\text{,}\) and \(F_0\)) will there be no practical resonance (that is, for what values of \(m\) is there no maximum of \(C(\omega)\) for \(\omega > 0\))?
Answer.
As was done in the previous exercise, we find \(m \leq \frac{c^2}{2k}\text{.}\)

Exercise 2.6.6.

A mass of 4 kg on a spring with \(k=\unitfrac[4]{N}{m}\) and a damping constant \(c=\unitfrac[1]{Ns}{m}\text{.}\) Suppose that \(F_0 = \unit[2]{N}\text{.}\) Using forcing function \(F_0 \cos (\omega t)\text{,}\) find the \(\omega\) that causes practical resonance and find the amplitude.
Answer.
\(\omega = \frac{\sqrt{31}}{4\sqrt{2}} \approx 0.984\)     \(C(\omega) = \frac{16}{3\sqrt{7}} \approx 2.016\)

Exercise 2.6.7.

A water tower in an earthquake acts as a mass-spring system. Assume that the container on top is full and the water does not move around. The container then acts as the mass and the support acts as the spring, where the induced vibrations are horizontal. The container with water has a mass of \(m=\unit[10,000]{kg}\text{.}\) It takes a force of 1000 newtons to displace the container 1 meter. For simplicity assume no friction. When the earthquake hits the water tower is at rest (it is not moving). The earthquake induces an external force \(F(t) = m A \omega^2 \cos (\omega t)\text{.}\)
  1. What is the natural frequency of the water tower?
  2. If \(\omega\) is not the natural frequency, find a formula for the maximal amplitude of the resulting oscillations of the water container (the maximal deviation from the rest position). The motion will be a high frequency wave modulated by a low frequency wave, so simply find the constant in front of the sines.
  3. Suppose \(A = 1\) and an earthquake with frequency 0.5 cycles per second comes. What is the amplitude of the oscillations? Suppose that if the water tower moves more than 1.5 meter from the rest position, the tower collapses. Will the tower collapse?
Solution.
No friction means that \(c=0\text{,}\) and the differential equation is
\begin{equation*} \begin{aligned} x''+\omega_0^2x=\frac{F_0}{m}\cos(\omega t) \end{aligned} \end{equation*}
The initial conditions are \(x(0)=0,\ x'(0)=0\text{,}\) and the applied force is \(F(t)=F_0\cos(\omega t),\ F_0=mA\omega^2\text{.}\)
a) \(\omega_0=\sqrt{k/m}=\sqrt{1/10}\text{.}\)
b) For undamped, forced motion with \(\omega\neq \omega_0\) we know that
\begin{equation*} \begin{aligned} x(t)=C_1\cos(\omega_0 t)+C_2\sin(\omega_0 t)+\frac{A\omega^2}{(\omega_0^2-\omega^2)}\cos(\omega t) \end{aligned} \end{equation*}
Applying the initial conditions, we find that
\begin{equation*} \begin{aligned} C_1=\frac{-A\omega^2}{\omega_0^2-\omega^2},\quad C_2=0 \end{aligned} \end{equation*}
Giving
\begin{equation*} \begin{aligned} x&=\frac{A\omega^2}{\omega_0^2-\omega}\left(\cos(\omega t)-\cos(\omega_0 t)\right) \\ &= \underbrace{\frac{2A\omega^2}{\omega_0^2-\omega}}_{C(\omega)} \sin\left( \frac{\omega+\omega_0}{2}t\right)\sin\left( \frac{\omega-\omega_0}{2}t\right) \end{aligned} \end{equation*}

c) Plugging the given values into the answer from (b), we find
\begin{equation*} \begin{aligned} C(\omega)=\frac{2\times 1\times \pi^2}{\pi^2-0.1^2}\approx 2\ \textrm{m} \end{aligned} \end{equation*}
So the tower collapses.

Exercise 2.6.8.

Suppose there is no damping in a mass and spring system with \(m = 5\text{,}\) \(k= 20\text{,}\) and \(F_0 = 5\text{.}\) Suppose \(\omega\) is chosen to be precisely the resonance frequency.
  1. Find \(\omega\text{.}\)
  2. Find the amplitude of the oscillations at time \(t=100\text{,}\) given the system is at rest at \(t=0\text{.}\)
Answer.
a) \(\omega = 2\)     b) \(25\)
K. Billah and R. Scanlan, Resonance, Tacoma Narrows Bridge Failure, and Undergraduate Physics Textbooks, American Journal of Physics, 59(2), 1991, 118–124, http://www.ketchum.org/billah/Billah-Scanlan.pdf 2 
www.ketchum.org/billah/Billah-Scanlan.pdf
www.geogebra.org/m/nhsahq8e
For a higher quality printout use the PDF version: https://www.jirka.org/diffyqs/diffyqs.pdf