Differential Equations for EngineersXYZ Homework Edition

⇩ Download ▾

3.6 Second order systems and applications

Undamped Mass-Spring Systems

While we did say that we will usually only look at first order systems, it is sometimes more convenient to study the system in the way it arises naturally. For example, suppose we have 3 masses connected by springs between two walls. We could pick any higher number, and the math would be essentially the same, but for simplicity we pick 3 right now. Let us also assume no friction, that is, the system is undamped. The masses are m1,m2m_1, m_2, and m3m_3 and the spring constants are k1,k2,k3k_1, k_2, k_3, and k4k_4. Let x1x_1 be the displacement from rest position of the first mass, and x2x_2 and x3x_3 the displacement of the second and third mass. We will make, as usual, positive values go right (as x1x_1 grows, the first mass is moving right). See Figure 1\PageIndex{1}.

Diagram of three wheeled carts m sub 1, m sub 2 and m sub 3 in a row, linked in series by four coil springs k sub 1 through k sub 4 that also anchor the chain to a wall at each end.

Figure 1\PageIndex{1}: System of masses and springs.

This simple system turns up in unexpected places. For example, our world really consists of many small particles of matter interacting together. When we try the above system with many more masses, we obtain a good approximation to how an elastic material behaves. By somehow taking a limit of the number of masses going to infinity, we obtain the continuous one dimensional wave equation (that we study in Section 4.7). But we digress.

Let us set up the equations for the three mass system. By Hooke’s law we have that the force acting on the mass equals the spring compression times the spring constant. By Newton’s second law we have that force is mass times acceleration. So if we sum the forces acting on each mass and put the right sign in front of each term, depending on the direction in which it is acting, we end up with the desired system of equations.

m 1 x 1 = k 1 x 1 + k 2 ( x 2 x 1 ) = ( k 1 + k 2 ) x 1 + k 2 x 2 , m 2 x 2 = k 2 ( x 2 x 1 ) + k 3 ( x 3 x 2 ) = k 2 x 1 ( k 2 + k 3 ) x 2 + k 3 x 3 , m 3 x 3 = k 3 ( x 3 x 2 ) k 4 x 3 = k 3 x 2 ( k 3 + k 4 ) x 3 . \begin{array}{lll}{m_1 x''_1=-k_1 x_1+k_2(x_2-x_1)}&{\quad}&{=-(k_1+k_2)x_1+k_2x_2,} \\ {m_2x''_2=-k_2(x_2-x_1)+k_3(x_3-x_2)}&{\quad}&{=k_2x_1-(k_2+k_3)x_2+k_3x_3,} \\ {m_3x''_3=-k_3(x_3-x_2)-k_4x_3}&{\quad}&{=k_3x_2-(k_3+k_4)x_3.}\end{array} \nonumber

We define the matrices

M = [ m 1 0 0 0 m 2 0 0 0 m 3 ] a n d K = [ ( k 1 + k 2 ) k 2 0 k 2 ( k 2 + k 3 ) k 3 0 k 3 ( k 3 + k 4 ) ] . M = \left[ \begin{array}{ccc} m_1 & 0 & 0 \\ 0 & m_2 & 0 \\ 0 & 0 & m_3 \end{array} \right] ~~~{\rm{~and~}}~~~ K= \left[ \begin{array}{ccc} -(k_1+k_2) & k_2 & 0 \\ k_2 & -(k_2+k_3) & k_3 \\ 0 & k_3 & -(k_3+k_4) \end{array} \right]. \nonumber

We write the equation simply as

M x = K x . M \vec{x}'' = K \vec{x}. \nonumber

At this point we could introduce 3 new variables and write out a system of 6 first order equations. We claim this simple setup is easier to handle as a second order system. We call x\vec{x} the displacement vector, MM the mass matrix, and KK the stiffness matrix.

As with a single equation we want to “divide by MM.” This means computing the inverse of MM. The masses are all nonzero and MMis a diagonal matrix, so comping the inverse is easy:

M 1 = [ 1 m 1 0 0 0 1 m 2 0 0 0 1 m 3 ] . M^{-1}= \left[ \begin{array}{ccc} \frac{1}{m_1} & 0 & 0 \\ 0 & \frac{1}{m_2} & 0 \\ 0 & 0 &\frac{1}{m_3}\end{array} \right]. \nonumber

This fact follows readily by how we multiply diagonal matrices. As an exercise, you should verify that MM1=M1M=I.MM^{-1}=M^{-1}M=I.

Let A=M1KA = M^{-1}K. We look at the system x=M1Kx\vec{x}'' = M^{-1} K \vec{x}, or

x = A x . \vec{x}'' = A \vec{x}. \nonumber

Many real world systems can be modeled by this equation. For simplicity, we will only talk about the given masses-and-springs problem. We try a solution of the form

x = v e α t . \vec{x} =\vec{v} e^{\alpha t}. \nonumber

We compute that for this guess, x=α2veαt.\vec{x}'' =\alpha^2 \vec{v} e^{\alpha t}. We plug our guess into the equation and get

α 2 v e α t = A v e α t . \alpha^2 \vec{v} e^{\alpha t} =A \vec{v} e^{\alpha t}. \nonumber

We divide by eαte^{\alpha t} to arrive at α2v=Av\alpha^2 \vec{v}=A\vec{v}. Hence if α2\alpha^2 is an eigenvalue of AA and v\vec{v} is a corresponding eigenvector, we have found a solution.

In our example, and in other common applications, AA has only real negative eigenvalues (and possibly a zero eigenvalue). So we study only this case. When an eigenvalue λ\lambda is negative, it means that α2=λ\alpha^2= \lambda is negative. Hence there is some real number ω\omega such that ω2=λ-\omega^2 = \lambda. Then α=±iω\alpha= \pm i \omega. The solution we guessed was

x = v ( cos ( ω t ) + i sin ( ω t ) ) . \vec{x}= \vec{v}(\cos(\omega t) +i \sin(\omega t)). \nonumber

By taking the real and imaginary parts (note that v\vec{v} is real), we find that vcos(ωt)\vec{v}\cos(\omega t) and vsin(ωt)\vec{v}\sin(\omega t) are linearly independent solutions.

If an eigenvalue is zero, it turns out that both v\vec{v} and vt\vec{v}t are solutions, where v\vec{v} is an eigenvector corresponding to the eigenvalue 0.

We use this solution and the setup from the introduction of this section even when some of the masses and springs are missing. For example, when there are only 2 masses and only 2 springs, simply take only the equations for the two masses and set all the spring constants for the springs that are missing to zero.

x = [ 1 2 ] c 1 cos ( t α 1 ) + [ 1 1 ] c 2 cos ( 2 t α 2 ) . \vec{x}= \left[ \begin{array}{c} 1 \\ 2 \end{array} \right] c_1 \cos (t-\alpha_1) + \left[ \begin{array}{c} 1 \\ -1 \end{array} \right] c_2 \cos (2t- \alpha_2). \nonumber

The first term,

[ 1 2 ] c 1 cos ( t α 1 ) = [ c 1 cos ( t α 1 ) 2 c 1 cos ( t α 1 ) ] , \left[ \begin{array}{c} 1 \\ 2 \end{array} \right] c_1 \cos (t-\alpha_1) = \left[ \begin{array}{c} c_1 \cos (t-\alpha_1) \\ 2c_1 \cos (t-\alpha_1) \end{array} \right], \nonumber

corresponds to the mode where the masses move synchronously in the same direction.

The second term,

[ 1 1 ] c 2 cos ( 2 t α 2 ) = [ c 2 cos ( 2 t α 2 ) c 2 cos ( 2 t α 2 ) ] , \left[ \begin{array}{c} 1 \\ -1 \end{array} \right] c_2 \cos (2t-\alpha_2) = \left[ \begin{array}{c} c_2 \cos (2t-\alpha_2) \\ -c_2 \cos (2t-\alpha_2) \end{array} \right], \nonumber

corresponds to the mode where the masses move synchronously but in opposite directions.

The general solution is a combination of the two modes. That is, the initial conditions determine the amplitude and phase shift of each mode. As an example, suppose we have initial conditions

x ( 0 ) = [ 1 1 ] , x ( 0 ) = [ 0 6 ] . \vec{x}(0) = \begin{bmatrix} 1 \\ -1 \end{bmatrix}, \qquad \vec{x}'(0) = \begin{bmatrix} 0 \\ 6 \end{bmatrix}. \nonumber

We use the aj,bja_j, b_j constants to solve for initial conditions. First

[ 1 1 ] = x ( 0 ) = [ 1 2 ] a 1 + [ 1 1 ] a 2 = [ a 1 + a 2 2 a 1 a 2 ] . \begin{bmatrix} 1 \\ -1 \end{bmatrix} = \vec{x}(0) = \begin{bmatrix} 1 \\ 2 \end{bmatrix} a_1 + \begin{bmatrix} 1 \\ -1 \end{bmatrix} a_2 = \begin{bmatrix} a_1+a_2 \\2a_1 - a_2 \end{bmatrix}. \nonumber

We solve (exercise) to find a1=0a_1 = 0, a2=1a_2 = 1. To find the b1b_1 and b2b_2, we differentiate first:

x=[12](a1sin(t)+b1cos(t))+[11](2a2sin(2t)+2b2cos(2t)).{\vec{x}}' = \begin{bmatrix} 1 \\ 2 \end{bmatrix} \bigl( - a_1 \sin (t) + b_1 \cos (t) \bigr) + \begin{bmatrix} 1 \\ -1 \end{bmatrix} \bigl( - 2a_2 \sin (2t) + 2 b_2 \cos (2t) \bigr). \nonumber

Now we solve:

[06]=x(0)=[12]b1+[11]2b2=[b1+2b22b12b2].\begin{bmatrix} 0 \\ 6 \end{bmatrix} = {\vec{x}}'(0) = \begin{bmatrix} 1 \\ 2 \end{bmatrix} b_1 + \begin{bmatrix} 1 \\ -1 \end{bmatrix} 2 b_2 = \begin{bmatrix} b_1+2b_2 \\ 2b_1-2b_2 \end{bmatrix}. \nonumber

Again solve (exercise) to find b1=2b_1 = 2, b2=1b_2 = -1. So our solution is

x=[12]2sin(t)+[11](cos(2t)sin(2t))=[2sin(t)+cos(2t)sin(2t)4sin(t)cos(2t)+sin(2t)].\vec{x} = \begin{bmatrix} 1 \\ 2 \end{bmatrix} 2 \sin (t) + \begin{bmatrix} 1 \\ -1 \end{bmatrix} \bigl( \cos (2t) - \sin (2t) \bigr) = \begin{bmatrix} 2 \sin (t) + \cos(2t)- \sin(2t) \\ 4 \sin (t) - \cos(2t) + \sin(2t) \end{bmatrix}. \nonumber

The graphs of the two displacements, x1x_1 and x2x_2 of the two carts is in Figure 4\PageIndex{4}.

Plot of two displacement curves against t from 0 to 10: the green one is larger, ranging from about 5.2 down to -3.9, while the smaller blue one is visibly rippled by the faster mode.
Figure 4\PageIndex{4}: Superposition of the two modes given the initial conditions.

OK, let us first set the system up. Let t=0t=0 be the time when the two cars link up. Let x1x_1 be the displacement of the first car from the position at t=0t=0, and let x2x_2 be the displacement of the second car from its original location. Then the time when x2(t)=10x_2(t)=10 is exactly the time when impact with wall occurs. For this t,x2(t)t,x'_2(t) is the speed at impact. This system acts just like the system of the previous example but without k1k_1. Hence the equation is

[ 2 0 0 1 ] x = [ 2 2 2 2 ] x . \left[ \begin{array}{cc} 2 & 0 \\ 0 & 1 \end{array} \right] \vec{x}''= \left[ \begin{array}{cc} -2 & 2 \\ 2 & -2 \end{array} \right] \vec{x}. \nonumber

or

x = [ 1 1 2 2 ] x . \vec{x}''= \left[ \begin{array}{cc} -1 & 1 \\ 2 & -2 \end{array} \right] \vec{x}. \nonumber

We compute the eigenvalues of AA. It is not hard to see that the eigenvalues are 00 and 3-3 (exercise). Furthermore, eigenvectors are [11]\left[ \begin{array}{c} 1 \\ 1 \end{array} \right] and [12]\left[ \begin{array}{c} 1 \\ -2 \end{array} \right] respectively (exercise). Then ω2=3\omega_2= \sqrt{3} and by the second part of the theorem we find our general solution to be

x = [ 1 1 ] ( a 1 + b 1 t ) + [ 1 2 ] ( a 2 cos ( 3 t ) + b 2 sin ( 3 t ) ) = [ a 1 + b 1 t + a 2 cos ( 3 t ) + b 2 sin ( 3 t ) a 1 + b 1 t 2 a 2 cos ( 3 t ) 2 b 2 sin ( 3 t ) ] \begin{align}\begin{aligned} \vec{x}&= \left[ \begin{array}{c} 1 \\ 1 \end{array} \right](a_1+b_1t) + \left[ \begin{array}{c} 1 \\ -2 \end{array} \right](a_2 \cos(\sqrt3 t)+ b_2 \sin(\sqrt3 t)) \\ &=\left[ \begin{array}{c} a_1 + b_1t+a_2 \cos(\sqrt3 t)+ b_2 \sin(\sqrt3 t) \\ a_1+ b_1t- 2a_2 \cos(\sqrt3 t)- 2b_2 \sin(\sqrt3 t) \end{array} \right]\end{aligned}\end{align} \nonumber

We now apply the initial conditions. First the cars start at position 00 so x1(0)=0x_1(0)=0 and x2(0)=0x_2(0)=0. The first car is traveling at 3ms3\:\frac{\text{m}}{\text{s}}, so x1(0)=3x'_1(0)=3 and the second car starts at rest, so x2(0)=0x'_2(0)=0. The first conditions says

0 = x ( 0 ) = [ a 1 + a 2 a 1 2 a 2 ] . \vec{0} = \vec{x}(0)= \left[ \begin{array}{c} a_1+ a_2 \\ a_1-2a_2 \end{array} \right]. \nonumber

It is not hard to see that a1=a2=0a_1=a_2=0. We set a1=0a_1=0 and a2=0a_2=0 in x(t)\vec{x}(t) and differentiate to get

x ( t ) = [ b 1 + 3 b 2 cos ( 3 t ) b 1 2 3 b 2 cos ( 3 t ) ] . \vec{x}'(t)= \left[ \begin{array}{c} b_1+ \sqrt3 b_2 \cos(\sqrt3 t)\\ b_1- 2 \sqrt3 b_2 \cos(\sqrt3 t) \end{array} \right]. \nonumber

So

[ 3 0 ] = x ( 0 ) = [ b 1 + 3 b 2 b 1 2 3 b 2 ] . \left[ \begin{array}{c} 3\\ 0 \end{array} \right] = \vec{x}'(0)= \left[ \begin{array}{c} b_1+ \sqrt3 b_2 \\ b_1- 2 \sqrt3 b_2 \end{array} \right]. \nonumber

Solving these two equations we find b1=2b_1 = 2 and b2=13b_2 = \frac{1}{\sqrt3}. Hence the position of our cars is (until the impact with the wall)

x = [ 2 t + 1 3 sin ( 3 t ) 2 t 2 3 sin ( 3 t ) ] . \vec{x}= \left[ \begin{array}{c} 2t +\frac{1}{\sqrt3} \sin(\sqrt3 t) \\ 2t -\frac{2}{\sqrt3} \sin(\sqrt3 t)\end{array} \right]. \nonumber

Note how the presence of the zero eigenvalue resulted in a term containing tt. This means that the carts will be traveling in the positive direction as time grows, which is what we expect.

What we are really interested in is the second expression, the one for x2x_2. We have x2(t)=2t23sin(3t)x_2(t)= 2t -\frac{2}{\sqrt3} \sin(\sqrt3 t). See Figure 6\PageIndex{6} for the plot of x2x_2 versus time.

Plot of the second car's position against time, t from 0 to 6, rising like a staircase: it climbs to a flat shelf near 7.3, pauses, then surges to about 12.9, never turning back.
Figure 6\PageIndex{6}: Position of the second car in time (ignoring the wall).

Just from the graph we can see that time of impact will be a little more than 5 seconds from time zero. For this we have to solve the equation 10=x2(t)=2t23sin(3t)10=x_{2}(t)=2t-\frac{2}{\sqrt{3}}\sin (\sqrt{3}t). Using a computer (or even a graphing calculator) we find that timpact5.22t_{\text{impact}}\approx 5.22 seconds.

As for the speed we note that x2=22cos(3t)x'_2 = 2-2 \cos(\sqrt3 t). At time of impact (5.225.22 seconds from t=0t=0) we get that x2(timpact)3.85x'_2(t_{\rm{impact}}) \approx 3.85.

The maximum speed is the maximum of 22cos(3t)2-2 \cos(\sqrt3 t), which is 44. We are traveling at almost the maximum speed when we hit the wall.

Suppose that Bob is a tiny person sitting on car 2. Bob has a Martini in his hand and would like not to spill it. Let us suppose Bob would not spill his Martini when the first car links up with car 2, but if car 2 hits the wall at any speed greater than zero, Bob will spill his drink. Suppose Bob can move car 2 a few meters towards or away from the wall (he cannot go all the way to the wall, nor can he get out of the way of the first car). Is there a “safe” distance for him to be at? A distance such that the impact with the wall is at zero speed?

Interactive figureHow far from the wall can Bob keep his Martini?Drag the Distance from car 2 to the wall slider from 5 to 11.
A solid curve climbs from the origin, steepening and then briefly levelling off before steepening again; a horizontal dashed line marks the wall, and a dotted curve below traces the speed, dropping to touch zero exactly beneath each level stretch. Where solid curve meets dashed line is the impact. Lowering the wall slides that crossing left onto the flat stretch until the line merely grazes it — the one distance at which the car arrives with no speed at all. Adjustable parameter: Distance from car 2 to the wall (d) = 10 m. Viewing window: x from -6 to 13, y from -0.28 to 11.48.
XYZ Graph · viewer build 5edf91b
Car 2's position x₂(t) = 2t − (2/√3)·sin(√3·t) (solid), the wall as a movable line at distance d (dashed), and the impact speed x₂′(t) = 2 − 2·cos(√3·t) (dotted); the horizontal axis is time t in seconds. At Lebl's d = 10 m the wall cuts the curve on a steep stretch — impact at t ≈ 5.22 s, speed ≈ 3.85 m/s. Drag d down and the crossing slides left and flattens until, at d ≈ 7.26 m (= 4π/√3), the line is tangent to the plateau at t = 2π/√3, right where the dotted speed curve touches zero. That tangency is Bob's safe distance; drag past it and the curve steepens again.

The answer is yes. Looking at Figure 6\PageIndex{6}, we note the “plateau” between t=3t=3 and t=4t=4. There is a point where the speed is zero. To find it we need to solve x2(t)=0x'_2(t)=0. This is when cos(3t)=1\cos(\sqrt3 t) = 1 or in other words when t=2π3,4π3,t=\frac{2 \pi}{\sqrt3}, \frac{4 \pi}{\sqrt3}, \ldots and so on. We plug in the first value to obtain x2(2π3)=4π37.26x_2 \left( \frac{2 \pi}{\sqrt3} \right)=\frac{4 \pi}{\sqrt3} \approx 7.26. So a “safe” distance is about 7 and a quarter meters from the wall.

Alternatively Bob could move away from the wall towards the incoming car 2 where another safe distance is 8π314.51\frac{8 \pi}{\sqrt3 } \approx 14.51 and so on, using all the different tt such that x2(t)=0x'_2(t)=0. Of course t=0t=0 is always a solution here, corresponding to x2=0x_2=0, but that means standing right at the wall.

Forced Oscillations

Finally we move to forced oscillations. Suppose that now our system is

x = A x + F cos ( ω t ) . \label{eq:30} \vec{x}'' = A \vec{x} + \vec{F} \cos(\omega t).

(3.6.1)

That is, we are adding periodic forcing to the system in the direction of the vector F\vec{F}.

As before, this system just requires us to find one particular solution xp\vec{x}_p, add it to the general solution of the associated homogeneous system xc\vec{x}_c, and we will have the general solution to (3.6.1). Let us suppose that ω\omega is not one of the natural frequencies of x=Ax\vec{x}'' = A \vec{x}, then we can guess

x p = c cos ( ω t ) , \vec{x}_p= \vec{c} \cos(\omega t), \nonumber

where c\vec{c} is an unknown constant vector. Note that we do not need to use sine since there are only second derivatives. We solve for c\vec{c} to find xp\vec{x}_p. This is really just the method of undetermined coefficients for systems. Let us differentiate xp\vec{x}_p twice to get

x p = ω 2 c cos ( ω t ) . \vec{x}''_p= -\omega^2 \vec{c} \cos(\omega t). \nonumber

Plug xp\vec{x}_p and xp\vec{x}''_p into the equation (3.6.1):

ω 2 c cos ( ω t ) x p = A c cos ( ω t ) A x p + F cos ( ω t ) . \overbrace{ -\omega^2 \vec{c} \cos (\omega t) }^{{\vec{x}_p}''} = \overbrace{ A \vec{c} \cos (\omega t) }^{A \vec{x}_p} + \vec{F} \cos (\omega t). \nonumber

We cancel out the cosine and rearrange the equation to obtain

( A + ω 2 I ) c = F . (A + \omega^2 I) \vec{c} = - \vec{F}. \nonumber

So

c = ( A + ω 2 I ) 1 ( F ) . \vec{c} = (A + \omega^2 I)^{-1} (- \vec{F}). \nonumber

Of course this is possible only if (A+ω2I)=(A(ω2)I)(A + \omega^2 I)= (A -(- \omega^2) I) is invertible. That matrix is invertible if and only if ω2- \omega^2 is not an eigenvalue of AA. That is true if and only if ω\omega is not a natural frequency of the system.

We simplified things a little bit. If we wish to have the forcing term to be in the units of force, say Newtons, then we must write

Mx=Kx+Gcos(ωt).M \vec{x}'' = K \vec{x} + \vec{G} \cos(\omega t). \nonumber

If we then write things in terms of A=M1KA = M^{-1} K, we have

x=M1Kx+M1Gcos(ωt)orx=Ax+Fcos(ωt),\vec{x}'' = M^{-1}K \vec{x} + M^{-1} \vec{G} \cos(\omega t) \qquad \text{or} \qquad \vec{x}'' = A \vec{x} + \vec{F} \cos(\omega t), \nonumber

where F=M1G\vec{F} = M^{-1} \vec{G}.

Note that given force f\vec{f}, we write the equation as Mx=Kx+fM\vec{x}''=K\vec{x}+\vec{f} to get the units right. Then we write x=M1Kx+M1f\vec{x}''=M^{-1}K\vec{x}+M^{-1}\vec{f}. The term g=M1f\vec{g}=M^{-1}\vec{f} in x=Ax+g\vec{x}''=A\vec{x}+\vec{g} is in units of force per unit mass.

If ω\omega is a natural frequency of the system resonance occurs because we will have to try a particular solution of the form

x p = c t sin ( ω t ) + d t cos ( ω t ) . \vec{x}_p = \vec{c}t \sin(\omega t)+ \vec{d}t \cos(\omega t). \nonumber

That is assuming that the eigenvalues of the coefficient matrix are distinct. Next, note that the amplitude of this solution grows without bound as tt grows.

Adapted from Differential Equations for Engineers by Jiří Lebl (https://www.jirka.org/diffyqs/), © Jiří Lebl, licensed under CC BY-SA 4.0. Changes were made. License: CC-BY-SA-4.0.

These eBooks are a prerelease and are not yet certified conformant with WCAG 2.1 AA or ADA Title II. Every page is built against an automated accessibility gate, and the published editions will meet ADA Title II requirements when they release in late September 2026. If something is unusable, please tell us.