We simulate the two dimensional stochastic Allen-Cahn equation by using a finite element method on the spatial domain $[0,1]^2$ and a semi-implicit Euler method on the temporal domain $[0,T]$.

To begin with, we take a family of triangulations $(\mathcal{T}_h)_h$ of $[0,1]^2$ and equipartitions $(\mathcal{I}_k)_k$ of $[0,T]$. We use an approximation $\overline{\xi}$ of the space time white noise on these families, defined as follows:

$\displaystyle \overline{\xi}(\mathbf{x},t) := \sum_{n=1}^N\sum_{m=1}^M \overline{\eta}_{m,n}\mathbf{1}_{K_m}(\mathbf{x})\mathbf{1}_{I_n}(t),\;\;\;K_m \in \mathcal{T}_h, I_n \in \mathcal{I}_h$

where the $\overline{\eta}_{m,n}$ are IID Gaussian random variables defined by

$ \displaystyle\overline{\eta}_{m,n} := \frac{\sqrt{2}}{kh^{3/2}}\int_{[0,T]}\int_{[0,1]^2} \mathbf{1}_{K_m}(\mathbf{x})\mathbf{1}_{I_n}(t)\,\xi(\mathrm{d}\mathbf{x},\mathrm{d}t) \sim N\left(0,\frac{1}{kh}\right)$

This approximations converges simultaneously with the grid. The numerical scheme then reads

$ \displaystyle\left\langle \frac{u^n - u^{n-1}}{k},\varphi\right\rangle + \left\langle \nabla u^n,\nabla \varphi\right\rangle + \left\langle f_{\varepsilon}(u^n),\varphi\right\rangle = \varepsilon^{\gamma}\left\langle \overline{\xi},\varphi\right\rangle\;\;\;\forall\varphi \in V_h $

where $V_h$ is a Galerkin space of continuous piecewise linear functions associated with $\mathcal{T}_h$ and $\gamma \in \mathbb{R}$ is a scaling parameter. By a Taylor approximation of $f_{\varepsilon}(u^n)$ we can rewrite this as

$\displaystyle\left\langle \frac{u^n - u^{n-1}}{k},\varphi\right\rangle + \left\langle \nabla u^n,\nabla \varphi\right\rangle + \left\langle f_{\varepsilon}'(u^{n-1})u^n,\varphi\right\rangle = \left\langle f_{\varepsilon}'(u^{n-1})u^{n-1} - f_{\varepsilon}(u^{n-1}),\varphi \right\rangle + \varepsilon^{\gamma}\left\langle \overline{\xi},\varphi\right\rangle\;\;\;\forall\varphi \in V_h $

which has the advantage of only needing to solve a linear system at each timestep: in matrix form this becomes

$\displaystyle \left[\frac{1}{k}\mathbf{M} +\mathbf{A} + \mathbf{N}(\mathbf{u}^{n-1})\right]\mathbf{u}^n = \mathbf{g}(\mathbf{u}^{n-1}) + \frac{1}{k}\mathbf{M}\mathbf{u}^{n-1} + \varepsilon^{\gamma}\mathbf{w} $

where the matrix terms $\mathbf{M}$, $\mathbf{A}$ and $\mathbf{N}(\mathbf{u}^{n-1})$ are all sparse. The entries of the vector $\mathbf{w}$ are scaled sums of a small number of IID $N(0,1)$ random variables.

The scheme was implemented in C++ using the DUNE numerics environment. You may now want to see some numerical experiments obtained with our code.