Reduction to a First-Order System
The second-order linear ordinary differential equation with constant coefficients represents a cornerstone of mathematical physics, describing a vast array of phenomena from mechanical vibrations to electrical circuits.
The classical model for such systems is the simple harmonic oscillator, which we will analyze using several analytical techniques.
Subsequently, we will extend this analysis to the more general case of the damped harmonic oscillator.
A simple harmonic oscillator describes a system where a mass m is subject to a restoring force proportional to its displacement x from an equilibrium position.
This is an idealized system without any dissipative (damping) or external (driving) forces. The physical law governing this system is Hooke’s Law:
F = -kx
where k is the spring constant.
Using Newton’s second law of motion, F = ma, we establish the equation of motion:
m \frac{\mathrm{d}^2 x}{\mathrm{d}t^2} = -k x
This equation is conventionally rearranged into the standard form of a second-order linear homogeneous ODE:
\frac{\mathrm{d}^2 x}{\mathrm{d}t^2} + \frac{k}{m} x(t) = 0
To specify a unique trajectory, two initial conditions are required, typically the initial position x(0) = x_0 and the initial velocity \dot{x}(0) = v_0.
We define the natural angular frequency \omega_0^2 = k/m, rewriting the equation as:
\frac{\mathrm{d}^2 x}{\mathrm{d}t^2} + \omega_0^2 x(t) = 0
We will now demonstrate three distinct but related methods for obtaining the solution.
A systematic approach for solving linear homogeneous ODEs with constant coefficients is to propose a solution of the form x(t) = e^{\lambda t}.
The rationale is that the exponential function retains its form under differentiation. Substituting into the equation yields:
\begin{aligned} & \frac{\mathrm{d}}{\mathrm{d}t}(e^{\lambda t}) = \lambda e^{\lambda t} \\ & \frac{\mathrm{d}^2}{\mathrm{d}t^2}(e^{\lambda t}) = \lambda^2 e^{\lambda t} \\ & \lambda^2 e^{\lambda t} + \omega_0^2 e^{\lambda t} = 0 \end{aligned}
Since e^{\lambda t} is never zero, we can divide by it to obtain the characteristic equation:
\begin{aligned} & \lambda^2 + \omega_0^2 = 0 \\ & \lambda = \pm i \omega_0 \end{aligned}
The two roots are purely imaginary and form a complex conjugate pair.
The principle of superposition for linear homogeneous equations states that the general solution is a linear combination of the two independent solutions corresponding to these roots:
x(t) = c_1 e^{i \omega_0 t} + c_2 e^{-i \omega_0 t}
Here, c_1 and c_2 are complex constants. While this is a complete solution, a real-valued physical system should have a real-valued solution.
Using Euler’s formula, e^{i\theta} = \cos(\theta) + i\sin(\theta), we can express the solution in terms of real-valued trigonometric functions:
\begin{aligned} x(t) & = c_1(\cos(\omega_0 t) + i\sin(\omega_0 t)) + c_2(\cos(\omega_0 t) - i\sin(\omega_0 t)) \\ & = (c_1 + c_2)\cos(\omega_0 t) + i(c_1 - c_2)\sin(\omega_0 t) \end{aligned}
By defining new constants A = c_1 + c_2 and B = i(c_1 - c_2), the solution takes the familiar form:
x(t) = A\cos(\omega_0 t) + B\sin(\omega_0 t)
For x(t) to be real, A and B must be real constants. We determine these constants by applying the initial conditions:
\begin{aligned} & x(0) = A\cos(0) + B\sin(0) = A = x_0 \\ & \dot{x}(t) = -A\omega_0\sin(\omega_0 t) + B\omega_0\cos(\omega_0 t) \\ & \dot{x}(0) = -A\omega_0\sin(0) + B\omega_0\cos(0) = B\omega_0 = v_0 \\ & B = \frac{v_0}{\omega_0} & \end{aligned}
The particular solution is:
x(t) = x_0 \cos(\omega_0 t) + \frac{v_0}{\omega_0} \sin(\omega_0 t)
Another technique involves seeking a solution as a power series expansion around t=0. We assume a solution of the form:
x(t) = \sum_{n=0}^\infty C_n t^n
Differentiating term-by-term and substituting into the ODE \ddot{x} + \omega_0^2 x = 0:
\sum_{n=2}^\infty n(n-1)C_n t^{n-2} + \omega_0^2 \sum_{n=0}^\infty C_n t^n = 0
To combine these sums, we shift the index of the first sum by letting k = n-2, which gives n = k+2:
\begin{aligned} & \sum_{k=0}^\infty (k+2)(k+1)C_{k+2} t^k + \sum_{k=0}^\infty \omega_0^2 C_k t^k = 0 \\ & \sum_{k=0}^\infty \left[ (k+2)(k+1)C_{k+2} + \omega_0^2 C_k \right] t^k = 0 \end{aligned}
For this identity to hold for all t, the coefficient of each power of t must be zero. This yields the recurrence relation:
C_{k+2} = -\frac{\omega_0^2}{(k+2)(k+1)} C_k
The initial conditions provide the first two coefficients: C_0 = x(0) = x_0 and C_1 = \dot{x}(0) = v_0. The recurrence relation couples coefficients with indices that differ by two. This separates the series into even and odd terms, determined by C_0 and C_1 respectively.
For the even coefficients:
\begin{aligned} & C_2 = -\frac{\omega_0^2}{2 \cdot 1} C_0 = -\frac{\omega_0^2}{2!} x_0 \\ & C_4 = -\frac{\omega_0^2}{4 \cdot 3} C_2 = (-1)^2 \frac{\omega_0^4}{4!} x_0 \\ & C_{2m} = (-1)^m \frac{\omega_0^{2m}}{(2m)!} x_0 \end{aligned}
For the odd coefficients:
\begin{aligned} & C_3 = -\frac{\omega_0^2}{3 \cdot 2} C_1 = -\frac{\omega_0^2}{3!} v_0 \\ & C_5 = -\frac{\omega_0^2}{5 \cdot 4} C_3 = (-1)^2 \frac{\omega_0^4}{5!} v_0 \\ & C_{2m+1} = (-1)^m \frac{\omega_0^{2m}}{(2m+1)!} v_0 \end{aligned}
The full solution is the sum of the even and odd series:
x(t) = x_0 \sum_{m=0}^\infty (-1)^m \frac{(\omega_0 t)^{2m}}{(2m)!} + \frac{v_0}{\omega_0} \sum_{m=0}^\infty (-1)^m \frac{(\omega_0 t)^{2m+1}}{(2m+1)!}
Recognizing the Maclaurin series for \cos(z) = \sum_{m=0}^\infty (-1)^m \frac{z^{2m}}{(2m)!} and \sin(z) = \sum_{m=0}^\infty (-1)^m \frac{z^{2m+1}}{(2m+1)!}, we recover the identical solution:
x(t) = x_0 \cos(\omega_0 t) + \frac{v_0}{\omega_0} \sin(\omega_0 t)
Any n^{th} order ODE can be converted into a system of n first-order ODEs. For the harmonic oscillator, we define a state vector \mathbf{v}(t) composed of position and velocity:
\mathbf{v}(t) = \begin{bmatrix} x(t) \\ \dot{x}(t) \end{bmatrix}
The time evolution of this state vector is given by a matrix equation:
\frac{\mathrm{d}}{\mathrm{d}t} \mathbf{v}(t) = \begin{bmatrix} \dot{x}(t) \\ \ddot{x}(t) \end{bmatrix} = \begin{bmatrix} \dot{x}(t) \\ -\omega_0^2 x(t) \end{bmatrix}
This can be written in the compact linear form \dot{\mathbf{v}} = \mathbf{A} \mathbf{v}:
\frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} x \\ \dot{x} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\omega_0^2 & 0 \end{bmatrix} \begin{bmatrix} x \\ \dot{x} \end{bmatrix}
The solution to this system is given by the matrix exponential \mathbf{v}(t) = e^{\mathbf{A}t} \mathbf{v}(0).
The eigenvalues of the matrix \mathbf{A} are found by solving:
\det(\mathbf{A} - \lambda\mathbf{I}) = 0
which yields:
\lambda^2 + \omega_0^2 = 0
The eigenvalues are:
\lambda = \pm i\omega_0
precisely the roots of the characteristic equation found with the exponential guess, demonstrating the consistency between these methods.
A more realistic physical model includes a damping force, typically proportional to the velocity of the mass, F_d = -d\dot{x}, where d > 0 is the damping coefficient.
The balance of forces now becomes:
m \ddot{x} = -k x - d \dot{x}
Rearranging into standard form gives the equation for the damped harmonic oscillator:
m \ddot{x} + d \dot{x} + k x = 0
We again use the exponential x(t) = e^{\lambda t}. Substituting this into the equation results in a new characteristic equation:
\begin{aligned} & m\lambda^2 e^{\lambda t} + d\lambda e^{\lambda t} + k e^{\lambda t} = 0 \\ m \lambda^2 + d \lambda + k = 0 \end{aligned}
The roots of this quadratic equation determine the behavior of the system:
\lambda_{1,2} = \frac{-d \pm \sqrt{d^2 - 4mk}}{2m}
The general solution is:
x(t) = c_1 e^{\lambda_1 t} + c_2 e^{\lambda_2 t}
The nature of the solution depends on the sign of the discriminant:
\Delta = d^2 - 4mk
Overdamped Case (\Delta > 0): If d^2 > 4mk, there are two distinct, real, negative roots:
\lambda_1, \lambda_2 \in \mathbb R > 0
The solution is a sum of two decaying exponential functions. The system returns to its equilibrium position without oscillating.
Critically Damped Case (\Delta = 0): If d^2 = 4mk, there is one repeated real negative root:
\lambda = -\frac{d}{2m}
The two linearly independent solutions are e^{\lambda t} and t e^{\lambda t}. The general solution is x(t) = (c_1 + c_2 t)e^{\lambda t}. This represents the fastest possible return to equilibrium without any oscillation.
Underdamped Case (\Delta < 0): If d^2 < 4mk, the roots are a complex conjugate pair. Let \alpha = d/2m be the decay constant and \omega_d = \frac{\sqrt{4mk - d^2}}{2m} be the damped frequency.
The roots are:
\lambda_{1,2} = -\alpha \pm i\omega_d
The general real-valued solution is:
x(t) = e^{-\alpha t} \left( A\cos(\omega_d t) + B\sin(\omega_d t) \right)
This describes an oscillation with a decaying amplitude, where the system oscillates around the equilibrium point with progressively smaller swings.