Higher-Order Linear Differential Equations

A linear n-th order equation with constant coefficients is

aₙy⁽ⁿ⁾ + … + a₁y′ + a₀y = f(x).

Its general solution is y = y_c + y_p: the complementary function (the general solution of the homogeneous equation, with n constants) plus any particular integral.

The complementary function

Substitute y = e^mx to get the auxiliary equation aₙmⁿ + … + a₀ = 0.

Example. y″ − 5y′ + 6y = 0: m² − 5m + 6 = 0 gives m = 2, 3, so y = c₁e^2x + c₂e^3x. For y″ + 4y = 0, m = ±2i and y = c₁ cos 2x + c₂ sin 2x.

Particular integrals

Write the equation as F(D)y = f(x), with D = d/dx.

• f = e^ax: y_p = e^ax/F(a) if F(a) ≠ 0. If a is a root of multiplicity k, y_p = xᵏe^ax/F⁽ᵏ⁾(a).
• f = sin ax or cos ax: replace D² by −a² in F(D) (if that does not give zero).
• f = polynomial: expand 1/F(D) as a power series in D.
• f = e^axV(x): shift: 1/F(D) [e^axV] = e^ax · 1/F(D + a) [V].

Variation of parameters works for any f: with y_c = c₁y₁ + c₂y₂ and Wronskian W = y₁y₂′ − y₂y₁′, y_p = −y₁∫(y₂f/W) dx + y₂∫(y₁f/W) dx.

Euler-Cauchy equations

x²y″ + axy′ + by = 0. Substitute y = xᵐ: m(m − 1) + am + b = 0. Equivalently, put x = eᵗ to get a constant-coefficient equation in t.

Example. x²y″ − 2xy′ + 2y = 0: m² − 3m + 2 = 0, so m = 1, 2 and y = c₁x + c₂x².

Initial and boundary value problems

An initial value problem fixes y and its derivatives at one point; a boundary value problem fixes values at two points (for example y(0) and y(L)). A boundary value problem may have one, none or infinitely many solutions: y″ + y = 0 with y(0) = 0 and y(π) = 0 has infinitely many (y = c sin x).

Worked example

Solve y″ + 3y′ + 2y = 0 with y(0) = 1 and y′(0) = 0, and check y(1) with a numerical solver. Then find the particular integral of y″ − 3y′ + 2y = e^3x.

import math

# Roots -1 and -2: y = c1 e^-x + c2 e^-2x, with c1 + c2 = 1 and -c1 - 2 c2 = 0 -> c1 = 2, c2 = -1
exact = 2 * math.exp(-1) - math.exp(-2)

def f(y, v):
    return v, -3 * v - 2 * y

y, v, h = 1.0, 0.0, 1e-3
for _ in range(1000):                       # RK4 steps up to x = 1
    k1 = f(y, v)
    k2 = f(y + h / 2 * k1[0], v + h / 2 * k1[1])
    k3 = f(y + h / 2 * k2[0], v + h / 2 * k2[1])
    k4 = f(y + h * k3[0], v + h * k3[1])
    y += h / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0])
    v += h / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])
print(f"y(1): numerical {y:.8f}, exact 2/e - 1/e^2 = {exact:.8f}")

# PI of y'' - 3y' + 2y = e^{3x}: F(3) = 9 - 9 + 2 = 2, so y_p = e^{3x}/2
F = lambda a: a * a - 3 * a + 2
print("F(3) =", F(3), "-> y_p = e^(3x) /", F(3))

y(1): numerical 0.60042360, exact 2/e - 1/e^2 = 0.60042360
F(3) = 2 -> y_p = e^(3x) / 2

Practice questions