Ordinary Differential Equations (2): First-Order Methods

A bank account, a drug clearing the bloodstream, a tank of brine, a charging capacitor — they all obey the same kind of equation: a first-order ODE. The trick is recognising which of four shapes you are looking at, because each shape has a closed-form move that solves it cleanly. By the end of this chapter you will pattern-match an unfamiliar first-order equation in seconds and know exactly which lever to pull.

What You Will Learn#

• How to spot a separable equation and integrate it directly.
• The integrating factor that turns a linear equation into a perfect derivative.
• The exactness test and what it means geometrically: solutions are level curves of a potential.
• The Bernoulli substitution $v = y^{1-n}$ that linearises a whole family of nonlinear equations.
• Five real applications: drug half-life, RC charging, salt mixing, logistic growth, exponential interest — each one a small case study you can reuse.

Prerequisites#

• Chapter 1 ideas: what an ODE is, what an initial value problem is, how a slope field looks.
• Standard integration: substitution, integration by parts, partial fractions.

Separable equations: the most natural form#

The shape#

The constant $C$ is the family parameter — every choice picks out one solution curve. In the figure below, separating $\frac{dy}{dx} = -x/y$ gives $y\,dy + x\,dx = 0$ , which integrates to $x^2 + y^2 = C$ : a one-parameter family of circles. The slope-field arrows are tangent to these circles everywhere — exactly what “solution” means.

Worked example: exponential growth and decay#

Solve $\dfrac{dy}{dx} = ky$ for constant $k$ .

1. Separate. $\;\dfrac{dy}{y} = k\,dx$ .
2. Integrate both sides. $\;\ln|y| = kx + C_1$ .
3. Exponentiate. $\;y = A e^{kx}$ , where $A = \pm e^{C_1}$ absorbs the sign and the constant.

What just happened? The growth rate $k$ controls the time scale, not the shape. Two interpretations:

• $k > 0$ : doubling every $\ln 2 / k$ — populations with no resource limit, continuously compounded interest, runaway nuclear chain reaction.
• $k < 0$ : half-life of $\ln 2 / |k|$ — radioactive decay, a charged capacitor draining, a drug being cleared by the liver.

Application: drug metabolism and the half-life rule#

For ibuprofen $t_{1/2} \approx 2$ h. Starting from a 400 mg dose:

That’s why ibuprofen is dosed every 4–6 h: any longer and the level falls below the therapeutic window; any shorter and it accumulates. The “5 half-lives to clear” rule of thumb (about 99% gone) drops out of the same equation.

Application: the logistic equation#

Three structural facts fall out of the equation itself, before you ever solve it:

• $P \ll K$ : the bracket is $\approx 1$ , so growth is nearly exponential.
• $P = K/2$ : the right-hand side is $rK/4$ , the maximum possible rate.
• $P \to K$ : the bracket tends to zero — population saturates.

The right-hand panel below visualises this: the growth rate $\dot P$ as a function of $P$ is a downward parabola with peak at $K/2$ and roots at $P=0$ and $P=K$ .

First-order linear equations: the integrating factor#

Standard form#

This is the workhorse of applied math. Whenever a quantity changes in proportion to itself plus an external forcing term, you get this shape.

The trick#

That’s the entire method. Why does it work? $\mu$ is engineered so that $\mu' = \mu P$ , exactly what the product rule demands.

Worked example: $y' + 2xy = x$ #

As $x \to \pm\infty$ the term $C e^{-x^2}$ vanishes, so every solution funnels into the equilibrium $y = 1/2$ . The constant of integration only matters in a transient near the origin.

What the integrating factor does visually#

Three side-by-side panels make it obvious. We pick $y' + 2y = 4$ (so $P = 2$ , $\mu = e^{2x}$ , equilibrium $y=2$ ).

Application: charging an RC circuit#

The time constant $\tau$ is what every electrical engineer memorises:

Same formula, with a sign flip and $V_s = 0$ , governs discharging — and explains why your camera flash takes a noticeable second to recharge between shots.

Exact equations#

The geometric idea#

The solution curves are the level sets of $F$ . That is the whole picture, and it is profound: solving the ODE has been replaced by recovering a scalar function whose contour map is the solution family.

The exactness test#

If this fails, the equation is not exact (yet — see §3.4).

Solving an exact equation, step by step#

Take $(2x + y)\,dx + (x + 2y)\,dy = 0$ . Check exactness: $M_y = 1 = N_x$ . Good. Now reconstruct $F$ :

1. Integrate $M$ with respect to $x$ : $\;F = x^2 + xy + g(y)$ .
2. Differentiate with respect to $y$ and match $N$ : $\;x + g'(y) = x + 2y \Rightarrow g'(y) = 2y \Rightarrow g(y) = y^2$ .
3. Therefore $F(x,y) = x^2 + xy + y^2$ , and the solution family is $x^2 + xy + y^2 = C$ .

Geometrically these are tilted ellipses. The figure below overlays the gradient field $\nabla F = (M, N)$ on the level curves; the gradient is everywhere perpendicular to the level set, which is the formal way of saying that $\nabla F$ tells you “uphill” while the solution curve stays at the same height.

Salvage: making a non-exact equation exact#

If $M_y \neq N_x$ , multiply by $\mu(x,y)$ and look for a $\mu$ that fixes things. Two common cases admit closed forms:

This is the same construction as the linear-equation integrating factor — in fact, the linear method is a special case of this one.

Bernoulli equations#

The shape#

When $n = 0$ it’s already linear; when $n = 1$ it’s separable. The interesting range is everything else, where the $y^n$ term ruins linearity.

The substitution#

That is linear in $v$ . Solve it with the integrating factor of §2, then convert back via $y = v^{1/(1-n)}$ .

What the substitution does geometrically#

Take $y' + y = y^2$ (so $n = 2$ , $v = 1/y$ ). The substitution turns it into the linear equation $v' - v = -1$ , whose solutions are exponentials around the equilibrium $v = 1$ . In $v$ -space the family is straight and orderly; in the original $y$ -space it is curved and contains a finite-time blow-up.

The same substitution trick recovers the logistic equation: it is Bernoulli with $n = 2$ , $P(t) = -r$ , $Q(t) = -r/K$ . So the messy partial-fraction integration of §1.4 is really just an integrating-factor calculation in disguise.

A worked application: salt in a tank#

This is the classic “compartment model” — the same skeleton you will meet again in epidemiology, pharmacokinetics and chemical engineering.

Setup. A 1000 L tank starts with pure water. Brine with concentration $2$ g/L flows in at $5$ L/min. The tank is well-stirred and drains at the same rate, so the volume stays constant. Find the salt mass $Q(t)$ .

Read the answer. The equilibrium is $Q^\ast = 2000$ g, i.e. concentration $2$ g/L — exactly the inflow concentration. The time constant $\tau = 200$ min controls how fast we get there. The same equation with $Q(0) > Q^\ast$ describes a salty tank being washed clean — the curves below show both regimes converging to the same equilibrium.

Numerical methods preview#

Closed-form tricks fail the moment $f(x,y)$ is messy enough — and most real-world equations are. The answer is to step along the slope field.

Euler’s method#

• Local truncation error per step: $O(h^2)$ .
• Global error after $1/h$ steps: $O(h)$ . Hence “first-order method”.

Cheap, easy, often inadequate.

Runge–Kutta 4 (RK4)#

Global error $O(h^4)$ — four orders of magnitude better than Euler at the same step size, for the cost of three extra evaluations.

Slope fields tell you what the integrator will see#

Before you reach for a numerical solver, look at the right-hand side. Two superficially similar equations can have radically different long-term behaviour:

SciPy in three lines#

Newton’s law of cooling: $T' = -k(T - T_{\text{env}})$ with $T(0) = 90$ °C, $T_{\text{env}} = 20$ °C, $k = 0.1$ /min.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

from scipy.integrate import solve_ivp
import numpy as np

def cooling(t, T, k=0.1, T_env=20):
    return -k * (T - T_env)

sol = solve_ivp(cooling, [0, 60], [90.0], method="RK45", dense_output=True)

t = np.linspace(0, 60, 300)
T = sol.sol(t)[0]
print(f"At t=10 min: {T[50]:.1f} C")
print(f"At t=30 min: {T[150]:.1f} C")

solve_ivp adapts the step size automatically. For most problems you do not need to write your own RK4; you need to know when the standard solver will struggle (stiff systems — see Chapter 11 ).

Summary: the four-question filter#

When you meet a first-order equation, run it through these in order — the first “yes” tells you the method.

The four methods are not independent. The integrating factor for a linear equation is the special case of the integrating factor that fixes a non-exact equation. The Bernoulli substitution reduces nonlinear to linear. So learning these four well is really learning one idea — find the change of variable that turns the equation into something that integrates by inspection — applied four ways.

Exercises#

Basic.

1. Solve $y' = y^2 \sin x$ with $y(0) = 1$ . Where does the solution blow up?
2. Solve $y' + y = e^{-x}$ . Identify the transient and steady-state parts.
3. Verify that $(2xy + 3)\,dx + (x^2 + 1)\,dy = 0$ is exact and solve it.
4. A 1000 L tank holds 50 kg of dissolved salt. Pure water enters at 10 L/min; the well-mixed solution leaves at the same rate. Find $Q(t)$ and the time to halve the salt mass.

Advanced.

5. Solve the Bernoulli equation $y' + y = y^3$ . Sketch the phase line and identify the equilibria.
6. For $y' = 3y^{2/3}$ with $y(0) = 0$ , exhibit at least two distinct solutions on $[0, \infty)$ . Why does this not contradict the existence-uniqueness theorem?
7. A bacterial population obeys $P' = 0.5\,P(1 - P/1000)$ with $P(0) = 50$ . Find $P(t)$ explicitly and compute the time at which the growth rate is maximal.

Programming.

8. Implement Euler and RK4 from scratch. Compare their global error on $y' = -y$ , $y(0) = 1$ over $[0, 5]$ at step sizes $h = 0.5, 0.1, 0.01$ . Plot $\log(\text{error})$ versus $\log(h)$ and verify the predicted slopes 1 and 4.
9. Reproduce the slope-field comparison of §6.3 yourself: use numpy.meshgrid and matplotlib.pyplot.quiver for the field, then scipy.integrate.solve_ivp for the trajectories. Try $y' = y - y^5$ and identify the equilibria.

References#

• Boyce & DiPrima, Elementary Differential Equations and Boundary Value Problems, Wiley (2012). The standard undergraduate reference; Chapters 2–3 cover everything above with classical rigour.
• Zill, A First Course in Differential Equations with Modeling Applications, Cengage (2017). Heavier on applications; the mixing and circuit examples are particularly clean.
• Strogatz, Nonlinear Dynamics and Chaos, Westview (2014). Chapter 2’s phase-line treatment of the logistic equation is essential reading after this chapter.
• SciPy documentation: scipy.integrate ↗ . The reference for solve_ivp, including its choice of methods (RK45, Radau, BDF) and event handling.

This is Part 2 of the 18-part series on Ordinary Differential Equations.

• Part 1: Origins and Intuition
• Part 2: First-Order Methods (current)
• Part 3: Higher-Order Linear Theory
• Part 4: The Laplace Transform
• Part 5: Power Series and Special Functions
• Part 6: Linear Systems and the Matrix Exponential
• Part 7: Stability Theory
• Part 8: Nonlinear Systems and Phase Portraits
• Part 9: Chaos Theory and the Lorenz System
• Part 10: Bifurcation Theory
• Part 11: Numerical Methods
• Part 12: Boundary Value Problems
• Part 13: Introduction to PDEs
• Part 14: Epidemic Models
• Part 15: Population Dynamics
• Part 16: Fundamentals of Control Theory
• Part 17: Physics and Engineering Applications
• Part 18: Frontiers and Series Finale