Fitzhugh-Nagumo --- Labyrinthine Patterns

Simulation of the Fitzhugh-Nagumo equations $$ u_t = u - u^3 - v + \nabla^2 u\\ v_t = \epsilon(u - a_1 v - a_0) + \delta \nabla^2 v, $$ using the semi-spectral time integration method.

This simultation was heavily inspired by Aric Hagberg's simulation in "From Labyrinthine Patterns to Spiral Turbulence", PRL 1994.

The code below provides 3 initial conditions, "squiggle, blocks, and random". For time integration, besides the spectral method, we also provide the Euler method. Details about the semi-spectral method can be found after the code.

Parameters: $\epsilon=0.3$, $\delta=2.0$, $a_1=1.4$, and $a_0=0$.

Other simulations and Python examples can be found on my website: yairmau.com.

In [1]:
from IPython.display import HTML

HTML('<iframe width="560" height="315" src="https://www.youtube.com/embed/5au-G5FuI_A?rel=0" frameborder="0" allow="autoplay; encrypted-media" allowfullscreen></iframe>')
Out[1]:
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import os
from matplotlib import rcParams
rcParams['font.family'] = 'monospace'
In [ ]:
class FitzHughNagumo(object):
    def __init__(self, epsilon=0.3, delta=2.0, a1=1.4, a0=0.0,
                 n=(256, 256), l=(400, 400),
                 start=0.0, step=1.0, finish=2000.0, dt=0.1,
                 integration_type="spectral"):
        self.epsilon = epsilon
        self.delta = delta
        self.a1 = a1
        self.a0 = a0
        self.n = n
        self.l = l
        self.start = start
        self.step = step
        self.finish = finish
        self.dt = dt
        self.integration_type = integration_type
        self.rhs_a = np.zeros((2, self.n[0], self.n[1]))

    def spectral_multiplier(self):
        dx = float(self.l[0]) / self.n[0]
        dy = float(self.l[1]) / self.n[1]
        # wave numbers
        fx = 2.0 * np.pi * np.fft.fftfreq(self.n[0], dx)
        fy = 2.0 * np.pi * np.fft.fftfreq(self.n[1], dy)
        kx = np.outer(fx, np.ones(self.n[0]))
        ky = np.outer(np.ones(self.n[1]), fy)
        # multiplier
        mult_a = np.zeros((2, self.n[0], self.n[1]))
        mult_a[0] = np.exp(-(kx ** 2 + ky ** 2) * self.dt)  # u
        mult_a[1] = np.exp(-self.delta * (kx ** 2 + ky ** 2) * self.dt)  # v
        return mult_a

    def rhs_reaction(self, a):
        u = a[0]  # alias
        v = a[1]  # alias
        # FHN right hand side
        self.rhs_a[0] = u - u ** 3 - v
        self.rhs_a[1] = self.epsilon * (u - self.a1 * v - self.a0)
        return self.rhs_a

    def rhs_euler(self, a):
        # boundary conditions in laplacian
        laplacian = self.periodic_laplacian
        u = a[0]  # alias
        v = a[1]  # alias
        dx = float(self.l[0]) / self.n[0]
        # FHN right hand side
        self.rhs_a[0] = u - u ** 3 - v + laplacian(u, dx=dx)
        self.rhs_a[1] = self.epsilon * (u - self.a1 * v - self.a0) + \
                        self.delta * laplacian(v, dx=dx)
        return self.rhs_a

    def draw(self, a, t):
        u = a[0]
        self.im = plt.imshow(u.real, cmap="Greys_r", origin='lower',
                             vmin=-0.534522, vmax=0.534522,
                             interpolation="gaussian")
        self.title = plt.title('time = {:>4.0f}'.format(0))
        plt.xticks([])
        plt.yticks([])
        self.im.figure.canvas.draw()

    def draw_update(self, a, t):
        u = a[0]
        self.title.set_text('time = {:>4.0f}'.format(t))
        self.im.set_data(u.real)
        self.im.figure.canvas.draw()

    def save_frame(self, i):
        fname = "_tmp{:05d}.png".format(i)
        self.frame_names.append(fname)
        self.fig.savefig(fname, bbox_inches='tight', dpi=300)

    def periodic_laplacian(self, u, dx=1):
        """Return finite difference Laplacian approximation of 2d array.
        Uses periodic boundary conditions and a 2nd order approximation."""
        laplacian = (np.roll(u, -1, axis=0) +
                     np.roll(u, +1, axis=0) +
                     np.roll(u, -1, axis=1) +
                     np.roll(u, +1, axis=1) - 4.0 * u) / (dx ** 2)
        return laplacian

    def random_ic(self):
        return 0.5 * (np.random.random((2, self.n[0], self.n[1])) - 0.5)

    def blocks_ic(self):
        a = np.ones((2, self.n[0], self.n[1]))
        a[0] = 0.534522
        a[1] = 0.381802
        n = self.n
        p = n[0] / 8
        a[0][3 * p - 4:3 * p + 4, 5 * p - 4:5 * p + 4] = -0.534522
        a[0][6 * p - 4:6 * p + 4, 3 * p - 4:3 * p + 4] = -0.534522
        return a

    def squiggle_ic(self):
        a = np.ones((2, self.n[0], self.n[1]))
        l = self.l
        uplus = 0.534522
        vplus = 0.381802
        uminus = -uplus
        X, Y = np.meshgrid(np.linspace(0, self.l[0], self.n[0]),
                           np.linspace(0, self.l[0], self.n[0]))
        cos_term = 0.05 * l[0] * np.sin(10 * (2 * np.pi) * Y / l[1] + np.pi * 0.3)
        exp_term = np.exp(-((Y - l[1] / 2) / (0.1 * l[1])) ** 2)
        width = 0.05 * l[0]
        Z = np.exp(-((X - l[0] / 2 + cos_term * exp_term) / width) ** 2)
        a[0] = uplus
        a[1] = vplus
        a[0][Z > 0.8] = uminus
        return a
In [ ]:
plt.ion()
plt.clf()
foo = FitzHughNagumo()
foo.fig = plt.figure(1)
ax = foo.fig.add_subplot(111)
a = foo.squiggle_ic()
mult_a = foo.spectral_multiplier()
fft_a = np.fft.fftn(a, axes=(1, 2))
t = foo.start
foo.draw(a, t)
foo.frame_names = []
foo.save_frame(0)
for i, tout in enumerate(np.arange(foo.start + foo.step,
                                   foo.finish + foo.step,
                                   foo.step)):
    while t < tout:
        if foo.integration_type == "spectral":
            rhs_a = foo.rhs_reaction(a)
            fft_a = mult_a * (fft_a + foo.dt * np.fft.fftn(rhs_a, axes=(1, 2)))
            a = np.fft.ifftn(fft_a, axes=(1, 2))
        if foo.integration_type == "euler":
            a = a + foo.dt * foo.rhs_euler(a)
        t += foo.dt
    foo.draw_update(a, t)
    foo.save_frame(i + 1)
In [ ]:
fps = 24
frames = "_tmp%5d.png"
movie_command = "ffmpeg -y -r {:} -i {:} fhn.mp4".format(fps, frames)
os.system(movie_command)
for fname in foo.frame_names:
    os.remove(fname)

The semi-spectral method

The explanation below was taken from my thesis: "Pattern Formation in Spatially Forced Systems: Application to Vegetation Restoration".

The semi-spectral method is extremely useful when working with reaction-diffusion systems, and with parabolic PDEs in general. This was the method used to run all the simulations of the Swift-Hohenberg model in this thesis, and it proved to be reliable and fast. The explanation below is a summary of "Spectral algorithms for reaction-diffusion equations", by Richard V. Craster and Roberto Sassi, with a step by step recipe, so the reader can easily apply the method to any suitable problem.

The method

The semi-spectral transform method is very useful when we have to integrate a system that evolves really slowly. Let us say we have a (parabolic) system of the form: $$ u_t=\epsilon u + f(u)+D\nabla^2u, \qquad\qquad (1) $$ where $f(u)$ is a nonlinear function. First, we compute the Fourier transform of (1): $$ \hat{u}_t=\epsilon \hat{u} + \hat{f}(u)-k^2D\hat{u}, \qquad\qquad (2) $$ where the hat denotes the Fourier transform.

We rearrange (2) in the following way: $$ \hat{u}_t+a\hat{u}=\hat{f}(u), \qquad\qquad (3) $$ where $a=-\epsilon +k^2D$, and now we make a variable substitution $$ \hat{v}(k,t)=\;\hat{u}(k,t)\,e^{at} \qquad\qquad (4a)\\ \hat{v}_t=\;\hat{u}_te^{at}+a\hat{u}\,e^{at}. \qquad\qquad (4b) $$

We multiply (3) by $e^{at}$ and we finally get $$ \hat{v}_t=e^{at}\hat{f}(u). \qquad\qquad (5) $$

We can now advance $\hat{v}$ in time using a simple Euler step $$ \hat{v}^{t_{n+1}}=\hat{v}^{t_n}+\Delta t\left( e^{at_n}\hat{f}(u) \right). \qquad\qquad (6) $$

What we really want is $\hat{u}$, which, according to (4a), is given by $$ \begin{align} \hat{u}^{t_{n+1}}=&\;\hat{v}^{t_{n+1}}e^{-at_{n+1}}\\ =&\;\hat{v}^{t_{n+1}}e^{-at_{n}}e^{-a\Delta t}\\ =&\;\left(\hat{v}^{t_n}+\Delta t\; e^{a t_n}\hat{f}(u)\right)e^{-at_{n}}e^{-a\Delta t}\\ =&\;\left(\hat{v}^{t_n}e^{-at_{n}}+\Delta t \;{e^{a t_n}}\hat{f}(u) {e^{-at_{n}}}\right)e^{-a\Delta t}\\ =&\;\left( \hat{u}^{t_n}+\Delta t \hat{f}(u) \right)e^{-a\Delta t}. \end{align} $$

There is actually no need to use the variable substitution in (4). We now have an expression for $\hat{u}^{t_{n+1}}$:

$$ \hat{u}^{t_{n+1}}=\left( \hat{u}^{t_n}+\Delta t \hat{f}(u) \right)e^{-a\Delta t}. $$

Now it is time to go back from the Fourier space to the real space, and for that we use an inverse Fourier transform $$ u^{t_{n+1}}=\mathcal{F}^{-1}[\hat{u}^{t_{n+1}}]. $$

Step by step

To implement this technique, one just has to follow the steps below:

  1. Calculate the Fourier transform of $u$: $\hat{u}=\mathcal{F}[u]$.
  2. Have $f(u)$ calculated and then take its Fourier transform: $\hat{f}(u)=\mathcal{F}[f(u)]$.
  3. For a given lattice with $N$ points, and $\delta x$ being the distance between them, make the frequency bin vector (matrix) $k$ for your one (two) dimensional system. In python the command would be numpy.fft.fftfreq(N, dx). The frequency bin vector $k$ looks like: $$ \begin{align} k&=2\pi\cdot\left[ 0,1,\cdots,\tfrac{N}{2}-1,-\tfrac{N}{2},\cdots,-1 \right]/(N\, \delta x), \qquad \mbox{if N is even;}\\ k&=2\pi\cdot\left[ 0,1,\cdots,\tfrac{N-1}{2},-\tfrac{N-1}{2},\cdots,-1 \right]/(N\, \delta x), \qquad \mbox{if N is odd.} \end{align} $$ Remember that the domain size is given by $L=N\, \delta x$, which means that the denominator in the expressions above can be written simply as $L$. It is clear from that fact that $\delta k$, the tiniest slice of the Fourier space is $\delta k=2\pi/L$. Corollary: if you want to divide the Fourier space into very many parts, simply have a huge domain.
    If the system is two-dimensional, then have $k_x$ and $k_y$ calculated separately. The domain might not be square ($L_x\neq L_y$), and you might want to divide the domain into a different number of points ($Nx\neq Ny$). Anyway, prepare one-dimensional arrays of $k_x$ and $k_y$ as explained above, and then make an outer product of these arrays with a ones array of length $N$, as following: $$ k_{x,2d} = \begin{pmatrix}

          1 \\ 1 \\ \vdots \\ 1
         \end{pmatrix} 
         \begin{pmatrix}
          k_{x1} & k_{x2} & ... & k_{xN}
         \end{pmatrix}
         =
         \begin{pmatrix}
          k_{x1} & k_{x2} & ... & k_{xN}\\
          k_{x1} & k_{x2} & ... & k_{xN}\\
                 & \vdots &        &       \\
          k_{x1} & k_{x2} & ... & k_{xN}
         \end{pmatrix}
    

    $$ and $$ k_{y,2d} = \begin{pmatrix}

          k_{x1} \\ k_{x2} \\ \vdots \\ k_{xN}
         \end{pmatrix} 
         \begin{pmatrix}
          1 & 1 & ... & 1
         \end{pmatrix}
         =
         \begin{pmatrix}
          k_{y1} & k_{y1} &        & k_{y1}\\
          k_{y2} & k_{y2} & ... & k_{y2}\\
          \vdots & \vdots &        & \vdots\\
          k_{yN} & k_{yN} &        & k_{yN}
         \end{pmatrix}.
    

    $$ Then factor $e^{-a\Delta t}$ equals $$ e^{-a\Delta t}= e^{\left[\epsilon-D(k_x^2+k_y^2)\right]\Delta t}, $$ where $kx^2$ is the element-wise exponentiation of the 2d array $k{x,2d}$.

  4. Now that we have all the factors we need, we simply calculate $$ \hat{u}^{t_{n+1}}=\left( \hat{u}^{t_n}+\Delta t \hat{f}(u) \right)e^{\left[\epsilon-D(k_x^2+k_y^2)\right]\Delta t}. $$

  5. We finally go back to the real space by applying the inverse Fourier transform: $u^{t_{n+1}}=\mathcal{F}^{-1}[\hat{u}^{t_{n+1}}]$.

Example

For the parametrically forced Swift-Hohenberg equation $$ \frac{\partial u}{\partial t} = [\epsilon + \gamma \cos(k_f x)]u - u^3 -(\nabla^2+k_0^2)^2 u, $$ we have $$ f(u)= -u^3 + \gamma u \cos(k_f x),\qquad a = \epsilon - \left(k_0- k_x^2 - k_y^2 \right)^2. $$