Ordinary Differential Equations (3): Higher-Order Linear Equations

When first-order equations aren't enough to describe a system, we need higher-order differential equations. Spring oscillations, bridge swaying, circuit resonance — all these phenomena require second-order or higher ODEs for modeling. This chapter systematically covers the theory and solution methods for higher-order linear ODEs.

Starting with Physical Intuition: Spring-Mass Systems

The Simplest Oscillation Model

Imagine an object of mass hanging from a spring. According to Hooke's Law, the restoring force is proportional to displacement:whereis the spring constant, and the negative sign indicates the force opposes displacement.

According to Newton's second law:Rearranging:whereis the natural angular frequency.

This is a second-order homogeneous linear ODE with constant coefficients!

Intuitive Understanding

Question: What's the solution to this equation?

Guess: What function's second derivative equals itself multiplied by a negative number?

Answer: Sine and cosine functions!

General solution:Or equivalently:whereis the amplitude andis the initial phase.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False

def simple_harmonic():
    omega0 = 2 * np.pi  # Frequency 1Hz
    t = np.linspace(0, 3, 500)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Different initial conditions
    A, phi = 1.0, 0
    x = A * np.cos(omega0 * t + phi)
    v = -A * omega0 * np.sin(omega0 * t + phi)
    
    # Displacement-time graph
    axes[0, 0].plot(t, x, 'b-', linewidth=2.5)
    axes[0, 0].set_xlabel('Time (s)', fontsize=12)
    axes[0, 0].set_ylabel('Displacement x (m)', fontsize=12)
    axes[0, 0].set_title('Simple Harmonic Motion', fontsize=14, fontweight='bold')
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].axhline(y=0, color='k', linewidth=0.5)
    
    # Velocity-time graph
    axes[0, 1].plot(t, v, 'r-', linewidth=2.5)
    axes[0, 1].set_xlabel('Time (s)', fontsize=12)
    axes[0, 1].set_ylabel('Velocity v (m/s)', fontsize=12)
    axes[0, 1].set_title('Velocity vs Time', fontsize=14, fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].axhline(y=0, color='k', linewidth=0.5)
    
    # Phase space trajectory
    axes[1, 0].plot(x, v, 'g-', linewidth=2.5)
    axes[1, 0].plot(x[0], v[0], 'ro', markersize=10, label='Start')
    axes[1, 0].set_xlabel('Displacement x (m)', fontsize=12)
    axes[1, 0].set_ylabel('Velocity v (m/s)', fontsize=12)
    axes[1, 0].set_title('Phase Portrait (Ellipse)', fontsize=14, fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].set_aspect('equal')
    axes[1, 0].legend()
    
    # Energy conservation
    KE = 0.5 * 1.0 * v**2  # Kinetic energy, assuming m=1
    PE = 0.5 * (omega0**2) * x**2  # Potential energy
    TE = KE + PE
    
    axes[1, 1].plot(t, KE, 'r-', linewidth=2, label='Kinetic Energy')
    axes[1, 1].plot(t, PE, 'b-', linewidth=2, label='Potential Energy')
    axes[1, 1].plot(t, TE, 'k--', linewidth=2, label='Total Energy')
    axes[1, 1].set_xlabel('Time (s)', fontsize=12)
    axes[1, 1].set_ylabel('Energy (J)', fontsize=12)
    axes[1, 1].set_title('Energy Conservation', fontsize=14, fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].legend()
    
    plt.tight_layout()
    plt.savefig('fig1_shm.png', dpi=150, bbox_inches='tight')
    plt.close()

simple_harmonic()

General Theory of Higher-Order Linear ODEs

Definition and Standard Form

th-order linear ODE:

Standard form (leading coefficient equals 1):

Homogeneous equation: Non-homogeneous equation:

Structure of Solutions

Theorem 1 (Superposition Principle): Ifandare solutions to the homogeneous equation, thenis also a solution.

Theorem 2 (General Solution Structure): The general solution of anth-order homogeneous linear ODE consists oflinearly independent particular solutions:

Theorem 3 (Non-homogeneous Equation): General solution of non-homogeneous = General solution of homogeneous + One particular solution

Linear Independence and the Wronskian

Definition: Functionsare linearly independent on intervalifimplies.

Wronskian determinant:

Theorem: Ifat some point, thenare linearly independent.

Example: Verify thatandare linearly independentSo they are linearly independent.

Constant-Coefficient Homogeneous Equations: The Characteristic Equation Method

Core Idea

Consider anth-order homogeneous linear ODE with constant coefficients:

Key insight: Try solutions of the formSubstituting into the equation: Since, we get the characteristic equation:

Three Cases

Case 1:Distinct Real Roots

If the characteristic equation hasdistinct real roots:

Example 1: SolveCharacteristic equation: , givingGeneral solution:

Case 2: Repeated Roots

Ifis a root of multiplicity, the corresponding solutions are:

Example 2: SolveCharacteristic equation: ,is a double root

General solution:

Case 3: Complex Roots

If there are complex roots(appearing in pairs), the corresponding real-valued solutions are:

Example 3: SolveCharacteristic equation: General solution:

def characteristic_equation_cases():
    x = np.linspace(0, 5, 300)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Case 1: Different real roots
    # y'' - 5y' + 6y = 0, y = c1*e^2x + c2*e^3x
    y1 = np.exp(2*x)
    y2 = np.exp(3*x)
    axes[0, 0].plot(x, y1, 'b-', linewidth=2, label='$e^{2x}$')
    axes[0, 0].plot(x, y2, 'r-', linewidth=2, label='$e^{3x}$')
    axes[0, 0].plot(x, y1 - y2, 'g--', linewidth=2, label='$e^{2x} - e^{3x}$')
    axes[0, 0].set_xlabel('x', fontsize=12)
    axes[0, 0].set_ylabel('y', fontsize=12)
    axes[0, 0].set_title("Case 1: Distinct Real Roots\n$y'' - 5y' + 6y = 0$", fontsize=12, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].set_ylim(-10, 50)
    
    # Case 2: Repeated roots
    # y'' - 4y' + 4y = 0, y = (c1 + c2*x)e^2x
    y1 = np.exp(2*x)
    y2 = x * np.exp(2*x)
    axes[0, 1].plot(x, y1, 'b-', linewidth=2, label='$e^{2x}$')
    axes[0, 1].plot(x, y2, 'r-', linewidth=2, label='$x e^{2x}$')
    axes[0, 1].plot(x, (1 + x) * np.exp(2*x), 'g--', linewidth=2, label='$(1+x)e^{2x}$')
    axes[0, 1].set_xlabel('x', fontsize=12)
    axes[0, 1].set_ylabel('y', fontsize=12)
    axes[0, 1].set_title("Case 2: Repeated Root\n$y'' - 4y' + 4y = 0$", fontsize=12, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].set_ylim(-10, 100)
    
    # Case 3: Complex roots - underdamped
    # y'' + 2y' + 5y = 0, r = -1 ± 2i
    y1 = np.exp(-x) * np.cos(2*x)
    y2 = np.exp(-x) * np.sin(2*x)
    env = np.exp(-x)
    
    axes[1, 0].plot(x, y1, 'b-', linewidth=2, label='$e^{-x}\\cos 2x$')
    axes[1, 0].plot(x, y2, 'r-', linewidth=2, label='$e^{-x}\\sin 2x$')
    axes[1, 0].plot(x, env, 'k--', linewidth=1.5, alpha=0.5, label='Envelope$e^{-x}$')
    axes[1, 0].plot(x, -env, 'k--', linewidth=1.5, alpha=0.5)
    axes[1, 0].set_xlabel('x', fontsize=12)
    axes[1, 0].set_ylabel('y', fontsize=12)
    axes[1, 0].set_title("Case 3: Complex Roots (Damped Oscillation)\n$y'' + 2y' + 5y = 0$", fontsize=12, fontweight='bold')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    # Case 3b: Pure imaginary roots - undamped
    # y'' + 4y = 0, r = ±2i
    y1 = np.cos(2*x)
    y2 = np.sin(2*x)
    
    axes[1, 1].plot(x, y1, 'b-', linewidth=2, label='$\\cos 2x$')
    axes[1, 1].plot(x, y2, 'r-', linewidth=2, label='$\\sin 2x$')
    axes[1, 1].plot(x, y1 + 0.5*y2, 'g--', linewidth=2, label='$\\cos 2x + 0.5\\sin 2x$')
    axes[1, 1].set_xlabel('x', fontsize=12)
    axes[1, 1].set_ylabel('y', fontsize=12)
    axes[1, 1].set_title("Case 3b: Pure Imaginary Roots (Undamped)\n$y'' + 4y = 0$", fontsize=12, fontweight='bold')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig2_characteristic.png', dpi=150, bbox_inches='tight')
    plt.close()

characteristic_equation_cases()

Damped Oscillation: Complete Analysis

Physical Model

Consider a spring system with damping: -: mass -: damping coefficient -: spring constant

Normalized form:where: -: natural frequency -: damping ratio

Characteristic Equation Analysis

Characteristic equation: Based on the value of, we have three cases:

Underdamped:

, whereis the damped frequency

Characteristic: Oscillates while gradually decaying

Critically Damped:

(repeated root)

Characteristic: Fastest return to equilibrium without oscillation

Overdamped:

(two negative real roots)

Characteristic: Slow return to equilibrium without oscillation

def damping_comparison():
    omega0 = 2 * np.pi  # Natural frequency
    t = np.linspace(0, 4, 500)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Different damping ratios
    zeta_values = [0.1, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(zeta_values)))
    
    for zeta, color in zip(zeta_values, colors):
        if zeta < 1:  # Underdamped
            omega_d = omega0 * np.sqrt(1 - zeta**2)
            x = np.exp(-zeta * omega0 * t) * np.cos(omega_d * t)
            label = f'ζ = {zeta} (underdamped)'
        elif zeta == 1:  # Critically damped
            x = (1 + omega0 * t) * np.exp(-omega0 * t)
            label = f'ζ = {zeta} (critical)'
        else:  # Overdamped
            r1 = -zeta * omega0 + omega0 * np.sqrt(zeta**2 - 1)
            r2 = -zeta * omega0 - omega0 * np.sqrt(zeta**2 - 1)
            # Initial conditions x(0)=1, x'(0)=0
            A = r2 / (r2 - r1)
            B = -r1 / (r2 - r1)
            x = A * np.exp(r1 * t) + B * np.exp(r2 * t)
            label = f'ζ = {zeta} (overdamped)'
        
        axes[0].plot(t, x, color=color, linewidth=2, label=label)
    
    axes[0].axhline(y=0, color='k', linewidth=0.5)
    axes[0].set_xlabel('Time (s)', fontsize=12)
    axes[0].set_ylabel('Displacement x', fontsize=12)
    axes[0].set_title('Effect of Damping Ratio ζ', fontsize=14, fontweight='bold')
    axes[0].legend(fontsize=9, loc='upper right')
    axes[0].grid(True, alpha=0.3)
    axes[0].set_xlim(0, 4)
    
    # Phase space comparison
    for zeta, color in zip([0.2, 1.0, 2.0], ['blue', 'green', 'red']):
        if zeta < 1:
            omega_d = omega0 * np.sqrt(1 - zeta**2)
            x = np.exp(-zeta * omega0 * t) * np.cos(omega_d * t)
            v = -np.exp(-zeta * omega0 * t) * (zeta * omega0 * np.cos(omega_d * t) + omega_d * np.sin(omega_d * t))
            label = f'ζ = {zeta} (underdamped)'
        elif zeta == 1:
            x = (1 + omega0 * t) * np.exp(-omega0 * t)
            v = -omega0**2 * t * np.exp(-omega0 * t)
            label = f'ζ = {zeta} (critical)'
        else:
            r1 = -zeta * omega0 + omega0 * np.sqrt(zeta**2 - 1)
            r2 = -zeta * omega0 - omega0 * np.sqrt(zeta**2 - 1)
            A = r2 / (r2 - r1)
            B = -r1 / (r2 - r1)
            x = A * np.exp(r1 * t) + B * np.exp(r2 * t)
            v = A * r1 * np.exp(r1 * t) + B * r2 * np.exp(r2 * t)
            label = f'ζ = {zeta} (overdamped)'
        
        axes[1].plot(x, v, color=color, linewidth=2, label=label)
        axes[1].plot(x[0], v[0], 'o', color=color, markersize=8)
    
    axes[1].plot(0, 0, 'k*', markersize=15, label='Equilibrium')
    axes[1].set_xlabel('Displacement x', fontsize=12)
    axes[1].set_ylabel('Velocity v', fontsize=12)
    axes[1].set_title('Phase Portrait', fontsize=14, fontweight='bold')
    axes[1].legend(fontsize=10)
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig3_damping.png', dpi=150, bbox_inches='tight')
    plt.close()

damping_comparison()

Engineering Application: Car Suspension

Goal: Design a suspension system that allows the car to smoothly and quickly recover after hitting a bump.

Optimal choice: Slightly underdamped ()

Too small: Too much oscillation, uncomfortable ride -: Theoretically fastest, but sensitive to parameter variations
Too large: Recovery too slow

Non-homogeneous Equations: Method of Undetermined Coefficients

Method Overview

For the non-homogeneous equation:

Steps: 1. Solve the homogeneous equation to get$y_hy_pf(x)y = y_h + y_p$

Guessing Rules

Form of	Guessed

or
(polynomial)

Important correction: If a term in the guessedis already a solution of, multiply by(or, etc.) to correct.

Detailed Example

Example 4: Solve Step 1: Solve homogeneous equationCharacteristic equation: Step 2: Guess particular solution, by the rule we guessButis already part ofWe need to correct toSubstituting into original equation:Substituting:So Step 3: General solution

Resonance Phenomenon

Example 5: SolveHomogeneous solution:(because)

Note:is already part ofGuessAfter calculation (details omitted):

Key observation: The particular solution contains a factor of, meaning the amplitude grows linearly with time!

This is the resonance phenomenon — when the driving frequency equals the natural frequency, amplitude increases without bound (in the undamped case).

Variation of Parameters

Why Do We Need It?

The method of undetermined coefficients only works for specific forms of(exponential, trigonometric, polynomial). For more general cases, we need variation of parameters.

Formula for Second-Order Equations

ForIf we know homogeneous solutions, the particular solution is:whereis the Wronskian.

Example

Example 6: SolveHomogeneous solutions:

RLC Circuits: A Classic Application of Higher-Order ODEs

Circuit Equations

Consider a series RLC circuit connected to a voltage source.

Kirchhoff's voltage law:Differentiating with respect to time (assumingis differentiable):This is a second-order ODE for currentOr using capacitor voltage:

Characteristic Equation Analysis

Characteristic equation for the homogeneous equation: Define: -: resonant frequency -: damping coefficient -: quality factor

Root classification: -: underdamped (oscillatory) -: critically damped -: overdamped

Converting to First-Order Systems

Transformation

Anyth-order ODE can be converted to a system offirst-order ODEs.

ForLet:We get:Matrix form:

Example: Second-Order Equation in System FormLet

The eigenvalues of matrix are the characteristic roots of the original equation!

Summary

In this chapter, we deeply studied the theory and applications of higher-order linear differential equations:

Core Methods

Equation Type	Solution Method
Constant-coefficient homogeneous	Characteristic equation
Constant-coefficient non-homogeneous	Undetermined coefficients + Characteristic equation
Variable-coefficient non-homogeneous	Variation of parameters
Complex equations	Convert to system + Numerical methods

Relationship Between Characteristic Roots and Solutions

Root Type	Corresponding Solution
Distinct real root
-fold real root
Complex roots

Physical Intuition

Damping ratio: Determines oscillation characteristics -: Underdamped (decaying oscillation) -: Critical damping (fastest non-oscillatory return) -: Overdamped (slow return)
Resonance: When driving frequency ≈ natural frequency, amplitude increases dramatically

Engineering Applications

Mechanical vibration and shock absorber design
Circuit analysis (RLC circuits)
Structural engineering (bridges, buildings)
Control system design

Exercises

Basic Problems

Find the general solution of.
Solve(Hint: repeated root).
Solve,,.
Solveand plot the solution.
Use undetermined coefficients to find a particular solution of.

Advanced Problems

Solveand discuss the resonance phenomenon.
Use variation of parameters to find the general solution of.
RLC circuit with,,:
- Calculate natural frequency and damping ratio
- Determine if underdamped or overdamped
- Plot the step response curve
Prove: For, when, the period of zero crossings is.
For the spring-mass-damper system, find the amplitude and phase of the steady-state response.

Programming Problems

Use Python to plot unit step response curves for different damping ratios.
Implement a general algorithm for variation of parameters and test with several examples.
Simulate a double pendulum system (this is a chaotic system) and observe sensitivity to initial conditions.
Design a simple PID controller and analyze its effect on second-order system stability.

References

Boyce, W. E., & DiPrima, R. C. (2012). Elementary Differential Equations. Wiley.
Kreyszig, E. (2011). Advanced Engineering Mathematics. Wiley.
Nagle, R. K., Saff, E. B., & Snider, A. D. (2017). Fundamentals of Differential Equations. Pearson.
MIT OCW 18.03 - Differential Equations (Video Lectures)

Next Chapter Preview: The Laplace Transform — Turning differential equations into algebraic equations magically