Differential Equations — Theory, Methods, and Applications

1. ODE Basics, Classification, Order, Linearity

A differential equation is an equation involving an unknown function and its derivatives. If the unknown depends on one variable, we call it an ordinary differential equation (ODE). If it depends on several variables, we have a partial differential equation (PDE). Differential equations are the language of change — they appear in mechanics, circuits, biology, economics, climate, and nearly every quantitative science.

Order, degree, and normal form

The order is the highest derivative present. y'' + 3y' + 2y = sin x is second-order. The general explicit first-order ODE is y' = f(x, y); the implicit form is F(x, y, y') = 0. An n-th order ODE can be written as y^(n) = f(x, y, y', ..., y^(n-1)) when solvable for the highest derivative.

Linear vs nonlinear

An ODE is linear if the unknown and its derivatives appear to the first power, with no products or nonlinear functions of y. The standard linear n-th order form is a_n(x) y^(n) + ... + a_1(x) y' + a_0(x) y = g(x). If g(x)=0 it's homogeneous; otherwise nonhomogeneous. Linearity unlocks superposition: if y1 and y2 solve the homogeneous equation, so does c1 y1 + c2 y2. Nonlinear equations like y' = y^2 - x, y'' + sin y = 0 (pendulum), or y' = |y| do not satisfy superposition, and their solution structure can be far richer.

Initial and boundary value problems

An initial value problem (IVP) prescribes y and its derivatives at a single point: y(x0)=y0, y'(x0)=y1. A boundary value problem (BVP) prescribes conditions at different points, e.g., y(0)=0, y(L)=0 for a vibrating string. IVPs typically model evolution in time; BVPs model steady states and spatial profiles.

Existence and uniqueness

For y' = f(x, y) with y(x0)=y0, the Picard–Lindelöf theorem states: if f and ∂f/∂y are continuous in a rectangle around (x0,y0), then a unique solution exists locally. Failure of Lipschitz continuity yields non-uniqueness: y' = 3 y^{2/3}, y(0)=0 has infinitely many solutions (including y=0 and y=x^3). Geometric intuition helps: the slope field must not cross itself for uniqueness to hold.

Geometric viewpoint

The equation y' = f(x,y) assigns a slope to every point (x,y). The collection of tiny line segments is the slope field. A solution is a curve everywhere tangent to the field. This picture predicts behavior without formulas: equilibria where f=0, nullclines where slope is zero or infinite, regions of increase/decrease.

2. First-Order ODEs: Separable, Linear, Exact, Bernoulli

First-order equations dominate modeling because many processes are specified by a rate.

Separable equations

If dy/dx = g(x) h(y), separate: ∫ dy/h(y) = ∫ g(x) dx + C. Example: logistic growth y' = r y (1 - y/K). Separate: dy/[y(1-y/K)] = r dx. Partial fractions give y = K / (1 + A e^{-rt}). Note singular solution y=0.

Linear first-order and integrating factors

Standard form: y' + P(x) y = Q(x). Multiply by μ(x)=exp(∫P dx). Then (μ y)' = μ Q, so y = μ^{-1}(∫ μ Q dx + C). Example: y' + 2x y = x, μ = e^{x^2}, solution y = 1/2 + C e^{-x^2}. Integrating factors turn non-exact derivatives into exact ones — the core trick of linear theory.

Exact equations

Write M(x,y) dx + N(x,y) dy = 0. It's exact if there exists ψ with ψ_x=M, ψ_y=N, which requires M_y = N_x on a simply connected domain. Then ψ(x,y)=C implicitly solves the ODE. If not exact, seek integrating factor μ(x) or μ(y): for μ(x), condition (M_y - N_x)/N depends only on x.

Homogeneous and Bernoulli

Homogeneous (in the scaling sense): y' = F(y/x). Substitute v=y/x ⇒ y = v x ⇒ y' = v + x v'. Reduces to separable in v. Bernoulli: y' + P y = Q y^n (n≠0,1). Let v = y^{1-n}, then v' + (1-n)P v = (1-n)Q, linear in v. Example: y' + y = x y^3, set v=y^{-2} ⇒ v' -2v = -2x.

Autonomous equations and phase line

For y' = f(y), no explicit x. Equilibria are roots of f(y)=0. Sign of f determines flow direction on the y-axis. Stability: if f' < 0 at equilibrium, it's asymptotically stable (flow toward); if f'>0, unstable. This 1D picture is the foundation for higher-dimensional phase portraits.

Qualitative methods often beat formulas. Even when closed-form solutions exist (involving special functions), plotting the slope field reveals basins of attraction, blow-up times, and separatrices that formulas hide.

3. Second-Order Linear with Constant Coefficients

The workhorse model across physics: a y'' + b y' + c y = g(t). Homogeneous part a y'' + b y' + c y = 0 yields characteristic polynomial a r^2 + b r + c = 0.

• Distinct real roots r1≠r2: y = C1 e^{r1 t} + C2 e^{r2 t}. Overdamped oscillator.
• Repeated root r: y = (C1 + C2 t) e^{r t}. Critical damping.
• Complex r = α ± iβ: y = e^{α t}(C1 cos βt + C2 sin βt). Underdamped oscillations.

Wronskian and fundamental sets

For y1, y2 solutions, W = y1 y2' - y1' y2. For linear ODE y''+p y'+q y=0, Abel's identity gives W(t)=W0 exp(-∫p dt). If W≠0, {y1,y2} spans the solution space.

Nonhomogeneous: undetermined coefficients

For constant-coefficient RHS like polynomials, exponentials, sines/cosines, guess a particular form of same family. If overlap with homogeneous solution occurs, multiply by t (or t^s). Example: y''+y = sin t. Homogeneous cos t, sin t. Guess yp = t(A cos t + B sin t) due to resonance. Substituting yields yp = -½ t cos t.

Variation of parameters

For general g(t): seek yp = u1 y1 + u2 y2 with u1' y1 + u2' y2 =0. Solving gives u1' = -y2 g/W, u2' = y1 g/W. Works for variable coefficients too.

Resonance and damping

For forced oscillator m x'' + c x' + k x = F0 cos ωt, steady-state amplitude A(ω)=F0 / sqrt((k - m ω^2)^2 + (c ω)^2). Resonance peak near ω0 = sqrt(k/m). With c→0, amplitude grows linearly in time when ω=ω0 (pure resonance). Damping limits amplitude and shifts peak.

Physical interpretations

Mechanical: m x'' + c x' + k x = 0 models mass-spring-damper. Electrical: L Q'' + R Q' + Q/C = E(t) for RLC circuit. The math is identical — a triumph of linear theory.

4. First-Order Systems, Phase Portraits, Eigenvalues

Any n-th order ODE becomes a first-order system. For y'' + p y' + q y = g, set x1 = y, x2 = y' ⇒ x' = [0 1; -q -p] x + [0; g]. Systems unify theory.

Linear systems x' = A x

Solution: x(t)=e^{At} x0. For 2×2 with eigenvalues λ, eigenvectors v, fundamental solutions e^{λ t} v. Classification by trace τ = tr A, determinant Δ=det A:

The τ-Δ plane organizes all possibilities. Repeated eigenvalues yield defective nodes (improper). Nonlinear systems locally resemble their linearization if Re(λ)≠0 (Hartman–Grobman).

Nonhomogeneous and fundamental matrix

For x' = A x + f(t), solution x = Φ(t) C + Φ(t) ∫ Φ^{-1}(s) f(s) ds, where Φ is fundamental matrix (columns independent solutions). Matrix exponential e^{At} = Σ (At)^k/k! computes directly for constant A.

5. Laplace Transforms for IVPs

Define L{f}(s)=F(s)=∫₀^∞ e^{-st} f(t) dt. Key properties:

• Linearity: L{af+bg}=aF+bG
• Derivative: L{f'}=sF - f(0), L{f''}=s²F - s f(0)-f'(0)
• Shift: L{e^{at}f}=F(s-a)
• Convolution: L{f∗g}=F G where (f∗g)=∫₀^t f(τ)g(t-τ)dτ
• Step: L{u(t-a)}=e^{-as}/s ; Delta: L{δ(t-a)}=e^{-as}

Procedure: transform ODE → algebraic equation in s → solve for Y(s) → partial fractions → invert using table. Laplace automatically incorporates initial conditions and handles discontinuous forcing cleanly.

Example: y''+4y=0, y(0)=1, y'(0)=0 ⇒ (s²Y - s)+4Y=0 ⇒ Y=s/(s²+4) ⇒ y=cos 2t.

For y''+2y'+2y = u(t-π) sin(t-π), transform yields Y = e^{-π s}/[(s²+1)((s+1)²+1)], invert via convolution or residues.

6. PDEs: Heat, Wave, Laplace, Separation of Variables, Fourier Series

PDEs govern fields: temperature u(x,t), displacement u(x,t), potential φ(x,y). Second-order linear PDE in two variables: A u_xx + 2B u_xy + C u_yy + ... =0. Type: elliptic if B²-AC<0 (Laplace), parabolic if =0 (heat), hyperbolic if >0 (wave).

Heat equation u_t = α u_xx

Fourier's law yields diffusion. On 0u= Σ b_n sin(nπx/L) e^{-α(nπ/L)² t}, where b_n are Fourier sine coefficients of initial data f(x): b_n = 2/L ∫₀^L f(x) sin(nπx/L) dx. High modes decay faster; solutions smooth instantly.

Wave equation u_tt = c² u_xx

Models string, acoustics. With fixed ends, separation yields X_n=sin(nπx/L), T_n=A_n cos ω_n t + B_n sin ω_n t, ω_n = c nπ/L. Solution is superposition of standing waves. d'Alembert form on infinite line: u = [f(x-ct)+f(x+ct)]/2 + (1/2c)∫_{x-ct}^{x+ct} g(s) ds for initial displacement f and velocity g.

Laplace equation u_xx + u_yy =0

Steady-state heat, electrostatics. Solutions are harmonic; satisfy maximum principle (no interior maxima). On rectangle, separation gives products of sines/hyperbolic sines. On disk, polar coordinates yield Fourier series in θ and powers r^n.

Fourier series

Periodic extensions decompose into sines/cosines. For period 2L, f∼ a0/2 + Σ [an cos(nπx/L)+bn sin(nπx/L)]. Convergence: piecewise smooth ⇒ pointwise to midpoint at jumps; Gibbs overshoot ≈9% near discontinuities. Fourier series are eigenfunctions of d²/dx² with periodic BCs — spectral basis for separation.

7. Nonlinear Dynamics and Chaos

Linear systems cannot produce isolated periodic orbits or chaos. Nonlinearity can.

Logistic growth and bifurcations

Continuous logistic: y' = r y(1 - y/K) → stable equilibrium at K. Discrete logistic map x_{n+1}=r x_n(1-x_n) undergoes period-doubling cascade to chaos as r increases, paradigmatic route to chaos.

Lotka–Volterra predator-prey

x' = ax - bxy, y' = -cy + dxy. Equilibria (0,0) saddle, (c/d, a/b) center for linearization; nonlinear system has closed orbits with conserved quantity V = d x - c ln x + b y - a ln y. Populations oscillate quarter-phase apart.

Van der Pol and limit cycles

x'' - μ(1 - x²) x' + x =0. For μ>0, trajectories approach isolated periodic orbit (limit cycle), modeling heartbeat and vacuum tubes.

Lorenz attractor

Simplified convection: x' = σ(y-x), y' = x(ρ - z)-y, z' = xy - βz. For σ=10, β=8/3, ρ=28, solutions never repeat, stay bounded, and exhibit sensitive dependence — the butterfly effect. The attractor has fractal structure (~2.06 dimension). Chaos requires ≥3 dimensions for continuous autonomous flows.

8. Numerical Methods: Euler, RK4, Stability, Stiffness

Most ODEs lack elementary closed forms. Numerical integrators approximate y(t_{n+1}) from y(t_n).

Euler method

y_{n+1}= y_n + h f(t_n,y_n). Local error O(h²), global O(h). Explicit, simple, but stability region small: for test equation y' = λ y, stable iff |1 + hλ| ≤1.

Improved Euler (Heun) and midpoint

Heun: predictor y* = y_n + h f_n, corrector y_{n+1}= y_n + h/2 (f_n + f(t_{n+1},y*)). Second-order. Midpoint: y_{n+1}= y_n + h f(t_n+h/2, y_n + h/2 f_n).

Classical RK4

k1=f(tn,yn), k2=f(tn+h/2, yn+h k1/2), k3=f(tn+h/2, yn+h k2/2), k4=f(tn+h, yn+h k3), yn+1=yn+h(k1+2k2+2k3+k4)/6. Fourth-order, excellent balance of accuracy and cost. Workhorse for non-stiff problems.

Stability and stiffness

Stiffness occurs when components decay at vastly different rates (eigenvalues with large negative real parts). Explicit methods require tiny h for stability, not accuracy. Implicit Euler y_{n+1}= y_n + h f(t_{n+1},y_{n+1}) is A-stable (stable for all Re(hλ)<0). BDF methods and Radau collocation handle stiff chemistry and circuits.

Adaptive step size

Embedded pairs (e.g., RK45 Fehlberg, Dormand-Prince) estimate local error e ≈ |y_{high} - y_{low}|, adjust h to keep e within tolerance. This is what ode45 implements.

Conservation and symplectic integrators

For Hamiltonian systems (planetary orbits), standard RK drifts energy. Symplectic methods (e.g., Verlet, implicit midpoint) preserve phase-space area and long-term energy behavior.

Octave / MATLAB Example Scripts (Self-Contained)

These six scripts run offline in Octave. Copy into files or the REPL.

1) First-order logistic growth with ode45

% logistic_ode.m
r = 1.5; K = 100; y0 = 10; tspan = [0 10];
f = @(t,y) r*y.*(1 - y/K);
[t,y] = ode45(f, tspan, y0);
plot(t,y,'LineWidth',2); grid on;
xlabel('t'); ylabel('y'); title('Logistic Growth');

2) Damped harmonic oscillator (second-order)

% damped_osc.m
w0=2; zeta=0.1;
osc = @(t,Y) [Y(2); -2*zeta*w0*Y(2) - w0^2*Y(1)];
[t,Y] = ode45(osc, [0 20], [1; 0]);
plot(t, Y(:,1)); grid on;
xlabel('t'); ylabel('x(t)'); title('Underdamped Oscillator');

3) Phase portrait of linear system

% phase_portrait.m
A = [ -0.1 1 ; -1 -0.1 ];
[X,Y] = meshgrid(linspace(-3,3,21));
U = A(1,1)*X + A(1,2)*Y; V = A(2,1)*X + A(2,2)*Y;
quiver(X,Y,U,V); axis equal; hold on;
for theta = linspace(0,2*pi,8)
  x0 = 2.5*[cos(theta); sin(theta)];
  [t,s] = ode45(@(t,s) A*s, [0 20], x0);
  plot(s(:,1), s(:,2), 'LineWidth',1.5);
end
title('Stable Spiral'); grid on;

4) Laplace transform solution (symbolic)

% laplace_demo.m
pkg load symbolic
syms t s Y
y0=1; dy0=0;
F = 1/(s^2+1); % forcing sin t
 eqn = s^2*Y - s*y0 - dy0 + 2*(s*Y - y0) + 2*Y - F;
 Ysol = solve(eqn, Y);
 y = ilaplace(Ysol, s, t);
 pretty(y)

5) 1D heat equation with pdepe

% heat_pdepe.m
m=0; x=linspace(0,1,50); t=linspace(0,0.5,100); alpha=0.01;
pdefun = @(x,t,u,DuDx) deal(alpha*DuDx, u, 1);
icfun = @(x) sin(pi*x) + 0.5*sin(pi*3*x);
bcfun = @(xl,ul,xr,ur,t) deal(ul, ur, 0, 0);
sol = pdepe(m,pdefun,icfun,bcfun,x,t);
surf(x,t,sol); xlabel('x'); ylabel('t'); zlabel('u'); title('Heat Equation');

6) Lorenz attractor with ode45

% lorenz.m
sigma=10; beta=8/3; rho=28;
lor = @(t,s) [sigma*(s(2)-s(1)); s(1)*(rho-s(3))-s(2); s(1)*s(2)-beta*s(3)];
[t,s] = ode45(lor, [0 50], [1;0;0]);
plot3(s(:,1), s(:,2), s(:,3)); grid on;
xlabel('x'); ylabel('y'); zlabel('z'); title('Lorenz Attractor');