ODEs and the Euler method

An ordinary differential equation (ODE) specifies the rate of change of a quantity. The derivative and an initial condition determine the quantity itself. The simplest form is

\frac{d y}{d t} = f(y, t), \qquad y(0) = y_0.

At each time

t

, the function

f

gives the slope of

y

at the current state

(y, t)

. Solving the ODE means finding the function

y(t)

that satisfies this relation.

Exponential growth

Consider the ODE

\frac{d y}{d t} = y, \qquad y(0) = 1.

The rate of change of

y

equals

y

itself. This law describes populations, compound interest, and unstable physical systems in some regimes. The exponential has this property and satisfies the initial condition, so

y(t) = e^t.

This is the analytical solution. Most ODEs encountered in practice have no analytical solution, so you must compute

y(t)

numerically. The simplest numerical method is the Euler method.

The Euler method

The Euler method approximates the solution by taking repeated short steps in the direction of the derivative. Choose a small step size

\Delta t

. Starting from

(t_0, y_0) = (0, 1)

, replace the continuous ODE with the discrete update

y_{n+1} = y_n + \Delta t \cdot f(y_n, t_n), \qquad t_{n+1} = t_n + \Delta t.

At each step, evaluate the slope at the current point and move a distance

\Delta t

in that direction. After

N

steps, the result approximates

y

at time

t = N \Delta t

. The following Python function implements this update.

import numpy as np
import matplotlib.pyplot as plt


def euler(f, y0, t0, t_end, dt):
    """Euler method: integrate dy/dt = f(y, t) from t0 to t_end.

    Returns arrays (t_grid, y_trajectory).
    """
    n_steps = int(np.round((t_end - t0) / dt))
    t = np.empty(n_steps + 1)
    y = np.empty(n_steps + 1)
    t[0], y[0] = t0, y0
    for n in range(n_steps):
        y[n + 1] = y[n] + dt * f(y[n], t[n])
        t[n + 1] = t[n] + dt
    return t, y


# Exponential growth: dy/dt = y
f_exp = lambda y, t: y

t_euler, y_euler = euler(f_exp, y0=1.0, t0=0.0, t_end=2.0, dt=0.2)
t_dense = np.linspace(0, 2, 200)
y_exact = np.exp(t_dense)

print(f"Euler at t=2.0 (dt=0.2):  {y_euler[-1]:.4f}")
print(f"Exact at t=2.0:           {np.exp(2.0):.4f}")
print(f"Relative error:           {abs(y_euler[-1] - np.exp(2.0)) / np.exp(2.0):.2%}")

Euler at t=2.0 (dt=0.2):  6.1917
Exact at t=2.0:           7.3891
Relative error:           16.20%

Show plotting code

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(t_dense, y_exact, "k-", lw=2, label="exact: $y(t)=e^t$")
ax.plot(t_euler, y_euler, "o-", color="C0", label="Euler, $\\Delta t=0.2$")
ax.set_xlabel("t"); ax.set_ylabel("y")
ax.set_title("Euler vs exact solution of $dy/dt = y$")
ax.legend()
ax.grid(alpha=0.3)
plt.show()

Euler undershoots the exact curve. With

\Delta t = 0.2

, each step uses the slope at the start of the interval. The true slope grows over the interval, so the discrete update consistently lags behind. The next experiment shows how shrinking

\Delta t

changes the error and how quickly the error approaches zero as

\Delta t \to 0

Step size and error

The Euler method is first-order accurate. The global error at a fixed final time scales linearly with

\Delta t

. Halving

\Delta t

roughly halves the error. You can verify this by running Euler at several step sizes and plotting the error at

t = 2

on a log-log axis. A first-order method produces a straight line with slope 1.

dts = np.array([0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005])
errors = []
for dt in dts:
    _, y = euler(f_exp, y0=1.0, t0=0.0, t_end=2.0, dt=dt)
    errors.append(abs(y[-1] - np.exp(2.0)))
errors = np.array(errors)

for dt, err in zip(dts, errors):
    print(f"dt={dt:6.3f}   |error| at t=2.0 = {err:.4e}")

dt= 0.500   |error| at t=2.0 = 2.3266e+00
dt= 0.200   |error| at t=2.0 = 1.1973e+00
dt= 0.100   |error| at t=2.0 = 6.6156e-01
dt= 0.050   |error| at t=2.0 = 3.4907e-01
dt= 0.020   |error| at t=2.0 = 1.4441e-01
dt= 0.010   |error| at t=2.0 = 7.3038e-02
dt= 0.005   |error| at t=2.0 = 3.6731e-02

Show plotting code

fig, ax = plt.subplots(figsize=(6, 5))
ax.loglog(dts, errors, "o-", color="C0", label="Euler error at $t=2$")
ax.loglog(dts, errors[0] * (dts / dts[0]), "k--", alpha=0.6, label="slope $=1$ reference")
ax.set_xlabel("step size $\\Delta t$"); ax.set_ylabel("|error|")
ax.set_title("Euler is first-order: error $\\propto \\Delta t$")
ax.legend()
ax.grid(alpha=0.3, which="both")
plt.show()

The data points lie on the reference line with slope

1

, as expected for a first-order accurate method. Higher-order methods such as Heun (slope

2

) and RK4 (slope

4

) give more accuracy per step. Like Euler, they advance the solution using the derivative. The Euler step is sufficient for the examples that follow. Better integrators solve the same ODE with greater accuracy.

Rotation in two dimensions

An ODE can also describe a vector state

\mathbf{x} \in \mathbb{R}^d

. Its right-hand side is then a vector-valued function

\mathbf{v}(\mathbf{x}, t) \in \mathbb{R}^d

\frac{d \mathbf{x}}{d t} = \mathbf{v}(\mathbf{x}, t), \qquad \mathbf{x}(0) = \mathbf{x}_0.

The vector field

\mathbf{v}

assigns a velocity to every point in space. A solution follows this velocity. The Euler update has the same form:

\mathbf{x}_{n+1} = \mathbf{x}_n + \Delta t \cdot \mathbf{v}(\mathbf{x}_n, t_n).

Consider the two-dimensional vector field for pure rotation:

\mathbf{v}(x, y) = \begin{bmatrix} -y \\ x \end{bmatrix}.

Each particle at position

(x, y)

has a velocity perpendicular to its position vector. Its exact trajectory is a circle around the origin, traced counterclockwise.

def euler_vec(v, x0, t0, t_end, dt):
    """Euler method for vector ODE dx/dt = v(x, t). x is a 1-D array."""
    n_steps = int(np.round((t_end - t0) / dt))
    xs = np.empty((n_steps + 1, len(x0)))
    xs[0] = x0
    t = t0
    for n in range(n_steps):
        xs[n + 1] = xs[n] + dt * v(xs[n], t)
        t += dt
    return xs


# Pure rotation
v_rot = lambda x, t: np.array([-x[1], x[0]])

# Release a few particles, integrate one full revolution (2*pi)
starts = [(1.0, 0.0), (2.0, 0.0), (0.0, 1.5), (-1.0, -1.0)]
trajectories = [euler_vec(v_rot, np.array(s), t0=0.0, t_end=2 * np.pi, dt=0.02)
                for s in starts]

Show plotting code

fig, ax = plt.subplots(figsize=(6, 6))

# Quiver of the field
xx, yy = np.meshgrid(np.linspace(-3, 3, 18), np.linspace(-3, 3, 18))
ax.quiver(xx, yy, -yy, xx, color="lightgray", alpha=0.7, scale=60)

# Euler trajectories
for s, traj in zip(starts, trajectories):
    ax.plot(traj[:, 0], traj[:, 1], lw=1.5)
    ax.scatter(*s, s=40, zorder=3, edgecolor="white", linewidth=1.0)
    ax.scatter(traj[-1, 0], traj[-1, 1], marker="x", s=60, color="black", zorder=3)

ax.set_aspect("equal")
ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
ax.set_xlabel("x"); ax.set_ylabel("y")
ax.set_title("Euler trajectories under pure rotation $v(x,y)=(-y,\\,x)$")
ax.grid(alpha=0.3)
plt.show()

The black ✕ marks show the position of each particle after

t = 2\pi

. The Euler endpoints drift slightly outward because of the same first-order error seen in one dimension. Shrinking

\Delta t

moves the endpoints toward their starting circles. A higher-order solver reduces the drift faster than decreasing

\Delta t

with Euler.

Generative models

The same construction appears in modern generative models. Replace the hand-written rotation field with a learned vector field

\mathbf{v}_\theta(\mathbf{x}, t)

parameterized by a neural network. Initialize the particles with Gaussian noise. Then apply the same Euler loop.

If $\mathbf{v}_\theta$ is trained to match a target velocity field that connects Gaussian noise to the data distribution, the result is flow matching.
If $\mathbf{v}_\theta$ encodes the score of a noise-perturbed data distribution, the result is a diffusion model in its probability-flow ODE form.

Both models use an ODE with a learned vector field. The vector field

\mathbf{v}_\theta

, the probability path, and the regression target used to train the network still have to be chosen.

References

L. C. Evans, Partial Differential Equations (Chapter 1), for ODE basics and existence/uniqueness.
E. Hairer, S. P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I, the standard reference on Euler and its successors.
Khan Academy, Differential equations, for an introductory video course.

Edit this page on GitHub or file an issue.

Prerequisites

Probability

Linear Algebra

Calculus

Programming

ODEs and the Euler method

Exponential growth

The Euler method

Step size and error

Rotation in two dimensions

Generative models

References

​Exponential growth

​The Euler method

​Step size and error

​Rotation in two dimensions

​Generative models

​References

Exponential growth

The Euler method

Step size and error

Rotation in two dimensions

Generative models

References