Chapter 21
18 min read
Section 179 of 353

Substitution Methods

First-Order Differential Equations

Learning Objectives

By the end of this section, you will be able to:

  1. Recognize when a first-order ODE can be tamed by a change of variable, even though it is neither linear nor exact.
  2. Apply the homogeneous substitution v=y/xv = y/x to any ODE of the form dy/dx=F(y/x)dy/dx = F(y/x).
  3. Linearize a Bernoulli equation y+P(x)y=Q(x)yny' + P(x)y = Q(x)y^n via u=y1nu = y^{1-n}.
  4. Handle equations of the form dy/dx=f(ax+by+c)dy/dx = f(ax + by + c) with the linear shift v=ax+by+cv = ax + by + c.
  5. Verify the algebra numerically using Python and PyTorch's autograd.
  6. Connect these substitutions to real-world models: population growth, fluid drag, mixing, and beyond.

The Big Idea: Why Substitute?

"We do not solve hard problems. We trade them for easy ones."

In Sections 21.1 and 21.2 you mastered two structurally easy equations: linear and exact. The integrating factor handled the first, the potential function the second. But the wild ODE you meet in a physics lab rarely walks in already wearing one of those costumes. Consider three real specimens:

Not linear
dy/dx = (x² + y²) / (x y)
Not exact
dy/dx = r y (1 − y/K)
Not separable
dy/dx = (x + y + 1)²

Each of these resists every technique you know. And yet — with the right change of variable, all three collapse into separable or linear equations you can solve in your sleep.

The substitution mindset

Substitution is the universal escape hatch. Whenever a first-order ODE doesn't fit any standard mould, ask: is there a single quantity inside this equation that, if I gave it its own name, everything would simplify around? That quantity is your vv (or uu), and naming it is half the battle.

The three patterns below come up so often that they have dedicated names. In every case the recipe is the same four-step dance:

  1. Diagnose the structure on the right side.
  2. Introduce a new variable that captures that structure.
  3. Rewrite the ODE in the new variable — it will be separable or linear.
  4. Solve and undo the substitution to get y(x)y(x) back.

Diagnosing the Right Substitution

Pattern matching the right side of dy/dx=f(x,y)dy/dx = f(x,y) tells you which substitution to reach for:

Right side depends on…Call itSubstitutionEquation becomes
only the ratio y/xHomogeneousv = y/xseparable in v, x
y^n with a linear y elsewhereBernoulliu = y^(1−n)linear in u
only ax + by + cLinear shiftv = ax + by + cseparable in v, x

Spotting homogeneity at a glance

A right side f(x,y)f(x,y) is homogeneous of degree zero exactly when f(λx,λy)=f(x,y)f(\lambda x, \lambda y) = f(x,y) for every λ0\lambda \neq 0. Try λ=1/x\lambda = 1/x: you get f(1,y/x)f(1, y/x), a function of y/xy/x alone.


Homogeneous Equations: v = y/x

A first-order ODE is homogeneous of degree zero when its right side depends only on the ratio y/xy/x:

dydx=F ⁣(yx)\frac{dy}{dx} = F\!\left(\frac{y}{x}\right)

Set v=y/xv = y/x, so y=vxy = v\,x. Differentiating with the product rule gives y=v+xvy' = v + x\,v'. Substituting:

v+xdvdx=F(v)dvdx=F(v)vxv + x\,\frac{dv}{dx} = F(v) \quad\Longrightarrow\quad \frac{dv}{dx} = \frac{F(v) - v}{x}

The right side is a product of a function of vv and a function of xx — it is separable:

dvF(v)v=dxx=lnx+C\int \frac{dv}{F(v) - v} = \int \frac{dx}{x} = \ln|x| + C

What is the geometric meaning?

Slopes of a homogeneous ODE are constant along rays from the origin. Walk out along the ray y=mxy = mxand the slope F(m)F(m) never changes. Move the slider in the demo below to feel this.

Loading homogeneous substitution visualizer…

The Three-Line Recipe

  1. Write the ODE so its right side is a function of y/xy/x alone, call it F(v)F(v).
  2. Replace with v=y/x,  y=vx,  y=v+xvv = y/x,\; y = vx,\; y' = v + xv' and rearrange to dvF(v)v=dxx\frac{dv}{F(v) - v} = \frac{dx}{x}.
  3. Integrate both sides, then back-substitute v=y/xv = y/x.

Worked Example — Homogeneous

Solve dydx=x+yxy,y(1)=0\dfrac{dy}{dx} = \dfrac{x + y}{x - y}, \quad y(1) = 0.

Click to expand step-by-step solution (try it yourself first!)

Step 1 — Confirm homogeneity. Divide numerator and denominator by xx:

dydx=1+y/x1y/x=F(v)   with v=y/x.\frac{dy}{dx} = \frac{1 + y/x}{1 - y/x} = F(v) \;\text{ with } v = y/x.

Step 2 — Substitute. Let y=vxy = vx, so y=v+xvy' = v + xv':

v+xv=1+v1v    xv=1+v1vv=1+vv+v21v=1+v21v.v + x\,v' = \frac{1+v}{1-v} \;\Longrightarrow\; x\,v' = \frac{1+v}{1-v} - v = \frac{1+v - v + v^2}{1-v} = \frac{1+v^2}{1-v}.

Step 3 — Separate variables.

1v1+v2dv=dxx.\frac{1-v}{1+v^2}\,dv = \frac{dx}{x}.

Step 4 — Integrate. Split the left integral into two pieces:

dv1+v2vdv1+v2=dxx\int \frac{dv}{1+v^2} - \int \frac{v\,dv}{1+v^2} = \int \frac{dx}{x}
arctanv12ln(1+v2)=lnx+C.\arctan v - \tfrac12 \ln(1+v^2) = \ln|x| + C.

Step 5 — Back-substitute v=y/xv = y/x:

arctan ⁣(yx)12ln ⁣(1+y2x2)=lnx+C.\arctan\!\left(\frac{y}{x}\right) - \tfrac12 \ln\!\left(1 + \frac{y^2}{x^2}\right) = \ln|x| + C.

Step 6 — Apply the initial condition y(1)=0y(1) = 0: the left side becomes arctan012ln1=0\arctan 0 - \tfrac12 \ln 1 = 0, and the right side is ln1+C=C\ln 1 + C = C, so C=0C = 0. The implicit solution is:

arctan ⁣(yx)=lnx+12ln ⁣(1+y2x2).\arctan\!\left(\frac{y}{x}\right) = \ln|x| + \tfrac12 \ln\!\left(1 + \frac{y^2}{x^2}\right).

Step 7 — Sanity check. Differentiate implicitly and simplify — the slope at (1,0)(1, 0) should equal (1+0)/(10)=1(1 + 0)/(1 - 0) = 1, and indeed it does. Plot the curve in the demo above to see it threads through (1,0)(1, 0).


Bernoulli Equations: u = y^(1−n)

A Bernoulli equation is anything of the form

dydx+P(x)y=Q(x)yn,n0,1.\frac{dy}{dx} + P(x)\,y = Q(x)\,y^n, \qquad n \neq 0, 1.

It is one tiny step removed from a linear equation — only that troublesome yny^n on the right makes it non-linear. The magic substitution is u=y1nu = y^{1 - n}.

Why this exponent?

We want a quantity whose derivative absorbs the yny^n. Differentiate u=y1nu = y^{1-n}:

dudx=(1n)yndydx.\frac{du}{dx} = (1-n)\,y^{-n}\,\frac{dy}{dx}.

Now multiply the original ODE by (1n)yn(1-n)\,y^{-n}:

(1n)yny+(1n)P(x)y1n=(1n)Q(x).(1-n)\,y^{-n}\,y' + (1-n)\,P(x)\,y^{\,1-n} = (1-n)\,Q(x).

The first term is exactly du/dxdu/dx, and y1n=uy^{1-n} = u. So:

  dudx+(1n)P(x)u=(1n)Q(x)  \boxed{\;\frac{du}{dx} + (1-n)\,P(x)\,u = (1-n)\,Q(x)\;}

A linear first-order ODE in uu! Solve it with the integrating factor method from Section 21.1, then convert back via y=u1/(1n)y = u^{1/(1-n)}.

The famous case n = 2 — the logistic equation

With n=2n = 2 the substitution becomes u=y1=1/yu = y^{-1} = 1/y, and the equation y=ry(1y/K)y' = ry\bigl(1 - y/K\bigr) rewritten as yry=(r/K)y2y' - ry = -(r/K)\,y^2 turns into the linear equation u+ru=r/Ku' + ru = r/K. Play with rr, KK, and y0y_0 below to see how the S-curve emerges.

Loading Bernoulli logistic demo…

Worked Example — Bernoulli (Logistic)

Solve the logistic ODE dydx=ry(1yK),y(0)=y0,\dfrac{dy}{dx} = ry\bigl(1 - \dfrac{y}{K}\bigr), \quad y(0) = y_0, with r,K,y0>0r, K, y_0 > 0.

Click to expand step-by-step solution

Step 1 — Rewrite in Bernoulli standard form.

y=ryrKy2    yry=rKy2.y' = ry - \tfrac{r}{K}\,y^2 \;\Longrightarrow\; y' - ry = -\tfrac{r}{K}\,y^2.

Here P(x)=r,  Q(x)=r/K,  n=2P(x) = -r,\; Q(x) = -r/K,\; n = 2.

Step 2 — Substitute u=y1n=y1u = y^{1-n} = y^{-1}:

dudx=y2dydx,y=1/u.\frac{du}{dx} = -y^{-2}\,\frac{dy}{dx}, \quad y = 1/u.

Step 3 — Multiply the ODE by (1n)yn=y2(1-n)\,y^{-n} = -y^{-2}:

y2y+ry1=rK.-y^{-2}y' + r\,y^{-1} = \tfrac{r}{K}.

The left side is exactly u+ruu' + ru, giving the linear equation:

dudx+ru=rK.\frac{du}{dx} + r\,u = \frac{r}{K}.

Step 4 — Integrating factor μ=erx\mu = e^{rx}:

ddx(erxu)=rKerx.\frac{d}{dx}\bigl(e^{rx}\,u\bigr) = \frac{r}{K}\,e^{rx}.

Step 5 — Integrate both sides:

erxu=1Kerx+C    u=1K+Cerx.e^{rx}\,u = \frac{1}{K}\,e^{rx} + C \;\Longrightarrow\; u = \frac{1}{K} + C\,e^{-rx}.

Step 6 — Apply u(0)=1/y0u(0) = 1/y_0:

1y0=1K+C    C=1y01K=Ky0Ky0.\frac{1}{y_0} = \frac{1}{K} + C \;\Longrightarrow\; C = \frac{1}{y_0} - \frac{1}{K} = \frac{K - y_0}{K\,y_0}.

Step 7 — Convert back to y=1/uy = 1/u and clean up:

y(x)=K1+Ky0y0erx.y(x) = \frac{K}{1 + \dfrac{K - y_0}{y_0}\,e^{-rx}}.

Step 8 — Read the physics off the formula.

  • x → 0: the exponential is 1, so denominator = K/y₀, giving y(0) = y₀. ✓
  • x → ∞: exponential → 0, so y → K (the carrying capacity).
  • Inflection point: y = K/2 — the population grows fastest when half of capacity is filled. This is why S-curves show up in advertising spend, virus spread, and tech adoption.

Linear-Shift Substitution: v = ax + by + c

If the entire right side depends on a single linear combination,

dydx=f(ax+by+c),\frac{dy}{dx} = f(ax + by + c),

the right substitution is to give that combination a name. Set v=ax+by+cv = ax + by + c. Differentiating:

dvdx=a+bdydx=a+bf(v).\frac{dv}{dx} = a + b\,\frac{dy}{dx} = a + b\,f(v).

The new equation is separable: dva+bf(v)=dx\dfrac{dv}{a + b\,f(v)} = dx. Integrate, then undo the substitution.

Loading linear-shift demo…

Worked micro-example

For y=(x+y+1)2y' = (x + y + 1)^2, set v=x+y+1v = x + y + 1. Then v=1+v2v' = 1 + v^2, separable to arctanv=x+C\arctan v = x + C. Undoing: x+y+1=tan(x+C)x + y + 1 = \tan(x + C), so y=tan(x+C)x1y = \tan(x + C) - x - 1.


Python: Hand-rolled Substitution Solver

Let's see the homogeneous substitution work end-to-end on a computer. We'll start from the original ODE, apply the substitution by hand to derive the simplified equation, then ask SciPy to solve that simplified equation — and finally un-substitute to recover y(x)y(x). As a last step we will plug the answer back into the original ODE and measure the residual — if our algebra was right, the residual should be tiny.

Solving a homogeneous ODE by substitution
🐍substitution_solver.py
1NumPy + SciPy imports

NumPy gives us array math; scipy.integrate.odeint is the workhorse ODE solver (LSODA under the hood). The substitution trick is independent of these libraries — it only changes what equation we hand to the solver.

9Pin down the ORIGINAL equation

The equation dy/dx = (x + y)/(x − y) is homogeneous: divide top and bottom by x and the right side becomes (1 + y/x)/(1 − y/x), a function of the single ratio v = y/x. That observation is the whole reason this section exists.

14F(v) — what the right side becomes after substitution

Plugging v = y/x into the original right side cancels every loose x. F encapsulates that simplified expression so the rest of the code never has to know what the original equation looked like.

17dv/dx = (F(v) − v) / x — the substituted ODE

Differentiating y = v·x with the product rule gives y' = v + x·v'. Setting y' equal to F(v) and solving for v' yields v' = (F(v) − v)/x. Notice it is now separable: dv/(F(v) − v) = dx/x. The if-guard avoids dividing by zero exactly at x = 0.

22Pick a point on the solution curve

We start at (x₀, y₀) = (1.0, 0.5). The very first thing the substitution requires is converting the initial condition: v₀ = y₀/x₀ = 0.5. We solve in v-coordinates from now on.

26Integrate in v-space

odeint marches v from x = 1 to x = 4 using the v-equation. Internally it uses adaptive step control, but conceptually each step is Euler-like: v_{n+1} = v_n + h · (F(v_n) − v_n)/x_n.

27Undo the substitution: y = v · x

We only solved for v — but the user wants y. Multiplying element-wise rebuilds the original curve. This 'forward map then inverse map' pattern is the universal recipe of all substitution methods.

30Numerical derivative for verification

np.gradient uses central differences to approximate dy/dx from the y values we just produced. If our substitution work is correct, this should match the original ODE's right-hand side closely.

31Plug y back into the ORIGINAL right side

We compute (x + y)/(x − y) at every grid point and compare. The maximum absolute difference is our residual — it should be on the order of the solver tolerance, around 1e-4 or smaller.

36Read the printout

Expected output: v moves from 0.5 toward a value where the algebraic solution implicitly defines the curve. The residual line is the trust check — if it were 1.0 you would know the substitution had a bug.

25 lines without explanation
1import numpy as np
2from scipy.integrate import odeint
3import matplotlib.pyplot as plt
4
5# ------------------------------------------------------------------
6# Original homogeneous ODE:   dy/dx = (x + y) / (x - y)
7# Substitution:               v = y / x   =>   y = v x
8# Derived ODE in v:           dv/dx = (F(v) - v) / x
9# where F(v) = (1 + v) / (1 - v).
10# ------------------------------------------------------------------
11
12def F(v):
13    return (1.0 + v) / (1.0 - v)
14
15def dv_dx(v, x):
16    if abs(x) < 1e-9:
17        return 0.0
18    return (F(v) - v) / x
19
20# Initial point in (x, y), then convert to v.
21x0, y0 = 1.0, 0.5
22v0 = y0 / x0                 # = 0.5
23
24x_grid = np.linspace(1.0, 4.0, 200)
25v_sol  = odeint(dv_dx, v0, x_grid).flatten()
26y_sol  = v_sol * x_grid      # convert v back to y
27
28# Verify by plugging back into the ORIGINAL ODE.
29dy_numeric = np.gradient(y_sol, x_grid)
30dy_formula = (x_grid + y_sol) / (x_grid - y_sol)
31max_residual = np.max(np.abs(dy_numeric - dy_formula))
32
33print(f"v(x0) = {v0:.4f}, y(x0) = {y0:.4f}")
34print(f"v(x_end) = {v_sol[-1]:.4f}, y(x_end) = {y_sol[-1]:.4f}")
35print(f"max |dy/dx_numeric - F(x,y)| = {max_residual:.4e}")

PyTorch: Verifying with Autograd

Once we hand-derive a closed-form solution from a substitution, how do we know we didn't flip a sign or forget a chain-rule factor? PyTorch's autograd engine gives us a one-line answer: differentiate the candidate solution and compare with the right side of the original ODE. If they match at every point, the derivation is correct.

Verifying the logistic solution with PyTorch autograd
🐍autograd_verify.py
1Why PyTorch at all?

We are not training a model. We are using PyTorch's automatic differentiation engine as a derivative oracle: it computes dy/dx exactly from any expression we build, so it can verify that a closed-form solution actually satisfies its ODE.

14Logistic parameters

Growth rate r = 1.2 (units 1/time), carrying capacity K = 10 (max value y can approach), initial population y₀ = 1.0. With y₀ < K the solution starts below capacity and grows toward K.

19x as a leaf tensor with requires_grad=True

Setting requires_grad turns x into a leaf node in the autograd graph. Anything we compute downstream remembers how to send a gradient back to x. Without this flag the next call would error out.

23A = (K − y₀)/y₀ — the integration constant

When we solved u = 1/y → du/dx + ru = r/K and applied the initial condition u(0) = 1/y₀, the resulting integration constant in y-form is exactly A. This single number is what 'memory' of the initial condition becomes.

24y(x) = K / (1 + A · e^(−r x))

This is the famous logistic curve. As x → ∞ the exponential decays to 0, so y → K (saturation). At x = 0 the denominator is 1 + A = K/y₀, giving y(0) = K · y₀/K = y₀. ✓

27Why grad_outputs=ones?

torch.autograd.grad computes the sum Σᵢ gᵢ · dyᵢ/dxⱼ where gᵢ is the i-th entry of grad_outputs. Passing a vector of 1's gives the diagonal — i.e. dyᵢ/dxᵢ at each grid point — which is exactly the pointwise derivative we want.

28create_graph=False

We only need a first derivative, so we tell autograd to discard intermediate computation graphs after use. Setting it to True would let us differentiate dy/dx again (useful for second-order ODE checks).

31Recompute the RHS of the ORIGINAL ODE

r · y · (1 − y/K) is what dy/dx is supposed to equal. If our algebra was right, autograd's answer and this expression must agree everywhere.

34L∞ error — the trust check

We take the absolute difference at every point and grab the maximum. A correct derivation produces an error near machine epsilon (1e-6 to 1e-7 in float32). A non-zero residual would mean a sign error or a missed term.

37Reading the printout

The first line proves the substitution was algebraically correct. The next two confirm the boundary conditions: y(0) hits y₀, y(end) approaches K. Three checks, all green ⇒ the closed-form is bulletproof.

27 lines without explanation
1import torch
2
3# ------------------------------------------------------------------
4# Bernoulli case: logistic equation
5#     dy/dx = r * y * (1 - y/K)
6# Substitution u = 1/y turns it into the LINEAR equation
7#     du/dx + r*u = r/K
8# with closed-form solution
9#     u(x) = 1/K + (1/y0 - 1/K) * exp(-r*x)
10# Therefore
11#     y(x) = K / (1 + ((K - y0)/y0) * exp(-r*x))
12# ------------------------------------------------------------------
13
14r  = 1.2
15K  = 10.0
16y0 = 1.0
17
18# 1. Build x as a PyTorch tensor that tracks gradients.
19x = torch.linspace(0.0, 6.0, 200, requires_grad=True)
20
21# 2. Evaluate the closed-form y(x) we derived from the substitution.
22A = (K - y0) / y0
23y = K / (1 + A * torch.exp(-r * x))
24
25# 3. Ask autograd for dy/dx at every grid point in one shot.
26ones = torch.ones_like(y)
27dy_dx, = torch.autograd.grad(y, x, grad_outputs=ones, create_graph=False)
28
29# 4. Compute what dy/dx SHOULD be from the original ODE.
30rhs = r * y * (1 - y / K)
31
32# 5. The substitution is correct iff autograd matches the ODE rhs.
33err = (dy_dx - rhs).abs().max()
34
35print(f"max |autograd(y') - r*y*(1 - y/K)| = {err.item():.3e}")
36print(f"y(0)   = {y[0].item():.4f}    (expected {y0})")
37print(f"y(end) = {y[-1].item():.4f}   (expected {K})")

Real-World Applications

🌱 Population Biology (Bernoulli)

The logistic ODE is the workhorse of ecology — every bounded population, from yeast in a flask to wolves in Yellowstone, settles onto an S-curve toward its carrying capacity KK.

💨 Fluid Drag (Bernoulli, n = 2)

Newtonian drag gives mdv/dt=mgcv2m\,dv/dt = mg - c\,v^2. Dividing by v2v^2 and substituting u=v1u = v^{-1} turns terminal-velocity problems into linear ODEs.

🔬 Reaction Kinetics (Homogeneous)

Self-catalytic reactions A+B2BA + B \to 2B obey dB/dt=kABdB/dt = kAB. In dimensionless form (concentrations scaled by the total) the rate law is homogeneous of degree zero, solvable by v=B/Av = B/A.

🧠 Machine Learning (Bernoulli)

The replicator equation for evolutionary game theory, x˙i=xi(fi(x)fˉ(x))\dot x_i = x_i\bigl(f_i(x) - \bar f(x)\bigr), is a multi-species Bernoulli system and underlies modern reinforcement-learning dynamics.


Common Pitfalls

Forgetting the product rule

When you substitute y=vxy = vx, the derivative is y=v+xvy' = v + xv', not just vv'. Missing the vv term is the single most common mistake in homogeneous substitutions.

Bernoulli with n = 0 or n = 1

The substitution u=y1nu = y^{1-n} is degenerate when n=1n = 1 (gives u=y0=1u = y^0 = 1) or trivial when n=0n = 0. But in those cases the original equation is already linear — solve it directly with the integrating factor method.

Don't forget to back-substitute

After integrating the simplified equation in vv or uu, you MUST express the answer back in terms of the original variables yy and xx. A solution "in v" is not yet a solution to the original problem.

Signs in the linear-shift case

For v=ax+by+cv = ax + by + c with b=0b = 0, the substitution degenerates — the right side is then a function of xx alone and the ODE is already directly integrable.


Test Your Understanding


Summary

Substitution is the most strategic tool in the first-order ODE toolbox. With three patterns you can crack equations that look forbidding at first glance:

PatternSubstitutionBecomesHallmark
dy/dx = F(y/x)v = y/xseparable in v, xright side depends only on the ratio y/x
y' + Py = Qy^nu = y^(1−n)linear in urogue y^n term ruins linearity
dy/dx = f(ax+by+c)v = ax+by+cseparable in v, xright side groups x and y as a single block

Key Takeaways

  1. Substitution trades a difficult ODE for an easier one in a new variable — the rest of the work is back-substitution.
  2. Homogeneous: slopes are constant along rays from the origin → v=y/xv = y/x always works.
  3. Bernoulli: a single yny^n spoils linearity → kill it with u=y1nu = y^{1-n}.
  4. Linear shift: when ax+by+cax + by + c appears as a unit, name it vv.
  5. Numerical verification with SciPy and analytic verification with PyTorch autograd are cheap insurance against algebra mistakes — use them.
  6. Real-world S-curves, drag, and reaction kinetics are all logistic / Bernoulli equations in disguise.
The Substitution Principle:
"Don't solve the problem you have. Rename it so it becomes a problem you've already solved."
Coming Next: Section 21.4 turns these tools toward their most famous applications — exponential growth and decay. We'll meet radioactive isotopes, compound interest, and the half-life formula as one unified story.
Loading comments...