Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Recognize when a first-order DE is separable, i.e., can be written as $\frac{dy}{dx} = f(x)\,g(y)$ .
Execute the four-step procedure — separate, integrate, combine constants, apply the initial condition — to obtain a particular solution.
Explain why moving the differentials $dy$ and $dx$ across the equality is legitimate (the chain-rule justification, not just \"multiply both sides\").
Detect singular and lost solutions that the routine division by $g(y)$ can quietly throw away.
Verify a separable-equation solution by independent numerical integration (RK4) and by symbolic computation (SymPy / PyTorch autograd).

The Big Picture

We met direction fields in the previous section: a forest of tiny tangent segments that hint at every solution at once. Beautiful, but qualitative. Now we want actual closed-form solutions — formulas like $y(x) = y_0\,e^{-kx}$ that you can evaluate at any $x$ with a calculator. The first technique that achieves this is also the simplest, the oldest, and by far the most common in physics and engineering: separation of variables.

The idea is embarrassingly clean: if every appearance of $y$ and $dy$ can be herded onto one side, and every appearance of $x$ and $dx$ onto the other, then both sides are now ordinary single-variable integrals. Calculus we already know solves the rest.

One-line summary: Separation of variables turns the differential equation $\frac{dy}{dx} = f(x)\,g(y)$ into the integral equation $\int \frac{dy}{g(y)} = \int f(x)\,dx$ . After that, it's just antiderivatives.

Where the Trick Comes From

The technique is essentially as old as calculus itself. Leibniz, who invented the notation $dy/dx$ , deliberately designed it so the symbols on top and bottom behave as if they were ordinary numbers in many contexts. He used separation in 1691 to solve the catenary problem (the shape of a hanging chain). Johann Bernoulli, his student, used it to solve the brachistochrone problem in 1696 (the curve of fastest descent). Both problems reduce to a separable ODE after a clever substitution.

For about two centuries afterward, "solve the differential equation" mostly meant "try to massage it into a separable form." Even today, when you see an exponential, a Gaussian, a logistic curve, or Newton's law of cooling, you are looking at the output of a separation of variables that someone, somewhere, did once and recorded.

What Does "Separable" Mean?

A first-order ODE is separable if its right-hand side factors as a product of a function of $x$ alone and a function of $y$ alone:

\frac{dy}{dx} = f(x)\,g(y).

The factorisation is the whole game. Look at three quick examples:

Equation	f(x)	g(y)	Separable?
dy/dx = x·y	x	y	Yes
dy/dx = sin(x)/y	sin(x)	1/y	Yes
dy/dx = y² − x²	—	—	No (subtraction, not product)
dy/dx = e^{x+y}	e^x	e^y	Yes (because e^{x+y}=e^x·e^y)
dy/dx = x + y	—	—	No (sum of two terms, neither factoring)

Trick to memorize: A sum on the right-hand side is the death of separability. A product is the friend of separability. When in doubt, try to write the RHS as a single multiplicative chunk.

The Separation Procedure

Once the equation is in the form $\frac{dy}{dx} = f(x)\,g(y)$ , you carry out exactly four moves.

Step 1 — Separate

Divide both sides by $g(y)$ , multiply both sides by $dx$ :

\frac{1}{g(y)}\,dy = f(x)\,dx.

Step 2 — Integrate both sides

The left side is a $y$ -integral; the right is an $x$ -integral. Crucially they are now independent:

\int \frac{1}{g(y)}\,dy = \int f(x)\,dx.

Step 3 — Combine the constants

Each integral introduces its own constant; merge them into a single $C$ . After integration you have an implicit equation of the form

G(y) = F(x) + C,

where $G$ is an antiderivative of $1/g$ and $F$ is an antiderivative of $f$ . Sometimes this implicit form is the end of the road (e.g., when $G^{-1}$ has no elementary expression). Often, though, you can finish:

Step 4 — Solve for y, then apply the initial condition

Algebraically invert to get $y(x) = G^{-1}(F(x) + C)$ . Then substitute the initial value $y(x_0) = y_0$ and solve for the single unknown $C$ . The result is the unique solution that passes through the prescribed point.

Mental rhythm: separate, integrate, combine, pin down. Four words, four steps. Memorise them and 80 % of first-order ODEs in any physics textbook will fall.

Why Moving Differentials Like Symbols Is Allowed

At step 1 you wrote $\frac{1}{g(y)}\,dy = f(x)\,dx$ by treating $dy/dx$ as a quotient of differentials. Is that rigorous, or just a happy notation? Let's check it cleanly.

Suppose $y(x)$ is a solution. Define $H(y) = \int \frac{1}{g(y)}\,dy$ , so that $H'(y) = 1/g(y)$ . Then by the chain rule:

\frac{d}{dx}\,H(y(x)) \;=\; H'(y(x))\,\frac{dy}{dx} \;=\; \frac{1}{g(y)}\cdot f(x)\,g(y) \;=\; f(x).

So $H(y(x))$ has derivative $f(x)$ with respect to $x$ , which means it is an antiderivative of $f$ . That is exactly the conclusion you get from $\int H'(y)\,dy = \int f(x)\,dx$ . The Leibniz manipulation is a shortcut for this chain-rule argument — never a sleight of hand.

Read this twice. The reason every calculus textbook tells you to "move

dy

across the equals sign" without explanation is that the rigorous justification is a one-line chain-rule calculation. Once you've seen it, the move stops feeling like cheating.

Interactive Lab: Direction Fields and Solutions

Pick a separable DE from the panel, drag the parameter slider, and click anywhere on the plane to plant an initial condition. The thin coloured segments are the slope field (tangents to every solution); the thick coloured curve is the actual solution computed by a black-box RK4 integrator. The bracketed steps on the right show the symbolic separation evolving in real time.

Loading interactive separable-equations lab…

What to notice: with

dy/dx = ky

, every click produces a curve that runs perfectly along the slope arrows. The arrows are drawn from the equation; the curve is drawn from the numerical integrator; the formula on the right is drawn from the closed form. All three agree because the separation algebra is correct.

Interactive: Watch a Separation Step by Step

Below, each press of the play button advances the symbolic manipulation by exactly one move. Use the per-step pills above to jump around, or step manually with the arrow buttons.

Loading step-by-step animator…

Worked Example: Newton's Cooling by Hand

A barista pulls an espresso shot at $T_0 = 95^{\circ}\text{C}$ into a room at $T_{\mathrm{env}} = 20^{\circ}\text{C}$ . Newton's law of cooling models the drink as

\frac{dT}{dt} = -k\,(T - T_{\mathrm{env}}),\quad k = 0.05\,\mathrm{min}^{-1}.

Find $T(t)$ and the time at which the coffee reaches a drinkable $30^{\circ}\text{C}$ .

Show full hand-computation (try it yourself first, then peek)

Step 1 — substitute $u(t) = T(t) - T_{\mathrm{env}}$ . Since $T_{\mathrm{env}}$ is constant, $du/dt = dT/dt$ , and the equation becomes

\frac{du}{dt} = -k\,u.

Step 2 — separate. Divide by $u$ , multiply by $dt$ :

\frac{1}{u}\,du = -k\,dt.

Step 3 — integrate both sides.

\int \frac{1}{u}\,du = \int -k\,dt \;\Longrightarrow\; \ln|u| = -k\,t + C.

Step 4 — exponentiate and absorb signs into $A = \pm e^{C}$ :

u(t) = A\,e^{-k\,t}.

Step 5 — apply the initial condition. At $t=0$ , $u(0) = T_0 - T_{\mathrm{env}} = 95 - 20 = 75$ . So $A = 75$ .

u(t) = 75\,e^{-0.05\,t} \;\Longrightarrow\; T(t) = 20 + 75\,e^{-0.05\,t}.

Step 6 — when does T = 30 °C?

30 = 20 + 75\,e^{-0.05\,t} \;\Longrightarrow\; e^{-0.05\,t} = \tfrac{10}{75} \;\Longrightarrow\; t = -\tfrac{1}{0.05}\,\ln\tfrac{10}{75}.

Compute it numerically: $t \approx 40.30\ \text{min}$ .

Sanity table:

t (min)	T(t) (°C)	Plain English
0	95.00	Just poured
10	65.49	Still scalding
30	36.73	Comfortable sip
40.30	30.00	Target reached
60	23.73	Lukewarm
120	20.19	Essentially room temperature

Run the Python block below: every row in this table reproduces to four decimal places.

A Subtler Example: Logistic Growth

Some separable equations require partial fractions before they can be integrated. The classic example is the logistic equation, which models any population whose growth slows as it approaches a carrying capacity $K$ :

\frac{dy}{dx} = y\left(1 - \frac{y}{K}\right).

Show the full derivation (we'll use it in the visualizer above)

Step 1 — separate.

\frac{K}{y(K - y)}\,dy = dx.

Step 2 — partial fractions. Write $\frac{K}{y(K - y)} = \frac{A}{y} + \frac{B}{K - y}$ and solve: $A = 1$ , $B = 1$ . So

\frac{1}{y}\,dy + \frac{1}{K - y}\,dy = dx.

Step 3 — integrate.

\ln|y| - \ln|K - y| = x + C \;\Longrightarrow\; \ln\left|\tfrac{y}{K - y}\right| = x + C.

Step 4 — exponentiate. With $A = \pm e^{C}$ :

\frac{y}{K - y} = A\,e^{x}.

Step 5 — solve for y. Cross-multiply, expand, and collect $y$ :

y = (K - y)\,A\,e^{x} \;\Longrightarrow\; y\,(1 + A\,e^{x}) = K\,A\,e^{x}.

y(x) = \frac{K\,A\,e^{x}}{1 + A\,e^{x}} = \frac{K}{1 + A^{-1}\,e^{-x}}.

Let $B = A^{-1}$ . Then $y(x) = \dfrac{K}{1 + B\,e^{-x}}$ .

Step 6 — initial condition. At $x = 0$ : $y_0 = K/(1 + B) \Rightarrow B = (K - y_0)/y_0$ .

With $K = 1000$ and $y_0 = 10$ , $B = 99$ and the population takes about $x = \ln 99 \approx 4.60$ units to reach the inflection point $K/2 = 500$ .

x	y(x)	Behaviour
0	10.00	Initial seed
1	26.72	Still slow
3	168.66	Take-off
5	599.86	Past the inflection
7	917.20	Saturating
10	995.53	Practically at K

A Catalog of Common Separable DEs

These five appear in almost every applied-mathematics course. Pattern-recognise them and you save a lot of derivation time.

Equation	Separated form	Closed-form solution	Where it appears
dy/dx = ky	dy/y = k dx	y = y₀ e^{k(x-x₀)}	Exponential growth/decay, RC circuits, half-lives
dT/dt = -k(T - Tenv)	du/u = -k dt (u = T-Tenv)	T = Tenv + (T₀ - Tenv)e^{-kt}	Newton's cooling, heat-exchangers
dy/dx = y(1 - y/K)	K dy/[y(K-y)] = dx	y = K / (1 + ((K-y₀)/y₀)e^{-x})	Population dynamics, learning curves, S-shaped adoption
dy/dx = -x/y	y dy = -x dx	x² + y² = R²	Circular orbits, perpendicular trajectories
dy/dx = ky(1 - y)·... etc.	(partial fractions)	various	Chemical kinetics, epidemic SIR

Real-World Applications

Radioactive decay. $dN/dt = -\lambda N$ gives $N(t) = N_0\,e^{-\lambda t}$ ; the half-life is $t_{1/2} = \ln 2 / \lambda$ .
RC circuit discharge. $dV/dt = -V/(RC)$ separates the same way; the time constant is $\tau = RC$ .
Chemical kinetics. First-order reactions $d[A]/dt = -k[A]$ . Second-order reactions like $d[A]/dt = -k[A]^{2}$ also separate, giving $1/[A] = kt + 1/[A]_0$ .
Mixing problems. Concentration in a stirred tank with inflow and outflow: $dC/dt = (c_{\mathrm{in}} - C)\,q/V$ separates after substituting $u = c_{\mathrm{in}} - C$ .
Logistic populations & epidemics. The SI model and the Verhulst growth model are both separable after partial fractions.

Connection to Machine Learning

Separable equations sneak into machine learning more often than people realise:

Exponential learning-rate schedules. The schedule $d\eta/dt = -\lambda\,\eta$ gives $\eta(t) = \eta_0\,e^{-\lambda t}$ , which is the continuous-time analogue of multiplying $\eta$ by a fixed factor each epoch.
Diffusion models. The forward process is governed by an SDE whose deterministic skeleton is separable in $t$ ; the closed-form marginals come from this separation.
Population-style RL. Replicator dynamics in evolutionary algorithms look like $d p_i/dt = p_i (f_i - \bar f)$ — separable, with logistic-style solutions.
Continuous-time normalising flows. Whenever the ODE is $dz/dt = a(t)\,z$ , separation gives $z(t) = z_0\,\exp\!\int_0^{t} a(s)\,ds$ — the basis of time-conditional Gaussian flows.

Python: Solving a Separable DE Numerically

Code is the second proof. Below, plain Python implements the same Newton's- cooling problem two independent ways: directly from the closed-form solution $y_0\,e^{-k x}$ , and via a hand-written RK4 integrator that knows nothing about the closed form. We then compute the maximum disagreement. If the separation algebra is correct, the gap is tiny round-off.

Closed-form vs. RK4: verifying the separation algebra

🐍python

Explanation(21)

Code(46)

1Import math

We need math.exp for the closed-form curve. We deliberately avoid numpy/scipy so every step of the integrator is visible as plain Python.

EXAMPLE

math.exp(0.6)  → 1.8221

3Function signature

All five inputs are the things you would write on a board when stating the IVP: initial value y0, decay rate k, the x-interval [x_start, x_end], and how many integration steps to use.

EXAMPLE

solve_separable_decay(95.0, 0.05, 0.0, 60.0, 600)

4Docstring

Spells out exactly what "two ways" means: closed-form (from separation) and a black-box RK4 that does not know the answer. The point of the function is to verify the separation algebra by independent computation.

11Step size h

Total interval length divided by number of steps. For x in [0, 60] with n_steps=600 we get h = 0.1. Smaller h → smaller numerical error.

EXAMPLE

h = (60 - 0) / 600 = 0.1

12x grid

Build the list of x values we are going to evaluate: x_start, x_start+h, x_start+2h, …, x_end. Length is n_steps+1.

EXAMPLE

xs[:3] = [0.0, 0.1, 0.2]

15Closed-form curve

This is the formula we derived by separation: y(x) = y0 · exp(-k·(x - x_start)). For each x in our grid we evaluate it directly. No iteration.

EXAMPLE

y0=95, k=0.05, x=10  →  95 · e^{-0.5} = 57.6149 if x_start=0

18Define the right-hand side f(x, y)

f returns dy/dx at the point (x, y). For our DE the answer is -k·y; x is unused but we keep the signature general so the same RK4 loop works for any first-order DE.

EXAMPLE

f(0, 95) = -0.05 * 95 = -4.75

21Initialize the numerical trajectory

Start the integrator at the initial condition. y will be updated step by step; numerical will collect the running value.

EXAMPLE

numerical = [95.0]

22Set the starting point

Two scalar trackers: y is the current value of the solution, x is the current independent variable.

24RK4 stage k1

k1 is the slope at the start of the step — exactly what Euler's method would use. Numerically: f(x, y) = -0.05 · 95 = -4.75 at the very first iteration.

EXAMPLE

k1 = -4.75

25RK4 stage k2

k2 is the slope at the midpoint of the step, using a half-step Euler prediction from k1. This is what makes RK4 fourth-order accurate instead of first-order.

EXAMPLE

k2 = f(0.05, 95 + 0.05·(-4.75)) = f(0.05, 94.7625) = -4.7381

26RK4 stage k3

Another midpoint slope, but using k2 as the predictor. RK4 averages two different midpoint estimates to cancel error.

EXAMPLE

k3 ≈ -4.7382

27RK4 stage k4

Slope at the end of the step, predicted from k3. Together (k1, k2, k3, k4) form Simpson's rule on the slope.

EXAMPLE

k4 ≈ -4.7263

28Combine the four slopes

The classical RK4 update: y_new = y + (h/6)·(k1 + 2k2 + 2k3 + k4). The coefficients (1, 2, 2, 1)/6 are Simpson's-rule weights.

EXAMPLE

y_new ≈ 95 + (0.1/6) · (-4.75 + 2·(-4.7381) + 2·(-4.7382) + (-4.7263)) ≈ 94.5262

29Advance x

Move the x cursor by one step h. Now we are ready to compute the next slope at the new (x, y).

30Record the new value

Append the freshly computed y so that the returned list lines up index-for-index with xs.

33Maximum absolute error

Over every grid point, compare the closed-form value to the RK4 value. If the separation algebra was correct AND RK4 is accurate, this number is tiny.

EXAMPLE

Typically max_err < 1e-8 for h=0.1 over [0, 60]

34Return everything

Return the x grid, both y arrays, and the max error so the caller can plot or print them.

37Driver block

Run a concrete instance: a cup of coffee at 95 °C cooling toward an environment at 0 °C (the substitution u = T − T_env reduces Newton's cooling to plain exponential decay). k = 0.05 /min, run for an hour, 600 steps.

41Spot-check printouts

Print closed-form and RK4 side by side at a few representative times. They should agree to ~4 decimal places for h=0.1.

EXAMPLE

y(0)  exact=95.0000   rk4=95.0000
y(10) exact=57.6150   rk4=57.6150
y(30) exact=21.1939   rk4=21.1939
y(60) exact=4.7269    rk4=4.7269

45Confirm the error bound

The final line prints max_err. For these settings you should see something like 8e-11 — well below any practical concern.

25 lines without explanation

1import math
2
3def solve_separable_decay(y0, k, x_start, x_end, n_steps):
4    """
5    Solve dy/dx = -k * y on [x_start, x_end] two ways:
6      1. The closed-form formula y(x) = y0 * exp(-k * (x - x_start))
7      2. A hand-rolled numerical integrator (RK4) that knows nothing
8         about the closed form -- so any agreement is honest evidence
9         that the separation algebra was correct.
10    """
11    h = (x_end - x_start) / n_steps     # step size
12    xs = [x_start + i * h for i in range(n_steps + 1)]
13
14    # --- exact (from separation of variables) -------------------------------
15    exact = [y0 * math.exp(-k * (x - x_start)) for x in xs]
16
17    # --- numerical (RK4) ----------------------------------------------------
18    def f(x, y):                        # the right-hand side of the DE
19        return -k * y                   # x is unused but kept for generality
20
21    numerical = [y0]
22    y = y0
23    x = x_start
24    for _ in range(n_steps):
25        k1 = f(x,           y)
26        k2 = f(x + h / 2,   y + h * k1 / 2)
27        k3 = f(x + h / 2,   y + h * k2 / 2)
28        k4 = f(x + h,       y + h * k3)
29        y  = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
30        x  = x + h
31        numerical.append(y)
32
33    # --- compare ------------------------------------------------------------
34    max_err = max(abs(a - b) for a, b in zip(exact, numerical))
35    return xs, exact, numerical, max_err
36
37
38if __name__ == "__main__":
39    xs, exact, num, err = solve_separable_decay(
40        y0=95.0, k=0.05, x_start=0.0, x_end=60.0, n_steps=600,
41    )
42    print(f"y(0)  exact={exact[0]:.4f}   rk4={num[0]:.4f}")
43    print(f"y(10) exact={exact[100]:.4f}   rk4={num[100]:.4f}")
44    print(f"y(30) exact={exact[300]:.4f}   rk4={num[300]:.4f}")
45    print(f"y(60) exact={exact[600]:.4f}   rk4={num[600]:.4f}")
46    print(f"max absolute error over the run: {err:.2e}")

PyTorch: Symbolic vs. Autograd Verification

SymPy can carry out separation of variables symbolically, and PyTorch's autograd can independently differentiate the closed-form expression to confirm it satisfies the DE. The two together are about as airtight as a one-page proof.

SymPy solves the ODE; PyTorch autograd verifies the result

🐍python

Explanation(23)

Code(34)

1SymPy import

SymPy is the canonical Python library for symbolic mathematics. It can perform the same separation-of-variables steps we did by hand and return a sympy expression.

2PyTorch import

We use PyTorch's autograd engine purely as an independent way to differentiate the closed-form solution and confirm it satisfies the DE.

3math import

Used only for tiny numeric sanity checks if you extend the script.

6Symbol x

Declare x as a symbolic real variable. SymPy treats it as a formal symbol, not a number, so derivatives and integrals stay exact.

EXAMPLE

type(x)  →  <class 'sympy.Symbol'>

7Symbolic function y(x)

y is declared as a function (not just a symbol) because the DE contains the derivative dy/dx, and we need SymPy to know that y depends on x.

8Symbol k

Declare the rate constant k as positive — this restricts SymPy's case-analysis when it integrates 1/y (the sign of y matters for ln|y|).

10Build the ODE

sp.Eq writes the equation dy/dx = -k·y exactly as you would on a board. y(x).diff(x) is SymPy's symbolic derivative.

EXAMPLE

ode  →  Eq(Derivative(y(x), x), -k*y(x))

11Solve with initial condition y(0) = 95

sp.dsolve internally performs separation of variables on this equation: it splits the variables, integrates both sides, combines the constants, exponentiates, and plugs in the initial condition. The ics keyword does the final "apply IC" step automatically.

EXAMPLE

solution.rhs  →  95*exp(-k*x)

12Print the symbolic solution

solution.rhs is the right-hand side y(x). SymPy returns exactly what we derived by hand: 95·exp(-k·x). This is independent confirmation of the algebra.

15Pick a numeric k

Switch to numerics for the PyTorch half. k = 0.05 matches the Newton's-cooling example used in the Python code block above.

16Tensor of x values, requires_grad=True

Build a tensor of 7 evenly spaced x values from 0 to 60. requires_grad=True tells autograd to track operations on x so we can differentiate w.r.t. it later.

EXAMPLE

x_t  →  tensor([ 0., 10., 20., 30., 40., 50., 60.], requires_grad=True)

17Apply the closed-form formula

y_t is computed from x_t by the formula 95·exp(-k·x). Because x_t has requires_grad=True and we only use differentiable ops, autograd has built a graph from x_t to y_t.

EXAMPLE

y_t  →  tensor([95.0000, 57.6150, 34.9342, 21.1862, 12.8462, 7.7888, 4.7237], grad_fn=…)

20Call torch.autograd.grad

Ask autograd: "what is the derivative of (sum of y_t) with respect to x_t?" Since each y_t[i] depends only on x_t[i], the resulting vector is exactly dy/dx evaluated at each x[i].

21outputs = y_t.sum()

torch.autograd.grad needs a scalar output. Summing y_t gives a scalar whose gradient w.r.t. x_t equals the element-wise dy/dx (because the i-th term contributes only to x_t[i]).

22inputs = x_t

The variable we are differentiating with respect to. autograd will return a tensor of the same shape as x_t.

23create_graph = False

We do NOT need to differentiate the result again, so do not build a second-order graph. Setting True here would let you compute higher-order derivatives at extra cost.

24Indexing [0]

torch.autograd.grad returns a tuple (one entry per input). We have one input, so we unpack the first element.

27Compute the residual of the DE

The original DE is dy/dx = -k·y, equivalently dy/dx + k·y = 0. If the closed-form is correct, this expression must be (numerically) zero at every x.

EXAMPLE

residual  →  tensor([0., 0., 0., 0., 0., 0., 0.])   (up to float precision)

29Print x grid

Detach before listing — we only want the values, not the autograd graph. Output: [0.0, 10.0, …, 60.0].

30Print y values

Closed-form temperatures at each time. Same numbers we expect from the Newton's-cooling worked example.

EXAMPLE

[95.0, 57.6150, 34.9342, 21.1862, 12.8462, 7.7888, 4.7237]

31Print dy/dx values

Autograd's verdict for dy/dx. Compare against -k·y manually if you wish.

EXAMPLE

[-4.75, -2.8808, -1.7467, -1.0593, -0.6423, -0.3894, -0.2362]

32Print residual

Each entry should be on the order of 1e-7 or smaller. If you see anything large, the closed-form is wrong, OR the substitution into the DE is wrong.

34Print max |residual|

A single number summarizing the worst-case violation across all sampled x values. For correctly separated DEs this prints something like 1.9e-6 (pure float round-off).

EXAMPLE

max |residual| = 1.9073e-06

11 lines without explanation

1import sympy as sp
2import torch
3import math
4
5# --- 1. SYMBOLIC: let SymPy solve the same DE by separation -----------------
6x = sp.Symbol("x", real=True)
7y = sp.Function("y")
8k = sp.Symbol("k", positive=True)
9
10ode      = sp.Eq(y(x).diff(x), -k * y(x))                # dy/dx = -k y
11solution = sp.dsolve(ode, y(x), ics={y(0): 95})           # apply IC at x=0
12print("symbolic solution:", solution.rhs)
13
14# --- 2. NUMERICAL VERIFICATION using PyTorch + autograd ---------------------
15k_val = 0.05
16x_t   = torch.linspace(0.0, 60.0, 7, requires_grad=True) # 0, 10, 20, …, 60
17y_t   = 95.0 * torch.exp(-k_val * x_t)                   # plug into formula
18
19# Use autograd to compute dy/dx at the same x's
20dy_dx = torch.autograd.grad(
21    outputs      = y_t.sum(),    # sum to get a scalar
22    inputs       = x_t,
23    create_graph = False,
24)[0]
25
26# The DE says dy/dx should equal -k * y.  Check it:
27residual = dy_dx + k_val * y_t   # should be ≈ 0 elementwise
28
29print("x        ", x_t.detach().tolist())
30print("y(x)     ", [round(v, 4) for v in y_t.detach().tolist()])
31print("dy/dx    ", [round(v, 4) for v in dy_dx.detach().tolist()])
32print("residual ", [round(v, 1e-12) if False else round(v, 8)
33                    for v in residual.detach().tolist()])
34print("max |residual| =", residual.abs().max().item())

What this buys you in research: any time you have a candidate closed-form solution to a first-order DE — whether you derived it by separation, guessed it, or pulled it from a paper — you can verify it in three lines of PyTorch. The residual

y' - f(x, y)

should be machine-zero at every sampled point.

Common Pitfalls

Pitfall 1 — Dividing by zero (lost equilibrium solutions). Step 1 requires dividing by

g(y)

. If

g(y_*) = 0

for some constant

y_*

, then

y(x) \equiv y_*

is also a solution (the derivative is zero on both sides) — but you cannot recover it by plugging values of

C

into your formula. For the logistic equation, this means

y \equiv 0

and

y \equiv K

are real equilibrium solutions that the separation formula does not produce.

Pitfall 2 — Sign of the absolute value. Step 3 produces

\ln|y|

, not

\ln y

. After exponentiating, write

A = \pm e^{C}

, then let the initial condition pick the sign.

Pitfall 3 — Domain of validity. The implicit solution

G(y) = F(x) + C

is only locally a function. The hyperbolas example,

dy/dx = -x/y

, gives circles

x^2 + y^2 = R^2

; selecting the upper or lower semicircle requires committing to a sign of

y_0

Pitfall 4 — Forgetting that "separable" is a structural property, not a calculation. Many DEs become separable only after a substitution (e.g.

v = y/x

turns any homogeneous DE into a separable one in

v

). If a problem does not separate at first glance, consider whether a substitution opens it up.

Sanity protocol you can apply to any candidate solution: plug it back into the original DE and check that both sides agree. The Python and PyTorch blocks above are the automated form of this check.

Summary

A first-order DE is separable iff its right-hand side factors as $f(x)\,g(y)$ — a product, never a sum.
The four-step recipe — separate, integrate, combine constants, apply IC — turns the DE into two ordinary integrals.
The Leibniz move "multiply both sides by $dx$ " is rigorously justified by the chain rule, not by ambiguity of notation.
Watch for lost equilibrium solutions when $g(y_*) = 0$ , and for domain restrictions when inverting an implicit solution.
Numerical verification with RK4 and symbolic verification with SymPy / PyTorch autograd give you an independent, automatic sanity check that any candidate solution truly satisfies the original equation.

Looking ahead: Separation is the workhorse but it only handles first-order equations whose right-hand side factors. The next section introduces Euler's method: a numerical scheme that handles arbitrary first-order DEs — even non-separable ones — by following the direction field one tiny step at a time.