Analytical Solutions to 2nd-order ODEs

Next, we’ll focus on analytical solutions for 2nd-order ODEs, and specifically initial-value problems (IVPs). As before, let’s categorize problems based on their solution approach.

1. Solution by direct integration

If you have a 2nd-order ODE of this form:

\begin{equation} \frac{d^2 y}{dx^2} = f(x) \end{equation}

then you can solve by direct integration.

For example, let’s say we are trying to solve for the deflection of a cantilever beam \(y(x)\) with a force \(P\) at the end, where \(E\) is the modulus and \(I\) is the moment of intertia, and the initial conditions are \(y(0)=0\) and \(y^{\prime}(0) = 0\):

\begin{align} \frac{d^2 y}{dx^2} &= \frac{-P (L-x)}{EI} \\ \frac{d}{dx} \left(\frac{dy}{dx}\right) &= \frac{-P}{EI} (L-x) \\ \int d \left(\frac{dy}{dx}\right) &= \frac{-P}{EI} \int (L-x) dx \\ y^{\prime} = \frac{dy}{dx} &= \frac{-P}{EI} \left(Lx - \frac{x^2}{2}\right) + C_1 \\ \int dy &= \int \left( \frac{-P}{EI} \left(Lx - \frac{x^2}{2}\right) + C_1 \right) dx \\ y(x) &= \frac{-P}{EI} \left(\frac{L}{2} x^2 - \frac{1}{6} x^3\right) + C_1 x + C_2 \end{align}

That is our general solution; we can obtain the specific solution by applying our two initial conditions:

\begin{align} y(0) &= 0 = C_2 \\ y^{\prime}(0) &= 0 = C_1 \\ \therefore y(x) &= \frac{P}{EI} \left( \frac{x^3}{6} - \frac{L x^2}{2} \right) \end{align}

2. Solution by substitution

If we have a 2nd-order ODE of this form:

\begin{equation} \frac{d^2 y}{dx^2} = f(x, y^{\prime}) \end{equation}

then we can solve by substitution, meaning by substituting a new variable for \(y^{\prime}\). (Notice that \(y\) itself does not show up in the ODE.)

Let’s substitute \(u\) for \(y^{\prime}\) in the above ODE:

\begin{align} u &= y^{\prime} \\ u^{\prime} &= y^{\prime\prime} \\ \rightarrow u^{\prime} &= f(f, u) \end{align}

Now we have a 1st-order ODE! Then, we can apply the methods previously discussed to solve this; once we find \(u(x)\), we can integrate that once more to get \(y(x)\).

Example: falling object

For example, consider a falling mass where we are solving for the downward distance as a function of time, \(y(t)\), that is experiencing the force of gravity downward and a drag force upward. It starts at some reference point so \(y(0) = 0\), and has a zero initial downward velocity: \(y^{\prime}(0) = 0\). The governing equation is:

\begin{equation} m \frac{d^2 y}{dt^2} = mg - c \left( \frac{dy}{dt} \right)^2 \end{equation}

where \(m\) is the mass, \(g\) is acceleration due to gravity, and \(c\) is the drag proportionality constant. We can substitute \(V\) for \(y^{\prime}\), which gives us a first-order ODE:

\begin{align} \text{let} \quad \frac{dy}{dt} = V \\ \rightarrow m \frac{dV}{dt} &= mg - c V^2 \\ \end{align}

Then, we can solve this for \(V(t)\) using our initial condition for velocity \(V(0) = 0\). Once we have that, we can integrate once more:

\begin{equation} y(t) = \int V(t) dt \end{equation}

and apply our initial condition for position, \(y(0) = 0\), to obtain \(y(t)\).

Here is the full process:

\begin{align} \frac{dV}{dt} &= g - \frac{c}{m} V^2 \\ \frac{dV}{g - \frac{c}{m} V^2} &= dt \\ \frac{m}{c} \int \frac{dV}{a^2 - V^2} &= \int dt = t + \bar{c}, \quad \text{where} \quad a = \sqrt{\frac{mg}{c}} \\ \frac{m}{c} \frac{1}{a} \tanh^{-1} \left(\frac{V}{a}\right) &= t + c_1 \\ V &= a \tanh \left( \frac{a c}{m} t + c_1 \right) \\ \therefore V(t) &= \sqrt{\frac{mg}{c}} \tanh \left(\sqrt{\frac{gc}{m}} t + c_1\right) \end{align}

Applying the initial condition for velocity, \(V(0) = 0\):

\begin{align} V(0) &= 0 = \sqrt{\frac{mg}{c}} \tanh \left(0 + c_1\right) \\ \therefore c_1 &= 0 \\ V(t) &= \sqrt{\frac{mg}{c}} \tanh \left(\sqrt{\frac{gc}{m}} t\right) \end{align}

Then, to get \(y(t)\), we just need to integrate once more:

\begin{align} \frac{dy}{dt} = V(t) &= \sqrt{\frac{mg}{c}} \tanh \left(\sqrt{\frac{gc}{m}} t\right) \\ \int dy &= \sqrt{\frac{mg}{c}} \int \tanh \left(\sqrt{\frac{gc}{m}} t\right) dt \\ y(t) &= \sqrt{\frac{mg}{c}} \sqrt{\frac{m}{gc}} \log\left(\cosh\left(\sqrt{\frac{gc}{m}} t\right)\right) + c_2 \\ \rightarrow y(t) &= \frac{m}{c} \log\left(\cosh\left(\sqrt{\frac{gc}{m}} t\right)\right) + c_2 \end{align}

Finally, we can apply the initial condition for position, \(y(0) = 0\), to get our solution:

\begin{align} y(0) &= 0 = \frac{m}{c} \log\left(\cosh\left(0\right)\right) + c_2 = c_2 \\ \rightarrow c_2 &= 0 \\ y(t) &= \frac{m}{c} \log\left(\cosh\left(\sqrt{\frac{gc}{m}} t\right)\right) \end{align}

Example: catenary problem

The catenary problem describes the shape of a hanging chain or rope fixed between two points. (It was also a favorite of one of my professors, Joe Prahl, and I like to teach it as an example in his honor.) The downward displacement of the hanging string/chain/rope as a function of horizontal position, \(y(x)\), is governed by the equation

\begin{equation} y^{\prime\prime} = \sqrt{1 + (y^{\prime})^2} \end{equation}

This is actually a boundary value problem, with the boundary conditions for the displacement at one side \(y(0) = 0\) and that the slope is zero in the middle: \(\frac{dy}{dx}\left(\frac{L}{2}\right) = 0\). (Please note that I have skipped the derivation of the governing equation, and left some details out.)

We can solve this via substitution, by letting a new variable \(u = y^{\prime}\); then, \(u^{\prime} = \frac{du}{dx} = y^{\prime\prime}\). This gives is a first-order ODE, which we can integrate:

\begin{align} \frac{du}{dx} &= \sqrt{1 + u^2} \\ \int \frac{du}{\sqrt{1 + u^2}} &= \int dx \\ \sinh^{-1}(u) &= x + c_1, \quad \text{where } \sinh(x) = \frac{e^x - e^{-x}}{2} \\ u(x) &= \sinh(x + c_1) \end{align}

Then, we can integrate once again to get \(y(x)\):

\begin{align} \frac{dy}{dx} &= u(x) = \sinh(x + c_1) \\ \int dy &= \int \sinh(x + c_1) dx = \int \left(\sinh(x)\cosh(c_1) + \cosh(x)\sinh(c_1)\right)dx \\ y(x) &= \cosh(x)\cosh(c_1) + \sinh(x)\sinh(c_1) + c_2 \\ \rightarrow y(x) &= \cosh(x + c_1) + c_2 \end{align}

This is the general solution to the catenary problem, and applies to any boundary conditions.

For our specific case, we can apply the boundary conditions and find the particular solution, though it involves some algebra…:

\begin{align} y(0) &= 0 = \cosh(c_1) + c_2 \\ \frac{dy}{dx}\left(\frac{L}{2}\right) &= u(0) = \sinh \left(\frac{L}{2} + c_1\right) \\ \rightarrow c_1 &= -\frac{L}{2} \\ 0 &= \cosh \left( -\frac{L}{2} \right) + c_2 \\ \rightarrow c_2 &= -\cosh\left( -\frac{L}{2} \right) = -\cosh\left(\frac{L}{2} \right) \end{align}

So, the overall solution for the catenary problem with the given boundary conditions is

\begin{equation} y(x) = \cosh \left(x - \frac{L}{2}\right) - \cosh\left( \frac{L}{2} \right) \end{equation}

Let’s see what this looks like:

L = 1.0;
x = linspace(0, 1);
y = cosh(x - (L/2.)) - cosh(L/2.);
plot(x, y)

Please note that I’ve made some simplifications in the above work, and skipped the details of how the ODE is derived. In general, the solution for the shape is

\begin{equation} y(x) = C \cosh \frac{x + c_1}{C} + c_2 \end{equation}

where you would solve for the constants \(C\), \(c_1\), and \(c_2\) using the constraints:

\begin{align} \int_{x_a}^{x_b} \sqrt{1 + (y^{\prime})^2} dx &= L \\ y(x_a) &= y_a \\ y(x_b) &= y_b \;, \end{align}

where \(L\) is the length of the rope/chain.

You can read more about the catenary problem here (for example): http://euclid.trentu.ca/aejm/V4N1/Chatterjee.V4N1.pdf

3. Homogeneous 2nd-order ODEs

An important category of 2nd-order ODEs are those that look like

\begin{equation} y^{\prime\prime} + p(x) y^{\prime} + q(x) y = 0 \end{equation}

“Homogeneous” means that the ODE is unforced; that is, the right-hand side is zero.

Depending on what \(p(x)\) and \(q(x)\) look like, we have a few different solution approaches:

• constant coefficients: \(y^{\prime\prime} + a y^{\prime} + by = 0\)
• Euler-Cauchy equations: \(x^2 y^{\prime\prime} + axy^{\prime} + by = 0\)
• Series solutions

First, let’s talk about the characteristics of linear, homogeneous 2nd-order ODEs:

Solutions have two parts: \(y(x) = c_1 y_1 + c_2 y_2\), where \(y_1\) and \(y_2\) are each a basis of the solution.

The two parts of the solution \(y_1\) and \(y_2\) are linearly independent.

One way of defining this is that \(a_1 y_1 + a_2 y_2 = 0\) only has the trivial solution \(a_1=0\) and \(a_2=0\).

Another way of thinking about this is that \(y_1\) and \(y_2\) are linearly dependent if one is a multiple of the other, like \(y_1 = x\) and \(y_2 = 5x\). This cannot be solutions to a linear, homogeneous ODE.

\(y_1\) and \(y_2\) each satisfy the ODE. Meaning, you can plug each of them into the ODE for \(y\) and obtain 0.

However, we need both parts together to fully solve the ODE.

If \(y_1\) is known, we can get \(y_2\) by reduction of order. Let \(y_2 = u y_1\), where \(u\) is some unknown function of \(x\). Then, put \(y_2\) into the ODE \(y^{\prime}{\prime} + p(x) y^{\prime} + q(x) y = 0\):

\begin{align} y_2 &= u y_1 \\ y_2^{\prime} &= u y_1^{\prime} + u^{\prime} y_1 \\ y_2^{\prime\prime} &= 2 u^{\prime} y_1^{\prime} + u^{\prime\prime} y_1 + u y_1^{\prime\prime} \\ \rightarrow u^{\prime\prime} &= - \left[ p(x) + \left(\frac{2 y_1^{\prime}}{y_1}\right) \right] u^{\prime} \text{or, } u^{\prime\prime} &= - \left( g(x) \right) u^{\prime} \end{align}

Now, we have an ODE with only \(u^{\prime\prime}\), \(u^{\prime}\), and some function \(g(x)\)—so we can solve by substitution! Let \(u^{\prime} = v\), and then we have \(v^{\prime} = -g(x) v\):

\begin{align} \frac{dv}{dx} &= - \left( p(x) + \frac{2 y_1^{\prime}}{y_1} \right) v \\ \int \frac{dv}{v} &= - \int \left(p(x) + \frac{2 y_1^{\prime}}{y_1} \right) dx \\ \text{Recall } 2 \frac{d}{dx} \left( \ln y_1 \right) &= 2 \frac{y_1^{\prime}}{y_1} \\ \therefore \int \frac{dv}{v} &= - \int \left(p(x) + 2 \frac{d}{dx} \left( \ln y_1 \right) \right) dx \\ \ln v &= -\int p(x) dx - 2 \ln y_1 \\ \rightarrow v &= \frac{\exp\left( -\int p(x)dx \right)}{y_1^2} \end{align}

So, the actual solution procedure is then:

1. Solve for \(v\): \(v = \frac{\exp\left( -\int p(x)dx \right)}{y_1^2}\)
2. Solve for \(u\): \(u = \int v dx\)
3. Get \(y_2\): \(y_2 = u y_1\)

Here’s an example, where we know one part of the solution \(y_1 = e^{-x}\):

\begin{align} y^{\prime\prime} + 2 y^{\prime} + y &= 0 \\ \text{Step 1:} \quad v = \frac{\exp \left( -\int 2dx \right)}{ \left(e^{-x}\right)^2} = \frac{e^{-2x}}{e^{-2x}} &= 1 \\ \text{Step 2:} \quad u = \int v dx = \int 1 dx &= x \\ \text{Step 3:} \quad y_2 &= x e^{-x} \end{align}

Then, the general solution to the ODE is \(y(x) = c_1 e^{-x} + c_2 x e^{-x}\).

3a. Equations with constant coefficents

A common category of 2nd-order homogeneous ODEs are equations with constant coefficients, of the form:

\begin{equation} y^{\prime\prime} + a y^{\prime} + by = 0 \end{equation}

Note that these are unforced, and the right-hand side is zero.

Solutions to these equations take the form \(y(x) = e^{\lambda x}\), and inserting this into the ODE gives us the characteristic equation

\begin{equation} \lambda^2 + a \lambda + b = 0 \end{equation}

which we can solve to find the solution for given coefficients \(a\) and \(b\) and initial conditions. Depending on those coefficients and the solution to the characteristic equation, our solution can fall into one of three cases:

• Real roots: \(\lambda_1\) and \(\lambda_2\). This is an overdamped system and the full solution takes the form

  \begin{equation} y(x) = c_1 e^{\lambda_1 x} + c_2 e^{\lambda_2 x} \end{equation}

• Repeated roots: \(\lambda_1 = \lambda_2 = \lambda\). This is a critically damped system and the full solution is

  \begin{equation} y(x) = c_1 e^{\lambda x} + c_2 x e^{\lambda x} \end{equation}

  (Where does that second part come from, you might ask? Well, we know that \(y_1\) is \(e^{\lambda x}\), but the second part cannot also be \(e^{\lambda x}\) because those are linearly dependent. So, we use reduction of order to find \(y_2\), which is \(x e^{\lambda x}\).

• Imaginary roots: \(\lambda = \frac{-a}{2} \pm \beta i\), where \(\beta = \frac{1}{2} \sqrt{4b - a^2}\). This is an underdamped system and the solution takes the form

  \begin{equation} y(x) = e^{-ax/2} \left( c_1 \sin \beta x + c_2 \cos \beta x \right) \end{equation}

Some examples:

1. \(y^{\prime\prime} + 3 y^{\prime} + 2y = 0\)

   \begin{align} \rightarrow \lambda^2 + 3\lambda + 2 &= 0 \\ (\lambda + 2)(\lambda + 1) &= 0 \\ \lambda &= -2, -1 \\ y(x) &= c_1 e^{-x} + c_2 e^{-2x} \end{align}

   Then, we would use the initial conditions given for \(y(0)\) and \(y^{\prime}(0)\) to find \(c_1\) and \(c_2\).

2. \(y^{\prime\prime} + 6 y^{\prime} + 9y = 0\)

   \begin{align} \rightarrow \lambda^2 + 6\lambda + 9 &= 0 \\ (\lambda + 3)(\lambda + 3) &= 0 \\ \lambda &= -3 \\ y(x) &= c_1 e^{-3x} + c_2 x e^{-3x} \end{align}

3. \(y^{\prime\prime} + 6 y^{\prime} + 25 y = 0\)

   \begin{align} \rightarrow \lambda^2 + 6\lambda + 25 &= 0 \\ \lambda &= -3 \pm 4i \\ y(x) &= e^{-3x} \left( c_1 \sin 4x + c_2 \cos 4x \right) \end{align}

3b. Euler-Cauchy equations

Euler-Cauchy equations are of the form

\begin{equation} x^2 y^{\prime\prime} + axy^{\prime} + by = 0 \end{equation}

Solutions take the form \(y = x^m\), which when plugged into the ODE leads to a different characterisic equation to find \(m\):

\begin{align} y &= x^m \\ y^{\prime} &= m x^{m-1} \\ y^{\prime\prime} &= m (m-1) x^{m-2} \\ \rightarrow x^2 m (m-1) x^{m-2} + axmx^{m-1} + bx^m &= 0 \\ m^2 + (a-1)m + b &= 0 \end{align}

This is our new characteristic formula for these problems, and solving for the roots of this equation gives us \(m\) and thus our general solution.

Like equations with constant coefficients, we have three solution forms depending on the roots of the characteristic equation:

• Real roots: \(y(x) = c_1 x^{m_1} + c_2 x^{m_2}\)
• Repeated roots: \(y(x) = c_1 x^m + c_2 x^m \ln x\)
• Imaginary roots: \(m = \alpha \pm \beta i\), and \(y(x) = x^{\alpha} \left[c_1 \cos (\beta \ln x) + c_2 \sin (\beta \ln x)\right]\)

4. Inhomogeneous 2nd-order ODEs

Inhomogeneous, or forced, 2nd-order ODEs with constant coefficients take the form

\begin{equation} y^{\prime\prime} + a y^{\prime} + by = F(t) \end{equation}

with initial conditions \(y(0) = y_0\) and \(y^{\prime}(0) = y_0^{\prime}\). Depending on the form of the forcing function \(F(t)\), we can solve with techniques such as

• the method of undetermined coefficients
• variation of parameters
• LaPlace transforms

The solution in general to inhomogeneous ODEs includes two parts:

\begin{equation} y(t) = y_{\text{H}} + y_{\text{IH}} = c_1 y_1 + c_2 y_2 + y_{\text{IH}} \;, \end{equation}

where \(y_{\text{H}}\) is the solution from the equivalent homogeneous ODE \(y^{\prime\prime} + a y^{\prime} + b y = 0\).

The forcing function \(F(t)\) may be

• continuous
• periodic
• aperiodic/discontinuous

4a. Continuous \(F(t)\): method of undetermined coefficients

For continuous forcing functions, we have two solution methods: the method of undetermined coefficients, and variation of parameters.

Generally you’ll want to use the method of undetermined coefficients when possible, which depends on if \(F(t)\) matches one of a set of functions. In that case, the form of the inhomogeneous solution \(y_{\text{IH}}(t)\) follows that of the forcing function \(F(t)\), with one or more unknown constants:

\(F(t)\)	\(y_{\text{IH}}(t)\)
constant	\(K\)
\(\cos \omega t\)	\(K_1 \cos \omega t + K_2 \sin \omega t\)
\(\sin \omega t\)	\(K_1 \cos \omega t + K_2 \sin \omega t\)
\(e^{-at}\)	\(K e^{-at}\)
\((A) t\)	\(K_0 + K_1 t\)
\(t^n\)	\(K_0 + K_1 t + K_2 t^2 + \ldots + K_n t^n\)

For combinations of these functions, we can combine functions; for example, given

\begin{align} F(t) &= e^{-at} \cos \omega t \quad \text{or} e^{-at} \sin \omega t \\ y_{\text{IH}} &= K_1 e^{-at} \cos \omega t + K_2 e^{-at} \sin \omega t \end{align}

(Note how in all the above cases how the inhomogeneous solution follows the functional form of the forcing function; for example, the exponential decay rate \(a\) or the sinusoidal frequency \(\omega\) match.

The method of undetermined coefficients works by plugging the candidate inhomogeneous solutionn \(y_{\text{IH}}\) into the full ODE, and solving for the constants (e.g., \(K\))—but not from the initial conditions.

For example, let’s solve

\begin{equation} y^{\prime\prime} + 2y^{\prime} + y = e^{-x} \end{equation}

with initial conditions \(y(0) = y^{\prime}(0) = 0\). First, we should find the solution to the homogeneous equation

\begin{equation} y^{\prime\prime} + 2y^{\prime} + y = 0 \;. \end{equation}

We can do this by using the associated characteristic formula

\begin{align} \lambda^2 + 2 \lambda + 1 &= 0 \\ (\lambda + 1)(\lambda + 1) &= 0 \\ \rightarrow y_{\text{H}} &= c_1 e^{-x} + c_2 x e^{-x} \end{align}

To find the inhomogeneous solution, we would look at the table above to find what matches the forcing function \(e^{-x}\). Normally, we’d grab \(K e^{-x}\), but that would not be linearly independent from the first part of the homogeneous solution \(y_{\text{H}}\). The same is true for \(K x e^{-x}\), which is linearly dependent with the second part of \(y_{\text{H}}\), but \(K x^2 e^{-x}\) works! Then, we just need to find \(K\) by plugging this into the ODE:

\begin{align} y_{\text{IH}} &= K x^2 e^{-x} \\ y^{\prime} &= K e^{-x} (2x - x^2) \\ y^{\prime\prime} &= K e^{-x} (x^2 - 4x + 2) \\ 2 K &= 1 \\ \rightarrow K &= \frac{1}{2} \\ y_{\text{IH}} &= \frac{1}{2} x^2 e^{-x} \end{align}

Thus, the overall general solution is

\begin{equation} y(x) = c_1 e^{-x} + c_2 x e^{-x} + \frac{1}{2} x^2 e^{-x} \end{equation}

and we would solve for the integration constants \(c_1\) and \(c_2\) using the initial conditions.

Important points to remember:

• The constants of the inhomogeneous solution \(y_{\text{IH}}\) come from the ODE, not the initial conditions.
• Only solve for the integration constants \(c_1\) and \(c_2\) (part of the homogeneous solution) once you have the full general solution \(y = c_1 y_1 + c_2 y_2 + y_{\text{IH}}\).

4b. Continuous \(F(t)\): variation of parameters

We have the variation of parameters approach to solve for inhomogeneous 2nd-order ODEs that are more general:

\begin{equation} y^{\prime\prime} + p(x) y^{\prime} + q(x) y = r(x) \end{equation}

In this case, we can assume a solution \(y(x) = y_1 u_1 + y_2 u_2\).

The solution procedure is:

1. Obtain \(y_1\) and \(y_2\) by solving the homogeneous equation: \(y^{\prime\prime} + p(x) y^{\prime} + q(x) y = 0\)

2. Solve for \(u_1\) and \(u_2\):

   \begin{align} u_1 &= - \int \frac{y_2 r(x)}{W} dx + c_1 \\ u_2 &= \int \frac{y_1 r(x)}{W} dx + c_2 \\ W &= \begin{vmatrix} y_1 & y_2\\ y_1^{\prime} & y_2^{\prime}\\ \end{vmatrix} = y_1 y_2^{\prime} - y_2 y_1^{\prime} \;, \end{align}

   where \(W\) is the Wronksian.

3. Then, we have the general solution:

   \begin{align} y &= u_1 y_1 + u_2 y_2 \\ &= \left( -\int \frac{y_2 r(x)}{W} dx + c_1 \right) y_1 + \left( \int \frac{y_1 r(x)}{W} dx + c_2 \right) y_2 \;, \end{align}

   where we solve for \(c_1\) and \(c_2\) using the two initial conditions.