Appendix A: Crash Course on Differential Equations

Differential equations (DEs) are fundamental tools for modeling dynamic systems and can be broadly categorized into ordinary differential equations (ODEs), stochastic differential equations (SDEs), and partial differential equations (PDEs).

ODEs describe how a system’s state changes over time according to a precise rule, so that knowing the starting point determines the future path exactly. SDEs add randomness to this evolution, modeling how noise or uncertainty influences the system’s behavior, making the outcome probabilistic rather than fixed. PDEs explain how functions depending on several variables, such as time and space, evolve together, capturing phenomena like heat spreading, waves moving, or the time evolution of probability densities in stochastic systems (Spoiler: Fokker-Planck equation). These types of differential equations form a fundamental language for understanding how systems evolve over time and space under both deterministic and random influences.

In this chapter, we provide essential prerequisites on differential equations.

Foundation of Ordinary Differential Equations

This section introduces the fundamental theory of ODEs, emphasizing the uniqueness of solutions given an initial condition. It also covers practical methods for solving ODEs using numerical solvers.

Intuition of Ordinary Differential Equation

The deterministic process is called an ordinary differential equation (ODE). In the multivariate case, we consider systems of the form: \[\begin{equation} \label{eq:general_ode} \frac{\mathop{}\!\mathrm{d}\mathbf{x}(t)}{\mathop{}\!\mathrm{d}t} = \mathbf{v}(\mathbf{x}(t), t), \end{equation}\] where \(\mathbf{x}(t) \in \mathbb{R}^D\) is a vector-valued function representing the state of the system at time \(t\), and \(\mathbf{v}: \mathbb{R}^D \times \mathbb{R} \to \mathbb{R}^D\) is a vector field specifying the direction and magnitude of change at each point in space and time.

High-Level Intuition for Solving ODEs.

To build intuition, imagine the vector field \(\mathbf{v}(\mathbf{x}, t)\) as a dynamic landscape of arrows that tells you how a point \(\mathbf{x}\) should move at any given time \(t\). Solving the differential equation means tracing out a curve \(\mathbf{x}(t)\) through this field such that the tangent (i.e., the instantaneous velocity) of the curve at any point aligns with the vector given by \(\mathbf{v}(\mathbf{x}(t), t)\).

• Vector Field Perspective: The function \(\mathbf{v}(\mathbf{x}, t)\) defines how things should move: it gives the local “instructions” for motion or change.

• Trajectory Perspective: The solution \(\mathbf{x}(t)\) is a path that a particle would follow if it obeys the rule set by the vector field \(\mathbf{v}\) at every instant.

Thus, solving an ODE is like placing a particle in a flow field and observing where it goes over time.

Existence and Uniqueness of Ordinary Differential Equations

So far, we have seen that solving an ODE means finding a path that follows the directions given by the vector field at every point. Intuitively, this is like tracing the trajectory of a particle as it moves along the flow defined by the velocities.

But this picture leads to an important question:

Answering these questions is essential because it tells us whether the system’s behavior can be reliably predicted from its starting position. The Existence and Uniqueness Theorem provides conditions on the vector field that guarantee exactly one path starting from any given initial point. This ensures the solution behaves consistently and forms a cornerstone of the theory of ODEs.

Local (in Time) Existence and Uniqueness Theorem.

Below, we state a local version of the theorem, which asserts existence and uniqueness of a solution in a neighborhood of the initial time for a given initial condition.

The essence of the proof is rooted in the Picard–Lindelöf iteration method, whose core idea is also leveraged in Section ((Optional) ParaDiGMs(Optional) ParaDiGMs) in Chapter 9 to accelerate the sampling process of diffusion models.

Global (in Time) Existence and Uniqueness Theorem.

While the Local Existence and Uniqueness Theorem guarantees the existence of solutions on a small time interval, the “global (in time) existence and uniqueness theorem” extends this result to the entire interval \([t_0, T]\) under additional regularity conditions. A well-known result in this category is the Carathéodory theorem, which ensures the global existence and uniqueness of solutions to ODEs under two key assumptions: local Lipschitz continuity in the state variable and a linear growth bound.

1. Local Lipschitz condition in \(\mathbf{x}\): There exists a function \(\mathrm{Lip}(t)\), integrable on \([0, T]\), such that for all \(\mathbf{x}_1, \mathbf{x}_2 \in \mathbb{R}^D\), \[\|{\mathbf{v}}(\mathbf{x}_1, t) - {\mathbf{v}}(\mathbf{x}_2, t)\| \leq \mathrm{Lip}(t) \|\mathbf{x}_1 - \mathbf{x}_2\|.\]

2. Linear growth condition: There exists a function \(M(t)\), integrable on \([0, T]\), such that for all \(\mathbf{x} \in \mathbb{R}^D\), \[\|{\mathbf{v}}(\mathbf{x}, t)\| \leq M(t)(1 + \|\mathbf{x}\|).\]

We refer the reader to (Reid 1971) for a comprehensive discussion of the assumptions, formal statement, and detailed proof of the theorem.

In summary, when an initial condition is given to an ODE defined by a time-dependent velocity field, the trajectory of the particle flow is uniquely determined.

Uniqueness Implies Non-Intersection of Solutions

The uniqueness of solutions in ODEs, as guaranteed by the Local Existence and Uniqueness Theorem, implies a fundamental property: two different solution trajectories, starting from different initial conditions, cannot cross each other. This reflects the deterministic nature of ODEs, ensuring that each state evolves along a unique path. The following corollary formalizes this result.

By guaranteeing non-intersecting solution paths, this theorem offers hidden yet crucial support for the flow map model (see Chapter 10 and Chapter 11).

Exponential Integration Factor

Even an ODE determined by a general time-varying velocity \({\mathbf{v}}\) does not admit a closed-form solution, in some special cases, we can solve them analytically or by reducing its formulation to a better structural one.

An Illustrative Example.

Consider the following linear scalar ODE: \[\frac{\mathop{}\!\mathrm{d}{\mathbf{x}}(t)}{\mathop{}\!\mathrm{d}t} = L(t){\mathbf{x}}(t),\] where \(L(t) \in \mathbb{R}\) is a continuous function. This equation is solvable in closed form, and its solution is well known (for any \(s\) and \(t\)): \[{\mathbf{x}}(t) = {\mathbf{x}}(s) \cdot \exp\left( \int_s^t L(\tau) \,\mathop{}\!\mathrm{d}\tau \right).\] This formula demonstrates how the solution evolves according to an exponential factor that accumulates the effect of the time-dependent coefficient \(L(t)\). This motivates the use of exponential integration factors: \[\begin{align} \label{eq:exp-int-factor} \mathcal{E}(s \rightarrow t) := \exp\left( \int_s^t L(\tau)\mathop{}\!\mathrm{d}\tau \right), \end{align}\] especially in more general settings where the dynamics include both linear and nonlinear components.

Semilinear ODEs and Exponential Integration Factors.

We now consider a broader class of ODEs known as semilinear ODEs. These equations separate the dynamics into a linear part (in the state variable) and a nonlinear remainder: \[\begin{align} \label{eq:semilinear-ode} \frac{\mathop{}\!\mathrm{d}{\mathbf{x}}(t)}{\mathop{}\!\mathrm{d}t} = L(t){\mathbf{x}}(t) + {\mathbf{N}}({\mathbf{x}}(t), t), \end{align}\] where \({\mathbf{x}}(t) \in \mathbb{R}^D\) is the state vector, \(L(t)\) is a scalar-valued continuous function, and \({\mathbf{N}}: \mathbb{R}^D \times [0,T] \to \mathbb{R}^D\) is a nonlinear vector field. This semilinear structure arises naturally in many physical and engineering systems. In particular, it also appears in the probability flow ODE formulation of diffusion models (see Equation in Chapter 4). Recognizing this structure enables integrating-factor representations that handle the known linear dynamics analytically and leave the remaining term to be approximated numerically. This principle plays a central role in the diffusion ODE solvers studied in Chapter 9.

Step 1: Isolate the Non-Linear Term via an Integration Factor.

Observing that we can isolate the nonlinear part by subtracting the linear drift from the semiliner ODE in Equation \(\eqref{eq:semilinear-ode}\): \[\frac{\mathop{}\!\mathrm{d}{\mathbf{x}}(t)}{\mathop{}\!\mathrm{d}t} - L(t){\mathbf{x}}(t) = {\mathbf{N}}({\mathbf{x}}(t), t).\] To absorb the linear term, we multiply both sides by the inverse integration factor: \[\mathcal{E}^{-1}(s \rightarrow t) = \exp\left(-\int_s^t L(\tau)\mathop{}\!\mathrm{d}\tau \right).\] Now apply the product rule to the left-hand side: \[\begin{align*} \mathcal{E}^{-1}(s \rightarrow t) \left( \frac{\mathop{}\!\mathrm{d}{\mathbf{x}}(t)}{\mathop{}\!\mathrm{d}t} - L(t){\mathbf{x}}(t) \right) &= \frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t} \left[ \mathcal{E}^{-1}(s \rightarrow t) {\mathbf{x}}(t) \right]. \end{align*}\] Hence, the equation becomes: \[\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t} \left[ \mathcal{E}^{-1}(s \rightarrow t) {\mathbf{x}}(t) \right] = \mathcal{E}^{-1}(s \rightarrow t) {\mathbf{N}}({\mathbf{x}}(t), t).\]

This transformation simplifies the original equation by isolating the nonlinear component, allowing us to focus entirely on the nonlinear dynamics in a transformed coordinate system.

Step 2: Integrate Over Time.

We now integrate both sides from \(s\) to \(t\): \[\int_s^t \frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}\tau} \left[ \mathcal{E}^{-1}(s \rightarrow\tau) {\mathbf{x}}(\tau) \right] \mathop{}\!\mathrm{d}\tau = \int_s^t \mathcal{E}^{-1}(s \rightarrow\tau) {\mathbf{N}}({\mathbf{x}}(\tau), \tau)\mathop{}\!\mathrm{d}\tau.\] The left-hand side is simply the difference of the transformed variable evaluated at \(t\) and \(s\): \[\mathcal{E}^{-1}(s \rightarrow t) {\mathbf{x}}(t) - {\mathbf{x}}(s).\] Hence, we obtain: \[\mathcal{E}^{-1}(s \rightarrow t) {\mathbf{x}}(t) = {\mathbf{x}}(s) + \int_s^t \mathcal{E}^{-1}(s \rightarrow\tau) {\mathbf{N}}({\mathbf{x}}(\tau), \tau)\mathop{}\!\mathrm{d}\tau.\]

Step 3: Solve for \({\mathbf{x}}(t)\).

Multiplying both sides by the exponential flow \(\mathcal{E}(s \rightarrow t)\) gives the solution: \[\begin{align} \label{eq:exp-int-factor-integration} {\mathbf{x}}(t) = \underbrace{\mathcal{E}(s \rightarrow t) {\mathbf{x}}(s)}_{\text{linear part}} + \underbrace{\int_s^t \mathcal{E}(\tau \rightarrow t) {\mathbf{N}}({\mathbf{x}}(\tau), \tau)\mathop{}\!\mathrm{d}\tau}_{\text{nonlinear part}} . \end{align}\]

The solution naturally separates into a linear and a nonlinear component. Exponential integrators exploit this structure by solving the linear part in exactly closed form and discretizing only the nonlinear residual. This ensures that the step size is dictated by the nonlinear dynamics rather than by the potentially large linear coefficient (see, for instance, the comparison between the exponential Euler update and the vanilla Euler update in Section (Prologue) in Chapter 9).

Numerical Solvers of Ordinary Differential Equations

We consider the ODE in Equation \(\eqref{eq:general_ode}\), \[\frac{\mathop{}\!\mathrm{d}\mathbf{x}(t)}{\mathop{}\!\mathrm{d}t} = \mathbf{v}(\mathbf{x}(t),t),\] together with an initial condition. Solving the ODE means finding the continuous trajectory \(\mathbf{x}(t)\) whose instantaneous change agrees with the velocity field \(\mathbf{v}\) at every time.

Equivalently, the solution satisfies the integral equation \[\begin{equation} \label{eq:integral_form} \mathbf{x}(t) = \mathbf{x}(0) + \int_0^t \mathbf{v}(\mathbf{x}(\tau),\tau)\mathop{}\!\mathrm{d}\tau. \end{equation}\] If this integral cannot be evaluated in closed form, we approximate the trajectory numerically by evaluating the velocity field at finitely many points.

From a Continuous Trajectory to Numerical Updates.

A numerical solver first chooses a discrete time grid \[t_0,t_1,\ldots,t_M\] and produces corresponding approximations \[\mathbf{x}_0,\mathbf{x}_1,\ldots,\mathbf{x}_M\] to the continuous trajectory.

The gap \[h_i:=t_{i+1}-t_i\] is the step size. Smaller steps generally provide more local information about the trajectory but require more function evaluations. The purpose of a numerical solver is therefore to use the available evaluations of \(\mathbf{v}\) efficiently so that accurate trajectories can be obtained with relatively few steps.

Different solvers mainly differ in where the information used for the next update comes from. This leads to the following useful organization.

A Structural View of Numerical Solvers.

For the methods relevant to this book, it is useful to distinguish two broad numerical viewpoints.

• Time-Stepping Methods: Starting from the current state, the solver advances the trajectory from one time point to the next. Euler, Runge–Kutta, and Adams–Bashforth methods are representative examples.

  Within time-stepping methods, we further distinguish:

  • a single-step method, which constructs the next state from the current state and any additional evaluations generated within the current interval; and

  • a multistep method, which additionally reuses information from previous completed steps.

• Parallel-in-Time Viewpoint: Instead of completing all time points strictly one after another, one may formulate the trajectory computation so that tentative states at several times can be updated together. Fixed-point constructions such as Picard iteration provide a basic mechanism for this idea.

A single-step method should not be confused with a single function evaluation. For example, Euler uses one velocity evaluation per step, whereas a Runge–Kutta method may evaluate the velocity several times at intermediate stages while still being a single-step method. The defining distinction is whether information from earlier completed steps is required.

The distinction is especially useful later in diffusion sampling: Runge–Kutta-type constructions obtain richer information by looking within the current step, whereas Adams–Bashforth-type constructions obtain it by looking backward at previously computed values.

Solving ODEs in Forward and Reverse Time.

The taxonomy above does not depend on the direction in which the ODE is integrated. In ordinary forward integration, \[t_{n+1}>t_n, \qquad h_n=t_{n+1}-t_n>0.\] Diffusion sampling instead typically starts from a noisy terminal state and uses a decreasing time grid, \[t_{n+1}<t_n, \qquad h_n=t_{n+1}-t_n<0.\]

The same numerical formulas apply in either direction when the step size is treated as a signed quantity. Thus, reverse-time integration of an ODE does not require a new family of numerical solvers; Euler, Runge–Kutta, multistep, and related methods can all be applied on a decreasing time grid.

This simple reversal should be distinguished from reversing an SDE, where the stochastic dynamics require additional analysis. We discuss that issue in the following section.

Representative Numerical Solvers

The taxonomy above tells us where a solver obtains information for its next update. We now make these ideas concrete through a few representative examples.

We deliberately keep the list short. Euler illustrates the simplest single-step update; Runge–Kutta methods illustrate the use of additional stages within the current step; Adams–Bashforth illustrates reuse of past history; and Picard iteration illustrates the fixed-point viewpoint behind parallel-in-time methods.

Throughout the examples below, let \[t_{n+1}=t_n+h.\] The step size \(h\) is signed: \(h>0\) for forward integration and \(h<0\) for reverse integration.

Euler’s Method.

Euler’s method is the simplest single-step method. It uses only the velocity at the beginning of the current interval: \[\mathbf{x}_{n+1} = \mathbf{x}_n + h\,\mathbf{v}(\mathbf{x}_n,t_n).\] Geometrically, Euler looks at the current direction of motion and follows that direction for one complete step.

Under standard smoothness assumptions, Euler has local truncation error \(\mathcal O(h^2)\) and global error \(\mathcal O(h)\). Its simplicity makes it a useful first-order baseline, but its accuracy deteriorates when the velocity changes substantially over one step.

Runge–Kutta Methods.

Euler uses only the velocity at the beginning of the interval. Runge–Kutta (RK) methods obtain more information by evaluating the velocity at additional stages generated within the same step.

They are therefore single-step, multi-stage methods: all information used to construct \(\mathbf{x}_{n+1}\) is generated from the current starting point \(\mathbf{x}_n\), but one step may contain several function evaluations.

We illustrate the basic RK principle using two second-order methods.

Explicit Midpoint Method.

The explicit midpoint method first uses the current velocity to predict the state halfway through the interval: \[\mathbf{x}_{\mathrm{mid}} = \mathbf{x}_n + \frac{h}{2} \mathbf{v}(\mathbf{x}_n,t_n).\] It then evaluates the velocity at this predicted midpoint and uses that direction for the complete update: \[\mathbf{x}_{n+1} = \mathbf{x}_n + h\, \mathbf{v} \left( \mathbf{x}_{\mathrm{mid}}, t_n+\frac{h}{2} \right).\]

Thus, instead of assuming that the starting velocity remains representative over the entire interval, the midpoint method “looks inside” the step before deciding how to complete it.

Heun’s Method (Improved Euler).

Heun’s method is another two-stage, second-order Runge–Kutta method. Instead of looking at the midpoint, it first uses Euler’s method to predict the endpoint: \[\mathbf{x}_{\mathrm{pred}} = \mathbf{x}_n + h\,\mathbf{v}(\mathbf{x}_n,t_n).\] It then evaluates the velocity at this predicted endpoint and averages the starting and predicted ending velocities: \[\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{h}{2} \left[ \mathbf{v}(\mathbf{x}_n,t_n) + \mathbf{v}(\mathbf{x}_{\mathrm{pred}},t_n+h) \right].\]

The explicit midpoint and Heun methods both use two stages and, under standard smoothness assumptions, have local truncation error \(\mathcal O(h^3)\) and global error \(\mathcal O(h^2)\). They differ only in how the additional within-step information is collected and combined.

Higher-order Runge–Kutta methods extend the same principle by using more stages; their standard formulas are omitted here. In diffusion sampling, Karras et al. (2022) use Heun’s method directly as an ODE solver.

Adams–Bashforth Method.

Runge–Kutta methods obtain additional information by looking within the current step. A multistep method instead looks backward and reuses information from previous completed steps.

For example, with a uniform step size \(h\), the second-order Adams–Bashforth method (AB2) is \[\mathbf{x}_{n+1} = \mathbf{x}_n + h \left[ \frac{3}{2} \mathbf{v}(\mathbf{x}_n,t_n) - \frac{1}{2} \mathbf{v}(\mathbf{x}_{n-1},t_{n-1}) \right].\] The second piece of information is therefore a stored velocity from the previous step rather than a new evaluation inside the current interval.

Under standard smoothness assumptions and sufficiently accurate initialization, AB2 has local truncation error \(\mathcal O(h^3)\) and global error \(\mathcal O(h^2)\). Because the first step has no previous history, a multistep method must first be initialized using another solver.

The formula above assumes a uniform grid. On nonuniform grids, the multistep coefficients depend on the neighboring step sizes. This history-reuse principle will later reappear in the diffusion-specific multistep solvers in Section (DEISDEIS) in Chapter 9 and Section (DPM-Solver) in Chapter 9++-multi.

Picard Iteration.

The methods above approximate the trajectory by constructing local time-stepping updates. Picard iteration follows a different idea: it repeatedly refines an entire candidate trajectory through the integral form of the ODE.

Starting from an initial guess \(\mathbf{x}^{(0)}(t)\) satisfying \(\mathbf{x}^{(0)}(0)=\mathbf{x}(0)\), Picard iteration defines \[\mathbf{x}^{(k+1)}(t) = \mathbf{x}(0) + \int_0^t \mathbf{v} \bigl( \mathbf{x}^{(k)}(s),s \bigr) \mathop{}\!\mathrm{d}s.\] Each iteration therefore evaluates the velocity along the previous candidate trajectory and uses those evaluations to construct a new one.

A fixed point of this iteration satisfies the original ODE integral equation. Under the standard conditions of the Picard–Lindelöf theorem, the iteration converges locally to the ODE solution.

For our later purposes, the important observation is computational: within one Picard iteration, the velocity at each time depends on the previous candidate trajectory. Evaluations at different time points can therefore be carried out independently, which opens the possibility of parallel-in-time computation. This viewpoint will be revisited in Section ((Optional) ParaDiGMs(Optional) ParaDiGMs) in Chapter 9.

Foundation of Stochastic Differential Equations

Stochastic Differential Equations (SDEs) are an extension of ordinary differential equations (ODEs) that incorporate randomness, providing a mathematical framework for modeling systems affected by uncertainty. This chapter introduces SDEs, beginning with the discretization of ODEs, extending to the discretization of SDEs, and culminating in a discussion of general SDEs, including Ito’s calculus and Ito’s formula.

From ODEs to SDEs: An Intuitive Introduction

Let us begin with a ODE describing the deterministic evolution of a state variable \({\mathbf{x}}(t) \in \mathbb{R}^D\): \[\begin{equation} \frac{\mathop{}\!\mathrm{d}{\mathbf{x}}(t)}{\mathop{}\!\mathrm{d}t} = {\mathbf{f}}({\mathbf{x}}(t), t), \quad {\mathbf{x}}(0) = {\mathbf{x}}_0. \label{eq:ode} \end{equation}\] Here, \({\mathbf{f}}: \mathbb{R}^D \times [0,T] \rightarrow \mathbb{R}^D\) is a time-dependent velocity field that governs the dynamics of \({\mathbf{x}}(t)\). The solution to this ODE is a smooth trajectory \(t \mapsto {\mathbf{x}}(t)\), fully determined by the initial condition \({\mathbf{x}}_0\).

Discretization Perspective.

To build intuition, consider an Euler discretization of Equation \(\eqref{eq:ode}\) over small time steps \(\Delta t\): \[\begin{equation*} {\mathbf{x}}_{t+\Delta t} = {\mathbf{x}}_t + {\mathbf{f}}({\mathbf{x}}_t, t) \Delta t. \end{equation*}\] This approximation becomes more accurate as \(\Delta t \to 0\), converging (under standard regularity conditions on \({\mathbf{f}}\)) to the exact solution of the ODE.

Introducing Randomness: From ODE to SDE.

In many real-world systems, perfect knowledge of the dynamics is unrealistic. Noise, uncertainty, or unmodeled interactions may affect the evolution. To incorporate such randomness, we augment the ODE with a stochastic term: \[\begin{equation} {\mathbf{x}}_{t+\Delta t} = {\mathbf{x}}_t + {\mathbf{f}}({\mathbf{x}}_t, t)\Delta t + g(t)\sqrt{\Delta t} \cdot \bm{\epsilon}_t, \label{eq:euler-sde} \end{equation}\] where

• \(g: [0,T] \to \mathbb{R}\) is a diffusion coefficient (possibly dependent on both state and time, though here assumed time-dependent only),

• \(\bm{\epsilon}_t \sim \mathcal{N}(\bm{0}, {\mathbf{I}}_D)\) are i.i.d. standard Gaussian vectors.

This modified update rule reflects not just deterministic drift, but also random perturbations scaled by \(\sqrt{\Delta t}\). The scaling ensures that the stochastic perturbation remains finite in the limit \(\Delta t \to 0\). Importantly, this formulation gives rise to a continuous-time stochastic process as \(\Delta t \to 0\), which leads us to the framework of SDE.

Stochastic Differential Equations.

Formally, the limit of the discrete update Equation \(\eqref{eq:euler-sde}\) as \(\Delta t \to 0\) defines the SDE: \[\begin{equation} \mathop{}\!\mathrm{d}{\mathbf{x}}(t) = {\mathbf{f}}({\mathbf{x}}(t), t)\mathop{}\!\mathrm{d}t + g(t) \mathop{}\!\mathrm{d}{\mathbf{w}}(t). \label{eq:sde} \end{equation}\] Here, \({\mathbf{w}}(t) \in \mathbb{R}^D\) is a Wiener process (standard Brownian motion), a continuous-time stochastic process characterized by:

• Initial State: \({\mathbf{w}}(0) = 0\) almost surely;

• Independent Increments: for \(0 \leq s < t\), the increment \({\mathbf{w}}(t) - {\mathbf{w}}(s)\) is independent of the past;

• Gaussian Increments: \[\begin{align} \label{eq:wiener-gaussian} {\mathbf{w}}(t) - {\mathbf{w}}(s) \sim \mathcal{N}\bigl(\bm{0}, (t - s){\mathbf{I}}_D \bigr) \end{align}\]

• Continuity: Sample paths \(t \mapsto {\mathbf{w}}(t)\) are almost surely continuous but nowhere differentiable.

In addition, the notation \[\mathop{}\!\mathrm{d}{\mathbf{w}}(t) := {\mathbf{w}}(t + \mathop{}\!\mathrm{d}t) - {\mathbf{w}}(t)\] is often used to denote the infinitesimal increment of the Wiener process.

While suggestive, this notation is heuristic and should not be interpreted as a classical differential (e.g., in the Riemann or Lebesgue sense), since Brownian paths are almost surely nowhere differentiable. Instead, it serves as a formal shorthand to express the Gaussian increments property: \[\mathop{}\!\mathrm{d}{\mathbf{w}}(t) \sim \mathcal{N}(0, \mathop{}\!\mathrm{d}t\, \mathrm{I}_D),\] meaning that over an infinitesimal time interval of length \(\mathop{}\!\mathrm{d}t\), the increment of the Wiener process behaves like a Gaussian random variable with zero mean and covariance \(\mathop{}\!\mathrm{d}t\, \mathrm{I}_D\).

Further Explanation of Equation \(\eqref{eq:sde}\)

The SDE in Equation \(\eqref{eq:sde}\) should be understood in its integral form: \[\begin{align} \label{eq:sde-integral-form} {\mathbf{x}}(t) = {\mathbf{x}}(0) + \int_0^t {\mathbf{f}}({\mathbf{x}}(s), s) \mathop{}\!\mathrm{d}s + \int_0^t g(s) \mathop{}\!\mathrm{d}{\mathbf{w}}(s), \end{align}\] interpreted in the Itô sense. Here, the first term is a classical (Riemann or Lebesgue) integral representing the accumulated deterministic drift, while the second term is an Itô stochastic integral, which integrates with respect to the Wiener process \({\mathbf{w}}(t)\). We do not provide a full rigorous construction of the Itô integral, but offer the following intuition.

Intuition for Itô Integration.

The Itô integral can be understood as the limit of discrete sums of the form¹ \[\sum_i g(t_i)\bigl(w(t_{i+1})-w(t_i)\bigr),\] where the integrand \(g(t)\) is evaluated at the left endpoint \(t_i\) of each subinterval. As the partition becomes finer, these sums converge to the Itô integral. The left-point evaluation rule is what distinguishes Itô integration from classical integrals, which may use midpoints or other evaluation rules. Because Brownian paths are continuous yet almost surely nowhere differentiable, classical integration does not directly apply. The Itô integral handles this irregularity and captures the cumulative effect of stochastic fluctuations over time.

Use of Differential Notation.

Expressions such as \(\mathop{}\!\mathrm{d}{\mathbf{x}}(t)\), \(\mathop{}\!\mathrm{d}t\), and \(\mathop{}\!\mathrm{d}{\mathbf{w}}(t)\) are not classical differentials. Instead, they are formal notations representing infinitesimal increments of the respective processes. While heuristic, they are widely used for their convenience in expressing SDEs analogously to ODEs and facilitate formal manipulations within Itô calculus.

How Itô calculus is applied in diffusion models will be explained in Appendix C.

Comparison with ODEs.

In ODEs, e.g., \[\frac{\mathop{}\!\mathrm{d}{\mathbf{x}}(t)}{\mathop{}\!\mathrm{d}t} = {\mathbf{f}}({\mathbf{x}}(t), t),\] the integral form \[{\mathbf{x}}(t) = {\mathbf{x}}(0) + \int_0^t {\mathbf{f}}({\mathbf{x}}(\tau), \tau) \mathop{}\!\mathrm{d}\tau\] is justified by the Fundamental Theorem of Calculus, which ensures that differentiable functions can be recovered from their derivatives.

By contrast, in SDEs such as Equation \(\eqref{eq:sde}\), there is no direct analog of this theorem because Brownian motion lacks differentiability, and stochastic integrals do not follow the classical chain rule. Instead, Itô calculus introduces alternative tools (e.g., Itô’s lemma) to analyze and manipulate stochastic dynamics.

Thus, while the differential notation for SDEs is compact and intuitive, a rigorous understanding depends on interpreting them via their integral formulation using Itô integrals.

A Numerical Solver for SDE.

Like ODEs, the SDE in Equation \(\eqref{eq:sde}\) admits a unique solution² if \(\mathbf{f}(\cdot, t)\) and \(g(\cdot)\) satisfy some smoothness conditions: \(\mathbf{f}(\cdot, t)\) is Lipschitz and of linear growth in \(\mathbf{x}\), and \(g(\cdot)\) is square integrable.

For general SDEs as in Equation \(\eqref{eq:sde}\), closed-form solutions are generally unavailable, so numerical methods are necessary. A common approach is the Euler–Maruyama method, which generalizes Euler’s method for ODEs and, indeed, we have already seen it in Equation \(\eqref{eq:euler-sde}\). It approximates the drift term \({\mathbf{f}}({\mathbf{x}}(t), t)\) over a time step \(\Delta t\) and simulates the stochastic noise \(g(t) \mathop{}\!\mathrm{d}{\mathbf{w}}(t)\) using Gaussian increments \(\sqrt{\Delta t} \, \boldsymbol{\epsilon}_t\) with \(\boldsymbol{\epsilon}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{I})\).

Later, in Section (Itô’s Formula Application: Closed-Form Solution of a Linear SDE) in Appendix C, we will see that a linear SDE admits a closed-form solution.