Bottle cap race

Imagine flicking a bottle cap (a chapa) across a table. Its movement isn't perfectly smooth. It's driven forward by your push, slowed down by friction, and knocked off course by tiny bumps and imperfections on the surface. This simulation captures that complex motion using a stochastic ordinary differential equation.

The model describes the chapa's state at any instant through its position ( $X_{t}$ ) and velocity ( $V_{t}$ ). It's defined by two simple-looking but profound equations:

\begin{aligned} d X_{t} & = V_{t} d t \\ d V_{t} & = (- γ V_{t} + μ) d t + σ d W_{t} \end{aligned}

1. Stratonovich Form and Analytical Solution

An SDE in Itô form, $d Y_{t} = a (Y_{t}) d t + b (Y_{t}) d W_{t}$ , can be converted to the Stratonovich form, $d Y_{t} = a^{*} (Y_{t}) d t + b (Y_{t}) \circ d W_{t}$ , using the correction term:

a^{*} (Y_{t}) = a (Y_{t}) - \frac{1}{2} b (Y_{t}) \frac{\partial b}{\partial y} (Y_{t})

See stochastic ordinary differential equation#Stratonovich form.
In the SDE for velocity, $V_{t}$ :

The drift is $a (V_{t}) = - γ V_{t} + μ$ .
The diffusion is $b (V_{t}) = σ$ .
Since the diffusion term $σ$ is a constant, its derivative with respect to $V_{t}$ is zero:

\frac{\partial b}{\partial V_{t}} = \frac{\partial σ}{\partial V_{t}} = 0

Therefore, the correction term is zero, and the drift terms are identical, $a^{*} (V_{t}) = a (V_{t})$ .
This means the Itô and Stratonovich forms of this SDE are the same. The system in the Stratonovich sense is:

\begin{aligned} d X_{t} & = V_{t} d t, \\ d V_{t} & = (- γ V_{t} + μ) d t + σ \circ d W_{t} . \end{aligned}

We can solve this system by first finding the solution for $V_{t}$ and then integrating to find $X_{t}$ .

1. Solving for Velocity $V_{t}$ :
The SDE for $V_{t}$ is a linear equation known as the Ornstein-Uhlenbeck process with a non-zero mean. We can solve it using an integrating factor, similar to a linear ordinary differential equation (see this). Let the integrating factor be $e^{γ t}$ .
Starting with $d V_{t} + γ V_{t} d t = μ d t + σ d W_{t}$ , we multiply by $e^{γ t}$ :

e^{γ t} d V_{t} + γ e^{γ t} V_{t} d t = μ e^{γ t} d t + σ e^{γ t} d W_{t}

The left side is the differential $d (e^{γ t} V_{t})$ according to Itô's product rule. Integrating from $0$ to $t$ :

\int_{0}^{t} d (e^{γ s} V_{s}) = \int_{0}^{t} μ e^{γ s} d s + \int_{0}^{t} σ e^{γ s} d W_{s}

e^{γ t} V_{t} - V_{0} = \frac{μ}{γ} (e^{γ t} - 1) + σ \int_{0}^{t} e^{γ s} d W_{s}

Isolating $V_{t}$ gives the final solution:

V_{t} = V_{0} e^{- γ t} + \frac{μ}{γ} (1 - e^{- γ t}) + σ \int_{0}^{t} e^{γ (s - t)} d W_{s}

2. Solving for Position $X_{t}$ :
The position is simply the integral of the velocity:

X_{t} = X_{0} + \int_{0}^{t} V_{s} d s

Substituting the expression for $V_{s}$ and integrating each term yields:

X_{t} = X_{0} + \frac{V_{0}}{γ} (1 - e^{- γ t}) + \frac{μ}{γ} t - \frac{μ}{γ^{2}} (1 - e^{- γ t}) + \frac{σ}{γ} \int_{0}^{t} (1 - e^{γ (s - t)}) d W_{s}

This is the exact analytical solution for the position $X_{t}$ .

2. Numerical Simulation Scheme

The Euler-Maruyama scheme is a simple and effective method for simulating this system. We discretize time into small steps of size $Δ t$ .

Let $t_{n} = n Δ t$ . The discrete updates for position $X_{n} \approx X_{t_{n}}$ and velocity $V_{n} \approx V_{t_{n}}$ are as follows:

Initialize: Start with the initial conditions $X_{0}$ and $V_{0}$ at time $t_{0} = 0$ . Choose a small time step $Δ t$ .
Iterate: For each step $n = 0, 1, 2, \dots$ :
- Generate a random number $Z_{n}$ from a standard normal distribution, $N (0, 1)$ .
- Update the velocity using the SDE's discretized form: $V_{n + 1} = V_{n} + (μ - γ V_{n}) Δ t + σ \sqrt{Δ t} Z_{n}$
- Update the position using the last step's velocity (this is the standard explicit Euler method): $X_{n + 1} = X_{n} + V_{n} Δ t$
Repeat: Continue the iteration until the desired final time is reached.

This scheme provides a sequence of points $(X_{n}, V_{n})$ that approximates a sample path of the solution.