Euler's method, step size and error, improved Euler, Runge-Kutta, numerical stability, and systems — the computational tools that make differential equations solvable in practice
This chapter introduces numerical methods for differential equations, a collection of techniques for approximating solutions when exact symbolic formulas are difficult or impossible to obtain. Students will learn Euler's method as the foundational algorithm, study how step size affects accuracy, build and read tables of numerical solutions, explore the improved Euler method and the Runge-Kutta method as more accurate alternatives, and understand what numerical stability means and why it matters.
By the end of this chapter, students should understand that numerical methods are not a fallback for students who cannot solve equations analytically. They are essential tools in modern science, engineering, computing, data modeling, biology, economics, and every field where differential equations describe real systems, because most real differential equations do not yield to simple symbolic techniques.
A population of bacteria grows according to the differential equation dP/dt = rP(1 − P/K) − H(t), where r is the intrinsic growth rate, K is the carrying capacity, and H(t) is a harvesting function that varies with time in a complicated way. A biologist wants to know how the population behaves over the next twelve months. There is no general formula that solves this equation exactly for an arbitrary harvesting function H(t). The equation cannot be separated in the usual way. An integrating factor will not help. Laplace transforms will not simplify it cleanly. What does the biologist do? This is where numerical methods begin — not as an admission of failure, but as the right tool for the problem.
Calculus gives us an extraordinarily powerful language for describing change. A differential equation captures a relationship involving a function and its rate of change. What calculus also teaches us, however, is that integration — the central tool for recovering functions from their rates of change — does not always produce a tidy formula.
Some integrals simply do not evaluate to elementary functions. Some differential equations involve functions of both the unknown and the independent variable in ways that prevent separation, linearization, or standard transformation. Some real systems change in ways that depend on external inputs, historical data, or environmental conditions that cannot be captured in a smooth symbolic expression.
In all of these situations, we cannot solve the differential equation by finding a formula for the solution. But we can still ask: at each step in time, what does the rate of change tell us about where the solution is going next? That is the intuition behind every numerical method for differential equations. Instead of finding the whole solution at once, we walk along the solution curve one small step at a time, using the differential equation itself to tell us the direction and slope at each point.
Suppose we have an initial value problem: dy/dx = f(x, y), y(x₀) = y₀. We know the solution passes through (x₀, y₀). We know the slope there is f(x₀, y₀). If we take a small step to x₁ = x₀ + h, we can estimate the solution at x₁ by moving along the tangent line. Step by step, we build an approximate portrait of the solution.
A numerical method does not produce a formula. It produces a sequence of points: (x₀, y₀), (x₁, y₁), (x₂, y₂), …, (xₙ, yₙ). These pairs approximate the solution at discrete points. We can tabulate them, plot them, and use them to answer specific questions about the system.
An exact solution to a differential equation is a function y = φ(x) that satisfies the equation identically for all values of x in some interval. When an exact solution exists, it gives complete information: we can evaluate it at any point, differentiate it, integrate it, and analyze its symbolic form.
Many equations that arise naturally in practice, however, do not yield exact solutions in elementary functions. Consider the equation dy/dx = e−x² + y². The function e−x² alone has no elementary antiderivative. Together with y², this equation resists every standard symbolic method.
Numerical methods allow us to compute approximate values of solutions even when no formula exists. The approximation is not a failure; it is a carefully controlled estimate with quantifiable error.
| Source of Error | Description | Main Control |
|---|---|---|
| Truncation error | Approximating the true curve with a simpler model (e.g., a tangent line instead of the actual curve) | Step size h; method order |
| Rounding error | Finite precision of computer arithmetic; accumulation across many operations | Usually small in practice; double-precision arithmetic |
A numerical solution that appears to produce numbers is not automatically reliable. We need to know whether those numbers are close to the true solution and how we can make them closer if needed. Understanding error is not optional — it is the core of responsible numerical computation.
Euler's method is the most fundamental numerical algorithm for ordinary differential equations. It is simple enough to carry out by hand for small examples, transparent enough to understand completely, and important enough to serve as the conceptual foundation for all more sophisticated methods.
Given dy/dx = f(x, y) with initial condition y(x₀) = y₀ and step size h, set xₙ = x₀ + nh and iterate:
y_{n+1} = yₙ + h · f(xₙ, yₙ)
At each step: multiply the current slope by the step size, then add the result to the current value. Repeat from the new point.
Euler's method is essentially a sequence of tangent-line approximations. At each point on the approximate solution path, we draw the tangent line given by the differential equation, move along it for a horizontal distance h, and use the resulting point as the starting position for the next step.
Draw a smooth curve representing the exact solution y = φ(x), starting at (x₀, y₀). At that point, draw the tangent line with slope f(x₀, y₀). Mark the point (x₁, y₁) on the tangent line, slightly above or below the true curve. From (x₁, y₁), draw a new tangent segment with slope f(x₁, y₁), leading to (x₂, y₂). Show how the Euler path — a sequence of line segments — follows but gradually diverges from the exact curve as x increases.
Step 1:
x₁ = 0.25
y₁ = 1 + 0.25 · f(0, 1) = 1 + 0.25 · (0 + 1) = 1.25
Step 2:
x₂ = 0.50
y₂ = 1.25 + 0.25 · f(0.25, 1.25) = 1.25 + 0.25 · 1.5 = 1.625
Step 3:
x₃ = 0.75
y₃ = 1.625 + 0.25 · f(0.5, 1.625) = 1.625 + 0.25 · 2.125 = 2.15625
Step 4:
x₄ = 1.00
y₄ = 2.15625 + 0.25 · f(0.75, 2.15625) = 2.15625 + 0.25 · 2.90625 = 2.88281
The choice of step size h is one of the most important decisions in any numerical computation. It determines the spacing between the x-values at which we compute approximate solution values, and it directly controls both the accuracy of the approximation and the computational effort required.
For Euler's method, the error at each individual step is proportional to h², and because we take approximately 1/h steps over a fixed interval, the total accumulated error over the interval is proportional to h. This means that halving the step size approximately halves the total error — a useful rule of thumb for checking whether a computation is behaving correctly.
There is a trade-off. Halving the step size doubles the number of steps. Modern solvers use adaptive step-size control, automatically making steps smaller in regions where the solution changes rapidly and larger in smoother regions.
| Situation | Guidance |
|---|---|
| Solution changes slowly on a long interval | A larger step size may be adequate for a quick check |
| Solution changes rapidly or is concave | Use a smaller step size; the tangent line drifts quickly |
| Checking reliability | Run with h and h/2; if results agree closely, step size is adequate |
| High-accuracy requirement | Use a higher-order method (improved Euler or RK4) rather than just shrinking h |
Do not assume that a numerical answer is correct simply because the computer produced it quickly and with many decimal places. More decimal places do not mean more accuracy. Only understanding the error behavior tells you whether a numerical result can be trusted.
A table is the natural way to record and communicate a numerical solution. Each row represents one step of the computation, showing the current x-value, the approximate y-value, and the slope used to advance to the next step. Including the slope column makes the computation transparent and checkable: if a row contains an unexpected value, it is easy to trace back to where the error occurred.
| n | xₙ | yₙ | f(xₙ, yₙ) | h · f |
|---|---|---|---|---|
| 0 | x₀ | y₀ | slope at step 0 | increment |
| 1 | x₁ | y₁ | slope at step 1 | increment |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| n | xₙ | yₙ | slope at step n | — |
n xₙ yₙ f(xₙ, yₙ) = 2xₙ − yₙ h·f
0 0.0 0.00000 0.00000 0.00000
1 0.1 0.00000 0.20000 0.02000
2 0.2 0.02000 0.38000 0.03800
3 0.3 0.05800 0.54200 0.05420
4 0.4 0.11220 0.68780 0.06878
5 0.5 0.18098 — —
A table of values gives numbers. A graph gives shape. After computing a table of numerical approximations, plotting the points (xₙ, yₙ) and connecting them with straight segments gives a visual picture of the approximate solution. This picture can reveal qualitative behavior that numbers alone obscure: whether the solution is increasing or decreasing, whether it is accelerating or leveling off, and whether it appears to be heading toward a finite limit or growing without bound.
Draw a coordinate plane with x running from 0 to 1 and y running from 0 to 4. Plot the Euler approximation from Example 11.3.1 as a connected sequence of straight segments through the five computed points. Then draw a smooth curve representing the exact solution y = 2eˣ − x − 1 on the same axes. Label both curves. Show how the Euler path hugs the exact curve at first, then gradually diverges below it as x increases, because the solution is concave up and Euler's method uses left-endpoint slopes.
Check that the graph agrees with the slope field of the equation. Verify that qualitative features — growth, decay, equilibrium — match expectations. Watch for wild oscillations: if the graph oscillates or grows unexpectedly for an equation that should be smooth, suspect numerical instability (Section 11.10), not genuine solution behavior.
When no exact solution is available, compare numerical solutions computed with different step sizes. If two numerical solutions computed with h and h/2 agree closely over the whole interval, we have good reason to trust the approximation. If they differ significantly, reduce the step size further before drawing conclusions about the solution's behavior.
Understanding error is the core of responsible numerical computation. A numerical result without an error estimate is incomplete information.
The local truncation error at a single step is the error introduced by taking one Euler step from a known exact point. Using the Taylor series expansion of y(xₙ₊₁), the Euler step uses only the first two terms of the series. The local truncation error is therefore approximately h²/2 · y″(xₙ), which is proportional to h².
The global error at a point xₙ is the total accumulated difference between the numerical approximation and the exact solution. Since we take approximately (xₙ − x₀)/h steps, each contributing a local error of order h², the total global error is proportional to h. Euler's method is therefore called a first-order method: halving the step size approximately halves the global error.
Local truncation error ≈ (h²/2) · y″(xₙ) → order h²
Global error ≈ C · h → order h (first-order method)
Halving h approximately halves the global error. To gain one more decimal place of accuracy requires approximately ten times as many steps.
h = 0.5 (two steps):
y₁ = 1 + 0.5·1 = 1.5
y₂ = 1.5 + 0.5·1.5 = 2.25
Error: 2.71828 − 2.25 = 0.46828
h = 0.25 (four steps):
y₁ = 1.25, y₂ = 1.5625, y₃ = 1.95313, y₄ = 2.44141
Error: 2.71828 − 2.44141 = 0.27687
h = 0.1 (ten steps, formula: y₁₀ = (1.1)¹⁰):
y₁₀ = 2.59374
Error: 2.71828 − 2.59374 = 0.12454
Step size h Approximation Error
0.5 2.25000 0.468
0.25 2.44141 0.277
0.1 2.59374 0.125
Error ratio from h = 0.5 to h = 0.25: ≈ 0.469/0.277 ≈ 1.69, close to 2. This confirms first-order global error behavior.Euler's method uses the slope at the left endpoint of each step to advance the solution. The improved Euler method, also called the Heun method or the predictor-corrector method, does better by using both the beginning and end of each step.
Instead of using only the slope at (xₙ, yₙ), we average two slopes: the slope at the beginning of the step, and the slope at the predicted end of the step. The two slopes are then averaged to get a better estimate of the direction across the interval.
Predictor: ỹ_{n+1} = yₙ + h · f(xₙ, yₙ)
Corrector: y_{n+1} = yₙ + (h/2)[f(xₙ, yₙ) + f(x_{n+1}, ỹ_{n+1})]
The predictor uses one Euler step to estimate the next location. The corrector evaluates the slope at that prediction and averages it with the slope at the current point. Global error: order h² (second-order method).
Step 1: x₀ = 0, y₀ = 1
k₁ = f(0, 1) = 1
ỹ₁ = 1 + 0.25·1 = 1.25
k₂ = f(0.25, 1.25) = 1.5
y₁ = 1 + (0.25/2)(1 + 1.5) = 1 + 0.3125 = 1.3125
Step 2: x₁ = 0.25, y₁ = 1.3125
k₁ = f(0.25, 1.3125) = 1.5625
ỹ₂ = 1.3125 + 0.25·1.5625 = 1.703125
k₂ = f(0.5, 1.703125) = 2.203125
y₂ = 1.3125 + (0.25/2)(1.5625 + 2.203125) = 1.783203
Step 3: x₂ = 0.50, y₂ = 1.783203
k₁ = 2.283203, ỹ₃ = 2.354004, k₂ = 3.104004
y₃ = 2.456604
Step 4: x₃ = 0.75, y₃ = 2.456604
k₁ = 3.206604, ỹ₄ = 3.258255, k₂ = 4.258255
y₄ = 3.389711
Method Approx. at x=1 Error
Euler (h=0.25) 2.88281 0.554
Improved Euler 3.38971 0.047
Exact 3.43656 —
The improved Euler method averages two slope estimates. The Runge-Kutta family of methods extends this idea by sampling the slope at four carefully chosen points within each step and combining them with particular weights. The most widely used is the fourth-order Runge-Kutta method, universally abbreviated RK4. It is the workhorse of scientific computing.
Given dy/dx = f(x, y) and current values (xₙ, yₙ), compute four slope estimates:
k₁ = f(xₙ, yₙ)
k₂ = f(xₙ + h/2, yₙ + (h/2)k₁)
k₃ = f(xₙ + h/2, yₙ + (h/2)k₂)
k₄ = f(xₙ + h, yₙ + h·k₃)
Then advance: y_{n+1} = yₙ + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)
Global error: order h⁴ (fourth-order method). Halving h reduces global error by a factor of 16.
The weights 1, 2, 2, 1 (summing to 6) are not arbitrary. They come from a careful derivation that matches the Taylor series of the exact solution up to terms involving h⁴. The first and fourth estimates bracket the interval. The second and third probe the midpoint. By weighting the midpoint estimates more heavily, RK4 corrects for curvature in a way that basic Euler and improved Euler cannot.
Step 1: x₀ = 0, y₀ = 1, h = 0.5
k₁ = f(0, 1) = 1
k₂ = f(0.25, 1.25) = 1.25
k₃ = f(0.25, 1.3125) = 1.3125
k₄ = f(0.5, 1.65625) = 1.65625
y₁ = 1 + (0.5/6)(1 + 2·1.25 + 2·1.3125 + 1.65625)
= 1 + (0.5/6)(7.78125)
= 1 + 0.64844 = 1.64844
Step 2: x₁ = 0.5, y₁ = 1.64844, h = 0.5
k₁ = 1.64844
k₂ = f(0.75, 2.06055) = 2.06055
k₃ = f(0.75, 2.16358) = 2.16358
k₄ = f(1.0, 2.73023) = 2.73023
y₂ = 1.64844 + (0.5/6)(1.64844 + 4.12110 + 4.32716 + 2.73023)
= 1.64844 + 1.06891 = 2.71735
Method h Approx. Error
Euler 0.5 2.250 0.468
Improved Euler 0.5 — better than Euler
RK4 0.5 2.717 0.001
Exact — 2.71828 —
| Method | Stages per step | Global error order | Effect of halving h |
|---|---|---|---|
| Euler | 1 | O(h) | Error halves |
| Improved Euler (Heun) | 2 | O(h²) | Error divides by 4 |
| RK4 | 4 | O(h⁴) | Error divides by 16 |
Accuracy and stability are different concepts, and confusing them is one of the most serious mistakes in numerical computation. A method can be accurate for small step sizes and still become violently unstable for larger step sizes. Stability is about whether errors introduced at early steps grow or shrink as the computation proceeds.
Suppose we are solving an equation whose exact solution is a smooth, slowly decaying function approaching zero as x increases. If the numerical method is unstable, the computed approximation may instead oscillate wildly or grow exponentially, even though the true solution is quiet and decaying.
Draw two graphs on the same axes. In the first, label it "exact solution" and draw a smooth, gently decaying curve approaching zero. In the second, label it "unstable numerical approximation" and draw a wildly oscillating curve that grows in amplitude rather than decaying. Mark the true equilibrium with a dashed horizontal line at y = 0. The contrast shows that instability produces behavior that looks nothing like the true solution — the numbers are not just inaccurate, they are qualitatively wrong.
For the test equation dy/dx = λy (real, negative λ), the Euler iteration is y_{n+1} = (1 + hλ)yₙ. For the approximation to remain bounded:
|1 + hλ| ≤ 1 → h ≤ 2/|λ|
If |λ| is large (rapidly decaying equation), the required step size is very small. This is the stability constraint for Euler's method.
A stiff differential equation is one that contains components decaying at very different rates. Some parts of the solution change rapidly while others change slowly. To capture the fast components stably, the step size must be very small, even when the behavior of primary interest is slow. Stiff equations arise frequently in chemistry, electrical circuits, and ecology. Specialized methods called implicit methods or stiff solvers (such as ode15s in MATLAB or Rodas4 in Julia) are designed to handle these cases.
If a numerical solution behaves wildly, try reducing the step size before concluding that the differential equation itself produces wild behavior. If reducing the step size makes the wild behavior disappear, it was numerical instability, not a property of the equation. If the step size required for stability is impractically small, the equation may be stiff and require a different solver.
The numerical methods we have studied for a single equation extend naturally to systems. This is one of the most important extensions, because most real-world models involve multiple interacting variables, each governed by its own differential equation.
For a system of two equations dx/dt = f(t, x, y) and dy/dt = g(t, x, y), the Euler update is applied simultaneously to both variables at each step:
x_{n+1} = xₙ + h · f(tₙ, xₙ, yₙ)
y_{n+1} = yₙ + h · g(tₙ, xₙ, yₙ)
Both updates use the current values (tₙ, xₙ, yₙ) — not the newly computed xₙ₊₁. Compute all increments first, then apply all updates together.
Step 1: t₀ = 0, x₀ = 1, y₀ = 1
f(0, 1, 1) = 2(1) − (1)(1) = 1
g(0, 1, 1) = −(1) + (1)(1) = 0
x₁ = 1 + 0.1·1 = 1.1
y₁ = 1 + 0.1·0 = 1.0
Step 2: t₁ = 0.1, x₁ = 1.1, y₁ = 1.0
f(0.1, 1.1, 1.0) = 2(1.1) − (1.1)(1.0) = 1.1
g(0.1, 1.1, 1.0) = −(1.0) + (1.1)(1.0) = 0.1
x₂ = 1.1 + 0.1·1.1 = 1.21
y₂ = 1.0 + 0.1·0.1 = 1.01
Every modern scientific computing environment provides high-quality numerical differential equation solvers implementing adaptive Runge-Kutta methods, stiff solvers, and many other algorithms automatically. Understanding what these solvers do is essential for using them intelligently. The solver is not a black box — it is an implementation of the ideas studied in this chapter.
In Python, scipy.integrate.solve_ivp solves initial value problems numerically. By default it uses an adaptive Runge-Kutta method of order 4 and 5 (RK45), which automatically adjusts the step size to achieve a specified tolerance.
MATLAB offers ode45 for non-stiff equations and ode15s for stiff equations. Julia's DifferentialEquations.jl package is one of the most advanced numerical libraries available, with hundreds of algorithms and automatic stiffness detection.
Technology solves the computational problem. It does not solve the modeling problem. A solver cannot tell you what equation to write, verify that your units are correct, or interpret the output for you. The student who understands Euler's method, RK4, error, and stability can use these tools wisely: checking whether output is plausible, understanding what error tolerance to request, and diagnosing instability when it appears. The student who treats a solver as a calculator, entering equations without understanding, is one modeling mistake away from a confidently computed wrong answer.
Correct: Accuracy depends on h, the method's order, and solution behavior — not on decimal display precision. Always accompany a numerical result with the step size and method used.
Correct: Compute with both h and h/2, then compare. If results agree closely, the step size is adequate. If they differ significantly, reduce h further.
Correct: Compute all slope values at the current state (xₙ, yₙ) before updating any variable. Store increments, then apply all updates together.
Correct: Reduce the step size. If the wild behavior disappears, it was numerical instability, not a property of the equation. Consider a stiff solver if instability persists with small h.
Correct: Before running a multi-step computation, verify f by substituting one or two known values and checking that the result matches expectations. If an exact solution exists, verify that f(x, y_exact) = y'_exact at one or two points.
Correct: Before computing, form expectations using the slope field, equilibrium analysis, or physical reasoning. After computing, compare the numerical result to these expectations and explain any discrepancy.
Correct: Recognize stiffness (large negative eigenvalue, multiple time scales) and switch to an appropriate stiff solver rather than trying to force a standard method to work with tiny step sizes.
Chapter 11 has introduced the computational side of differential equations. We have seen that when symbolic methods fail, numerical methods provide a way to compute approximate solutions step by step, to build tables, to produce graphs, and to study systems of equations.
But computation and approximation are not the only alternatives to symbolic solution. There is another approach which does not produce numerical values at all but instead describes the behavior of solutions: whether they increase or decrease, what equilibria exist, which equilibria are stable, and how long-term behavior depends on initial conditions. This approach is called qualitative analysis.
Chapter 12 will introduce qualitative methods for differential equations, focusing on autonomous equations, phase lines, equilibrium stability, and nonlinear systems. Qualitative methods ask not "what is the value of the solution at time t?" but rather "what does the solution do in the long run?"
Together, the numerical methods of Chapter 11 and the qualitative methods of Chapter 12 give us two powerful ways to study differential equations when exact symbolic formulas are unavailable. The boundary between these approaches is not sharp: numerical phase plane plots are among the most powerful tools for visualizing qualitative behavior, and qualitative stability theory explains why some numerical computations are reliable and others are not. As you move into Chapter 12, carry your numerical intuition with you. The two chapters form a unity.
Use Euler's method to approximate the solution of each initial value problem at the indicated point. Show all steps.
dy/dx = 2y, y(0) = 1, step h = 0.5. Approximate y(2).dy/dx = x − y, y(0) = 0, step h = 0.1. Approximate y(0.4).dy/dx = y², y(0) = 0.5, step h = 0.2. Approximate y(0.6).dy/dx = sin(x) + y, y(0) = 1, step h = 0.25. Approximate y(1).dy/dx = x/y, y(0) = 1, step h = 0.5. Approximate y(1.5).dP/dt = 0.3P, P(0) = 100, step h = 1. Approximate P(4). The exact solution is P = 100e^{0.3t}. Compare approximation with exact value and compute the error.For each problem, define variables and units, write the differential equation, apply the specified numerical method, and interpret your results.
Use Python, MATLAB, a graphing calculator, or any available numerical solver. Focus on interpreting results, not merely reporting numbers.
Each problem shows a student's work containing one or more errors. Identify the error, explain why it is wrong, and show the correct computation.
What did the student do wrong? What is the correct computation?
What is wrong with this claim? What would you need to do to actually assess accuracy?
What did the student do wrong when updating y₁? What is the correct computation?
Identify the error in this conclusion. What is the true long-term behavior, and what caused the numerical result?
Describe a strategy to catch this error before committing to a long computation. What check should the student perform first?
dy/dx = 1 − y, y(0) = 0, and complete the table for x = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0. The exact solution is y = 1 − e⁻ˣ. Compute the exact value at x = 1 and find the error of your Euler approximation.dP/dt = 0.4P, P(0) = 10, for two steps. What does each step represent in terms of population growth?dy/dx = 2 − y, y(0) = 0. Approximate y(1) using two steps. Show the predictor and corrector steps explicitly.dy/dx = x² + y, y(0) = 0. Show all four k-values and the final approximation for y(0.5).dy/dx = 3y, y(0) = 1, to y_{n+1} = y_n + h·3y_{n+1}, using the new value on the right side. Explain the conceptual error and write the correct Euler formula for this equation.Chapter 11 has given us the tools to compute approximate solutions step by step. Chapter 12 will give us the tools to see the structure of solutions without computing their values at all. Qualitative analysis asks questions about behavior rather than specific numbers: where are the equilibria? Which are stable? What does the solution do in the long run, regardless of initial conditions? These questions are answered by reading the structure of the differential equation through phase lines, stability analysis, and qualitative portraits of nonlinear systems. Numerical and qualitative methods complement each other naturally: compute a numerical solution to see approximate behavior, and use qualitative analysis to understand whether that behavior is typical, stable, or dependent on initial conditions. Together, the two chapters form the essential toolkit for studying differential equations when exact formulas are unavailable.