Imported from jeffreytse/grimoire-core (skills/science/mathematics/skills/apply-differential-equations/SKILL.md). Install upstream with npx skills add jeffreytse/grimoire-core --skill apply-differential-equations. Copyright stays with the author.

Apply Differential Equations

Classify, solve, and validate differential equations — choosing between analytical solutions (exact, separation of variables, integrating factor) and numerical methods (RK4, stiff solvers, finite difference/element) appropriate to the problem structure and stiffness.

Why This Is Best Practice

Why best: Matching solution method to equation structure (analytical form vs. stiffness-aware numerical solver) avoids both wasted derivation effort and numerically invalid results — the wrong solver on a stiff system produces exponentially growing errors, while an available closed form makes numerical approximation unnecessary. Adopted by: Differential equations model virtually all physical, biological, chemical, and engineering systems. Newton's laws, Maxwell's equations, Schrödinger equation, Navier-Stokes, the SIR epidemic model, Black-Scholes option pricing, and population dynamics are all differential equations. MATLAB ODE Suite, SciPy's solve_ivp, and Mathematica NDSolve are the standard numerical tools used in research and industry worldwide. Impact: Strogatz (2014) demonstrates that even qualitative analysis of ODEs (phase portraits, stability analysis) reveals system behavior without solving explicitly — a method used in epidemiology (epidemic thresholds), neuroscience (neural firing patterns), and climate modeling (tipping points). Hairer et al. (1993) showed that choosing the wrong solver for a stiff ODE (e.g., explicit RK4 on a stiff chemical kinetics problem) produces exponentially growing errors and requires time steps 1000× smaller than the implicit alternative.

Steps

1. Classify the differential equation

Before solving, identify the type — it determines the solution method:

Ordinary Differential Equations (ODE): one independent variable (usually time t)

• Order: highest derivative present (1st order: dy/dt = f(t,y); 2nd order: y'' + p(t)y' + q(t)y = g(t))
• Linearity: linear if y and all derivatives appear to the first power; nonlinear otherwise
• Autonomy: autonomous if f does not depend explicitly on t (dy/dt = f(y))

Partial Differential Equations (PDE): multiple independent variables

• Elliptic (Laplace, Poisson): equilibrium problems; no time dependence
• Parabolic (Heat equation): diffusion; marching in time
• Hyperbolic (Wave equation): wave propagation; second-order in both space and time

Stiffness: system is stiff if it contains time scales that differ by orders of magnitude; stiff systems require implicit solvers.

2. Apply analytical solution methods

1st order ODE (dy/dt = f(t,y)):

• Separable: if f(t,y) = g(t)h(y) → separate variables: dy/h(y) = g(t)dt; integrate both sides
• Linear 1st order: y' + P(t)y = Q(t) → use integrating factor μ(t) = e^∫P(t)dt
  • Solution: y = (1/μ) [∫μQ dt + C]
• Exact equations: M(x,y)dx + N(x,y)dy = 0 with ∂M/∂y = ∂N/∂x → exists F(x,y) = C

Linear 2nd order with constant coefficients (ay'' + by' + cy = g(t)):

• Characteristic equation: ar² + br + c = 0
• Roots:
  • Two real distinct: y = C₁e^(r₁t) + C₂e^(r₂t)
  • Repeated real root: y = (C₁ + C₂t)e^(rt)
  • Complex conjugate r = α ± βi: y = e^(αt)[C₁cos(βt) + C₂sin(βt)]
• Particular solution for non-homogeneous: method of undetermined coefficients (polynomial, exponential, sinusoidal forcing) or variation of parameters

Laplace transform (for IVP with forcing): transform to algebraic equation in s-domain; solve; inverse transform.

3. Select the numerical solver

For numerical solution, solver selection depends on problem type:

from scipy.integrate import solve_ivp

def f(t, y):
    return [-y[0] + y[1], -100*y[1]]  # stiff system

sol = solve_ivp(f, t_span=[0, 10], y0=[1, 0],
                method='Radau',  # or 'BDF' for stiff
                dense_output=True, rtol=1e-8, atol=1e-10)

4. Set tolerance and step size

• Relative tolerance (rtol): controls relative error; rtol=1e-6 means ~6 significant figures
• Absolute tolerance (atol): controls error for near-zero values; set to ~10× smaller than smallest meaningful value
• Never use rtol > 1e-3 for any published scientific computation
• Default tolerances (rtol=1e-3) are for quick checks only

Verify solution quality: halve the step size (or tighten tolerances) and check if solution changes meaningfully.

5. Conduct qualitative/stability analysis

For autonomous systems dy/dt = f(y), without solving explicitly:

• Fixed points: solve f(y*) = 0
• Linearization: compute Jacobian J = ∂f/∂y at each fixed point
• Stability: eigenvalues of J determine type:
  • All Re(λ) < 0: stable node/spiral (attracting)
  • Any Re(λ) > 0: unstable (repelling)
  • Re(λ) = 0: center or limit cycle (requires nonlinear analysis)

Phase portrait analysis reveals long-term behavior without numerical integration.

6. Validate the solution

For analytical solutions: verify by differentiating and substituting back into the ODE. Check initial/boundary conditions.

For numerical solutions:

• Compute at two different tolerances; compare solutions
• Compare with known analytical solution for simple limits
• Check conservation laws (energy, mass) if applicable
• Plot solution and residual; check for oscillations or drift

Common Mistakes

• Using explicit RK4 on a stiff system: Step size must be smaller than 1/max_eigenvalue(Jacobian). For stiff problems (eigenvalue ratio 1:10⁶), RK4 requires millions of tiny steps where an implicit solver takes hundreds.
• Not checking uniqueness: Without Lipschitz continuity of f, existence and uniqueness are not guaranteed. Blow-up solutions and non-uniqueness occur in finite time for some nonlinear ODEs.
• Confusing general solution with particular solution: The general solution contains integration constants (C₁, C₂). Apply initial conditions last, after finding the general form.

When NOT to Use

• Fractional-order differential equations (modeling anomalous diffusion, viscoelasticity): require specialized fractional calculus solvers; standard ODE integrators cannot handle fractional derivatives.