Numerical solutions of ordinary differential equations

State the model and the time grid

A supplied position model starts at 1 m and has rate x′(t) = (2 m/s²)t on [0 s, 2 s]. Can we advance it in short time steps? Here t is elapsed time, and x is a signed position in metres along a fixed axis, with the same tracked point and origin. The exact check solution is x(t) = 1 m + (1 m/s²)t². At 2 s it gives 5 m; we will see how two numerical step sizes compare with that value.

The prime means derivative with respect to time. More generally, a scalar initial-value problem supplies a rate y′(t) = F(t, y(t)) and a starting value y(t₀) = y₀. You need the meaning of a derivative and an initial condition; the published lessons Derivatives, Taylor series, and local approximation and Ordinary differential equations supply that preparation. The position model here is a stipulated test with a known solution, rather than a measured trajectory or verified force law. For other problems, assume a solution exists on the intended interval; a numerical method does not establish that existence.

Choose N positive integer steps, width h > 0, and times tₙ = t₀ + nh for n = 0, …, N. The subscript n counts steps. Set h = (t_N − t₀)/N so the final grid point is the chosen final time. A computed value Yₙ approximates y(tₙ); it is neither a new instrument reading nor an entire continuous solution between grid points.

Advance with explicit Euler

Explicit Euler uses Yₙ₊₁ = Yₙ + hF(tₙ, Yₙ), starting from Y₀ = y₀. It takes the rate at the beginning of the step as constant for that update. If y is a position, F has units m/s and hF has units m, so the sum has the right quantity kind and units.

The rate must be evaluated at the current numerical state, not at an unknown exact future state. For our position model F does not depend on position, so Xₙ₊₁ = Xₙ + h(2 m/s²)tₙ. With h = 1 s, X₀ = 1 m gives X₁ = 1 m and X₂ = 3 m. The exact positions at those times are 2 m and 5 m.

Following the update exactly does not make the numbers exact solutions of the differential equation. The first Euler update leaves position unchanged because its starting rate is zero, although the supplied rate becomes positive during that step. Between-point display interpolation is a separate choice; this lesson compares values at matching grid times.

Distinguish one-step and global errors

To isolate a step's approximation error, imagine starting that step from the exact y(tₙ), using exact arithmetic. Define dₙ = y(tₙ₊₁) − [y(tₙ) + hF(tₙ, y(tₙ))]. This is a one-step state discrepancy, with the same units as y. Terminology varies: some sources call dₙ/h the truncation or consistency error instead, so the formula and units matter.

If the exact solution has a continuous second derivative throughout the step, Taylor's theorem gives dₙ = (h²/2)y″(ξ) for some interior time ξ. An interval-wide bound |y″| ≤ M therefore gives |dₙ| ≤ Mh²/2. For a position, M has units m/s². This bound starts from the exact state; it does not include an error inherited from earlier numerical steps.

The global error at a grid point is eₙ = y(tₙ) − Yₙ, and its absolute size is Eₙ = |eₙ|. Actual Euler steps usually start from approximate Yₙ. A changed state can change the next rate, so global error requires an argument about both new step discrepancies and inherited errors. Rounding, inaccurate initial data and a wrong physical model are further contributions.

Propagate error with explicit conditions

Suppose |F(t,u) − F(t,v)| ≤ L|u − v| for the states being compared at each grid time, with L ≥ 0. This is a Lipschitz bound: it limits how much a rate changes when the state changes. L has units s⁻¹. Require the exact solution and computed states to stay in the domain where this bound holds, and retain a continuous second derivative with |y″| ≤ M on the whole computed interval.

Subtracting the Euler update from the exact step identity and taking absolute values gives Eₙ₊₁ ≤ (1 + hL)Eₙ + Mh²/2. The factor 1 + hL is dimensionless; every term has the units of y. This is an upper bound, not an assertion that errors always have the same sign or always grow. It assumes exact arithmetic and the stated initial data unless additional error terms are included.

For our position rate, L = 0 s⁻¹ and M = 2 m/s². With exact X₀, repeated use of the bound gives E_N ≤ N Mh²/2. At the fixed final time 2 s this is (1 m/s²)(2 s)h: 2 m for h = 1 s and 1 m for h = 0.5 s. This particular quadratic attains those bounds. General state-dependent rates require their own propagation bound; a smaller local discrepancy alone is not the full global-error analysis.

Evaluate two rates in an explicit step

A two-stage alternative first forms k₁ = F(tₙ, Yₙ) and a predictor P = Yₙ + hk₁. It then evaluates k₂ = F(tₙ + h, P) and accepts Yₙ₊₁ = Yₙ + (h/2)(k₁ + k₂). Both k values are rates. This predictor-and-average formula is an explicit version of improved Euler; P is a trial state, not the accepted next value.

For our position model at t₀ = 0 s with h = 1 s, k₁ = 0 m/s, P = 1 m and k₂ = 2 m/s. The accepted X₁ is 2 m. Repeating from 1 s averages 2 m/s and 4 m/s, giving X₂ = 5 m. The rates are affine functions of time and independent of the state, so their endpoint average gives the exact position increment for every step of this supplied model.

That exactness is specific to this rate rule. For a separate decay rate F(t,y) = −λy, the same two-stage formula has multiplier 1 − λh + (λh)²/2. At λh = 2 it is 1, so the update does not decay despite two rate evaluations. Additional stages, a method's name or exact arithmetic do not by themselves establish the required accuracy or stability.

Test decay stability and signs

Use a separate stipulated model: θ(t) is temperature excess in kelvin relative to a fixed 280 K reference. Supply θ′ = −λθ, λ = 2 s⁻¹ and θ(0 s) = 8 K on [0 s, 3 s]. The exact solution θ(t) = (8 K)e^(−λt) is positive and decreases there; the exponent is dimensionless. This is a teaching rate law, not verification of an actual cooling object.

Explicit Euler gives Θₙ₊₁ = rΘₙ with r = 1 − λh. A change in the initial value is multiplied by r at every step. For this linear test recurrence, |r| < 1 damps its magnitude, whereas |r| > 1 amplifies it. The condition for strict damping is 0 < λh < 2. At h = 1 s, r = −1: 8 K, −8 K, 8 K shows alternating signs without damping. At h = 1.5 s, r = −2 gives 8 K, −16 K, 32 K by 3 s.

Damping is different from sign preservation and accurate values. With h = 0.75 s, r = −0.5 damps the magnitude but the first value is −4 K, contrary to the exact positive excess. Nonnegative Euler values for every step from a positive datum require 0 ≤ λh ≤ 1; at λh = 1 they become zero immediately rather than match the positive exact solution. These tests concern this scalar recurrence, not a universal criterion for nonlinear or coupled systems.

Recognize an implicit update

Implicit Euler instead specifies Yₙ₊₁ = Yₙ + hF(tₙ₊₁, Yₙ₊₁). The unknown next state appears inside its own rate evaluation, so a step generally requires solving an equation. A nonlinear step may need further hypotheses and a solver; the label implicit does not guarantee that this equation has a unique solution.

For our linear decay problem, rearrangement gives Θₙ₊₁ = Θₙ/(1 + λh). Since λ > 0 and h > 0, the denominator is positive and the multiplier lies strictly between zero and one. This particular update preserves positive excess and damps initial-value perturbations for any positive step that remains within the declared calculation interval.

At h = 1 s the first implicit value is 8/3 K; the exact value is (8 K)e^(−2), which is different. Numerical damping therefore does not certify accuracy. This favorable multiplier is a property of the supplied linear test problem, not a claim that every implicit method is stable or physically faithful on every equation.

Verify the computation and the model

Record the equation, initial data, units, state domain, time interval, method, step widths and arithmetic precision. Check each rate evaluation and update independently, and compare different step widths at the same final time. For our position model, Euler gives 3 m, 4 m and 4.5 m at 2 s with widths 1 s, 0.5 s and 0.25 s; the exact check value is 5 m.

Refinement is a useful check, but agreement between two computations alone is not an error certificate. A bound needs its smoothness, state-sensitivity and domain conditions. More steps also mean more arithmetic operations, so reducing h need not reduce the total error in a finite-precision calculation. A zero discrepancy in the numerical update merely shows that the chosen recurrence was followed.

An accurate solution of a supplied equation remains a result within that model. Testing a real system also requires evidence for its rate law, parameters, initial data and operating conditions. Retain the declared intervals even when algebraic formulas can be evaluated elsewhere. The next lesson, Numerical solutions of partial differential equations, adds spatial variation and boundary conditions to numerical modelling.

Worked examples

Try it yourself

1. For the position model, take two explicit Euler steps with h = 0.5 s from 0 s to 1 s. Find X₁ and X₂, then compare X₂ with x(1 s).

   Show the answer and reasoning

   X₁ = 1 m and X₂ = 1.5 m. The exact check position is 2 m, so e₂ = +0.5 m.

   The starting rates are 0 m/s and 1 m/s, giving increments 0 m and 0.5 m. The first step's zero starting rate does not make the continuous rate zero over that interval. Both time points lie inside the stated domain.

2. Suppose a position solution has continuous second derivative and |y″| ≤ 6 m/s² on a 0.2 s step. Starting from its exact state in exact arithmetic, what bound applies to |dₙ|? Is it automatically the accumulated error of a long computation?

   Show the answer and reasoning

   The one-step state-error bound is 0.12 m. It is not automatically a global-error bound.

   Mh²/2 = (6 m/s²)(0.2 s)²/2 = 0.12 m. The hypothesis restarts the step from the exact state. A sequence starting from previous approximate values must also account for inherited error and state-dependent rate changes.

3. Use the supplied conditions of the Euler error recurrence with L = 3 s⁻¹, h = 0.2 s, M = 6 m/s² and Eₙ ≤ 0.1 m. What upper bound follows for Eₙ₊₁?

   Show the answer and reasoning

   Eₙ₊₁ ≤ 0.28 m, under the supplied regularity, state-domain and exact-arithmetic conditions.

   1 + hL = 1.6. The inherited-error term is at most 1.6(0.1 m) = 0.16 m, and the new one-step term is at most 0.12 m. Their sum is 0.28 m. This upper bound does not require the actual error to equal it or exceed the earlier one.

4. For λ = 2 s⁻¹ and initial excess 8 K, explicit Euler uses h = 0.75 s. Find its multiplier and first value. Does damping guarantee that the exact positive sign is preserved?

   Show the answer and reasoning

   r = −0.5 and Θ₁ = −4 K. Its magnitude is damped, but its sign is not preserved.

   λh = 1.5 lies inside the strict damping interval (0,2). The negative multiplier changes the sign. The exact excess (8 K)e^(−1.5) is positive. This is a temperature-excess comparison, not a claim that the absolute temperature is negative.

5. For the decay model, take the two-stage predictor-and-average update with h = 1 s. Find the predictor, accepted value and multiplier. Is strict decay obtained?

   Show the answer and reasoning

   P = −8 K, the accepted excess is 8 K and its multiplier is 1. Strict decay is not obtained.

   k₁ = −16 K/s gives P = 8 K + (1 s)(−16 K/s) = −8 K. Then k₂ = +16 K/s. The accepted update adds half a second times their zero sum, leaving 8 K. The multiplier formula at λh = 2 also gives 1. Two rate evaluations do not remove this stability limitation.

6. The h = 1 s position values 1 m, 1 m, 3 m exactly satisfy the Euler recurrence. Does that prove they are exact solutions of the original ODE or that a real object follows it? What check can distinguish the numerical and exact values here?

   Show the answer and reasoning

   Neither conclusion follows. Substitution into the known check solution shows errors of 1 m and 2 m at the later grid points.

   Zero update discrepancy verifies the specified numerical recurrence. The stipulated continuous solution gives 2 m and 5 m at 1 s and 2 s, instead of the computed 1 m and 3 m. Refinement at the same final time offers another numerical check; physical validity would require independent observations and model evidence.

What to carry forward

• An ODE solver approximates an initial-value solution at declared times; its values are distinct from new measurements.
• Explicit Euler uses the current rate, and every update must retain compatible quantity units.
• One-step state error assumes an exact starting state; global error also contains inherited effects.
• Error propagation requires smoothness, state-sensitivity and domain assumptions, with further terms for other error sources.
• Damping, positive signs and accurate values are different properties even for a simple linear decay model.
• Check equal-time refinement, exact reference cases and model evidence before interpreting a computed trajectory.