Advanced Differential Equations

1. First-Order Ordinary Differential Equations

A first-order ODE relates a function y(x) to its first derivative y'. The general form is F(x, y, y') = 0 or equivalently y' = f(x, y). Several structural types admit systematic solution methods.

Separable Equations

A separable equation has the form dy/dx = g(x) h(y). Dividing both sides by h(y) and integrating gives the implicit solution integral(dy / h(y)) = integral(g(x) dx) + C. This method works whenever variables can be cleanly separated to opposite sides of the equation.

Linear First-Order Equations and Integrating Factors

The standard form is y' + P(x)y = Q(x). The integrating factor method multiplies both sides by mu(x) = exp(integral of P(x) dx), which makes the left side an exact derivative d/dx[mu(x) y]. Integrating both sides yields the general solution directly.

Exact Equations

An equation M(x,y) dx + N(x,y) dy = 0 is exact if the partial derivative of M with respect to y equals the partial derivative of N with respect to x. When exact, there exists a potential function F(x,y) with F_x = M and F_y = N, and the solution is F(x,y) = C. If the equation is not exact, one may seek an integrating factor mu(x) or mu(y) that makes it exact.

Bernoulli Equations

The Bernoulli equation y' + P(x)y = Q(x) y^n (n not 0 or 1) is nonlinear but reducible. The substitution v = y^(1-n) transforms it into a linear equation v' + (1-n)P(x)v = (1-n)Q(x), which is solvable by the integrating factor method. This technique applies to many biological and physical growth models.

Riccati Equations

The Riccati equation y' = P(x) + Q(x)y + R(x) y^2 is generally not solvable in closed form. However, if one particular solution y1(x) is known, the substitution y = y1 + 1/v converts it into a linear equation for v. Finding a particular solution often requires inspection, series methods, or recognizing special structure.

2. Second-Order Linear ODEs

The general second-order linear ODE is y'' + P(x)y' + Q(x)y = f(x). The associated homogeneous equation (f = 0) has a two-dimensional solution space. The general solution to the nonhomogeneous equation is y = y_h + y_p, where y_h is the general homogeneous solution and y_p is any particular solution.

Constant Coefficient Equations: Characteristic Equation

For ay'' + by' + cy = 0, the trial solution y = e^(rx) leads to the characteristic equation ar^2 + br + c = 0. The nature of the roots determines the form of the general solution.

The Wronskian and Linear Independence

For two solutions y1 and y2 of the homogeneous equation, the Wronskian is W(x) = y1(x)y2'(x) - y2(x)y1'(x). The solutions form a fundamental set of solutions (are linearly independent) if and only if W is nonzero at some (hence every) point of the interval. Abel's theorem gives the Wronskian explicitly without solving the ODE: W(x) = W(x0) * exp(-integral from x0 to x of P(t) dt).

Method of Undetermined Coefficients

For ay'' + by' + cy = f(x) where f is a polynomial, exponential, sine, cosine, or a product of these, guess a particular solution of the same form. If the guess overlaps the homogeneous solution, multiply by x (or x^2 if necessary). This method is efficient when applicable but limited to constant-coefficient equations with special right-hand sides.

Variation of Parameters

Variation of parameters finds a particular solution for any continuous right-hand side f(x), even when undetermined coefficients does not apply. Given fundamental solutions y1 and y2, assume y_p = u1(x) y1 + u2(x) y2 and impose the condition u1' y1 + u2' y2 = 0. This yields the system for u1' and u2', which is solved using the Wronskian.

Cauchy-Euler Equations

The equation x^2 y'' + ax y' + by = 0 has non-constant coefficients but the substitution x = e^t (or equivalently y = x^r) reduces it to a constant-coefficient equation. The trial solution y = x^r leads to the indicial equation r(r-1) + ar + b = 0, with the same three cases (distinct real, repeated, complex roots) as constant- coefficient equations.

3. Systems of Ordinary Differential Equations

Many physical phenomena are described by coupled ODEs. A system of n first-order linear ODEs can be written in matrix form x' = Ax + g(t), where x is a vector of unknowns and A is an n-by-n matrix. Systems also arise naturally from converting higher-order ODEs into first-order systems.

Matrix Form and the Eigenvalue Method

For the homogeneous system x' = Ax, seek solutions of the form x = v e^(lambda t) where v is a constant vector. Substituting gives Av = lambda v, so lambda must be an eigenvalue of A and v the corresponding eigenvector. The general solution is a linear combination of all such solutions over all eigenvalues.

Phase Plane Analysis

The phase plane for a 2D autonomous system x' = f(x,y), y' = g(x,y) shows trajectories in the (x,y) plane parametrized by time. The direction field is tangent to each trajectory. Equilibrium points (where f = g = 0) are classified by linearization: the Jacobian at the equilibrium determines whether it is a stable node, unstable node, saddle, stable spiral, unstable spiral, or center.

Fundamental Matrix and Matrix Exponential

A fundamental matrix Phi(t) is a matrix whose columns are n linearly independent solutions of x' = Ax. The general solution is x = Phi(t) c for constant vector c. For constant-coefficient systems, the fundamental matrix equals the matrix exponential e^(At), which can be computed via eigendecomposition (if A is diagonalizable) or the Cayley-Hamilton theorem. The matrix exponential satisfies d/dt[e^(At)] = A e^(At) and e^(A*0) = I.

4. Laplace Transforms

The Laplace transform converts a function of time t into a function of complex frequency s, transforming differential equations into algebraic equations. It is especially powerful for initial value problems and for handling discontinuous or impulsive forcing functions.

Definition and Basic Properties

For a function f(t) defined on t greater than or equal to 0, the Laplace transform is L(f)(s) = integral from 0 to infinity of e^(-st) f(t) dt, provided the integral converges. The transform exists whenever f is piecewise continuous and of exponential order (|f(t)| less than or equal to M e^(at) for some M, a, and all large t).

Transform Properties for Derivatives

The power of Laplace transforms for ODEs comes from the derivative property: L(f')(s) = s L(f)(s) - f(0), and L(f'')(s) = s^2 L(f)(s) - s f(0) - f'(0). These formulas encode the initial conditions automatically, transforming the ODE into an algebraic equation in Y(s) = L(y)(s).

Convolution Theorem

The convolution of two functions is (f * g)(t) = integral from 0 to t of f(tau) g(t - tau) d tau. The convolution theorem states L(f * g)(s) = L(f)(s) * L(g)(s). This is fundamental for solving ODEs where the right-hand side is a product in the s-domain, expressing the solution as a convolution integral. It also gives the response of a linear system to an arbitrary input in terms of the impulse response.

Partial Fractions for Inverse Transforms

After transforming an IVP, Y(s) is typically a rational function of s. Partial fraction decomposition breaks Y(s) into simpler terms that match the Laplace table. For a root at s = a with multiplicity m, the partial fraction expansion includes terms A/(s-a), B/(s-a)^2, ..., up to A/(s-a)^m. Complex conjugate roots a +/- ib contribute terms of the form (As + B) / ((s-a)^2 + b^2).

Heaviside Step Function and Discontinuous Forcing

The Heaviside unit step u(t - a) equals 0 for t less than a and 1 for t greater than or equal to a. It models a forcing function that switches on at time a. The second shifting theorem gives L(u(t-a) f(t-a))(s) = e^(-as) F(s). This enables closed-form Laplace treatment of piecewise-defined forcing functions, such as a hammer blow or a switch being thrown at a specific time.

5. Power Series Solutions

When an ODE has non-constant coefficients that do not yield to elementary methods, power series solutions provide systematic access to the solution near a specific point. The key distinction is whether the expansion point is ordinary or singular.

Ordinary vs. Singular Points

For y'' + P(x) y' + Q(x) y = 0, a point x0 is ordinary if P and Q are analytic at x0 (representable as convergent power series). At an ordinary point, two independent power series solutions exist and converge in a disk at least as large as the nearest singularity of P or Q in the complex plane. A singular point is where P or Q fails to be analytic. A singular point x0 is regular if (x-x0)P(x) and (x-x0)^2 Q(x) are analytic at x0; otherwise it is irregular.

The Frobenius Method

At a regular singular point x0, assume a solution of the form y = (x-x0)^r * sum over n of a_n (x-x0)^n with a_0 not 0. Substituting into the ODE and equating powers of (x-x0) to zero, the lowest power gives the indicial equation, a quadratic in r. The two roots r1 and r2 (with r1 greater than or equal to r2) of the indicial equation determine the structure of the solutions:

Bessel's Equation

Bessel's equation x^2 y'' + x y' + (x^2 - nu^2) y = 0 arises in problems with circular or cylindrical symmetry (drum vibrations, heat flow in cylinders, electromagnetic waveguides). The origin x = 0 is a regular singular point with indicial equation r^2 - nu^2 = 0, giving roots r = +/- nu. The Frobenius series with r = nu is the Bessel function of the first kind J_nu(x). For integer nu = n, the second solution Y_n(x) (Bessel function of the second kind or Neumann function) involves a logarithmic term.

Legendre's Equation

Legendre's equation (1-x^2)y'' - 2xy' + n(n+1)y = 0 arises in problems with spherical symmetry (gravitational potentials, quantum mechanics of hydrogen). The points x = +/- 1 are regular singular points; x = 0 is ordinary. For non-negative integer n, one solution is the Legendre polynomial P_n(x) of degree n. The first few are: P_0 = 1, P_1 = x, P_2 = (3x^2 - 1)/2, P_3 = (5x^3 - 3x)/2. Legendre polynomials are orthogonal on [-1, 1] with weight 1.

6. Sturm-Liouville Theory

Sturm-Liouville theory provides a unified framework for the orthogonal function expansions that arise throughout applied mathematics. It generalizes the familiar Fourier series to a broad class of second-order eigenvalue problems.

The Sturm-Liouville Problem

The regular Sturm-Liouville problem consists of the ODE d/dx[p(x) y'] - q(x) y + lambda w(x) y = 0 on [a, b], with boundary conditions alpha1 y(a) + alpha2 y'(a) = 0 and beta1 y(b) + beta2 y'(b) = 0. Here p, q, w are given functions with p and w positive on [a, b]. The parameter lambda is the eigenvalue; nontrivial solutions y_n are eigenfunctions.

Self-Adjoint Operators and Orthogonality

The Sturm-Liouville operator L[y] = -(1/w)(d/dx[p y'] - q y) is self-adjoint (symmetric) with respect to the weighted inner product (f, g) = integral from a to b of f(x) g(x) w(x) dx. As a consequence, all eigenvalues are real, and eigenfunctions corresponding to distinct eigenvalues are orthogonal with respect to this weighted inner product.

Generalized Fourier Series

By completeness, any suitable function f on [a, b] can be written as f(x) = sum over n of c_n y_n(x), where the coefficients are c_n = integral from a to b of f(x) y_n(x) w(x) dx / integral from a to b of y_n(x)^2 w(x) dx. This is the Sturm-Liouville version of the Fourier series. Standard Fourier sine and cosine series, Bessel series, and Legendre series are all special cases.

7. Partial Differential Equations: Classification

A second-order PDE in two independent variables x and y has the general form Au_xx + Bu_xy + Cu_yy + Du_x + Eu_y + Fu = G. The discriminant B^2 - 4AC classifies the PDE, analogous to the classification of conic sections.

Method of Characteristics for First-Order PDEs

For a first-order PDE a(x,y,u) u_x + b(x,y,u) u_y = c(x,y,u), characteristics are curves in the (x,y) plane along which the PDE reduces to an ODE. The characteristic equations are dx/ds = a, dy/ds = b, du/ds = c. By solving these ODEs along each characteristic curve, one propagates the solution from initial data. The method extends to quasi-linear and fully nonlinear first-order PDEs via the theory of contact transformations.

8. The Heat Equation

The heat equation u_t = k u_xx describes diffusion of heat (or any diffusing quantity) along a rod. Here k is the thermal diffusivity. The equation is parabolic; solutions are smooth for t greater than 0 regardless of initial data.

Separation of Variables

Assume u(x,t) = X(x) T(t). Substituting into the heat equation and dividing by X T gives T'/(kT) = X''/X = -lambda. Both sides equal the same constant (the separation constant). This yields two ODEs: T' + k lambda T = 0 and X'' + lambda X = 0. The spatial ODE with boundary conditions determines the allowable values of lambda (eigenvalues) and eigenfunctions X_n.

Neumann and Mixed Boundary Conditions

Neumann boundary conditions u_x(0,t) = u_x(L,t) = 0 (insulated ends) give eigenfunctions cos(n pi x / L) with n = 0, 1, 2, ... The n = 0 term is a_0/2 (constant), representing the conserved average temperature. Insulated BCs preserve the total thermal energy: d/dt integral from 0 to L of u dx = 0.

The Fundamental Solution (Heat Kernel)

On the full real line, the heat equation with initial data u(x,0) = f(x) has solution u(x,t) = integral from -infinity to infinity of G(x-y, t) f(y) dy, where the heat kernel is G(x,t) = (1 / sqrt(4 pi k t)) exp(-x^2 / (4kt)). The heat kernel is itself a Gaussian that broadens over time, reflecting diffusive spreading. As t approaches 0 from above, G approaches the Dirac delta distribution.

9. The Wave Equation

The wave equation u_tt = c^2 u_xx models vibrating strings, acoustic pressure waves, and electromagnetic fields. It is hyperbolic; disturbances propagate at finite speed c with no smoothing effect.

D'Alembert's Solution

The general solution on the infinite line is u(x,t) = F(x + ct) + G(x - ct), a superposition of waves traveling left and right at speed c. With initial conditions u(x,0) = f(x) (initial displacement) and u_t(x,0) = g(x) (initial velocity), d'Alembert's formula gives the solution explicitly without series.

Standing Waves on a Bounded String

For a string fixed at both ends (u(0,t) = u(L,t) = 0), separation of variables gives standing wave solutions u_n(x,t) = sin(n pi x / L) [A_n cos(n pi c t / L) + B_n sin(n pi c t / L)]. The natural frequencies are omega_n = n pi c / L. The general solution is a superposition of all modes; the Fourier coefficients A_n and B_n are determined by initial displacement and velocity.

Energy Conservation

The total energy of the wave equation system is E(t) = (1/2) integral from 0 to L of [u_t^2 + c^2 u_x^2] dx. For homogeneous Dirichlet or Neumann boundary conditions, dE/dt = 0: the wave equation conserves energy. This contrasts sharply with the heat equation, which dissipates energy.

Domain of Dependence and Influence

The value u(x0, t0) depends only on initial data in the interval [x0 - ct0, x0 + ct0] — the domain of dependence. Conversely, the initial data at x0 influences the solution only in the cone |x - x0| less than or equal to c t — the domain of influence. This finite propagation speed (Huygens' principle in 1D) is a fundamental feature of hyperbolic equations.

10. The Laplace Equation and Harmonic Functions

The Laplace equation u_xx + u_yy = 0 (or Delta u = 0 in higher dimensions) governs electrostatic potentials, steady-state heat distribution, and incompressible irrotational fluid flow. Solutions are called harmonic functions and possess remarkable regularity properties.

Mean Value Property

If u is harmonic in a domain, then at any interior point the value u(x0, y0) equals the average of u over any circle centered at that point lying in the domain. In three dimensions, the average is over a sphere. This property characterizes harmonic functions and has profound consequences: a harmonic function cannot have an interior maximum or minimum unless it is constant.

Maximum Principle

A nonconstant harmonic function on a bounded domain attains its maximum and minimum values only on the boundary, not in the interior. The strong maximum principle states that if the maximum is attained at any interior point, the function is identically constant. This implies uniqueness for the Dirichlet problem: two harmonic functions with the same boundary data must be equal.

Separation of Variables in Polar Coordinates

On a disk of radius a, the Laplace equation in polar coordinates (r, theta) is u_rr + (1/r) u_r + (1/r^2) u_theta_theta = 0. Separation gives u = R(r) Theta(theta), leading to Euler equations for R and periodic Theta. The bounded solution is u(r, theta) = A_0/2 + sum over n of r^n [A_n cos(n theta) + B_n sin(n theta)], with coefficients determined by the boundary condition u(a, theta) = f(theta) via Fourier series in theta.

Green's Functions

The Green's function G(x, y; x0, y0) for a domain is the response to a unit point source at (x0, y0) with zero boundary data. Once known, the solution to Delta u = f with zero Dirichlet data is u(x0, y0) = double integral of G(x,y; x0,y0) f(x,y) dA. For the upper half-plane, the Green's function is constructed by the method of images. For more complex geometries, conformal mapping can transform the domain to one where the Green's function is known.

11. Nonlinear ODEs and Qualitative Analysis

Most real-world ODEs are nonlinear and have no closed-form solutions. Qualitative analysis — studying the long-term behavior without explicit formulas — is essential. Phase plane methods, nullclines, and stability theory provide rich information about the dynamics.

Nullclines and the Phase Plane

For the autonomous system x' = f(x, y), y' = g(x, y), the x-nullclines are curves where f = 0 (so x' = 0), and the y-nullclines are curves where g = 0 (so y' = 0). Their intersections are equilibrium points. The sign of f and g in each region between nullclines gives the direction of flow, allowing a qualitative sketch of trajectories without solving the system.

Linearization at Equilibria

Near an equilibrium (x*, y*), the nonlinear system is approximated by its linearization: the system with coefficient matrix equal to the Jacobian J evaluated at (x*, y*). The Hartman-Grobman theorem guarantees that the local phase portrait of the nonlinear system near a hyperbolic equilibrium (eigenvalues with nonzero real parts) is topologically equivalent to the phase portrait of the linearization.

Poincaré-Bendixson Theorem

The Poincaré-Bendixson theorem is the principal existence result for limit cycles in the plane. If a trajectory of a continuously differentiable autonomous system remains in a compact region that contains no equilibrium point, then the trajectory's omega-limit set is a closed orbit. This theorem establishes the existence of periodic solutions in many biological and mechanical models.

Limit Cycles

A limit cycle is an isolated closed trajectory in the phase plane. A stable limit cycle attracts nearby trajectories as t approaches infinity. The van der Pol oscillator x'' - mu(1 - x^2)x' + x = 0 (mu greater than 0) is the canonical example: it has exactly one stable limit cycle. For small mu the cycle is nearly circular; for large mu the oscillation becomes a relaxation oscillation with rapid transitions and slow drift.

12. Applications: Population Models, Vibrations, and Circuits

Differential equations are indispensable tools in biology, mechanical engineering, and electrical engineering. The following canonical models illustrate the power and range of the methods developed above.

Logistic Population Model

The logistic equation dP/dt = rP(1 - P/K) models population growth with carrying capacity K. Unlike pure exponential growth, the logistic model saturates: as P approaches K the growth rate approaches zero. Separating variables yields the explicit solution P(t) = K / (1 + ((K - P0)/P0) e^(-rt)), which is an S-shaped (sigmoid) curve. The equilibria P = 0 (unstable) and P = K (stable) are confirmed by phase line analysis.

Lotka-Volterra Predator-Prey System

The system x' = ax - bxy (prey), y' = -cy + dxy (predator) models the coupled dynamics of prey (population x) and predator (population y). The constants a, b, c, d are positive. There are two equilibria: the trivial equilibrium (0, 0) (saddle point) and the coexistence equilibrium (c/d, a/b) (center in the linearization). All trajectories in the positive quadrant are closed orbits encircling the coexistence equilibrium — the populations oscillate periodically without ever reaching a stable steady state.

Mechanical Vibrations: Mass-Spring-Damper

Newton's second law for a mass m on a spring (constant k) with damping (coefficient b) gives mx'' + bx' + kx = F(t). Dividing by m: x'' + 2 gamma x' + omega_0^2 x = f(t), where omega_0 = sqrt(k/m) is the natural frequency and gamma = b/(2m) is the damping ratio. The system's behavior depends on whether gamma^2 compared to omega_0^2:

RLC Circuit Analysis

A series RLC circuit with resistance R, inductance L, and capacitance C driven by voltage E(t) satisfies L Q'' + R Q' + Q/C = E(t), where Q is the charge on the capacitor. This is mathematically identical to the mass-spring equation (L plays the role of m, R the role of damping, 1/C the role of the spring constant). The natural frequency is omega_0 = 1/sqrt(LC), and the damping constant is gamma = R/(2L). Laplace transforms are particularly efficient for circuit analysis with switched sources and initial charge/current conditions.

13. Practice Problems with Solutions

Work through these problems, then expand the solutions to check your understanding.

Quick Reference: Key Formulas

Frequently Asked Questions