2  Second-order linear IVPs

Code 
using Plots, DifferentialEquations
default(label="", linewidth=3, markersize=4, size=(500,320))

Definition 2.1 (Second-order ODE) A second-order ODE is an equation of the form

\[ x'' = F(t,x,x'), \tag{2.1}\]

to be solved for \(x(t)\).

Note

Technically, this is not the most general form possible, since an equation in the implicit form \(F(t,x,x',x'')=0\) might be impossible to manipulate into an explicit expression for \(x''\). We won’t be considering such equations.

As with first-order problems, we can learn a little from the trivial case \(x''=0\). This ODE implies that \(x'\) is constant, so \(x(t)=c_1 t + c_2\) for arbitrary constants \(c_1,c_2\). While a generalization of this observation is far from obvious at this point, it is nonetheless true.

Important

The general solution of a second-order ODE has two integration constants.

As a consequence, we note that two initial conditions are now needed to specify the integration constants to get a unique solution.

Definition 2.2 (Second-order IVP) A second-order IVP (initial-value problem) is

\[ x'' = F(t,x,x'), \quad x(a) = x_0,\, x'(a) = v_0. \tag{2.2}\]

2.1 Linear problems

We will focus on linear problems.

Definition 2.3 A linear, 2nd-order ODE in standard form is \[ x'' + p(t)x' + q(t) x = f(t), \tag{2.3}\] where \(p\), \(q\), and \(f\) are specified functions of the independent variable.

For a linear initial-value problem, we would also specify two conditions \(x(a)=x_0\) and \(x'(a)=v_0\), like in Equation 2.2.

Equation Equation 2.3 is in the form \(\opA[x]=f\) for the linear operator defined by

\[ \opA[x] = \ddd{x}{t} + p \dd{x}{t} + q x. \]

Much of what we learned using the abstract language of operators we introduced for first-order problems applies here as well. In particular, we still have a superposition theorem like Theorem 1.1 and a solving strategy like Theorem 1.3:

  1. Find the general solution \(x_h\) of the associated homogeneous problem \(\opA[x]=0\).
  2. Find any particular solution \(x_p\) of the original \(\opA[x]=f\).
  3. Add them to get the general solution.
  4. If given, use initial conditions to solve for the integration constants.

Steps 1 and 2 above will require some new explanations, however, beginning with a major new wrinkle.

2.1.1 Independence

Consider again the trivial problem \(x''=0\). It’s certainly true that \(x=t\) and \(x=5t\) are different solutions of this equation. So it would be tempting to write the general solution of this homogeneous problem as

\[ c_1\cdot (t) + c_2\cdot (5t) \]

for arbitrary constants \(c_1\) and \(c_2\). While this expression does always give a solution, it doesn’t give all of them, such as the constant solution \(x\equiv 2\). The issue is that \(t\) and \(5t\) are not different enough, in a way we now make precise.

Definition 2.4 (Linear independence) Functions \(u_1(t),u_2(t),\ldots,u_n(t)\) are linearly dependent on interval \(I\) if there is a way to choose constants \(c_1,\ldots,c_n\), not all zero, such that

\[ c_1u_1(t) + c_2u_2(t) + \cdots c_n u_n(t) = 0 \]

for all \(t \in I\). If the functions are not linearly dependent, then they are linearly independent on \(I\).

Note

Linear independence is a notion that extends well beyond differential equations, as we will see.

Example 2.1 The functions \(u_1=t\) and \(u_2=5t\) are linearly dependent on any interval, because \(-5u_1 + u_2 = -5t+5t=0\) for all \(t\).

Example 2.2 The functions \(u_1=\cos(t)\) and \(u_2=\sin(t)\) are independent on the real line \((-\infty,\infty)\), because the equation

\[ c_1\cos(t) + c_2\sin(t) = 0 \]

is equivalent to \(\tan(t)=-c_1/c2\) if \(c_2\neq 0\), which is impossible for all \(t\). Hence \(c_2=0\), which implies that \(c_1 \cos(t)=0\) for all \(t\), and this requires \(c_1=0\) as well.

2.1.2 Wronskian

There is a mechanical way to determine the independence of given functions.

Definition 2.5 Suppose \(u_1,\ldots,u_n\) all have at least \(n-1\) continuous derivatives. Their Wronskian is

\[ W(u_1,\ldots,u_n) = \begin{vmatrix} u_1 & u_2 & \cdots & u_n \\ u_1' & u_2' & \cdots & u_n' \\ \vdots & \vdots & & \vdots \\ u_1^{(n-1)} & u_2^{(n-1)} & \cdots & u_n^{(n-1)} \end{vmatrix}. \tag{2.4}\]

Note

Equation 2.4 is a determinant, which we will not learn about in general for a while yet. In the \(2\times 2\) case, it’s

\[ \begin{vmatrix} a & b \\ c & d \end{vmatrix} = ad-bc. \]

In the \(3\times 3\) case, it’s the same process as finding cross products in vector calculus. You can look ahead at Example 5.1 for a refresher.

Here’s the punch line.

Theorem 2.1 If \(x_1,\ldots,x_n\) are solutions of a linear homogeneous ODE \(\opA[x]=0\), all defined for all \(t\in I\), then they are linearly independent on \(I\) if and only if their Wronskian is nonzero at all \(t\in I\).

The upshot is that we can check the Wronskian at any one value of \(t\) we want; zero means linear dependence, and nonzero means independence.

Example 2.3 You can check that \(\cos(t)\) and \(\sin(t)\) are solutions of \(x''+x=0\). Since

\[ W(t) = \twodet{\cos(t)}{\sin(t)}{-\sin(t)}{\cos(t)} = \cos^2(t) + \sin^2(t) = 1, \]

they are linearly independent.

2.1.3 General solution

Now we can state precisely what goes into a general solution in the linear, second-order case.

Theorem 2.2 Suppose \(x_1\) and \(x_2\) are linearly independent solutions of the second-order linear problem \(\opA[x]=0\). Then the general solution of \(\opA[x]=0\) is

\[ x_h(t) = c_1 x_1(t) + c_2 x_2(t), \tag{2.5}\]

where \(c_1\) and \(c_2\) are arbitrary constants. Furthermore, the general solution of the nonhomogeneous problem \(\opA[x]=f\) is

\[ x(t) = x_h(t) + x_p(t), \tag{2.6}\]

where \(x_h\) is the homogeneous solution Equation 2.5 and \(x_p\) is any one particular solution of \(\opA[x]=f\).

2.1.4 Constant coefficients

The general linear second-order problem \(\opA[x]=f\) has the structure described above, but only a few such problems are straightforward to solve. When it comes to actually finding solutions, we will limit ourselves to a friendly subset of problems:

Definition 2.6 The second-order, linear, constant-coefficient problem in standard form is \[ x'' + bx' + k x = f(t), \tag{2.7}\]

where \(b\) and \(k\) are constants (positive, negative, or zero). The function \(f(t)\) is the forcing function.

Note

In Equation 2.7 the coefficient \(b\) has units \(1/T\), the coefficient \(k\) has units \(1/T^2\), and \(f(t)\) has units \(X/T^2\), where \(X\) and \(T\) are the units of \(x\) and \(t\), respectively.

2.2 Homogeneous solutions

We start with solutions of the homogeneous equation

\[ x'' + b x' + k x = 0. \]

We will use \(x_h\) to name the general solution, since it might be part of the solution to a nonhomogeneous equation.

2.2.1 Characteristic polynomial

We employ the super-sophisticated method of guessing the answer and then proving we’re right. This starts with

\[ x_h(t) = e^{\lambda t} \]

for a to-be-determined value of \(\lambda\). We substitute \(x_h\) into the ODE and get

\[ \begin{split} 0 &= \lambda^2 e^{\lambda t} + b \bigl(\lambda e^{\lambda t}\bigr) + k \bigl(e^{\lambda t}\bigr)\\ &= e^{\lambda t} \bigl(\lambda^2 + b \lambda + k \bigr). \end{split} \]

We therefore know that \(x_h\) is a homogeneous solution provided that \(\lambda^2 + b \lambda + k =0\).

Definition 2.7 (Characteristic polynomial of a linear ODE) The characteristic polynomial of the linear operator \(\opA[x]=x''+bx'+kx\) is

\[ p(s) = s^2 + b s + k. \]

Its roots are called the characteristic values of \(\opA\). For brevity, we will often refer to these as simply the roots of \(\opA\), although this usage is not standard.

Theorem 2.3 (Homogeneous solution, distinct roots) Let \(\lambda_1,\lambda_2\) be the characteristic values of \(\opA\), where \(\opA[x]=x''+bx'+kx.\) If \(\lambda_1\neq \lambda_2\), then the general solution of \(\opA[x]= 0\) is

\[ x_h(t) = c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t}. \tag{2.8}\]

Part of what makes Equation 2.8 work is that \(e^{\lambda_1 t}\) and \(e^{\lambda_2 t}\) are independent solutions, because

\[ \twodet{e^{\lambda_1 t}}{e^{\lambda_2 t}}{\lambda_1 e^{\lambda_1 t}}{\lambda_2 e^{\lambda_2 t}} = (\lambda_2-\lambda_1)e^{(\lambda_1+\lambda_2)t}, \]

which is never zero.

We have to handle the case \(\lambda_1=\lambda_2\) separately.

Theorem 2.4 (Homogeneous solution, repeated root) Let \(\lambda\) be a double root of \(\opA\), where \(\opA[x]=x''+bx'+kx.\) The general solution of \(\opA[x]= 0\) is

\[ x_h(t) = e^{\lambda t} \bigl( c_1 t + c_2 \bigr). \tag{2.9}\]

Example 2.4 The trivial problem \(x''=0\) is linear, with \(b=k=0\). The characteristic polynomial is \(s^2\), which makes \(\lambda=0\) a double root. So the general solution is

\[ x_h(t) = e^{0t} (c_1 t + c_2) = c_1 t + c_2. \]

This was not exactly a mystery to us, but it’s a good idea to check new formulas on cases you already understand.

Example 2.5 Find the general solution of \(x''-x'-2x=0\).

Solution. The characteristic polynomial is \(s^2-s-2\), which has roots \(\lambda_1=-1\), \(\lambda_2=2\). This gives the general solution \(x_h(t)=c_1 e^{-t} + c_2 e^{2t}\).

Example 2.6 Solve the IVP \(x'' - 5 x = 0\), \(x(0)=6\), \(x'(0)=0\).

Solution. The roots of \(s^2-5\) are \(\lambda_1=\sqrt{5}\), \(\lambda_2=-\sqrt{5}\). The general solution is

\[ x(t) = c_1 e^{\sqrt{5}\, t} + c_2 e^{-\sqrt{5}\, t}. \]

The initial conditions lead to

\[\begin{split} 6 &= x(0) = c_1 e^0 + c_2 e^0 = c_1 + c_2, \\ 0 &= x'(0) = \sqrt{5} c_1 e^0 - \sqrt{5} c_2 e^0 = \sqrt{5}(c_1-c_2). \end{split}\]

It’s easy to conclude from here that \(c_1=c_2=3\). In general we might have to solve a \(2\times 2\) linear algebraic system for the constants.

All this seems simple enough. But remember that the roots of a real quadratic polynomial may come as a pair of complex conjugate numbers. How then do we interpret Equation 2.8 or Equation 2.9? The answer is what makes second-order problems really different from first-order ones, and we tackle it by taking a brief deep dive into complex numbers.

2.3 Complex characteristic values

Code 
using Plots,LaTeXStrings,Printf
default(linewidth=2,label="",size=(400,280))
using Logging
disable_logging(Logging.Info);

For a really smart and entertaining introduction to complex numbers, I recommend this video series.

There’s often a lot of wariness about complex numbers. Terminology is part of the reason. Using “real” and “imaginary” to label numbers is the residue of a value judgment that ruled mathematics for centuries. But complex numbers are actually just as “real” as so-called real numbers. If anything, they are arguably more fundamental.

As a practical matter, you can pretty much always replace a complex value with two real ones, and vice versa. But the complex point of view often provides a unifying perspective. Sometimes the formula manipulations are easier in the complex form, too; in particular, you may be able to replace trigonometry with algebra.

2.3.1 Basic arithmetic

We can write a complex number \(z\in \complex\) as \(z=x+iy\), where \(i^2=-1\) and \(x\) and \(y\) are real numbers known as the real part and imaginary part of \(z\),

\[ z = x + i y \quad \Leftrightarrow \quad x = \text{Re}(z), \quad y = \text{Im}(z). \]

Writing a complex number this way is equivalent to using rectangular or Cartesian coordinates in the plane to specify a point \((x,y)\). Thus complex numbers are on a number plane rather than a number line.

The generalization of absolute value to complex numbers is known as the modulus or magnitude, a real nonnegative quantity defined as

\[ |z|^2 = [\text{Re}(z)]^2 + [\text{Im}(z)]^2. \]

Like the absolute value, \(|z|\) is the distance from \(z\) to the origin, and \(|w-z|\) is the distance between two complex numbers. Here, though, the distances are in the plane rather than on a line.

An important operation on complex numbers that has no real-value counterpart is the conjugate,

\[ \bar{z} =\text{Re}(z) - i \,\text{Im}(z). \]

Geometrically the conjugate is the reflection across the real axis of the plane. No matter how complicated an expression is, you just replace \(i\) by \(-i\) everywhere to get the conjugate.

You add, subtract, and multiply complex numbers by applying the usual algebraic rules, applying \(i^2=-1\) as needed. They should give little trouble. Division can be a little trickier, even though the rules are always the same. One trick is to give a complex ratio a purely real denominator:

\[ \frac{w}{z} = \frac{w \bar{z}}{z \bar{z}} = \frac{w \bar{z}}{|z|^2}. \]

This is a lot like rationalizing a denominator with square roots.

Example 2.7 \[ \frac{2-i}{3+2i} = \frac{2-i}{3+2i}\cdot \frac{3-2i}{3-2i} = \frac{6-4i-3i+2i^2}{3^2+2^2} = \frac{4-7i}{13}. \]

Example 2.8 Suppose that \(|z|=1\). Then

\[ \frac{1}{z} = \frac{1}{z}\cdot \frac{\bar{z}}{\bar{z}} = \frac{\bar{z}}{|z|^2} = \bar{z}. \]

Tip

Memorize the special case \(\dfrac{1}{i} = -i\).

Here are some more simple rules to know:

Theorem 2.5 (Complex arithmetic) For complex numbers \(w\) and \(z\),

  1. \(|\bar{z}| = |z|\)
  2. \(|wz| = |w|\cdot |z|\)
  3. \(|w+z|\le |w| + |z|\) (triangle inequality)
  4. \(\left| \dfrac{1}{z} \right| = \dfrac{1}{|z|}\)
  5. \(\overline{wz}=\bar{w}\cdot \bar{z}\)
  6. \(\overline{w+z}=\bar{w}+\bar{z}\)
  7. \(\overline{\left(\dfrac{1}{z} \right)} = \dfrac{1}{\bar{z}}\)
  8. \(|z|^2 = z\cdot \bar{z}\)

2.3.2 Euler’s identity

For second-order equations we need to make sense of \(e^{\lambda t}\) when \(\lambda\) is complex. The key is one of the most famous equations in mathematics.

Theorem 2.6 (Euler’s identity) \[ :label: secondlin-euler e^{it} = \cos(t) + i \sin(t) \]

Note

Euler is a German name, and as such, it rhymes with “boiler,” not “Bueller.”

Euler’s identity shows that imaginary exponents produce oscillation, rather than the growth/decay of real exponents.

Code 
c = t -> real(exp(1im*t))
s = t -> imag(exp(1im*t))
an = @animate for t in range(0,4π,length=60)
    plot([c s],0,t,label=["Re part" "Im part"],xaxis=(L"t",(0,4π)),yaxis=(L"e^{it}",(-1.05,1.05)))
    title!("Complex exponential")
end
gif(an, "cpxexp.gif")

Alternatively, if we take \(x\) and \(y\) to be the real and imaginary parts of \(z=e^{it}\), then they parametrically describe a unit circle in the complex plane.

Code 
an = @animate for t in range(0,4π,length=60)
    plot(c,s,0,t,aspect_ratio=1,xaxis=("Re z",(-1.1,1.1)),yaxis=("Im z",(-1.1,1.1)))
    title!("Complex exponential = unit circle")
end
gif(an, "cpxexp-circle.gif")

2.3.3 Polar form

Suppose \(z\) is any nonzero complex value. We can write it in the form

\[ z = |z|\cdot \frac{z}{|z|}, \]

where \(|z|\) is the modulus (distance from origin in the plane). Since \(\frac{z}{|z|}\) has modulus equal to 1, it must lie on the unit circle. Hence there is a real \(\theta\) such that

\[ z = |z|\, e^{i\theta}. \tag{2.10}\]

We call Equation 2.10 the polar form of a complex number, because it expresses \(z\) as a distance from zero and an angle from the positive Re axis. Just as with any point in the plane, we can express a complex number either in Cartesian form using Re and Im parts, or in polar form using modulus and Euler’s identity.

Example 2.9 Suppose \(r\) is a positive real number. Then \(-r\) lies at a distance \(r\) from the origin along the negative real axis. Hence

\[ -r = |r|\, e^{i\pi}. \]

Supposing that we may take the log of both sides, we get

\[ \ln(-r) = \ln |r| + i\pi. \]

Using complex numbers, then, we can take the log of a negative number. You will find that this is the case in MATLAB.

2.3.4 Complex exponentials

Exponential functions still obey the properties you already know, even when the exponents are imaginary or complex numbers. This allows us to handle the exponential of any complex number. Writing \(a+i \omega\) for real \(a\) and \(\omega\), we have the function

Theorem 2.7 (Complex exponential function) \[ e^{(a+i \omega)t} = e^{at} \cdot e^{i\omega t} = e^{at} \bigl[ \cos(\omega t) + i \sin(\omega t)\bigr]. \tag{2.11}\]

In the function \(e^{\lambda t}\), the real part of \(\lambda\) controls growth or decay in time, as we are familiar with, and the imaginary part of \(\lambda\) controls a frequency of oscillation around a circle in the complex plane.

2.3.4.1 Growing

If \(\text{Re} \lambda > 0\), the magnitude of the function grows exponentially. The result is an outward spiral. The imaginary part of \(\lambda\) controls the frequency or pitch of the spiral. The real and imaginary parts of \(e^{\lambda t}\) oscillate between two growing real exponentials.

Code 
a = 0.025;  ω = 0.75;
t = range(0, 25, length=120)
z = @. exp((a+1im*ω)*t)

an = @animate for k in [eachindex(t);fill(length(t),5)]
    plot(t[1:k],real(z[1:k]),imag(z[1:k]),xaxis=("t",(0,30),:flip),yaxis=("Re z",(-2,2)),zaxis=("Im z",(-2,2)))
    plot!(t[1:k],real(z[1:k]),fill(-1.5,k))
    plot!(t[1:k],fill(2,k),imag(z[1:k]))
    plot!(fill(30,k),real(z[1:k]),imag(z[1:k]),l=:black)
    title!(@sprintf("Growth and oscillation: a = %.2f, ω = %.2f",a,ω))
end
gif(an, "cexp-growing.gif")

2.3.4.2 Neutral

When \(\lambda\) is purely imaginary, the function values stay on the unit circle in the complex plane. Taking the real and imaginary parts of \(e^{\lambda t}\) gives cosine and sine, respectively.

Code 
a = 0.;  ω = 0.75;
t = range(0, 25, length=120)
z = @. exp((a+1im*ω)*t)

an = @animate for k in [eachindex(t);fill(length(t),5)]
    plot(t[1:k],real(z[1:k]),imag(z[1:k]),xaxis=("t",(0,30),:flip),yaxis=("Re z",(-2,2)),zaxis=("Im z",(-2,2)))
    plot!(t[1:k],real(z[1:k]),fill(-1.5,k))
    plot!(t[1:k],fill(2,k),imag(z[1:k]))
    plot!(fill(30,k),real(z[1:k]),imag(z[1:k]),l=:black)
    title!(@sprintf("Pure oscillation: a = %.2f, ω = %.2f",a,ω))
end
gif(an, "cxep-neutral.gif")

2.3.4.3 Decaying

Finally, if \(\text{Re} \lambda < 0\), the spiral is a decaying one. The real and imaginary parts are attenuated oscillations.

Code 
a = -0.07;  om = 2.4;
t = range(0, 25, length=120)
z = @. exp((a+1im*ω)*t)

an = @animate for k in [eachindex(t);fill(length(t),5)]
    plot(t[1:k],real(z[1:k]),imag(z[1:k]),xaxis=("t",(0,30),:flip),yaxis=("Re z",(-2,2)),zaxis=("Im z",(-2,2)))
    plot!(t[1:k],real(z[1:k]),fill(-1.5,k))
    plot!(t[1:k],fill(2,k),imag(z[1:k]))
    plot!(fill(30,k),real(z[1:k]),imag(z[1:k]),l=:black)
    title!(@sprintf("Decay and oscillation: a = %.2f, ω = %.2f",a,ω))
end
gif(an, "cexp-decay.gif")

2.3.5 Homogeneous solutions

Now we can finish the story of computing solutions to homogeneous constant-coefficient ODEs. We can still use Theorem 2.3 and Theorem 2.4, now with the complex-valued interpretations of the exponentials.

Example 2.10 Solve the IVP \(x''+9x=0\), \(x(0)=2\), \(x'(0)=-12\).

Example 2.11 Solve the IVP \(x''+9x=0\), \(x(0)=2\), \(x'(0)=-12\).

Solution. The characteristic polynomial is \(s^2+9\), giving the roots \(\pm 3i\). Hence the general solution is \(x_h(t) = c_1e^{3it} + c_2e^{-3it}\).

The initial conditions require

\[\begin{split} 2 &= x(0) = c_1e^0 + c_2e^0 = c_1 +c_2,\\ -12 &= x'(0) = 3i c_1e^0 - 3i c_2e^0 = 3i c_1 - 3i c_2. \end{split}\]

This system is easy to solve for \(c_1\) and \(c_2\) to get \(c_1=1+2i\), \(c_2=1-2i\).

There are some helpful nuances to point out about the preceding example.

Theorem 2.8 If a real second-order IVP has complex conjugate characteristic values, then the integration constants satisfy \(c_2=\overline{c_1}\).

This can simplify the algebra a bit. If we set \(c_1=\alpha + i\beta\), then

\[ c_1 + c_2 = 2\alpha, \quad c_1-c_2 = 2i\beta. \tag{2.12}\]

In Example 2.11 this would have led us right to \(2\alpha=2\) and \(3i(2i\beta)=-12\).

Tip

When you have a complex expression “stuff \(+\) conjugate of that stuff”, the result has to be real, so you needn’t bother computing the imaginary part of that stuff.

Example 2.12 Solve the IVP \(x'' + 4x' + 5x=0\), \(x(0)=6\), \(x'(0)=0\).

Solution. The roots of \(s^2+4s+5\) are \(-2\pm i\). Thus the general solution is

\[ x(t) = c_1 e^{(-2+i)t} + \overline{c_1}e^{(-2-i)t}. \]

Note that

\[ x'(t) = c_1 (-2+i) e^{(-2+i)t} + \overline{c_1}(-2-i) e^{(-2-i)t}. \]

Writing \(c_1=\alpha + i\beta\), the initial condition on \(x\) requires that

\[ 6 = x(0) = \alpha + i\beta + \alpha - i\beta = 2\alpha, \]

so \(\alpha=3\). The other initial condition leads to

\[ \begin{split} 0 = x'(0) &= (3+i\beta) (-2+i)\cdot 1 + (3-i\beta) (-2-i)\cdot 1 \\ &= (-6-2i+3i-\beta) + (-6+2i-3i-\beta) = -12 - 2\beta, \end{split} \]

which yields \(\beta=-6\).

2.3.6 Fully real version

It might seem odd to use complex numbers to represent what we know must be a real-valued solution. There’s nothing wrong with doing so, but we can also convert it to an explicitly real form. The following is easily proved using Euler’s identity.

Theorem 2.9 The general homogeneous solution

\[ x_h(t) = c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t}, \]

where \(\lambda_{1,2} = a \pm i\omega\) are complex, is equivalent to the real expression

\[ x_h(t) = e^{at} \bigl[ b_1 \cos(\omega t) + b_2 \sin(\omega t) \bigr], \]

where \(b_1,b_2\) are real constants.

Example 2.13 In Example 2.11, the roots were \(0 \pm 3i\). Thus another expression for the general solution is

\[ x_h(t) = c_1\cos(3t) + c_2\sin(3t). \]

The initial conditions now yield

\[ \begin{split} 2 &= x(0) = c_1\cos(0) + c_2\sin(0) = c_1,\\ -12 &= x'(0) = -3 c_1 \sin(0) + 3 c_2 \cos(0) = 3c_2. \end{split} \]

Hence the IVP solution is \(2\cos(3t) - 4\sin(3t).\)

Example 2.14 Solve the IVP of Example 2.12 without using complex numbers to express the solution.

Solution. The roots of \(s^2+4s+5\) are \(-2\pm i\). Thus the general solution is

\[ x(t) = b_1 e^{-2t}\cos(t) + b_2 e^{-2t}\sin(t). \]

Note that

\[ x'(t) = b_1[ -2e^{-2t}\cos(t) - e^{-2t}\sin(t)] + b_2[ -2e^{-2t}\sin(t) + e^{-2t}\cos(t)]. \]

The initial condition on \(x\) requires that

\[ 6 = x(0) = b_1, \]

so \(b_1=6\). The other initial condition leads to

\[ 0 = x'(0) = 6(-2+0) + b_2(0+1) = -12 + b2, \]

which yields \(b_2=12\).

Since complex exponentials are much less familiar to you than sin and cos, the real form might seem more appealing to you. If you acquaint yourself with complex exponentials, though, there is nothing wrong with leaving things in that form, and you have less to memorize.

2.4 Variation of parameters

Section 2.2 and Section 2.3 explained how to find the general solution \(x_h(t)\) of a homogeneous linear system \(\opA[x]=0\). We need to perform that same step for a nonhomogeneous problem \(\opA[x]=f\), but then also find a particular solution of the original equation. One possibility is to adapt the variation of parameters technique used for first-order problems in Section 1.7.

The form of the homogeneous solution is

\[ x_h(t) = c_1 x_1(t) + c_2 x_2(t), \]

where \(x_1\) and \(x_2\) are independent homogeneous solutions. Now we replace the constants \(c_1\) and \(c_2\) with functions of \(t\):

\[ x_p(t) = u_1(t) x_1(t) + u_2(t) x_2(t). \tag{2.13}\]

It follows that

\[ x_p' = u_1'x_1 + u_1 x_1' + u_2'x_2 + u_2 x_2'. \]

We are trying to find a way to choose \(u_1\) and \(u_2\) that makes \(x_p\) a particular solution. Our lives will be a lot easier if we start with the constraint

\[ u_1'x_1 + u_2'x_2 = 0, \tag{2.14}\]

which implies \(x_p' = u_1 x_1' + u_2 x_2'\). Let’s see if we can get away with that. We differentiate again to get

\[ x_p'' = u_1'x_1' + u_1 x_1'' + u_2'x_2' + u_2 x_2''. \]

Since \(\opA[x] = x'' + bx' +k\), we get

\[ \opA[x_p] = u_1 (x_1'' + bx_1' + kx_1) + u_2 (x_2'' + bx_2' + kx_2) + ( u_1'x_1' + u_2'x_2'). \]

We are in great shape now. Since \(x_1\) and \(x_2\) are homogeneous solutions, we just have \(\opA[x_p] = u_1'x_1' + u_2'x_2'\). So we need only to introduce the additional requirement

\[ u_1'x_1' + u_2'x_2' = f, \tag{2.15}\]

and then \(x_p\) will be the desired particular solution.

Equation 2.14 and Equation 2.15 impose

\[ \begin{split} x_1 u_1' + x_2 u_2' &= 0, \\ x_1'u_1' + x_2'u_2' &= f. \end{split} \]

Writing \(u_1'=v_1\) and \(u_2'=v_2\), we see that this is a pair of linear equations for \(v_1\) and \(v_2\), which is routine to solve, leading to the following.

Theorem 2.10 (Variation of parameters, 2nd order) Let \(x_1\) and \(x_2\) be independent solutions of the homogeneous equation \(\opA[x]=0\), and let \(W(x_1,x_2)\) be their Wronskian. For any given forcing function \(f\), define

\[ u_1 = \int \frac{-x_2 f}{W(x_1,x_2)}\, dt, \quad u_2 = \int \frac{x_1 f}{W(x_1,x_2)}\, dt. \tag{2.16}\]

Then a particular solution of \(\opA[x]=f\) is

\[ x_p(t) = u_1(t) x_1(t) + u_2(t) x_2(t). \]

Important

The integration constants for \(u_1\) and \(u_2\) in Equation 2.16 can be taken to be zero.

Example 2.15 Find a particular solution of \(x'' - x = 4\).

Solution. There’s no getting around the need for the homogeneous solution here. The characteristic roots are \(\pm 1\), so we take \(x_1=e^t\), \(x_2=e^{-t}\). Then

\[ W(e^t, e^{-t}) =(e^t)(-e^{-t}) - (e^{-t})(e^t) = 2, \]

and

\[ u_1' = \frac{-4 e^{-t}}{-2}, \quad u_2' = \frac{4 e^{t}}{-2}. \]

From here we find \(u_1 = -2e^{-t}\) and \(u_2 = -2e^t\). Hence

\[ x_p = (-2e^{-t})(e^t) + (-2e^t)(e^{-t}) = -4. \]

You can instantly check that this result is correct.

Example 2.16 Find a particular solution of \(x'' - \lambda^2 x = 2 e^{\lambda t}\), where \(\lambda\neq 0\) is a constant parameter.

Solution. The homogeneous equation has roots \(\pm \lambda\). We choose \(x_1=e^{\lambda t}\), \(x_2=e^{-\lambda t}\). Now

\[ W(x_1, x_2) = \twodet{e^{\lambda t}}{e^{-\lambda t}}{\lambda e^{\lambda t}}{-\lambda e^{-\lambda t}} = -2\lambda. \]

This gives

\[ u_1' = \frac{ - 2 e^{-\lambda t} e^{\lambda t}}{-2\lambda}, \quad u_2' = \frac{ 2 e^{\lambda t} e^{\lambda t} }{-2\lambda}. \]

Thus,

\[ u_1 = \frac{1}{\lambda} \int dt = \frac{t}{\lambda}, \quad u_2 = -\frac{1}{\lambda} \int e^{2 \lambda t}\, dt = -\frac{1}{2 \lambda^2} e^{2 \lambda t}. \]

Finally,

\[ x_p(t) = \frac{1}{\lambda} t e^{\lambda t} - \frac{1}{2 \lambda^2} e^{\lambda t} \]

is a particular solution. Our recipe now gives us the general solution

\[ x(t) = x_h(t) + x_p(t) = c_1 e^{\lambda t} + c_2e^{-\lambda t} + \frac{1}{\lambda} t e^{\lambda t} - \frac{1}{2 \lambda^2} e^{\lambda t}. \]

However, that last term is redundant—all it does is change the value of \(c_1\). So it’s just as correct to write

\[ x(t) = c_1 e^{\lambda t} + c_2e^{-\lambda t} + \frac{1}{\lambda} t e^{\lambda t}. \]

As you can see, the VoP formula is pretty great if you happen to be a robot. [Note to editor: insert CAPTCHA here.] Us humans, though, could use something with a bit less algebra to manage, and that is coming next.

2.5 Undetermined coefficients

When it comes to finding particular solutions of nonhomogeneous problems, there is a simpler alternative to variation of parameters in the important special case when

  1. The coefficients in the ODE are constants, and
  2. The forcing function \(f(t)\) is a polynomial, exponential, sin, or cos, or a combination of these.

We have already been assuming point 1 above for a while. Now we consider

\[ x'' + bx' + kx = f(t), \tag{2.17}\]

where \(b\) and \(k\) are constants, and \(f\) appears in Table 2.1.

Table 2.1: Particular solutions for the method of undetermined coefficients
\(f(t)\) \(x_p(t)\)
\(b_n t^n + \cdots b_0\) \(B_n t^n + \cdots + B_0\)
\(e^{rt}(b_n t^n + \cdots b_0)\) \(e^{rt}(B_n t^n + \cdots B_0)\)
\(\cos(\omega t)\) \(A \cos(\omega t) + B \sin(\omega t)\)
\(\sin(\omega t)\) \(A \cos(\omega t) + B \sin(\omega t)\)

Corresponding to each \(f\) in Table 2.1 is a specific form for \(x_p\) containing one or more constants (given by capital letters) that must be deduced by substituting into the particular ODE under consideration. This procedure is what gives the method its name, the method of undetermined coefficients.

Example 2.17 Rework Example 2.15, finding a particular solution of \(x'' - x = 4\).

Solution. The forcing function is a polynomial of degree zero, so \(x_p\) should be as well. Let \(x_p=A.\) Then it is a solution if

\[ 4 = A'' - A = -A. \]

Hence \(A=-4\).

Caution

In the case where \(f\) is a polynomial that may have some zero coefficients, you cannot assume that \(x_p\) has zero coefficients in the same positions.

Example 2.18 Find a particular solution of \(x'' +4x'+4x=8t^2\).

Solution. The correct form of \(x_p\) is a quadratic polynomial, i.e.,

\[ x_p(t) = A + Bt + Ct^2. \]

Plugging that into the ODE yields

\[\begin{split} 2C + 4(B+2Ct) + 4(A+Bt+Ct^2) & = 8t^2,\\ 4C t^2 + (8C+4B) + (2C+4B+4A) & = 8t^2. \end{split}\]

This has to be an identity for all \(t\). Matching like powers, we get a \(3\times 3\) linear system for the coefficients. Fortunately, it can be solved easily if we go in order: first we read off \(C=2\), then \(B=-2C=-4\), and finally \(A=-B-C/2=3\). This provides us with

\[ x_p(t) =2t^2+4t+3. \]

Example 2.19 Find a particular solution of \(x'' - 2x'-3x=10e^{4t}\).

Solution. The proper choice is

\[ x_p(t) = Ae^{4t}. \]

Everything else is algebra.

\[\begin{split} 16A e^{4t} - 2 (4A e^{4t}) - 3 (Ae^{4t}) & =10 e^{4t},\\ 16A -8A -3A &= 10. \end{split}\]

.

From this it’s clear that \(A=2\).

Example 2.20 Find the general solution of \(x''+x'=\sin(2t)\).

Solution. The homogeneous problem \(x''+x'=0\) has roots \(0\) and \(-1\), hence

\[ x_h(t) = c_1 + c_2 e^{-t}. \]

For the particular solution we must choose

\[ x_p(t) = A\cos(2t) + B\sin(2t). \]

Inserting this into the original ODE leads to

\[ [-4A\cos(2t) - 4B\sin(2t)] + [-2A\sin(2t) + 2B\cos(2t)] = \sin(2t). \]

This is identically true for all time if and only if we match the coefficients of like trig functions:

\[ -4A + 2B = 0, \qquad -4B-2A = 1. \]

The solution of these equations is \(A=-1/10\), \(B=-1/5\). Thus

\[ x(t) = c_1 + c_2 e^{-t} - \frac{1}{10}\cos(2t) - \frac{1}{5}\sin(2t). \]

The examples above are the fundamental ones. If you have a forcing function that is a sum of terms from Table 2.1, then you can treat each term separately as above and add together the particular solutions found.

2.5.1 Limitations

There are rules for more complicated combinations of the same functions, but we won’t go into them. At some point the calculus/algebra becomes pretty intense. Furthermore, we are not going to cover a refinement of the rules necessary in certain circumstances where our simplified version fails.

Example 2.21 The equation \(x''+x=\cos(\omega t)\) suggests the particular solution \(x_p(t)=A\cos(\omega t)+B\sin(\omega t)\). Upon substitution,

\[ [-\omega^2 A\cos(\omega t) - \omega^2 B\sin(\omega t) ] + [ A\cos(\omega t) + B\sin(\omega t)] =\cos(\omega t), \]

which leads to the conclusion that \(B=0\) and, if \(\omega^2 \neq 1\), \(A=1/(1-\omega^2)\). However, if \(\omega = 1\), the substitution would leave us with \(0=\cos(t)\), which is impossible to satisfy for all \(t\).

The failure of Example 2.21 at \(\omega = 1\) was due to the fact that the \(x_p\) we tried is actually a homogeneous solution. This is not too hard to fix, but a general description is not worth the effort; we can always fall back to variation of parameters.

Example 2.22 We can make the tricky case \(x''+x=\cos(t)\) from Example 2.21 quite easy by relating it to complex exponentials. Since \(\cos(t)\) is the real part of \(e^{it}\), it follows that the solution we seek is the real part of any particular solution of

\[ x'' + x = e^{it}. \]

Referring to Example 2.16 with \(\lambda=i\), we use variation of parameters to find the particular solution

\[ \frac{1}{2i} t e^{i t} = -\frac{i}{2}t [\cos(t) + i\sin (t)]. \]

Upon taking the real part, we get

\[ x_p(t) = \frac{1}{2} t \sin(t). \]

2.6 Oscillators

Code 
using Plots,LaTeXStrings,Printf
default(linewidth=2,label="")

The distinguishing feature of second-order ODEs compared to the first-order case is that they permit oscillations as well as exponential growth and decay. These equations appear in models throughout engineering and science as a result.

For instance, an object attached to an ideal spring satisfies the ODE

\[ mx'' + b x' + kx = 0, \tag{2.18}\]

where \(x(t)\) is the displacement from the natural resting position, \(m\) is the mass of the object, \(b\) is an attenuation coefficient due to friction or mechanical damping, and \(k\) is a property of the spring called the spring constant. This constant is the ratio of the restoring force of the spring to the amount by which it is stretched or compressed.

All named parameters such as \(m\), \(k\), and \(b\) are assumed to be nonnegative unless stated otherwise.

Note

Equation 2.18 applies equally well to horizontal and vertical oscillations. Gravity does not appear in the vertical case because it is accounted for by measuring \(x\) from the static equilibrium position.

If given, initial values for Equation 2.18 supply the initial position and initial velocity, which determine the solution uniquely.

A pendulum is another type of mechanical oscillator. The proper ODE is

\[ \theta'' + \gamma \theta' + \frac{g}{L}\sin (\theta) = 0, \]

where \(\theta(t)\) is the angle made with vertically-down position, \(L\) is the arm length, and \(g\) is gravitational acceleration. This equation is nonlinear and difficult to analyze without developing new tools. But if the angle remains small, then a reasonable approximation is

\[ \theta'' + \gamma \theta' + \frac{g}{L}\theta = 0, \]

which is a linear oscillator with constant coefficients.

An AC circuit typically has elements of resistance, capacitance, and inductance. These analogize perfectly to friction/damping, spring constant, and mass. If these elements are wired in series, the governing ODE is

\[ LI'' + RI' + \frac{1}{C}I = 0, \]

where \(I(t)\) is the current flowing through the circuit, \(L\) is inductance, \(R\) is resistance, and \(C\) is capacitance.

These homogeneous models are examples of free oscillations. Of course, any of the models above might come with an external forcing function. For a pendulum, this could be someone pushing you on a swing, or for a circuit, it could be an AC generator. We call these forced oscillations.

2.6.1 Unifying notation

When you have many versions of the same fundamental problem, each using different symbols and units, you have three options.

  1. Solve each new problem from scratch.
  2. Derive custom formulas for each application.
  3. Find a minimal set of parameters and express the problem and solution in terms of them.

Option 1 is highly inefficient. Option 2 would be appropriate for an engineering course. Here, we take option 3 and develop generalized knowledge that can be reinterpreted for each new application.

The ODE of a free mass-spring oscillator, for example, has standard form

\[ x'' + \frac{b}{m} x' + \frac{k}{m} x = 0. \]

This suggests that we only need two parameters, not three, to express the full range of behavior. (An added convenience is that both \(b/m\) and \(\sqrt{k/m}\) have units of inverse time.) Accordingly, we introduce

\[ \omega_0 = \sqrt{\frac{k}{m}}, \qquad \zeta = \frac{b/m}{2\omega_0} = \frac{b}{\sqrt{4km}}. \]

The parameter \(\omega_0\) is known as the natural frequency, with units of inverse time, and \(\zeta\) is a dimensionless damping coefficient describing the relative intensity of the damping.

Important

In math we usually use frequency to mean the multiplier of \(t\) in a sin or cos function. That is our usage. In some fields this is called angular frequency, and frequency is used to mean the number of cycles per time unit, as in Hz.

If we include a forcing term, we arrive at the ODE

\[ x'' + 2 \zeta \omega_0\, x' + \omega_0^2\, x = f(t). \tag{2.19}\]

A similar derivation can be done starting from the pendulum or AC circuit equations. The only type of forcing we shall consider is harmonic forcing of the form

\[ f(t) = \cos(\omega t). \]

2.6.2 Classifications

The characteristic roots of this ODE are

\[ \lambda_{1,2} = -\zeta \omega_0 \pm \omega_0 \sqrt{\zeta^2-1}. \]

The discussion now splits into four cases, marked by increasing values of \(\zeta\).

2.6.2.1 Undamped oscillator, \(\zeta=0\)

The oscillator \(x'' + \omega_0^2 x\) has homogeneous solutions that can be expressed either as combinations of

\[ e^{i\omega_0 t},\, e^{-i\omega_0 t} \]

or of

\[ \cos(\omega_0 t),\, \sin(\omega_0 t). \]

Either way, these describe pure oscillation at frequency \(\omega_0\). This is known as simple harmonic motion. A particularly useful third form is the amplitude–phase form

\[ x_h(t) = R\cos(\omega_0 t+\phi), \tag{2.20}\]

where constants \(R\) and \(\phi\) can be used to satisfy initial conditions.

Example 2.23 When a 2 kg mass is hung on a spring, the spring stretches by 0.25 m. What is the natural frequency of the mass-spring system? Suppose the mass is pulled down 0.2 m past equilibrium and then thrown upward at 1 m/s. What is the amplitude of the motion?

Solution. Hooke’s Law for a spring states that \(F=k x\), so we find the spring constant from \(k=F/x=2g/0.25=8g\), where \(g=9.8\) m per second squared. The ODE for free motion of the system is thus \(2 x'' + 8g x = 0\), or

\[ x'' + 4g x = 0. \]

From this we identify the natural frequency

\[ \omega_0 = \sqrt{4g} = 2\sqrt{g} \approx 6.26 \text{s}^{-1}. \]

We can apply the initial conditions directly to the phase–amplitude form:

\[ \begin{split} -0.2 & = x(0) = R\cos(\phi) \\ 1 & = x'(0) = \omega_0 R\sin(\phi). \end{split} \]

Therefore,

\[ R^2 = [R\cos(\phi)]^2 + [R\sin(\phi)]^2 = 0.04 + \omega_0^{-2}, \]

which works out to \(R \approx 0.256 \text{m}\).

The particular solution for harmonic forcing \(f(t)=\cos(\omega t)\) is

\[ x_p(t) = \frac{1}{\omega_0^2-\omega^2} \cos(\omega t), \]

provided \(\omega\neq \omega_0\). Note that the harmonic amplitude grows without bound as \(\omega\to \omega_0\).

If the forcing is exactly at the natural frequency, then the situation is like that of Example 2.16:

\[ x_p(t) = \frac{1}{2 \omega_0} t \sin(\omega_0 t), \]

which is unbounded as \(t\to \infty\). This is known as resonance, one of the most important phenomena in physics.

2.6.2.2 Underdamped oscillator, \(0 < \zeta < 1\)

For \(0< \zeta < 1\) the roots of Equation 2.19 are complex, with negative real part:

\[ \lambda_{1,2} = -\zeta \omega_0 \pm i \omega_0 \sqrt{1-\zeta^2}. \]

Define the damped frequency

\[ \omega_d=\omega_0 \sqrt{1-\zeta^2}. \tag{2.21}\]

The unforced part of the solution is

\[ x_h(t) = e^{- \omega_0 \zeta t} [ c_1 \cos( \omega_d t) + c_2 \sin(\omega_d t) ]. \tag{2.22}\]

This solution is pseudoperiodic, combining oscillation at frequency \(\omega_d < \omega_0\) inside an exponential decay envelope. We call this an underdamped oscillator. The homogeneous solution Equation 2.22 is also called a transient solution, because it vanishes as \(t \to \infty\).

If harmonic forcing \(f(t)=\cos(\omega t)\) is added, then it determines the steady-state behavior \(x_s=A\cos(\omega t + \theta)\), where

\[ A = \frac{1}{\rule{0pt}{1em} \sqrt{\rule{0pt}{0.8em}(\omega_0^2-\omega^2)^2 + 4\omega_0^2\omega^2 \zeta^2 } } \tag{2.23}\]

The value \(A\) is known as the gain of the oscillator; it is the ratio of amplitudes of the steady response and the forcing. A little calculus shows that as a function of \(\omega\), the gain is maximized at

\[ \omega_\text{max} = \begin{cases} \omega_0 \sqrt{1-2\zeta^2}, & \text{ if } 0 < \zeta^2 < \frac{1}{2},\\ 0, & \text{ if } \frac{1}{2} \le \zeta^2. \end{cases} \tag{2.24}\]

Therefore, if the damping coefficient \(\zeta\) is small but finite, we have a pseudoresonance of finite amplitude at a frequency just a bit less than the natural frequency. The following figure shows the gain as a function of driving frequency and damping when \(\omega_0=1\). The black curve shows the maximal driving frequency at any given \(\zeta\).

Code 
ω = range(1e-3,2,length=90)
ζ = range(1e-3,1.2,length=90)
A = [ 1/(1-ω^2 + 2im*z*ω ) for ω in ω, z in ζ ]
log_g = @. log10(abs(A))
contour(ω,ζ,log_g',levels=20,l=(:bluesreds,2),clims=(-1.5,1.5))
ζ = range(1e-3,1/sqrt(2),length=100)
ωmax = @. real( sqrt(1+0im-2ζ^2) ) 
plot!(ωmax,ζ,l=(:black),
  xlabel="driving frequency ω", ylabel="damping coefficient ζ", title="Log gain for ω₀=1")

Here’s a short video showing how the gain of an oscillator varies with the forcing frequency \(\omega\). (I use the term “eigenvalues” to mean characteristic values here.)

2.6.2.3 Critically damped oscillator, \(\zeta = 1\)

At \(\zeta=1\) the complex roots coalesce into a double real root,

\[ \lambda_1 = \lambda_2 = -\omega_0, \]

with general homogeneous solution

\[ x_h(t) = e^{-\omega_0 t} (c_1 + c_2 t). \]

There is no longer any oscillation present, and we have a critically damped system. The linear growth of \(c_2 t\) doesn’t make much of a difference against the exponential decay, and \(x_h\to 0\) as \(t\to\infty\). Any steady response is due to the forcing term, but no resonance is possible.

2.6.2.4 Overdamped, \(\zeta >1\)

For \(\zeta > 1\) the roots are

\[ \lambda_{1,2} = -\omega_0 \zeta \pm \omega_0 \sqrt{\zeta^2-1}, \]

which are negative and real. This gives an exponentially decaying homogeneous solution. In this case we have an overdamped oscillator. Again, no resonance is possible.

2.6.3 Examples

Table 2.2: Damping coefficient and characteristic values
Damping coefficient Root property
\(\zeta=0\) imaginary
\(0 < \zeta < 1\) complex
\(\zeta=1\) real, negative, repeated
\(\zeta > 1\) real, negative, distinct

Example 2.24 A 5 kg mass is hung on a spring with constant \(11\) N per m and connected to a dashpot that provides 8 N-sec per meter of resistance. Is this system underdamped, overdamped, or critically damped?

Solution. The ODE for the mass-spring system is

\[ \begin{split} 5 x'' + 8x' + 11 x & = 0 ,\\ x'' + 1.6 x' + 2.2 x & = 0, \end{split} \]

from which we identify the natural frequency

\[ \omega_0 = \sqrt{2.2} \approx 1.483 \text{s}^{-1}. \]

The damping coefficient is therefore

\[ \zeta = \frac{1.6}{2\omega_0} \approx 0.539. \]

Since this value is less than one, the system is underdamped.

Example 2.25 Suppose the system from the previous example is initially at equilibrium when the mass is suddenly pushed downward at 0.5 m/sec. Find the motion of the mass.

Solution. We derived the governing ODE \(x'' + 1.6 x' + 2.2 x = 0\). The roots are the roots of \(\lambda^2 + 1.6\lambda + 2.2\), which are found numerically to be

\[ \lambda \approx -0.8000 \pm 1.2490i. \]

(The imaginary part is smaller than the natural frequency found in the last example, as it must be.) Choosing the sin-cos form of the general solution, we have

\[ x_h(t) = c_1 e^{-0.8 t} \cos(1.249 t) + c_2 e^{-0.8 t} \sin(1.249 t). \]

We apply the initial conditions \(x(0)=0\), \(x'(0)=-0.5\) to find

\[ \begin{split} 0 & = x_h(0) = c_1, \\ -0.5 & = x_h'(0) = c_1( -0.8 ) + c_2 (1.249 ), \end{split} \]

thus \(c_2 = -0.4003\). The motion is therefore given by \(x(t)=-0.4003\, e^{-0.8 t} \sin(1.249 t)\).