Euler's Method

Learning Objectives

By the end of this section you will be able to:

1. Derive Euler's update rule directly from the definition of the derivative, without memorization.
2. Visualize the method as walking along tangent lines of the slope field.
3. Execute the algorithm by hand on a small initial-value problem and check the answer against the exact solution.
4. Distinguish local truncation error $O(h^2)$ from global error $O(h)$ and explain why one is bigger than the other.
5. Empirically verify the first-order convergence rate using a log-log plot.
6. Implement Euler's method in plain Python and as a vectorized PyTorch routine that supports batching and autograd.
7. Recognize the stability barrier that limits how big $h$ can be on stiff problems.

Why We Need Numerical Methods

In Section 20.3 we solved differential equations by separation of variables. That trick is wonderful when it works — but it works rarely. Most equations that describe real physical systems have no elementary closed-form solution at all. Try to analytically solve

and you will quickly discover there is no antiderivative we can write down with familiar functions. Yet a physicist still needs to know what $y(2.5)$ is. The way out is to give up on a formula and ask only for numbers: a table of $(t_i, y_i)$ pairs that approximate the true curve as closely as we want.

The numerical pivot. Instead of asking "what is the function y(t)?" we ask "starting from y(0), where will y be a tiny step later, and a tiny step after that, and after that…?"

Euler's method, published by Leonhard Euler in 1768 in his Institutionum Calculi Integralis, is the oldest and simplest answer to that question. Every modern ODE solver — from Runge–Kutta to the implicit BDF integrators that underlie scientific Python — can be understood as a smarter cousin of Euler's idea. Once you grasp Euler, the rest of the family becomes natural.

The Core Idea: Walk Along the Tangent

Imagine you are standing on a hillside in dense fog. You cannot see the valley, but the differential equation acts like a compass that tells you, at every spot you stand on, the exact direction of steepest descent. What do you do? You take one careful step in that direction, look at the compass again, take another step, and so on. Each step is a straight line, but stitched together they trace out a path that closely follows the true downhill curve.

That is precisely Euler's method. The differential equation $dy/dt = f(t, y)$ gives the slope of the unknown solution at every point, even though we do not know the solution itself. We choose a small step size $h$ and we walk along the tangent for that length, recompute the slope at the new point, walk again, recompute, repeat.

Deriving the Update Rule

Start from the very definition of the derivative:

For a small but nonzero $h$ we drop the limit and accept an approximation:

Multiply through by $h$ and rearrange:

Now let $y_n$ stand for our running approximation to $y(t_n)$ and let $t_{n+1} = t_n + h$. Replacing the unknown true values by our approximations we obtain Euler's rule:

It is worth pausing on what just happened. We replaced an infinite limit — the derivative — by a finite difference. That single substitution is the bridge between calculus and computation. Every numerical ODE method does some version of it; what differs is how cleverly each method recovers the curvature that this naive truncation throws away.

We can also reach the same formula from the first-order Taylor expansion:

Euler keeps the linear term $h\,y'(t_n)$ and discards everything from $h^2$ onward. That discarded tail is the local error per step — we will quantify it in a moment.

Geometric Picture

Recall the slope field from Section 20.2: at every point $(t, y)$ we draw a tiny arrow whose direction is the slope $f(t, y)$. The true solution through $(t_0, y_0)$ is the curve tangent to every arrow it passes through.

Euler's method is the polygonal-line approximation to that curve: instead of following the field continuously, we read the arrow once, stride $h$ units along it, read the arrow at our new location, stride again, and so on. The result is a piecewise-linear path that secretly lives entirely on tangent lines of the true solution family — but never quite on the true curve we want.

The Algorithm in Plain Steps

1. Inputs. The slope function $f(t, y)$, the initial point $(t_0, y_0)$, the stopping time $t_{\text{end}}$, and a step size $h > 0$.
2. Setup. Set the step count $N = \lceil (t_{\text{end}} - t_0)/h \rceil$, set $t \leftarrow t_0$, $y \leftarrow y_0$, and record $(t, y)$.
3. Repeat $N$ times.
   • Evaluate the slope at the current point: $k = f(t, y)$.
   • Take a tangent step: $y \leftarrow y + h \cdot k$.
   • Advance time: $t \leftarrow t + h$.
   • Record the new $(t, y)$.
4. Output. The list of pairs $(t_0, y_0), (t_1, y_1), \dots, (t_N, y_N)$ which discretely approximate the solution curve.

Interactive Explorer

The visualization below lets you take Euler steps yourself on six different equations. The green curve is the true solution; the orange polygon is Euler's approximation. Pull the step size slider toward zero and watch the orange path snap onto the green one. Press Step to take a single Euler step at a time — the dashed red line is the tangent actually used to take the next step.

Worked Example (by hand)

We pick an equation whose true solution we happen to know, so we can score each Euler step directly.

Its analytic solution (verifiable by direct substitution) is $y(t) = (t + 1)^2 - \tfrac{1}{2} e^t$. We will compute four Euler steps to reach $t = 2$ and compare against the truth.

Here is a summary of the trace, computed and verified in Python:

How Wrong Is It? Local and Global Error

Two different errors are at play and it is essential to keep them apart.

Local truncation error

Suppose by some miracle we are exactly on the true curve at time $t_n$. How wrong is a single Euler step? From the Taylor expansion we already wrote,

So per step, the error is $O(h^2)$. Halve $h$ and a single step becomes four times more accurate.

Global error

But we take roughly $N = (t_{\text{end}} - t_0)/h$ steps, and each step injects its own error into the input of the next step. A careful argument (using the Lipschitz continuity of $f$) shows that all those $O(h^2)$ contributions add up to an accumulated error

where the constant $C$ depends on the problem (specifically, on the size of $y''$ and on the Lipschitz constant of $f$). This is why we say Euler is a first-order method: global error scales with the first power of $h$.

Convergence Demo: Order of Accuracy

The convergence claim above is testable. The widget below solves the chosen IVP with successively halved step sizes and plots $\log|\text{error}|$ against $\log h$. A first-order method must land on a line of slope 1.

The empirically fitted slope (green) sits right next to the reference slope-1 line (purple dashed) for every problem. This is experimental confirmation of the theory: Euler is first-order, no more and no less.

Stability: When Small Steps Are Still Not Small Enough

Accuracy is not the only concern. Some equations punish big $h$ not by being slightly wrong, but by producing solutions that explode to infinity even though the true solution is calm. The cleanest example is the linear decay equation

whose true solution $y(t) = y_0 e^{-\lambda t}$ decays smoothly to zero. Applying Euler:

So $y_n = (1 - \lambda h)^n y_0$. For the Euler solution to also decay, we need $|1 - \lambda h| < 1$, which requires

Push $h$ past this barrier and Euler oscillates with growing amplitude — a complete failure mode unrelated to accuracy. Equations with very large $\lambda$ (or, more generally, with very fast and very slow dynamics mixed together) are called stiff, and explicit Euler is a poor tool for them. Implicit methods (Section 20.x and beyond) cure this.

Python Implementation

Below is the cleanest possible Python implementation, with every line annotated. Read the explanations to the right — they explain not just what each line does but why the algorithm needs it.

t0=0, t_end=1, h=0.25 → n_steps = 4

If y' = y and we're at (t=0, y=1), then slope = 1.

y_new = 1 + 0.25 × 1 = 1.25

1# Plain-Python Euler's method
2# Solves the IVP:  dy/dt = f(t, y),  y(t0) = y0
3# Returns the (t, y) pairs produced by walking along the slope field.
4
5def euler(f, t0, y0, t_end, h):
6    """
7    f      : callable f(t, y) returning the slope
8    t0, y0 : initial condition  y(t0) = y0
9    t_end  : final time
10    h      : step size (must be > 0 and small enough)
11    """
12    n_steps = int(round((t_end - t0) / h))   # number of Euler steps
13    ts = [t0]                                # collected time grid
14    ys = [y0]                                # collected solution values
15
16    t = t0
17    y = y0
18    for _ in range(n_steps):
19        slope = f(t, y)        # 1. read the slope at the current point
20        y = y + h * slope      # 2. take a tangent step of length h
21        t = t + h              # 3. advance time
22        ts.append(t)
23        ys.append(y)
24
25    return ts, ys
26
27
28# --- demo: solve  y' = y,  y(0) = 1  on [0, 1]
29def f(t, y):
30    return y                  # exponential growth: slope = current value
31
32ts, ys = euler(f, t0=0.0, y0=1.0, t_end=1.0, h=0.25)
33for t, y in zip(ts, ys):
34    print(f"t={t:.2f}  y_euler={y:.6f}")
35
36# True solution is y(t) = exp(t), so y(1) = e = 2.71828...
37print("Euler at t=1 :", ys[-1])
38print("True  at t=1 :", 2.718281828459045)

And here is the experiment that confirms the first-order rate — the same one whose plot you played with above, but now as code you can actually run:

On y'=y from t=0..1, err/h converges to about 1.35 — that's the empirical C.

1# Measuring the error of Euler as h shrinks.
2# We solve  y' = y,  y(0) = 1  on [0, 1].
3# The true answer is  y(1) = e = 2.718281828459045
4
5import math
6
7def f(t, y):
8    return y
9
10def euler_final(f, t0, y0, t_end, n):
11    h = (t_end - t0) / n
12    t, y = t0, y0
13    for _ in range(n):
14        y = y + h * f(t, y)
15        t = t + h
16    return y
17
18print(f"{'n':>5} {'h':>10} {'y_Euler':>12} {'|error|':>12} {'err/h':>10}")
19print("-" * 54)
20
21true_value = math.e
22for n in [4, 8, 16, 32, 64, 128, 256]:
23    h        = 1.0 / n
24    y_hat    = euler_final(f, 0.0, 1.0, 1.0, n)
25    error    = abs(true_value - y_hat)
26    # Empirical "constant" in the bound  error ≈ C · h
27    constant = error / h
28    print(f"{n:>5d} {h:>10.5f} {y_hat:>12.6f} {error:>12.6e} {constant:>10.4f}")

When you run it you should see the $|\text{error}|$ column roughly halve every time $n$ doubles, and the rightmost column (the empirical $C = |\text{error}|/h$) converge to a single constant — the punchline of the whole error analysis.

PyTorch / Vectorized Implementation

The plain version is great for understanding but slow for large problems. The PyTorch version below is structurally identical — same update rule, same loop — but every value is a tensor. Three real advantages drop out for free:

• Batching. Solve a whole population of trajectories (different initial conditions, different parameters) in a single call.
• GPU. Move tensors to cuda and the same code runs there.
• Autograd. Because every step is a differentiable tensor op, we can differentiate the final state with respect to anything that fed into it — initial conditions, model parameters, the integration time. This is exactly the trick that makes Neural ODEs trainable.

If B=4, one tensor add replaces four scalar updates.

For y'=-0.5 y the true sensitivity is e^(-0.5·2) ≈ 0.3679; Euler reproduces it to four decimal places.

1# Vectorized Euler in PyTorch.
2# Same idea as the plain-Python version, but:
3#   • everything is a torch.Tensor, so it runs on GPU when you ask it to
4#   • we solve a batch of trajectories in parallel (different initial values)
5#   • autograd is preserved, so the *final state* can be differentiated
6#     with respect to the initial condition or any parameter of f.
7
8import torch
9
10def euler_batched(f, t0, y0, t_end, h):
11    """
12    f     : callable f(t: float, y: Tensor[B, D]) -> Tensor[B, D]
13    y0    : Tensor of shape [B, D]  — B independent initial conditions
14    Returns y(t_end) as a Tensor of shape [B, D].
15    """
16    n_steps = int(round((t_end - t0) / h))
17    t = torch.tensor(t0, dtype=y0.dtype, device=y0.device)
18    y = y0
19    for _ in range(n_steps):
20        y = y + h * f(t, y)
21        t = t + h
22    return y
23
24
25# --- demo: solve  y' = -lambda * y  (decay) for a batch of lambdas
26lambdas = torch.tensor([0.1, 0.5, 1.0, 2.0]).view(-1, 1)   # shape [4, 1]
27
28def f(_t, y):
29    # Element-wise: each row uses its own lambda
30    return -lambdas * y
31
32y0 = torch.ones(4, 1)               # all start at y=1
33y_T = euler_batched(f, t0=0.0, y0=y0, t_end=2.0, h=0.01)
34
35true = torch.exp(-lambdas * 2.0)    # closed form for comparison
36print("Euler  y(2) :", y_T.flatten().tolist())
37print("Truth  y(2) :", true.flatten().tolist())
38print("Errors      :", (y_T - true).abs().flatten().tolist())
39
40
41# --- bonus: differentiate the final state w.r.t. an initial condition
42y0 = torch.tensor([[1.0]], requires_grad=True)
43
44def f_single(_t, y):
45    return -0.5 * y
46
47y_T = euler_batched(f_single, 0.0, y0, 2.0, 0.01)   # y_T depends on y0
48y_T.sum().backward()                                # populate y0.grad
49print("dy(2)/dy(0) ≈", y0.grad.item())              # ≈ exp(-0.5*2) = 0.3679

Beyond Euler: A Glimpse at RK4 and Neural ODEs

Once you see that Euler is just "estimate the slope, take a step", every improvement is easy to anticipate:

Every one of these is a refinement on the same skeleton you just coded. The whole field of numerical analysis for ODEs is the search for ways to extract more accuracy or more stability per slope evaluation than Euler manages.

Common Pitfalls

Summary

Euler's method is the canonical bridge from calculus to computation. Its update rule

falls out directly from the definition of the derivative; geometrically it walks along tangent lines of the slope field. The cost of this simplicity is first-order accuracy — global error scales linearly with the step size — and a stability barrier on stiff equations. Yet every fancier ODE solver is a refinement of the same skeleton, and the PyTorch version is the first step toward Neural ODEs and the modern continuous-depth models built on top of them.

Take home. If you remember one thing, remember this rhythm: read the slope, walk a step, repeat. Everything beyond Euler is just how to read the slope more cleverly.