explainer
Ordinary differential equations: turn a rate law into a time course
Solve a cooling initial-value problem, compare its exact solution with Euler steps, and separate model behavior from numerical accuracy and stability.
What you will learn
- Interpret a differential equation as a rule for a rate of change.
- Use an initial condition to select a solution of a cooling model.
- Calculate forward Euler steps and compare them with the exact solution.
- Explain why step size affects accuracy and numerical stability differently.
Before you start
- Derivatives as rates and the chain rule
- Exponential functions and negative exponents
A cooling component does not usually lose the same number of degrees every minute. Its cooling rate changes as it approaches room temperature. An ordinary differential equation, or ODE, can express that changing rate before we know the full temperature curve.
This lesson solves one synthetic cooling model, then approximates it with forward Euler steps. The comparison shows how a numerical method can produce oscillations that the original model never predicts.
Read an ODE as a rule for a rate
An ODE relates an unknown function to derivatives with respect to one independent variable. That variable is often time. A first-order ODE uses a first derivative as its highest derivative:
dy/dt = f(t, y)
The function f specifies the current rate of change. It does not directly give y itself. If y is a position in meters, dy/dt is a velocity in meters per second; if y is temperature, the derivative has temperature-per-time units.
For example, dy/dt = −y says the rate depends on the current value. At y = 4 the rate is −4, while at y = 1 it is −1. The solution must change its slope as it moves.
The word ordinary distinguishes derivatives with respect to one variable from partial derivatives involving several independent variables. An ODE can still track several dependent quantities together, such as position and velocity. OpenStax defines differential equations, their order, and what counts as a solution.
Add the starting value
A rate law often describes a family of curves. Both y(t) = e⁻ᵗ and y(t) = 4e⁻ᵗ satisfy dy/dt = −y. They start at different values.
An initial condition states a value at a specified time. Pairing the equation with y(0) = 4 selects y(t) = 4e⁻ᵗ from this family. Together, the rate law and starting value form an initial-value problem.
To verify a proposed answer, check both requirements: differentiate the function and substitute it into the equation, then evaluate its initial value. Our linear cooling problem has a unique solution for each initial temperature. More general ODEs need appropriate existence and uniqueness conditions; supplying an initial value alone does not settle every problem.
Transfer functions express a linear system's input-output relation in the Laplace domain. That relation assumes zero initial conditions; the lesson separates the response to an applied torque from motion left over at the start.
Model one cooling robot component
Imagine a robot component cooling after its heat source switches off. Use a model with one uniform component temperature, a fixed room temperature, and a constant cooling coefficient. These assumptions simplify the physics; the numbers here do not come from a measured robot.
Let T(t) be the component temperature in °C, with time t in minutes. Room temperature is T_a = 20°C, and the positive coefficient k has units min⁻¹:
dT/dt = −k(T − T_a)
T(0) = T₀
This form of Newton's law of cooling makes the rate proportional to the temperature difference. If T exceeds room temperature, the derivative is negative. At T = T_a it is zero, giving an equilibrium. OpenStax develops this cooling model and its exponential solution.
For T₀ = 80°C and k = 1 min⁻¹, the initial rate is −1(80 − 20) = −60°C/min. That is an instantaneous slope. It does not mean the component must lose exactly 60 degrees over the next minute.
Temperature differences have the same numerical size in Celsius degrees and kelvins. The exponent we will use, kt, is dimensionless because minutes cancel inverse minutes.
Solve and verify the cooling equation
Subtract the constant ambient temperature: let z(t) = T(t) − T_a. Then z′ = −kz. Its solution is an exponentially shrinking temperature difference:
z(t) = z(0)e⁻ᵏᵗ
T(t) = T_a + (T₀ − T_a)e⁻ᵏᵗ
At t = 0, the exponential is 1, so the formula returns T₀. Differentiating gives T′ = −k(T₀ − T_a)e⁻ᵏᵗ = −k(T − T_a), which verifies the rate law. The chain-rule lesson explains the derivative of that exponential composition.
Our example becomes T(t) = 20 + 60e⁻ᵗ. After half a minute, the exact model gives about 56.392°C. After one minute, it gives 42.073°C, still well above room temperature.
The time constant τ = 1/k measures how quickly the initial difference shrinks. After one time constant, about 36.8% of that difference remains. For k = 1 min⁻¹, τ is one minute; that percentage applies to T − 20, not to the Celsius temperature itself.
Starting above ambient, this solution decreases toward 20°C without crossing it. Starting exactly at ambient gives a constant solution. A starting temperature below ambient would warm toward it under the same equation.
Take a forward Euler step
Many ODEs lack a convenient exact formula. Forward Euler uses the rate at the current numerical state to predict the next value:
t_(n+1) = t_n + h
T_(n+1) = T_n + h[−k(T_n − T_a)]
The step size h must use the same time unit as the rate. This is a repeated local linear prediction, connected to Taylor expansion and linearization. OpenStax describes the forward Euler construction.
Take h = 0.5 min. Starting at 80°C, the first prediction is 80 + 0.5(−60) = 50°C. The exact value at the same time is 56.392°C, so the first step cools too far.
For the next step, recompute the derivative at 50°C, the Euler value. It is −30°C/min, giving 50 + 0.5(−30) = 35°C at one minute. Euler carries its previous approximation forward; it does not restart each step from the exact curve.
Compare the exact curve with numerical samples
The Euler cooling experiment keeps room temperature fixed at 20°C and compares every run at six minutes. Changing T₀ or k changes the initial-value problem. Changing h changes the numerical method's grid while preserving that problem.
The default h = 0.5 min produces 12 steps. At six minutes, Euler gives about 20.015°C, while the exact solution gives 20.149°C. That small final gap hides an earlier sample error of about 7.073°C at one minute.
The maximum-error readout compares values at all Euler grid times. It does not claim to find the maximum error between samples. The table shows the initial point, first two steps, and final point; the plot includes every grid sample.
Try the numerical-oscillation presets while watching the unchanged exact curve. Select Already at ambient to see a constant solution. The plot expands its temperature range to keep inaccurate Euler samples visible.
Find the step-size stability boundary
Write the numerical deviation from ambient as z_n = T_n − T_a. The Euler recurrence becomes:
z_(n+1) = (1 − kh)z_n
z_n = (1 − kh)ⁿz₀
The multiplier r = 1 − kh determines what repeated steps do to a nonzero deviation or perturbation. Its magnitude controls growth; its sign controls alternation across ambient.
| Step condition | Euler deviation from ambient |
|---|---|
| 0 < kh < 1 | Shrinks without changing sign |
| kh = 1 | Reaches zero in one numerical step |
| 1 < kh < 2 | Shrinks while alternating sign |
| kh = 2 | Alternates without shrinking |
| kh > 2 | Alternates with growing magnitude |
Thus 0 < h < 2/k gives decaying deviations for this scalar equation. The boundary h = 2/k has multiplier −1, so it fails to reproduce decay. MIT derives this forward Euler stability restriction for exponential decay.
With k = 1 and h = 2, the temperatures are 80, −40, 80, −40°C. With h = 3, they are 80, −100, 260°C. These oscillations come from the numerical approximation; the exact model cools smoothly throughout.
At h = 1.5, deviations shrink with multiplier −0.5, yet Euler still crosses ambient. Decay of numerical disturbances does not guarantee an accurate or physically faithful trajectory. Even h = 1 sends Euler directly to ambient while the exact solution remains above it at every finite time.
If z₀ is exactly zero, multiplying it by any finite r keeps it zero. The behavior label still describes what would happen to a nonzero perturbation. A constant equilibrium run cannot establish that a step size handles nearby states well.
The ODE stability lesson examines the continuous equation itself: which equilibria attract nearby states, which repel them, and which starting states approach each equilibrium. Those questions concern the modeled dynamics before you choose a numerical method.
Separate accuracy, stability, and model quality
Euler omits the curvature terms in a Taylor expansion. For this smooth model, a single step starting from the exact state has an error of order h². Over a fixed time interval, repeated Euler steps generally accumulate an error of order h as h becomes small, ignoring roundoff. MIT distinguishes this local step error from the accumulated error.
That means halving a sufficiently small h roughly halves the fixed-time error, while doubling the step count. It does not imply that every finite step reduction halves every error measure. Compare approximations at the same times and check whether refinement changes the result enough to matter.
The numerical integration lesson compares Euler, symplectic Euler, and Runge–Kutta steps on an oscillator. It follows both trajectory error and energy drift to show what changing the integration method can improve.
Three separate questions remain:
- Accuracy: how close are the computed values to the solution of the chosen ODE?
- Numerical stability: do the numerical updates damp or amplify disturbances for the problem and step size?
- Model quality: does that ODE describe the physical system well enough for the task?
A smaller h cannot fix an incorrect cooling coefficient or missing heat generation. Multiple coupled temperatures may need a vector state with several rate equations. The Jacobian of that rate function describes how changing the state changes its rates; the scalar condition h < 2/k does not cover every such system.
Forward dynamics uses a vector state containing a robot arm's joint positions and velocities. Solving its mass-matrix equation gives acceleration, which supplies the velocity part of the state derivative. Integrating that derivative predicts motion under the chosen torques and model assumptions.
The same distinction appears in optimization. A continuous gradient-flow model uses θ′ = −∇L(θ); a forward Euler step gives the update θ_next = θ − h∇L(θ). This connects to gradient descent, whose step parameter belongs to the optimization algorithm.
Reproduce the comparison in Python
This Python 3 example uses only the standard library. Each listed step divides the six-minute horizon exactly. The method calculates every Euler value from the previous numerical value, then evaluates the exact solution at the same final time.
from math import exp
ambient, initial, k = 20.0, 80.0, 1.0
duration = 6.0
def exact(time):
return ambient + (initial - ambient) * exp(-k * time)
print(f"Initial rate: {-k * (initial - ambient):.3f} C/min")
print(f"Exact at 0.5 min: {exact(0.5):.6f} C")
print(f"Exact at 6 min: {exact(duration):.6f} C")
for h in (0.5, 1.0, 1.5, 2.0, 3.0):
steps = round(duration / h)
assert abs(steps * h - duration) < 1e-12
temperature = initial
for _ in range(steps):
derivative = -k * (temperature - ambient)
temperature += h * derivative
error = abs(temperature - exact(duration))
print(f"h={h:.1f}: steps={steps:2d}, Euler={temperature:.6f}, "
f"error={error:.6f}")
Expected output:
Initial rate: -60.000 C/min
Exact at 0.5 min: 56.391840 C
Exact at 6 min: 20.148725 C
h=0.5: steps=12, Euler=20.014648, error=0.134077
h=1.0: steps= 6, Euler=20.000000, error=0.148725
h=1.5: steps= 4, Euler=23.750000, error=3.601275
h=2.0: steps= 3, Euler=-40.000000, error=60.148725
h=3.0: steps= 2, Euler=260.000000, error=239.851275
This demonstration uses an exact solution as a reference. For an unfamiliar ODE, a solver usually needs error estimates, refinement checks, or another trusted reference. A smooth-looking numerical curve alone does not verify its accuracy.
Try it yourself
Exercise 1. A component starts at 60°C in a 20°C room, with k = 0.5 min⁻¹. Find its initial rate, exact temperature after one minute, and first Euler value with h = 1 min. Explain why the two temperatures differ.
Show the initial-value and first-step solution
The initial rate is −0.5(60 − 20) = −20°C/min. The exact solution is T(t) = 20 + 40 exp(−0.5t), so at one minute it gives 20 + 40 exp(−0.5) ≈ 44.261226°C.
Euler gives 60 + 1(−20) = 40°C. It extends the initial slope over the whole step. The actual cooling rate becomes less negative as the temperature falls, so the exact model stays warmer than this prediction.
Exercise 2. Return to T₀ = 80°C, T_a = 20°C, and k = 1 min⁻¹. For h = 2 min, compute the multiplier and first two Euler values. Does starting exactly at 20°C instead make this step size suitable for nearby initial temperatures?
Show the stability-boundary solution
The multiplier is 1 − kh = −1. Deviations from ambient go 60, −60, 60, giving temperatures 80, −40, 80°C at times 0, 2, and 4 minutes. The deviation does not decay.
Starting exactly at 20°C gives zero deviation forever, but any nonzero initial deviation alternates without shrinking. The exact model shrinks those deviations by e⁻² each two-minute interval. The constant equilibrium run therefore does not validate Euler's step size for neighboring states.
Sources and further study
- OpenStax, Calculus Volume 2: Basics of Differential Equations: rate equations, solution verification, order, and initial-value problems.
- OpenStax, Exponential Growth and Decay: Newton's cooling model and the exponential decay of the temperature difference.
- OpenStax, Direction Fields and Numerical Methods: forward Euler's repeated linear prediction.
- MIT, Forward and Backward Euler Methods: local and accumulated error, and the forward Euler step restriction for exponential decay.