Module 1.1
| Difficulty: | Beginner |
| Estimated Time: | 4–5 hours |
First-order ordinary differential equations are the simplest differential equations, yet they model an extraordinary range of phenomena. A single first-order ODE captures the essence of radioactive decay, population growth with carrying capacity, thermal equilibration, and chemical mixing—each governed by a rate of change that depends on the current state. Mastering separable, linear, and exact equations gives you an analytical toolkit that applies directly to engineering design, pharmacokinetics, ecology, and financial modeling. Perhaps more importantly, the existence-uniqueness theorem (Picard–Lindelöf) tells you when you can trust that a model has a well-defined solution and when multiple solutions may coexist, preventing you from blindly trusting a computer output that may be selecting one branch among many.
This module builds the algebraic and geometric foundations you will use in every subsequent module: integrating factors reappear in systems of ODEs, separation of variables returns in PDE methods, and the interplay between analytical and numerical solutions becomes a recurring theme through the entire curriculum.
After completing this module, you will be able to:
A first-order ODE is separable if it can be written in the form
so that the variables can be algebraically separated: \(\displaystyle \frac{dy}{g(y)} = f(x)\,dx\). Integration of both sides then yields the general solution (possibly in implicit form).
A first-order ODE is linear if it has the standard form
where \(p\) and \(g\) are continuous functions on an interval \(I\). The equation is homogeneous when \(g(t) \equiv 0\) and non-homogeneous otherwise.
For the linear ODE \(y' + p(t)\,y = g(t)\), the integrating factor is
Multiplying both sides by \(\mu(t)\) transforms the left side into the exact derivative \(\frac{d}{dt}[\mu(t)\,y]\), enabling direct integration.
An equation \(M(x,y)\,dx + N(x,y)\,dy = 0\) is exact on a simply connected region if there exists a function \(F(x,y)\) such that
A necessary and sufficient condition (when \(M\) and \(N\) are continuously differentiable) is \(\displaystyle \frac{\partial M}{\partial y} = \frac{\partial N}{\partial x}\).
An initial-value problem for a first-order ODE consists of the equation together with a prescribed value at a specific point:
The pair \((t_0, y_0)\) is the initial condition. A solution is a function \(y(t)\) satisfying both the ODE and the initial condition on some interval containing \(t_0\).
An equilibrium (or steady-state) solution of \(y' = f(t,y)\) is a constant function \(y(t) = y^*\) satisfying \(f(t, y^*) = 0\) for all \(t\). Equilibria are classified as stable (nearby solutions converge), unstable (nearby solutions diverge), or semi-stable.
Let \(f(t,y)\) be continuous on an open rectangle \(R = \{(t,y) : |t - t_0| < a,\; |y - y_0| < b\}\) and let \(f\) satisfy a Lipschitz condition in \(y\) on \(R\): there exists \(L > 0\) such that
for all \((t,y_1), (t,y_2) \in R\). Then the initial-value problem
has a unique solution \(y(t)\) on some interval \(|t - t_0| < h\) where \(h = \min\!\bigl(a,\, b/M\bigr)\) and \(M = \max_{R} |f|\).
If \(p(t)\) and \(g(t)\) are continuous on an open interval \(I\) containing \(t_0\), then the IVP
has the unique solution on \(I\):
In particular, because \(p\) and \(g\) are continuous on all of \(I\), the solution exists on the entire interval, not just locally. This is a stronger conclusion than Picard–Lindelöf provides for general nonlinear equations.
The ODE \(y' = F(x, y)\) is separable if and only if \(F\) can be factored as a product \(F(x,y) = f(x)\,g(y)\). A practical test: if \(\displaystyle \frac{\partial^2}{\partial x\,\partial y}\!\ln|F(x,y)| = 0\) wherever \(F \neq 0\), then \(F\) is separable.
Many students attempt to separate variables on equations like \(y' = x + y\). This equation is linear but not separable because the right side cannot be written as \(f(x)\,g(y)\). Forcing a separation leads to algebraic errors. Always classify before solving.
When solving \(\displaystyle \frac{dy}{y} = k\,dt\), the result is \(\ln|y| = kt + C\), not \(\ln y = kt + C\). Omitting the absolute value silently discards negative-valued solutions. In many applications the sign is determined by context (e.g., population is positive), but the mathematical step must be justified.
When computing \(\mu(t) = e^{\int p(t)\,dt}\), the arbitrary constant in the exponent is irrelevant because it produces a multiplicative constant that cancels. Students who carry it through create unnecessary complexity. Conversely, students sometimes omit the constant of integration in the final solution step, where it is essential.
The classic example \(y' = 3y^{2/3}\), \(y(0)=0\) has the trivial solution \(y=0\) and the family \(y = (t-c)^3\) for any \(c \ge 0\). Here \(\partial f/\partial y = 2y^{-1/3}\) is unbounded near \(y=0\), so the Lipschitz condition fails. Blindly applying a numerical solver returns one particular solution with no warning that others exist.
The test \(\partial M/\partial y = \partial N/\partial x\) determines whether the equation is already exact. If it fails, you need an integrating factor to make it exact. The integrating factor for exactness is in general different from the integrating factor for linear equations. Students often mix up the two procedures.
When separating \(\displaystyle \frac{dy}{g(y)} = f(x)\,dx\), any value \(y^*\) where \(g(y^*) = 0\) is a singular solution (an equilibrium). Dividing by \(g(y)\) silently excludes these constant solutions. Always check \(g(y)=0\) separately before declaring a general solution.
Problem. A 500-liter tank initially contains 100 liters of pure water. Brine with a salt concentration of 20 g/L flows in at 6 L/min. The well-stirred solution flows out at 4 L/min. Find the amount of salt \(Q(t)\) in the tank at time \(t\), and determine when the tank is full and the salt concentration at that moment.
Volume at time \(t\): \(V(t) = 100 + (6-4)t = 100 + 2t\) liters.
Rate in: \(6 \times 20 = 120\) g/min.
Rate out: \(\displaystyle 4 \times \frac{Q(t)}{V(t)} = \frac{4Q}{100+2t}\) g/min.
Balance:
Integrate both sides:
\(Q(0) = 0\): \(2500 \cdot 0 = 40 \cdot 125000 + C\), so \(C = -5{,}000{,}000\).
The tank is full when \(V(t) = 500\), i.e., \(t = 200\) min. At \(t=200\):
Concentration at \(t=200\): \(9920/500 = 19.84\) g/L, which is close to but slightly below the incoming concentration of 20 g/L, as expected.
As \(t \to \infty\) (if the tank could grow indefinitely), \(Q/V \to 40(50+t)/(50+t) = 40\)... but here the physical constraint \(V=500\) changes the regime at \(t=200\). After that, a new ODE with constant volume would apply.
Problem. A bacteria culture in a nutrient broth follows the logistic equation \(\displaystyle \frac{dP}{dt} = 0.5\,P\!\left(1 - \frac{P}{1000}\right)\), with initial population \(P(0) = 50\). Find \(P(t)\) and determine when the population reaches 900.
Use partial fractions on the left:
Solve for \(P\):
\(P(0)=50\): \(\displaystyle 50 = \frac{1000\,A}{1+A}\), so \(A = 50/950 = 1/19\).
The solution is a sigmoid curve. The population grows nearly exponentially when \(P \ll 1000\), inflects at \(P = 500\) (half the carrying capacity), and asymptotically approaches \(K=1000\). The inflection point occurs at \(t = 2\ln 19 \approx 5.89\), where the growth rate is maximal. After \(t \approx 10.3\), the population has consumed most of the available capacity and growth has slowed dramatically.
Problem. Solve \((2xy + 3)\,dx + (x^2 - 1)\,dy = 0\).
\(M = 2xy + 3\), \(N = x^2 - 1\).
Equal, so the equation is exact.
This can be solved explicitly: \(\displaystyle y = \frac{C - 3x}{x^2 - 1}\) for \(x \neq \pm 1\).
Below is a complete Python script (Pyodide-compatible) that solves the mixing problem from Example 1 using Euler's method and compares the result to the exact analytical solution. Copy this code into a Pyodide-enabled environment to run it.
import numpy as np
# -----------------------------------------------------------
# Mixing Problem: 500 L tank, starts at 100 L pure water
# Inflow: 6 L/min at 20 g/L
# Outflow: 4 L/min
# V(t) = 100 + 2t, tank full at t = 200 min
# -----------------------------------------------------------
def dQ_dt(t, Q):
"""Right-hand side of the mixing ODE: Q' = 120 - 2Q/(50+t)."""
return 120.0 - 2.0 * Q / (50.0 + t)
def Q_exact(t):
"""Exact solution derived in Example 1."""
return 40.0 * (50.0 + t) - 5_000_000.0 / (50.0 + t)**2
# ---- Euler's method ----
dt = 0.5 # step size (minutes)
t_end = 200.0 # tank full at t=200
N_steps = int(t_end / dt)
t_euler = np.zeros(N_steps + 1)
Q_euler = np.zeros(N_steps + 1)
t_euler[0] = 0.0
Q_euler[0] = 0.0 # initial condition: no salt
for n in range(N_steps):
Q_euler[n+1] = Q_euler[n] + dt * dQ_dt(t_euler[n], Q_euler[n])
t_euler[n+1] = t_euler[n] + dt
# ---- Evaluate exact solution on the same grid ----
Q_ex = Q_exact(t_euler)
# ---- Compute error metrics ----
abs_error = np.abs(Q_euler - Q_ex)
max_error = np.max(abs_error)
final_euler = Q_euler[-1]
final_exact = Q_ex[-1]
print("Mixing Problem — Euler vs. Exact")
print(f" Step size: {dt} min")
print(f" Number of steps: {N_steps}")
print(f" Q_euler(200): {final_euler:.4f} g")
print(f" Q_exact(200): {final_exact:.4f} g")
print(f" Max absolute error: {max_error:.4f} g")
print(f" Final concentration: {final_euler/500:.4f} g/L")
# ---- Optional: plot (works in Pyodide with matplotlib) ----
try:
import matplotlib
matplotlib.use('Agg') # non-interactive backend
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# Left: Q(t)
ax1.plot(t_euler, Q_ex, 'b-', linewidth=2, label='Exact')
ax1.plot(t_euler[::20], Q_euler[::20], 'ro', markersize=4, label=f'Euler (dt={dt})')
ax1.set_xlabel('Time (min)')
ax1.set_ylabel('Salt Q(t) (grams)')
ax1.set_title('Salt in Tank over Time')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Right: error
ax2.plot(t_euler, abs_error, 'r-', linewidth=1.5)
ax2.set_xlabel('Time (min)')
ax2.set_ylabel('|Euler - Exact| (grams)')
ax2.set_title('Absolute Error')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('mixing_problem_euler.png', dpi=120)
print(" Plot saved to mixing_problem_euler.png")
except ImportError:
print(" (matplotlib not available — skipping plot)")
Expected output (approximate):
Mixing Problem — Euler vs. Exact
Step size: 0.5 min
Number of steps: 400
Q_euler(200): 9919.4762 g
Q_exact(200): 9920.0000 g
Max absolute error: 1.0932 g
Final concentration: 19.8390 g/L
Experiment: Change dt to 2.0, 1.0, 0.1, and observe how the
maximum error scales roughly linearly with step size (first-order convergence). Then try
implementing a simple second-order method (Heun's method) and compare.
Viewing an ODE solver through the lens of reinforcement learning reveals a natural decision-making structure. Here we frame Euler's method (and more sophisticated solvers) as an agent that sequentially interacts with the differential equation.
At each step \(n\), the agent observes the state \(s_n = (t_n,\, y_n,\, \hat{e}_n)\) where \(t_n\) is the current time, \(y_n\) is the current numerical approximation, and \(\hat{e}_n\) is an estimate of the local truncation error (zero for fixed-step Euler, computed via step-doubling or embedded methods in adaptive solvers).
The agent selects an action \(a_n = (h_n, m_n)\): a step size \(h_n > 0\) and possibly
a method order \(m_n\) (e.g., Euler = order 1, RK4 = order 4). In fixed-step Euler the
action space is trivial (one constant), but adaptive solvers like scipy.integrate.solve_ivp
with method='RK45' make genuine step-size decisions at every stage.
The reward balances accuracy against computational cost:
The first term penalizes large local error; the second penalizes small steps (which cost more function evaluations). The constants \(\alpha, \beta > 0\) define the trade-off. A solver that takes too-large steps earns large negative rewards from error; one that takes too-small steps earns large negative rewards from cost.
Classical adaptive solvers implement a hand-crafted policy: if the estimated error exceeds a tolerance \(\varepsilon\), reject the step and halve \(h\); if the error is well below \(\varepsilon\), accept and grow \(h\) by a safety-factored ratio. This is analogous to a bang-bang controller. Modern neural-ODE solvers and learned step-size controllers replace this heuristic with a trained policy network, demonstrating that the agent-environment framing is not merely an analogy but a practical design paradigm.
In Module 1.4 (Numerical ODE Solvers) we will formalize the agent framework further, introducing the cost-of-computation perspective: given a fixed compute budget of \(N\) function evaluations, what is the optimal allocation of step sizes to minimize global error? This is a finite-horizon sequential optimization problem.
Exercise 1.1.1 (Separable). Solve the IVP \(\displaystyle \frac{dy}{dx} = \frac{x\,e^{-x^2}}{y}\), \(y(0) = 1\). Express the solution explicitly.
Separate: \(y\,dy = x\,e^{-x^2}\,dx\). Integrate:
Apply \(y(0)=1\): \(1/2 = -1/2 + C\), so \(C=1\).
Exercise 1.1.2 (Linear). Solve \(\displaystyle y' + \frac{3}{t}\,y = \frac{\sin t}{t^3}\) for \(t > 0\).
Integrating factor: \(\mu = e^{3\ln t} = t^3\).
Exercise 1.1.3 (Exact). Determine whether \((ye^{xy} + 2x)\,dx + (xe^{xy} - 1)\,dy = 0\) is exact, and if so, solve it.
\(M = ye^{xy} + 2x\), \(N = xe^{xy} - 1\).
\(M_y = e^{xy} + xye^{xy}\), \(N_x = e^{xy} + xye^{xy}\). Equal, so exact.
\(F = \int M\,dx = e^{xy} + x^2 + h(y)\). \(F_y = xe^{xy} + h'(y) = xe^{xy} - 1\), so \(h'(y) = -1\), \(h = -y\).
Exercise 1.1.4 (Existence-Uniqueness). Consider the IVP \(y' = \sqrt{|y|}\), \(y(0) = 0\). Does a solution exist? Is it unique? Find all solutions.
\(f(t,y) = \sqrt{|y|}\) is continuous everywhere, so Peano's theorem guarantees existence. However, \(\partial f/\partial y = 1/(2\sqrt{|y|})\) is unbounded near \(y=0\), so the Lipschitz condition fails. Indeed, uniqueness fails: \(y \equiv 0\) is a solution, and for any \(c \ge 0\), the function
is also a solution. There are infinitely many solutions.
Exercise 1.1.5. Implement Euler's method to solve the logistic equation \(P' = 0.5\,P(1 - P/1000)\), \(P(0)=50\), on \([0, 30]\) with step sizes \(h = 1, 0.5, 0.1\). Compute the error at \(t=10\) against the exact solution from Example 2 and verify the error is \(O(h)\).
Exercise 1.1.6. Write a Python function that takes a function \(f(t,y)\), an initial condition \((t_0, y_0)\), and an end time \(t_f\), and solves the IVP with a user-specified tolerance using adaptive step-size control: at each step compute both a full-step Euler result and two half-step results; if they disagree by more than the tolerance, halve \(h\); if they agree to within tolerance/10, double \(h\). Test on Newton's cooling: \(T' = -0.1(T - 20)\), \(T(0) = 90\).
import numpy as np
def adaptive_euler(f, t0, y0, tf, tol=1e-4, h_init=1.0):
"""Adaptive Euler method using step-doubling error estimation."""
t, y, h = t0, y0, h_init
ts, ys = [t], [y]
while t < tf - 1e-12:
h = min(h, tf - t)
# Full step
y_full = y + h * f(t, y)
# Two half steps
y_half = y + (h/2) * f(t, y)
y_half = y_half + (h/2) * f(t + h/2, y_half)
# Error estimate
err = abs(y_half - y_full)
if err > tol:
h /= 2.0
continue
# Accept step (use the more accurate half-step result)
t += h
y = y_half
ts.append(t)
ys.append(y)
if err < tol / 10:
h *= 2.0
return np.array(ts), np.array(ys)
# Test: Newton's cooling
f_cool = lambda t, T: -0.1 * (T - 20.0)
ts, Ts = adaptive_euler(f_cool, 0, 90.0, 50.0, tol=1e-3)
T_exact = lambda t: 20 + 70 * np.exp(-0.1 * t)
print(f"Steps taken: {len(ts)-1}")
print(f"T(50) adaptive: {Ts[-1]:.6f}")
print(f"T(50) exact: {T_exact(50):.6f}")
Exercise 1.1.7. For the radioactive decay equation \(N' = -\lambda N\) with \(\lambda = 0.0866\,\text{yr}^{-1}\) (corresponding to a half-life of 8 years) and \(N(0)=1000\), solve analytically and numerically on \([0, 40]\). Plot both on the same axes and add a horizontal line at \(N = 1000/2^5 = 31.25\) to visualize 5 half-lives.
Exercise 1.1.8.
Use scipy.integrate.solve_ivp to solve the mixing problem from Example 1 and
compare with the exact solution. Use both method='RK45' and method='Radau'.
Which uses fewer function evaluations? Why?
from scipy.integrate import solve_ivp
import numpy as np
def mixing_ode(t, Q):
return [120.0 - 2.0 * Q[0] / (50.0 + t)]
Q_exact = lambda t: 40.0*(50+t) - 5e6/(50+t)**2
sol_rk45 = solve_ivp(mixing_ode, [0, 200], [0.0], method='RK45',
dense_output=True, rtol=1e-8, atol=1e-10)
sol_radau = solve_ivp(mixing_ode, [0, 200], [0.0], method='Radau',
dense_output=True, rtol=1e-8, atol=1e-10)
t_eval = np.linspace(0, 200, 500)
print(f"RK45 steps: {sol_rk45.t.size}, nfev: {sol_rk45.nfev}")
print(f"Radau steps: {sol_radau.t.size}, nfev: {sol_radau.nfev}")
print(f"RK45 error at t=200: {abs(sol_rk45.sol(200)[0] - Q_exact(200)):.2e}")
print(f"Radau error at t=200: {abs(sol_radau.sol(200)[0] - Q_exact(200)):.2e}")
# RK45 typically uses fewer evaluations for this non-stiff problem.
Exercise 1.1.9 (Agent Design). Design a "step-size agent" for solving \(y' = -50y + 50\sin t\) on \([0, 2]\) with \(y(0)=0\). This equation is stiff: the transient \(e^{-50t}\) decays rapidly while the forced response \(\sin t\) varies slowly. Your agent should:
Report: Does your adaptive agent use fewer evaluations while achieving comparable accuracy? At what times does it choose large steps vs. small steps?
The stiff transient near \(t=0\) forces small steps (\(h \approx 0.001\) to maintain stability for explicit Euler). After \(t \approx 0.1\), the transient has decayed by \(e^{-5} \approx 0.007\) and the agent can safely increase \(h\) to \(\sim 0.05\) or more. The adaptive agent typically uses 200–300 evaluations compared to 2000 for the fixed agent, with comparable accuracy. The agent "learns" (via the error estimate) that the dynamics have become benign and increases its step size accordingly. Note: explicit Euler is unconditionally unstable for large \(h\) on stiff problems; if the agent grows \(h\) too aggressively, it must detect and reject bad steps.
Exercise 1.1.10 (Model Selection Agent). You are given noisy data from an unknown first-order process: time-concentration pairs measured every 2 minutes for 30 minutes. The agent must choose between three models:
Write a Python script that fits each model to the data (using least-squares), computes the AIC (Akaike Information Criterion) for each, and selects the best model. Test with synthetic data generated from exponential decay with added Gaussian noise (\(\sigma=0.5\)).
Which of the following is separable?
(a) \(y' = x + y\) (b) \(y' = xy^2\) (c) \(y' = x^2 + y^2\) (d) \(y' = e^{x+y}\)
The integrating factor for \(y' + 2ty = t\) is:
(a) \(e^{t^2}\) (b) \(e^{2t}\) (c) \(e^{t}\) (d) \(t^2\)
For the equation \(M\,dx + N\,dy = 0\) to be exact, we require:
(a) \(M = N\) (b) \(M_x = N_y\) (c) \(M_y = N_x\) (d) \(M_x = N_x\)
The Picard–Lindelöf theorem guarantees both existence and uniqueness when \(f(t,y)\) satisfies:
(a) Continuity in \(t\) only (b) Continuity in \(y\) only (c) A Lipschitz condition in \(y\) and continuity in both variables (d) Differentiability everywhere
How many solutions does the IVP \(y' = 3y^{2/3}\), \(y(0) = 0\) have?
(a) None (b) Exactly one (c) Exactly two (d) Infinitely many
In the logistic equation \(P' = rP(1 - P/K)\), the inflection point of \(P(t)\) occurs at:
(a) \(P = 0\) (b) \(P = K/4\) (c) \(P = K/2\) (d) \(P = K\)
Newton's law of cooling states \(T' = -k(T - T_{\text{env}})\). If \(T(0)=90\), \(T_{\text{env}}=20\), and \(T(10)=60\), then \(k\) equals:
(a) \(\tfrac{1}{10}\ln(7/4)\) (b) \(\tfrac{1}{10}\ln(7/3)\) (c) \(\tfrac{1}{10}\ln(9/4)\) (d) \(\tfrac{1}{10}\ln 2\)
Euler's method with step size \(h\) has a global truncation error of order:
(a) \(O(h^2)\) (b) \(O(h)\) (c) \(O(h^{1/2})\) (d) \(O(1)\)
In the mixing problem of Example 1, as \(t \to \infty\) (ignoring the tank-full constraint), the concentration \(Q/V\) approaches:
(a) 0 g/L (b) 10 g/L (c) 20 g/L (d) 40 g/L
Which statement about the solution to the linear ODE \(y' + p(t)y = g(t)\) is true?
(a) The solution exists only locally near \(t_0\) (b) The solution exists on any interval where \(p\) and \(g\) are continuous (c) Multiple solutions can pass through each initial point (d) The solution is always bounded
A drug is administered intravenously at a constant rate of \(D = 50\) mg/hr for the first 2 hours, then the infusion stops. The drug is eliminated from the bloodstream at a rate proportional to its current concentration, with elimination constant \(k_e = 0.3\,\text{hr}^{-1}\). The volume of distribution is \(V_d = 10\) L.
Let \(C(t)\) be the drug concentration (mg/L) in the blood. The governing ODE is:
with \(C(0) = 0\).
| Criterion | Points | Description |
|---|---|---|
| Correctness | 35 | Analytical solution is correct on both intervals with proper matching at \(t=2\). Numerical solutions converge to the analytical result. Clinical quantities are computed accurately. |
| Efficiency | 20 | Adaptive solver uses significantly fewer steps than fixed-step Euler for comparable accuracy. Convergence study uses a systematic sequence of step sizes (e.g., halving). |
| Robustness | 25 | Code handles the discontinuity at \(t=2\) correctly (e.g., event detection or explicit split). Edge cases (negative concentration, very large \(t\)) are addressed. Error estimation is included. |
| Reproducibility | 20 | Code is well-documented with clear variable names. Plots are labeled and captioned. A requirements specification or environment description is provided. Results can be reproduced by running the script. |
scipy.integrate.solve_ivp — https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html