Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Minor enhancements in Docs/Stochastic Forcing section #184

Merged
merged 6 commits into from
Jan 15, 2021
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 36 additions & 32 deletions docs/src/stochastic_forcing.md
Original file line number Diff line number Diff line change
Expand Up @@ -41,13 +41,21 @@ In a similar manner, a stochastic differential equation
\mathrm{d} x = f(x) \, \mathrm{d} t + g(x) \, \mathrm{d} W_t , \quad x(t_0) = 0 ,
```

with ``\mathrm{d} W_t`` a white-noise process, can be written in an integral form as:
with ``W_t`` a brownian motion or Wiener process.

!!! tip "Wiener process"
A Wiener process is a random variable ``W_t`` that depends continuously on ``t \ge 0`` and satisfies the following properties:
1. Independence. For ``0 \le s \le t`` the increment ``W_t - W_s`` is independent of any prior values, i.e., independent of all ``W_\tau``, ``\tau \le s``.
2. Stationarity. The statistical distribution of the increment ``W_{t+s} − W_s`` does not depend on ``s`` (and so is identical in distribution to ``W_t``).
3. Gaussianity. ``W_t`` is a Gaussian process with mean ``\langle W_t \rangle = 0`` and covariance ``\langle W_t W_s \rangle = \min(t, s)``.

The stochastic differential equation above can be written in an integral form as:

```math
x(t) = \int_{t_0}^{t} f(x(s)) \, \mathrm{d} s + \int_{t_0}^{t} g(x(s)) \, \mathrm{d} W_s .
```

Of course now, the last integral is a stochastic integral and there is not a single
Of course now, the last integral above is a stochastic integral and there is not a single
straight-forward way of computing it -- there are a lot of different ways we can
approximate it as a Riemannian sum and each of them leads to a different answer.
The two most popular ways for computing such stochastic integrals are:
Expand Down Expand Up @@ -120,9 +128,9 @@ The above is demonstrated by evaluating the simple stochastic integral:
```math
\begin{aligned}
{\color{Green} \text{Itô}}&: {\color{Green} \left \langle \int_{t_0}^{t} W_s \, \mathrm{d} W_s \right \rangle \approx \sum_{j} \left \langle W_j (W_{j+1} - W_j) \right \rangle} \\
&\hspace{7.3em} {\color{Green} = \sum_j \left \langle W_j W_{j+1} \right \rangle - \left \langle W_j W_j \right \rangle \sim \sum_{j} t_j - t_j = 0} , \\
& \hspace{7.3em} {\color{Green} = \sum_j \left \langle W_j W_{j+1} \right \rangle - \left \langle W_j W_j \right \rangle \sim \sum_{j} t_j - t_j = 0} , \\
{\color{Magenta}\text{Stratonovich}}&: {\color{Magenta}\left \langle \int_{t_0}^{t} W_s \circ \mathrm{d} W_s \right \rangle \approx \sum_j \left \langle \frac1{2}(W_j + W_{j+1}) (W_{j+1} - W_j)\right \rangle} \\
&\hspace{7.3em} {\color{Magenta} = \frac1{2} \sum_j \left \langle W_{j+1} W_{j+1} \right \rangle - \left \langle W_j W_j \right \rangle \sim \frac1{2} \sum_j t_{j+1} - t_j = \frac{t}{2}} .
& \hspace{7.3em} {\color{Magenta} = \frac1{2} \sum_j \left \langle W_{j+1} W_{j+1} \right \rangle - \left \langle W_j W_j \right \rangle \sim \frac1{2} \sum_j t_{j+1} - t_j = \frac{t}{2}} .
\end{aligned}
```

Expand Down Expand Up @@ -262,7 +270,7 @@ using Statistics: mean
using Random: randn, seed!
seed!(1234) # for reproducing the same plots

μ = 0.2 # drag
μ = 0.2
σ = 0.2 # noise strength
dt = 0.01 # timestep
nsteps = 2001 # total timesteps
Expand All @@ -278,29 +286,29 @@ E_ito = zeros(size(ΔW))
E_str = zeros(size(ΔW))
E_numerical = zeros(size(ΔW))

for j = 2:nsteps # time step the equation
for j = 2:nsteps # time step the equations

# time-step dxₜ = -μ xₜ + √σ dWₜ
@. x[j, :] = x[j-1, :] + -μ * x[j-1, :] * dt + sqrt(σ) * ΔW[j-1, :]
# time-step dx = - μ x dt + √σ dW
@. x[j, :] = x[j-1, :] - μ * x[j-1, :] * dt + sqrt(σ) * ΔW[j-1, :]

# time-step dEₜ = (-2μ Eₜ + ½σ) dt + √σ xₜ dWₜ
# time-step dE = (- E + ½σ) dt + √σ x dW
@. E_ito[j, :] = E_ito[j-1, :] + (-2μ * E_ito[j-1, :]
+ σ/2) * dt + sqrt(σ) * x[j-1, :] * ΔW[j-1, :]

# time-step dEₜ = -2μ Eₜ dt + √σ xₜdWₜ
xbar = @. x[j-1, :] - μ * x[j-1, :] * dt + sqrt(σ) * ΔW[j-1, :]
Ebar = @. E_str[j-1, :] + (-2μ * E_str[j-1, :]) * dt + sqrt(σ) * x[j-1, :] * ΔW[j-1, :]
# time-step dE = - E dt + √σ xdW
xbar = @. x[j-1, :] - μ * x[j-1, :] * dt + sqrt(σ) * ΔW[j-1, :]
Ebar = @. E_str[j-1, :] + (-2μ * E_str[j-1, :]) * dt + sqrt(σ) * x[j-1, :] * ΔW[j-1, :]
@. E_str[j, :] = E_str[j-1, :] + (-2μ * (E_str[j-1, :]
+ Ebar) / 2) * dt + sqrt(σ) * (x[j-1, :] + xbar) / 2 * ΔW[j-1, :]
end

# direct computation of E from xₜ
# direct computation of E from x
@. E_numerical = 0.5 * x^2

# compare the three E(t) solutions
plot(μ * t, [E_numerical[:, 1] E_ito[:, 1] E_str[:, 1]],
linewidth = [3 2 1],
label = ["½ xₜ²" "Eₜ (Ito)" "Eₜ (Stratonovich) "],
label = ["½ xₜ²" "Eₜ (Ito)" "Eₜ (Stratonovich)"],
linestyle = [:solid :dash :dashdot],
xlabel = "μ t",
ylabel = "E",
Expand Down Expand Up @@ -409,7 +417,7 @@ interpretations coincide.
The forcing ``\xi`` obeys:

```math
\langle \xi(\bm{x},t) \rangle = 0 \quad \text{and} \quad \langle \xi(\bm{x}, t) \xi(\bm{x}', t') \rangle= Q(\bm{x} - \bm{x}') \delta(t-t'),
\langle \xi(\bm{x}, t) \rangle = 0 \quad \text{and} \quad \langle \xi(\bm{x}, t) \xi(\bm{x}', t') \rangle = Q(\bm{x} - \bm{x}') \delta(t - t') ,
```

that is, the forcing is white in time but spatially correlated; its spatial correlation is
Expand All @@ -421,15 +429,15 @@ that is stochastically forced while dissipated by linear drag ``\mu``. The energ
fluid is:

```math
E = \tfrac{1}{2} \overline{|\bm{\nabla}\psi|^2}^{x,y} = -\tfrac{1}{2} \overline{\psi \nabla^2\psi}^{x,y},
E = \tfrac{1}{2} \overline{|\bm{\nabla} \psi|^2}^{x, y} = -\tfrac{1}{2} \overline{\psi \nabla^2 \psi}^{x, y} ,
```

where the overbar denotes average over ``x`` and ``y`` and an integration-by-parts was carried
through in the last equality. To obtain the energy equation we multiply (6) with ``-\psi`` and
average over the whole domain. Thus, the work done by the forcing is given by:

```math
P = - \, \overline{\psi \, \xi}^{x,y},
P = - \, \overline{\psi \, \xi}^{x, y} ,
```

but the above is a stochastic integral and it is meaningless without a rule for computing the stochastic integral.
Expand All @@ -455,7 +463,7 @@ But how much is the Itô drift term in this case? As in the previous section, th
*precisely* the ensemble mean of the Stratonovich work, i.e.:

```math
\textrm{Ito drift}= - \overline{\langle \underbrace{\psi(\bm{x}, t) \circ \xi(\bm{x}, t)}_{\textrm{Stratonovich}} \rangle}^{x,y} .
\textrm{Ito drift}= - \overline{\langle \underbrace{\psi(\bm{x}, t) \circ \xi(\bm{x}, t)}_{\textrm{Stratonovich}} \rangle}^{x, y} .
```

But again, the above can be computed using the "formal" solution of (6):
Expand All @@ -468,10 +476,10 @@ which implies

```math
\begin{aligned}
\text{drift} & = -\overline{e^{-\mu t} \underbrace{\left \langle \psi(\bm{x}, 0) \xi(\bm{x}, t) \right \rangle}_{=0}}^{x,y} - \int_0^t e^{- \mu (t - s)} \overline{\nabla^{-2} \left \langle \xi(\bm{x}, s) \xi(\bm{x}, t) \right\rangle}^{x,y} \, \mathrm{d} s \\
& = -\int_0^t e^{-\mu(t - s)} \overline{\underbrace{\left [ \nabla^{-2} Q (\bm{x}) \right ] \big|_{\bm{x}=0}}_{\text{independent of }x,y} \, \delta(t - s)}^{x,y} \, \mathrm{d} s \\
& = -\frac1{2} \nabla^{-2} Q(\bm{x}) \big|_{\bm{x}=0} \\
& = - \frac1{2} \left [ \nabla^{-2} \int \frac{\mathrm{d}^2 \bm{k}}{(2\pi)^2} \widehat{Q}(\bm{k}) \, e^{i \bm{k} \bm{\cdot}\bm{x}} \right ]_{\bm{x}=0} \\
\text{drift} & = -\overline{e^{- \mu t} \underbrace{\left \langle \psi(\bm{x}, 0) \xi(\bm{x}, t) \right \rangle}_{=0}}^{x, y} - \int_0^t e^{- \mu (t - s)} \overline{\nabla^{-2} \left \langle \xi(\bm{x}, s) \xi(\bm{x}, t) \right\rangle}^{x, y} \, \mathrm{d} s \\
& = - \int_0^t e^{-\mu(t - s)} \overline{\underbrace{\left [ \nabla^{-2} Q (\bm{x}) \right ] \big|_{\bm{x}=0}}_{\text{independent of }x, y} \, \delta(t - s)}^{x,y} \, \mathrm{d} s \\
& = - \frac1{2} \nabla^{-2} Q(\bm{x}) \big|_{\bm{x}=0} \\
& = - \frac1{2} \left [ \nabla^{-2} \int \frac{\mathrm{d}^2 \bm{k}}{(2\pi)^2} \widehat{Q}(\bm{k}) \, e^{i \bm{k} \bm{\cdot} \bm{x}} \right ]_{\bm{x}=0} \\
& = \int \frac{\mathrm{d}^2 \bm{k}}{(2\pi)^2} \frac{\widehat{Q}(\bm{k})}{2 |\bm{k}|^2} .
\end{aligned}
```
Expand All @@ -488,10 +496,10 @@ Using the above, the work for a single forcing realization at the ``j``-th times
computed as:

```math
{\color{Green} \text{Itô}} : {\color{Green} P_j = -\overline{\psi(\bm{x}, t_j) \xi(\bm{x}, t_{j+1})}^{x,y} + \varepsilon} , \tag{8}
{\color{Green} \text{Itô}} : {\color{Green} P_j = -\overline{\psi(\bm{x}, t_j) \xi(\bm{x}, t_{j+1})}^{x, y} + \varepsilon} , \tag{8}
```
```math
{\color{Magenta} \text{Stratonovich}} : {\color{Magenta} P_j = -\overline{\frac{\psi(\bm{x}, t_j) + \psi(\bm{x}, t_{j+1})}{2} \xi(\bm{x}, t_{j+1})}^{x,y}} . \tag{9}
{\color{Magenta} \text{Stratonovich}} : {\color{Magenta} P_j = -\overline{\frac{\psi(\bm{x}, t_j) + \psi(\bm{x}, t_{j+1})}{2} \xi(\bm{x}, t_{j+1})}^{x, y}} . \tag{9}
```

Remember, previously the work done by the stochastic forcing was:
Expand All @@ -505,11 +513,7 @@ and by sampling over various forcing realizations:

All modules in GeophysicalFlows.jl use Stratonovich calculus. For example, the domain-averaged
energy injected per unit time by the forcing in the `TwoDNavierStokes` module is computed
using (9) via

```@docs
GeophysicalFlows.TwoDNavierStokes.energy_work
```
using (9) via the [`energy_work`](@ref GeophysicalFlows.TwoDNavierStokes.energy_work) function.

## A bit more elaborate SPDE

Expand All @@ -525,12 +529,12 @@ stochastic forcing, remains the same. We can easily verify this from the "formal
of (10):

```math
\psi(\bm{x}, t) = e^{- \mu t} \psi(\bm{x}, 0) + \int_0^t e^{- \mu (t - s)} \nabla^{-2} \xi(\bm{x}, s) \, \mathrm{d} s - \int_0^t \nabla^{-2} \mathsf{J} \left ( \psi(\bm{x}, s), \nabla^2 \psi(\bm{x}, s) \right ) \mathrm{d} s ,
\psi(\bm{x}, t) = e^{- \mu t} \psi(\bm{x}, 0) + \int_0^t e^{- \mu (t - s)} \nabla^{-2} \xi(\bm{x}, s) \, \mathrm{d} s - \int_0^t \nabla^{-2} \mathsf{J} \left ( \psi(\bm{x}, s), \nabla^2 \psi(\bm{x}, s) \right ) \mathrm{d} s .
```

When multiplied with ``\xi(\bm{x}, t)`` the last term vanishes since its only non-zero
contribution comes from the point ``s = t``, which is of measure zero (in the integrated sense).

A demonstration of how the energy budgets can be computed for a case with stochastic forcing
can be found in an [example of the TwoDNavierStokes](../generated/twodnavierstokes_stochasticforcing_budgets/)
A demonstration of how the energy budgets can be computed when we have stochastic forcing is
illustrated in an [example of the TwoDNavierStokes](../generated/twodnavierstokes_stochasticforcing_budgets/)
module.