Boo-AI — Master Artificial Intelligence by Building from Scratch

The Question Behind the Equation

A 100-litre tank holds salt water. Fresh brine flows in at the top, and a stirrer keeps the contents perfectly mixed. A drain at the bottom releases the mixture at the same rate. How much salt is in the tank at every moment? Trivial-looking question — and yet it sits behind drug infusion, lake pollution, blood-oxygen dynamics, reservoir desalination, and every chemical reactor on the planet.

What makes mixing problems special is that the rate at which something leaves a container depends on how much is currently inside. That single feedback is what turns a one-line bookkeeping problem into a first-order linear differential equation — and what creates the universal $1 - e^{-t/\tau}$ curve you see in every textbook plot.

The mixing-tank ODE

Let $y(t)$ be the mass of dissolved substance at time $t$ . Inflow brings it in at concentration $c_{\text{in}}$ and rate $r_{\text{in}}$ . Outflow takes it out at rate $r_{\text{out}}$ at whatever concentration $y/V$ the tank currently has. Then:

\displaystyle \frac{dy}{dt} \;=\; r_{\text{in}}\,c_{\text{in}} \;-\; \frac{r_{\text{out}}}{V(t)}\,y(t)

When inflow and outflow rates are equal, $V$ stays constant and this becomes a linear ODE with constant coefficients — the kind we solved in $\S 21.1$ .

The Rate-In / Rate-Out Balance

Every conservation problem starts the same way. Pick a quantity. Account for everything that increases it. Subtract everything that decreases it. The rate of change of the quantity is the net of those two.

Quantity	Why it increases	Why it decreases
Salt in tank y(t)	Inflow carries brine of concentration c_in at rate r_in	Drain removes mixed solution of concentration y/V at rate r_out
Drug mass in blood	IV infusion at fixed mg/min	Kidneys clear at a rate proportional to current concentration
Pollutant in lake	Polluted stream feeding the lake	Outflow river carrying the well-mixed lake water

Notice the pattern: the in column is usually a constant per unit time, while the out column is proportional to the current amount. That asymmetry is the entire reason the solution involves an exponential — it is the same asymmetry that makes Newton's law of cooling, RC discharge, and radioactive decay all look alike.

Why outflow scales with y, not r_out alone

The drain doesn't care about the original salt: it pulls out $r_{\text{out}}$ litres per minute of whatever the tank currently is. Multiply by the current concentration $y/V$ to get the mass leaving per minute. That is the $(r_{\text{out}}/V)\,y$ term.

Building the ODE

Set up the conservation statement and translate it line-by-line into symbols. Assume for now that inflow and outflow rates are equal — call that common value $r$ — so the volume is constant at $V$ .

Rate in. Brine arrives at $r$ L/min carrying $c_{\text{in}}$ kg/L. Mass per minute entering: $r\,c_{\text{in}}$ .
Concentration in tank. Perfect mixing means the tank is uniform; concentration is $y(t)/V$ at every instant.
Rate out. Drain removes $r$ L/min at concentration $y/V$ . Mass per minute leaving: $(r/V)\,y$ .
Net rate. $dy/dt = r\,c_{\text{in}} - (r/V)\,y$ .

Rearranged into standard linear form,

\displaystyle \frac{dy}{dt} + \frac{r}{V}\,y \;=\; r\,c_{\text{in}}.

Compare with the canonical linear ODE $y' + p(t)\,y = q(t)$ : here $p = r/V$ is the constant decay rate and $q = r\,c_{\text{in}}$ is the constant forcing. Everything we did in $\S 21.1$ applies — and the integrating-factor recipe will give the answer in a few clean lines.

Solving the Linear ODE

The integrating factor is $\mu(t) = e^{(r/V)\,t}$ . Multiply both sides:

$\displaystyle \frac{d}{dt}\!\left[\,e^{(r/V)t}\,y\right] = r\,c_{\text{in}}\,e^{(r/V)t}.$
Integrate both sides in $t$ : $e^{(r/V)t}\,y = V\,c_{\text{in}}\,e^{(r/V)t} + C.$
Divide by $e^{(r/V)t}$ : $y(t) = V\,c_{\text{in}} + C\,e^{-(r/V)t}.$
Apply $y(0) = y_0$ : $y_0 = V\,c_{\text{in}} + C \;\Rightarrow\; C = y_0 - V\,c_{\text{in}}.$

Substitute the constant back and we have the full closed-form answer:

\displaystyle y(t) \;=\; V\,c_{\text{in}} \;+\; \bigl(y_0 - V\,c_{\text{in}}\bigr)\,e^{-(r/V)\,t}.

Sanity-check by differentiating

$y'(t) = -(r/V)\,(y_0 - V c_{\text{in}})\,e^{-(r/V)t}$ and $r c_{\text{in}} - (r/V) y(t) = r c_{\text{in}} - (r/V)\bigl[V c_{\text{in}} + (y_0 - V c_{\text{in}}) e^{-(r/V)t}\bigr] = -(r/V)(y_0 - V c_{\text{in}})e^{-(r/V)t}$ . The two match. ✓

Equilibrium and Time Constant

Two quantities deserve names because they describe the entire long-term behaviour of the tank.

Equilibrium

Set $dy/dt = 0$ in the ODE: $r\,c_{\text{in}} = (r/V)\,y_{\text{eq}}$ , which gives:

\displaystyle y_{\text{eq}} \;=\; V\,c_{\text{in}}.

At equilibrium, the tank's concentration matches the inflow concentration. That is the only steady state — and the only place the curve can flatten out.

Time Constant

The coefficient on $t$ in the exponential is $-r/V$ . Its reciprocal is the natural clock of the system:

\displaystyle \tau \;=\; \frac{V}{r}.

After one time constant the deviation from equilibrium has shrunk by a factor $1/e$ , i.e. you have closed $1 - 1/e \approx 63\%$ of the gap. After three time constants you are within $5\%$ . After five you are within $1\%$ — engineers casually call this "steady state".

Reading the solution. Rewrite $y(t) = y_{\text{eq}} + (y_0 - y_{\text{eq}})\,e^{-t/\tau}$ in words: start at $y_0$ , close the gap to $y_{\text{eq}}$ exponentially fast, with characteristic time $\tau$ . Every mixing-tank answer fits this template.

Interactive Tank Explorer

Drag the volume, flow rate, inflow concentration, and starting salt. The tank shades darker as concentration rises; the curve to the right plots $y(t)$ against time with the equilibrium line in red. The three coloured bars at the bottom show the rate balance at the current moment — note how the green (in) and red (out) bars start mismatched and gradually equalise as the curve approaches the red equilibrium line.

Loading tank mixing explorer…

What to play with first

Set $y_0 = 0$ and watch the tank charge up toward $V c_{\text{in}}$ . Now set $y_0 > V c_{\text{in}}$ and watch it discharge down. Both directions, same $\tau$ .
Triple $r$ . The equilibrium height doesn't change but the curve gets dramatically faster — because $\tau = V/r$ shrank.
Triple $V$ . Equilibrium triples (more tank = more salt at the same concentration) and $\tau$ triples (bigger tank reacts slowly).
Set $c_{\text{in}} = 0$ . Pure water in, drain pulls salt out. The tank empties exponentially: the classic decay equation falls out as a special case.

One Equilibrium, Many Starts

Different initial salt amounts all relax to the same long-time value. Below, the dashed red line is $y_{\text{eq}} = V c_{\text{in}}$ , and each coloured curve starts at a different $y_0$ . Watch how the spacing between curves shrinks exponentially — that is the geometric meaning of $e^{-t/\tau}$ .

Loading convergence visualizer…

A linear ODE's general solution always splits into steady-state plus transient. Here the steady state is $V c_{\text{in}}$ and the transient is $(y_0 - V c_{\text{in}})\,e^{-t/\tau}$ . The steady state forgets $y_0$ ; the transient remembers it. As $t \to \infty$ memory of the initial condition vanishes — that is what the exponential is doing.

When the Tank Is Filling or Draining

Drop the assumption $r_{\text{in}} = r_{\text{out}}$ and the volume itself becomes a function of time:

\displaystyle V(t) \;=\; V_0 + (r_{\text{in}} - r_{\text{out}})\,t.

The ODE keeps its linear form, but the coefficient on $y$ is no longer constant:

\displaystyle \frac{dy}{dt} + \frac{r_{\text{out}}}{V_0 + (r_{\text{in}} - r_{\text{out}})\,t}\,y \;=\; r_{\text{in}}\,c_{\text{in}}.

The recipe is still integrating factor — just integrate $r_{\text{out}}/V(t)$ instead of a constant. A clean example: tank starts with $V_0 = 100$ L at $y_0 = 10$ kg, brine of $c_{\text{in}} = 0.5$ kg/L flows in at $6$ L/min, drain runs at $4$ L/min. So $V(t) = 100 + 2t$ .

The integrating factor is $\mu(t) = \exp\!\left(\int \frac{4}{100 + 2t}\,dt\right) = (100 + 2t)^{2}$ . Multiplying through and integrating gives $(100+2t)^2\,y = (100+2t)^3/2 + C$ . Apply $y(0) = 10$ : $100^2 \cdot 10 = 100^3/2 + C \Rightarrow C = -400{,}000$ . The solution is:

\displaystyle y(t) \;=\; \frac{100 + 2t}{2} \;-\; \frac{400{,}000}{(100 + 2t)^{2}}.

What changes vs constant volume

With $V$ growing, the concentration $y/V$ still approaches $c_{\text{in}}$ , but $y$ itself can grow without bound because the tank keeps gaining volume. Conversely, if the tank drains ( $r_{\text{out}} > r_{\text{in}}$ ), there is a finite time at which the tank empties, and the ODE stops applying.

Worked Example: 100-Litre Brine Tank

Here is the canonical mixing problem. A tank holds $V = 100$ L of water with $y_0 = 10$ kg of dissolved salt. Brine of concentration $c_{\text{in}} = 0.5$ kg/L flows in at $r = 5$ L/min; the well-stirred mixture flows out at the same rate. Find $y(t)$ , the equilibrium, and the time at which the salt mass reaches 45 kg. Try it on paper first, then expand the panel.

Show step-by-step solution

Step 1. Write the ODE.

$dy/dt = 5 \cdot 0.5 - (5/100)\,y = 2.5 - 0.05\,y$ , with $y(0) = 10$ .

Step 2. Identify equilibrium and time constant.

Setting $dy/dt = 0$ : $y_{\text{eq}} = 2.5/0.05 = 50$ kg. $\tau = V/r = 100/5 = 20$ min.

Step 3. Write the closed-form solution.

$y(t) = 50 + (10 - 50)\,e^{-t/20} = 50 - 40\,e^{-t/20}$ .

Step 4. Evaluate at several times.

t (min)	Computation	y(t) (kg)
0	50 − 40·e^{0} = 50 − 40·1.0000	10.0000
10	50 − 40·e^{-0.5} = 50 − 40·0.6065	25.7388
20	50 − 40·e^{-1} = 50 − 40·0.3679	35.2848
40	50 − 40·e^{-2} = 50 − 40·0.1353	44.5866
60	50 − 40·e^{-3} = 50 − 40·0.0498	48.0085
100	50 − 40·e^{-5} = 50 − 40·0.0067	49.7305

Notice the entries don't step by a constant amount — the gap to 50 shrinks by a factor $e^{-0.5} \approx 0.6065$ every 10 minutes. After three time constants (60 min) the tank is within 2 kg of equilibrium.

Step 5. When does $y(t) = 45$ ?

Solve $45 = 50 - 40\,e^{-t/20}$ $\Rightarrow 40\,e^{-t/20} = 5 \Rightarrow e^{-t/20} = 0.125$ . Take logs: $-t/20 = \ln(0.125) = -2.0794$ , so $t = 41.589$ minutes.

That is about $2.08\,\tau$ — a useful rule for engineers: to close 90% of the gap to equilibrium takes a little over two time constants.

Step 6. Concentration in the tank.

At $t = 41.589$ minutes, the salt concentration is $y/V = 45/100 = 0.45$ kg/L. Compare with the inflow $c_{\text{in}} = 0.5$ kg/L — we are 90% of the way there.

Step 7. Verify on the explorer.

Set $V = 100$ , $r = 5$ , $c_{\text{in}} = 0.5$ , $y_0 = 10$ , marker at $t = 41.6$ min. The readout should show $y \approx 45.00$ kg and Rate IN $\approx$ 2.5 kg/min versus Rate OUT $\approx$ 2.25 kg/min — a small positive net flow still nudging the tank toward 50.

Application: IV Drug Infusion

Same equation, different costume. A drug is infused intravenously at constant rate $R_{\text{infuse}}$ (mg/min). The kidneys clear it at a rate proportional to the current blood concentration. Let $y(t)$ be the drug mass in the body and $V_d$ the apparent volume of distribution (a pharmacokinetic constant, typically tens of litres).

\displaystyle \frac{dy}{dt} \;=\; R_{\text{infuse}} \;-\; k_e\,y,

where $k_e$ is the elimination rate constant. The equilibrium drug mass is $y_{\text{eq}} = R_{\text{infuse}}/k_e$ and the steady-state plasma concentration is $C_{ss} = y_{\text{eq}}/V_d = R_{\text{infuse}}/(k_e V_d)$ . The time constant is $\tau = 1/k_e$ , and the elimination half-life is $t_{1/2} = (\ln 2)/k_e$ .

Clinical rule of thumb. A continuous IV reaches roughly 94% of its steady-state plasma concentration after $4\,t_{1/2}$ . For a drug with a 6-hour half-life, that is 24 hours of infusion before the plasma level stops climbing — which is exactly the kind of question this ODE answers.

Application: Pollutant in a Lake

A lake of volume $V$ is fed by an inflow stream and drained by an outflow river, both at the same volumetric rate $Q$ . Upstream activity contaminates the inflow at concentration $c_{\text{in}}$ . The ODE looks identical to the salt tank — only the names change:

\displaystyle \frac{dy}{dt} \;=\; Q\,c_{\text{in}} - \frac{Q}{V}\,y, \qquad \tau = \frac{V}{Q}.

Lake Erie has volume $\approx 480 \times 10^{9}$ m³ and an outflow of $\approx 5800$ m³/s, which works out to a time constant of about 2.6 years. So even if the polluting source is shut off completely, the lake takes roughly $3\tau \approx 8$ years to fall to 5% of its peak concentration — a sobering policy timeline that falls straight out of this one equation.

Python: Building the Simulator

Before trusting the closed form, let us reconstruct it from the ODE itself. The simulator below steps $dy = (r\,c_{\text{in}} - (r/V)\,y)\,dt$ forward many tiny steps and compares the result to the analytic answer at five checkpoints.

Hand-coded mixing simulator with ground-truth comparison

🐍mixing_scratch.py

Explanation(26)

Code(40)

1Import math

We need math.exp for the analytic closed-form solution. The simulator itself uses no special functions — just arithmetic — so this import is purely for the ground-truth comparison.

EXAMPLE

math.exp(-1) = 0.3678794411…

4Function mixing_step

One small Euler step of the ODE. Given the current salt mass y, tank volume V, flow rates and inlet concentration, it returns the updated (y, V) after a tiny time slice dt.

5Docstring locks the equation

The docstring spells out the exact ODE: rate-in minus rate-out, where rate-out is proportional to the current concentration y/V. Code without an explicit equation drifts; the docstring keeps the maths and the program aligned.

6Compute dy = (rate_in − rate_out) · dt

This single line is the entire physics. Salt enters at r_in·c_in kg/min. It leaves at (r_out/V)·y kg/min, because the outflow takes the same concentration y/V as the well-mixed tank. dy is the net change over one small interval.

EXAMPLE

y=10, V=100, r_in=r_out=5, c_in=0.5, dt=0.01 → dy = (5·0.5 − (5/100)·10)·0.01 = (2.5 − 0.5)·0.01 = 0.02 kg

7Compute dV = (r_in − r_out) · dt

If the inflow rate equals the outflow rate, dV is zero and the tank volume never changes. Keeping this term explicit lets us reuse the same code for filling tanks and draining tanks in the variable-volume case later.

EXAMPLE

r_in = r_out = 5 → dV = 0 every step, so V stays at 100 L.

8Return updated state

We don't mutate the inputs — we return fresh values. Pure functions like this are trivially testable: you can feed mixing_step any (y, V, …) and check the math by hand.

11Function simulate

Drives mixing_step in a loop. Each call advances time by dt. The output is a list of (t, y, V) triples — the discrete trajectory of the tank.

13Initial state

t, y, V are seeded with the initial conditions of the ODE. Without specifying these, the ODE has infinitely many solutions; with them, exactly one.

EXAMPLE

y0=10, V0=100 → we start the trajectory at (t=0, y=10 kg, V=100 L).

14history list

We keep every snapshot so plotting and inspection are easy later. Each entry is a tuple (t, y, V); we will index into this list to spot-check against the closed form.

15While loop with tiny tolerance

We march forward until t reaches t_final. The -1e-9 tolerance guards against floating-point drift: without it, adding 0.01 ten thousand times can leave t microscopically above 100.0 and skip the last write.

16One Euler step

Invokes the pure step function. y and V are simultaneously updated. The order matters subtly — we use the OLD (y, V) to compute the increments, then assign both new values. Modifying y first and then using the new y to compute dV would change the dynamics.

EXAMPLE

Step 1: y becomes 10.0200; V stays 100. Step 2: dy uses the new y=10.02, so dy = (2.5 − 0.05·10.02)·0.01 = 0.01999 kg.

17Advance time

t increases by dt after the state update. The first recorded pair is at t=0 (initial); the (i)-th update writes the state at t = i·dt.

18Append the new snapshot

The list grows by one row per step. With dt = 0.01 and t_final = 100, the loop runs 10 000 times and history ends with 10 001 entries.

19Return the full trajectory

Returning the history (not just the final state) costs almost nothing in memory and saves us from re-running the loop if we later want a plot or a different time slice.

22Function closed_form

The analytic solution we will derive in the next section: y(t) = y_eq + (y0 − y_eq) · exp(−k t). We use it as ground truth to grade our Euler steps.

24y_eq = V · c_in

The equilibrium salt mass: the steady-state value at which rate-in exactly matches rate-out. It does not depend on y0 — only on the inflow and tank size.

EXAMPLE

V=100, c_in=0.5 → y_eq = 50 kg.

25k = r / V

The decay constant of the relaxation toward equilibrium. r/V has units of 1/time. Its reciprocal V/r is the time constant τ — the natural clock of the tank.

EXAMPLE

r=5, V=100 → k = 0.05 per minute, τ = 20 min.

26Return the closed-form value

One line, no loop. The exponential captures the entire trajectory at once. For t = 0 it returns y0 + 0 = y0; as t → ∞ the exponential decays to 0 and the answer approaches y_eq. The (y0 − y_eq) factor is the 'distance to equilibrium' that shrinks over time.

EXAMPLE

y0=10, y_eq=50, k=0.05, t=20 → 50 + (10−50)·exp(−1) = 50 − 40·0.3679 = 35.2848 kg.

30Define the scenario

y0 = 10 kg salt initially. V0 = 100 L. Equal flow rates r = 5 L/min in and out keep the volume locked. The incoming brine carries c_in = 0.5 kg/L. Numbers chosen so the arithmetic in your head stays clean.

32Run the simulation

10 000 Euler steps of size dt = 0.01 min. Total simulated time: 100 minutes. Notice we pass r as both r_in and r_out — that is the constant-volume case the closed form covers.

34Header row for the table

f-strings with column-aligned widths produce a readable comparison. Right-aligning the numbers makes it obvious when Euler and the exact solution agree.

35Loop over check times

We sample five times: 0, 20, 40, 60, 100 — the start, three time constants, and well past equilibrium. They expose how Euler converges to the analytic curve.

36Index into trace by time

Since steps are uniform, the snapshot at t = t_check sits at position t_check/dt in the list. round() guards against fractional indices when t_check is not an exact multiple of dt.

EXAMPLE

t_check=20.0, dt=0.01 → idx = 2000. trace[2000] = (20.0, 35.2823, 100.0)

37Pull the Euler salt mass

history[idx] is a tuple (t, y, V). Index [1] grabs y. We ignore V because it is constant 100 for this run.

38Compute the exact reference

Calls closed_form with the same y0, V, r, c_in. This is the value Euler is trying to match.

EXAMPLE

t=20 → closed_form = 50 − 40·exp(−1) = 35.2848

39Print the row

When you run this you will see Euler within ~0.003 of exact at every t — proof that with a small enough dt the discrete update rule is faithful to the ODE. Bump dt to 5.0 and the discrepancy grows visibly.

EXAMPLE

Expected: t=20.0  Euler=35.2823  Exact=35.2848

14 lines without explanation

1import math
2
3
4def mixing_step(y, V, r_in, r_out, c_in, dt):
5    """One Euler step of  dy/dt = r_in * c_in - (r_out / V) * y."""
6    dy = (r_in * c_in - (r_out / V) * y) * dt
7    dV = (r_in - r_out) * dt
8    return y + dy, V + dV
9
10
11def simulate(y0, V0, r_in, r_out, c_in, t_final, dt):
12    """March the ODE forward in small Euler steps and record history."""
13    t, y, V = 0.0, y0, V0
14    history = [(t, y, V)]
15    while t < t_final - 1e-9:
16        y, V = mixing_step(y, V, r_in, r_out, c_in, dt)
17        t += dt
18        history.append((t, y, V))
19    return history
20
21
22def closed_form(y0, V, r, c_in, t):
23    """Exact answer when r_in = r_out = r (the volume never changes)."""
24    y_eq = V * c_in
25    k = r / V
26    return y_eq + (y0 - y_eq) * math.exp(-k * t)
27
28
29# A 100 L tank, 5 L/min flow each way, 0.5 kg/L brine,
30# starting with only 10 kg of salt dissolved.
31y0, V0, r, c_in = 10.0, 100.0, 5.0, 0.5
32
33trace = simulate(y0, V0, r, r, c_in, t_final=100.0, dt=0.01)
34
35print(f"{'t (min)':>8}  {'Euler y':>10}  {'Exact y':>10}")
36for t_check in (0.0, 20.0, 40.0, 60.0, 100.0):
37    idx = int(round(t_check / 0.01))
38    y_euler = trace[idx][1]
39    y_exact = closed_form(y0, V0, r, c_in, t_check)
40    print(f"{t_check:8.1f}  {y_euler:10.4f}  {y_exact:10.4f}")

Euler is the simplest viable stepper — and exactly what you would derive by replacing $dy/dt$ with $(y_{n+1} - y_n)/\Delta t$ . For production work you would reach for $\texttt{scipy.integrate.solve\_ivp}$ , but the ladder of understanding goes: ODE → discrete update → loop → closed form. Each rung confirms the next.

Python: Closed Form to Predict an Event

Once the analytic solution is in hand, every timing question becomes a one-line calculation — no loop required. Below we answer "at what time does the tank reach 90% of equilibrium?" for the canonical 100-litre setup.

Inverting the closed form to find a target time

🐍mixing_event_time.py

Explanation(14)

Code(26)

1Import math

We need math.log to invert the exponential. log here is the natural logarithm — the inverse of exp — and is the standard way to solve y_eq + Δ·exp(−kt) = target for t.

EXAMPLE

math.log(0.125) = -2.0794

4Tank parameters

Same brine tank as before. We will not run any simulation this time — the closed form lets us answer the timing question with three lines of arithmetic.

EXAMPLE

V=100 L, r=5 L/min, c_in=0.5 kg/L, y0=10 kg.

6y_eq = V · c_in

The steady-state salt mass. Driven only by the inflow concentration and the tank size; the initial salt drops out at long times.

EXAMPLE

y_eq = 100 · 0.5 = 50 kg.

7k = r / V

The decay rate of the deviation from equilibrium. r/V has units 1/min, so −k·t is dimensionless and exp(−kt) makes sense.

EXAMPLE

k = 5 / 100 = 0.05 per minute.

8τ = 1 / k = V / r

The time constant. After one τ, the gap to equilibrium has shrunk by a factor 1/e ≈ 0.368, i.e. you have closed 63% of the distance. After 3τ you are within 5% — a useful rule of thumb.

EXAMPLE

τ = V/r = 100/5 = 20 minutes.

11Choose the target fraction

We want the time at which the tank is at 90% of its long-term value. The choice of 0.90 is arbitrary — change it to 0.99 and ask 'how long until basically done?'.

12Convert fraction to a salt mass

target = 0.90 · y_eq = 0.90 · 50 = 45 kg. This is the salt mass we want to solve for.

19Compute the ratio inside the log

From y(t) = y_eq + (y0 − y_eq)·exp(−kt), we isolate the exponential. (target − y_eq) and (y0 − y_eq) are both negative when y0 < y_eq < target, so their ratio is positive and the log is well defined.

EXAMPLE

(target − y_eq)/(y0 − y_eq) = (45 − 50)/(10 − 50) = −5/−40 = 0.125

20Invert the exponential

Take the natural log, flip the sign, divide by k. The negative log of a number in (0,1) is positive, so t_hit comes out positive as expected.

EXAMPLE

t_hit = −ln(0.125)/0.05 = 2.0794/0.05 = 41.589 minutes.

22Print equilibrium

y_eq is what every long-time prediction will hover around. Reporting it first orients the reader: 'OK, the destination is 50 kg.'

EXAMPLE

Expected: equilibrium y_eq = 50.000 kg.

23Print time constant

20 minutes is the natural clock of this tank — every other timescale is reported as a multiple of it.

EXAMPLE

Expected: time constant tau = 20.000 min.

24Print target

Shows the absolute mass we are aiming for, not just the 0.90 fraction.

EXAMPLE

Expected: target (90% of y_eq) = 45.000 kg.

25Print the answer

41.589 minutes. Compare that to τ = 20: it takes a bit more than two time constants to reach 90% of equilibrium. The general rule is t/τ = ln(1/(1-frac)); plug frac=0.9 to get ln(10) ≈ 2.303 — but our numbers come out slightly less because we are not starting at zero.

EXAMPLE

Expected: time to reach target = 41.589 min.

26Report in units of τ

Dividing by τ makes the answer scale-free. Any tank with the same fractional-target question has the same t/τ — only the absolute time changes with V and r.

EXAMPLE

t_hit / tau = 41.589 / 20.000 = 2.079

12 lines without explanation

1import math
2
3# Constant-volume brine tank.
4V, r, c_in, y0 = 100.0, 5.0, 0.5, 10.0
5
6y_eq = V * c_in          # equilibrium salt mass (kg)
7k    = r / V             # relaxation rate (1/min)
8tau  = 1.0 / k           # time constant V/r (min)
9
10# Question: at what time does the tank reach 90% of equilibrium?
11frac   = 0.90
12target = frac * y_eq
13
14# y(t) = y_eq + (y0 - y_eq) * exp(-k t)
15# Solve for t:
16#   target - y_eq = (y0 - y_eq) * exp(-k t)
17#   exp(-k t) = (target - y_eq) / (y0 - y_eq)
18#   t = - ln( (target - y_eq) / (y0 - y_eq) ) / k
19arg   = (target - y_eq) / (y0 - y_eq)
20t_hit = -math.log(arg) / k
21
22print(f"equilibrium y_eq      = {y_eq:7.3f} kg")
23print(f"time constant tau     = {tau:7.3f} min")
24print(f"target (90% of y_eq)  = {target:7.3f} kg")
25print(f"time to reach target  = {t_hit:7.3f} min")
26print(f"that is {t_hit / tau:.3f} time constants")

The logarithm is the inverse of the exponential

Every timing question on this ODE reduces to the same template: $t = -\tau \ln\!\left(\dfrac{y_{\text{target}} - y_{\text{eq}}}{y_0 - y_{\text{eq}}}\right)$ . Plug in any target salt mass; if it lies between $y_0$ and $y_{\text{eq}}$ , the answer is a real positive number.

PyTorch: Recovering r/V from Noisy Data

In the lab you measure $y(t)$ with noise and want to fit the model. With a single ODE parameter the log trick works, but the moment the model gets more complex you want autograd to do the bookkeeping. Below: synthetic noisy samples, two free parameters ( $y_{\text{eq}}$ and $k = r/V$ ), MSE loss, gradient descent.

Fitting equilibrium and decay rate via Adam

🐍mixing_fit_pytorch.py

Explanation(21)

Code(35)

1Import torch

PyTorch handles the autograd machinery. We get tensors, gradient tracking, and an optimizer in a few lines — the heavy lifting of curve fitting becomes bookkeeping.

3Fix the random seed

torch.manual_seed(0) makes the noisy measurements reproducible. With the same seed, this script prints the same numbers on every machine.

6True parameters (the answer we are pretending not to know)

We construct synthetic data from V=100 L, r=5 L/min, c_in=0.5 kg/L, y0=10 kg. The fit will try to recover y_eq=50 and k=0.05 just from noisy y(t) samples.

7y_eq_true

V · c_in = 50 kg. We compute it once so the synthetic data is exactly consistent with a single linear-ODE solution.

8k_true

r/V = 0.05 per minute. This is the rate the fit must rediscover from noisy data.

1121 evenly spaced sample times

torch.linspace(0, 100, 21) gives t = 0, 5, 10, …, 100 — five-minute spacing over 100 minutes. Realistic for a lab measurement protocol.

EXAMPLE

t_data[0]=0, t_data[1]=5, …, t_data[20]=100

12Clean trajectory

Evaluate y(t) at every sample time using the exact closed form. This is the noiseless signal — what an infinitely precise instrument would measure.

EXAMPLE

y_clean[0] = 50 + (10−50)·exp(0) = 10. y_clean[4] = 50 + (10−40)·exp(−1) = 35.2848.

13Add Gaussian noise

torch.randn_like draws N(0,1) noise of the same shape as y_clean. Multiplying by 0.5 sets the standard deviation. The result mimics realistic instrument scatter.

EXAMPLE

y_noisy[4] ≈ 35.2848 + 0.5·(some small number) ≈ 35.5 or 34.9

16Initial guess for y_eq

We deliberately start at 80 — far from the truth (50). requires_grad=True tells PyTorch to track this tensor through every operation so we can ask for ∂loss/∂y_eq.

17Parametrize k via exp(log_k)

The decay rate must stay positive. Instead of constraining k directly, we let log_k be unconstrained and define k = exp(log_k). This is a standard trick that turns a constrained problem into an unconstrained one.

EXAMPLE

log_k=0.0 → k=exp(0.0)=1.0 initially. After training, log_k≈ln(0.05)=−3.0 → k≈0.05.

19Adam optimizer

Adam adapts the step size per parameter. lr=0.05 is generous because the loss landscape here is smooth and convex in y_eq. We pass both parameters in a list.

21Training loop: 801 steps

Each iteration of the loop does one forward pass, one loss computation, one backward pass, and one parameter update. 801 steps is enough to reach machine-precision agreement with the clean signal.

22Rebuild k from log_k

Inside the loop we re-evaluate k = exp(log_k). This stays inside the autograd graph so gradients flow from loss → k → log_k.

23Forward pass: the model prediction

Apply the exact-solution formula y_eq + (y0 − y_eq)·exp(−k·t) elementwise to t_data. PyTorch broadcasts over the 21-point tensor. This is our parametric model.

EXAMPLE

If y_eq=80, k=1.0 (initial), pred[0]=80+(10−80)·1=10; pred[5]=80+(−70)·exp(−25)≈80 — wildly wrong, exactly the kind of starting state Adam fixes.

24Mean-squared-error loss

Average of (pred − y_noisy)² across the 21 points. Squaring keeps the loss positive and penalises larger errors more. Mean (rather than sum) keeps lr scale-invariant in the number of samples.

26Zero the gradients

PyTorch ACCUMULATES gradients across backward passes. If we forget to zero them, the next loss.backward() would add to the previous gradients, and the optimizer would step on the wrong direction. zero_grad() is mandatory.

27Backward pass: compute ∂loss/∂y_eq and ∂loss/∂log_k

Autograd traces back through pred → exp(−k·t) → exp(log_k) and accumulates the partial derivatives into y_eq.grad and log_k.grad. We never write a derivative by hand.

28Optimizer step

Adam reads y_eq.grad and log_k.grad and updates the parameters. After ~200 steps y_eq is close to 50, after ~600 steps k is close to 0.05.

30Periodic progress print

Every 200 steps we print y_eq, k, and loss so we can watch convergence. Loss should drop monotonically (with small jitter) toward the irreducible noise floor of about (0.5)² ≈ 0.25.

31Format string

.item() unwraps the 0-dim tensor into a Python float. Width and decimal modifiers keep the columns aligned in the terminal.

35Final recovered values

After 800 steps, expect y_eq ≈ 49.9-50.1 and k ≈ 0.0498-0.0502 — both within noise. The exact closed form is a one-line miracle, but this gradient-descent fit works for any model where you can write the prediction in PyTorch — including ones with no closed form at all.

EXAMPLE

Expected: recovered y_eq ≈ 49.97  (true 50.0); recovered k ≈ 0.05002 (true 0.05).

14 lines without explanation

1import torch
2
3torch.manual_seed(0)
4
5# True (unknown to us) parameters
6V_true, r_true, c_in, y0 = 100.0, 5.0, 0.5, 10.0
7y_eq_true = V_true * c_in        # 50.0
8k_true    = r_true / V_true      # 0.05
9
10# Synthetic noisy measurements at 21 time points
11t_data  = torch.linspace(0.0, 100.0, 21)
12y_clean = y_eq_true + (y0 - y_eq_true) * torch.exp(-k_true * t_data)
13y_noisy = y_clean + 0.5 * torch.randn_like(y_clean)
14
15# Free parameters: y_eq (any real) and k (must stay positive -> use exp)
16y_eq   = torch.tensor(80.0, requires_grad=True)
17log_k  = torch.tensor(0.0,  requires_grad=True)
18
19optim = torch.optim.Adam([y_eq, log_k], lr=0.05)
20
21for step in range(801):
22    k    = torch.exp(log_k)
23    pred = y_eq + (y0 - y_eq) * torch.exp(-k * t_data)
24    loss = ((pred - y_noisy) ** 2).mean()
25
26    optim.zero_grad()
27    loss.backward()
28    optim.step()
29
30    if step % 200 == 0:
31        print(f"step {step:4d}  y_eq={y_eq.item():7.4f}  "
32              f"k={k.item():.5f}  loss={loss.item():.5f}")
33
34print(f"\nrecovered y_eq = {y_eq.item():.3f}  (true {y_eq_true})")
35print(f"recovered k    = {torch.exp(log_k).item():.5f}  (true {k_true})")

Why this generalises

Real reactors and biological compartments rarely fit a single-tank model. They might have time-varying inflow, multiple coupled compartments, or unknown delays. Closed-form solutions vanish, but the PyTorch fit template stays identical — define a model in forward(), pick a loss, call backward, step the optimizer. The mixing problem is the simplest possible instance of an ODE-based model-fitting pipeline that scales to entire pharmacokinetic and environmental models.

Common Pitfalls

Mixing units silently. If $r$ is in L/min and $c_{\text{in}}$ is in g/L, then $r\,c_{\text{in}}$ is in g/min — so $y$ must be in grams, not kilograms. Pick a unit system at the top of the problem and never change it.
Forgetting the well-mixed assumption. The whole derivation needs the outflow concentration to equal $y/V$ at every instant. Real tanks stratify; a layered tank obeys a different (much harder) PDE.
Confusing equilibrium with steady-state concentration. The equilibrium mass is $V c_{\text{in}}$ — proportional to volume. The equilibrium concentration is $c_{\text{in}}$ — independent of volume. Watch which one the question is asking about.
Treating variable-volume as constant-volume. When $r_{\text{in}} \neq r_{\text{out}}$ , the formula $y_{\text{eq}} = V c_{\text{in}}$ no longer applies — there is no equilibrium because $V$ never stops changing.
Sign errors when $y_0 > y_{\text{eq}}$ . The closed form still works — the transient is now negative and the curve falls toward equilibrium. The same $\tau$ sets the speed of falling and rising.

Summary

A well-mixed tank with rate-in $r_{\text{in}}\,c_{\text{in}}$ and rate-out $(r_{\text{out}}/V)\,y$ obeys the linear first-order ODE $dy/dt + (r_{\text{out}}/V)\,y = r_{\text{in}}\,c_{\text{in}}$ .
When inflow equals outflow, the closed-form solution is $y(t) = V c_{\text{in}} + (y_0 - V c_{\text{in}})\,e^{-t/\tau}$ , with time constant $\tau = V/r$ and equilibrium $y_{\text{eq}} = V c_{\text{in}}$ .
After one time constant the gap to equilibrium has shrunk by $1/e$ ; after three time constants the tank is within 5% of steady state.
When inflow and outflow rates differ, $V(t) = V_0 + (r_{\text{in}} - r_{\text{out}})\,t$ and the integrating factor is $V(t)^{r_{\text{out}}/(r_{\text{in}}-r_{\text{out}})}$ — still solvable in closed form for any constant flows.
The same equation governs IV drug infusion, lake pollution, continuous reactors, building HVAC loads, and many more. Only the names of $r$ , $V$ , and $c_{\text{in}}$ change.
Computationally, Euler gives the simplest discrete update; the closed form turns every timing question into a one-line log; and autograd plus gradient descent fits the model to real noisy measurements without ever writing a derivative by hand.