## Introduction

We simulate 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.

## The code

### import packages

import numpy as np
import matplotlib.pyplot as plt
import os
from matplotlib import rcParams
rcParams['font.family'] = 'monospace'


### define class with all the methods

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, self.n))

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

def rhs_reaction(self, a):
u = a  # alias
v = a  # alias
# FHN right hand side
self.rhs_a = u - u ** 3 - v
self.rhs_a = 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  # alias
v = a  # alias
dx = float(self.l) / self.n
# FHN right hand side
self.rhs_a = u - u ** 3 - v + laplacian(u, dx=dx)
self.rhs_a = 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
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
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, self.n)) - 0.5)

def blocks_ic(self):
a = np.ones((2, self.n, self.n))
a = 0.534522
a = 0.381802
n = self.n
p = n / 8
a[3 * p - 4:3 * p + 4, 5 * p - 4:5 * p + 4] = -0.534522
a[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, self.n))
l = self.l
uplus = 0.534522
vplus = 0.381802
uminus = -uplus
X, Y = np.meshgrid(np.linspace(0, self.l, self.n),
np.linspace(0, self.l, self.n))
cos_term = 0.05 * l * np.sin(10 * (2 * np.pi) * Y / l + np.pi * 0.3)
exp_term = np.exp(-((Y - l / 2) / (0.1 * l)) ** 2)
width = 0.05 * l
Z = np.exp(-((X - l / 2 + cos_term * exp_term) / width) ** 2)
a = uplus
a = vplus
a[Z > 0.8] = uminus
return a


### run simulation, save snapshots

plt.ion()
plt.clf()
foo = FitzHughNagumo()
foo.fig = plt.figure(1)
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)


### make movie, delete snapshots

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: $$\begin{equation*} u_t=\epsilon u + f(u)+D\nabla^2u, \label{eq:1} \tag{1} \end{equation*}$$

where $f(u)$ is a nonlinear function. First, we compute the Fourier transform of \eqref{eq:1}: $$\begin{equation*} \hat{u}_t=\epsilon \hat{u} + \hat{f}(u)-k^2D\hat{u}, \label{eq:2} \tag{2} \end{equation*}$$ where the hat denotes the Fourier transform.

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

We multiply \eqref{eq:3} by $e^{at}$ and we finally get $$\begin{equation*} \hat{v}_t=e^{at}\hat{f}(u). \label{eq:5} \tag{5} \end{equation*}$$

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

What we really want is $\hat{u}$, which, according to \eqref{eq:4a}, is given by

\begin{align*} \displaystyle\hat{u}^{t_{n+1}}=&\;\hat{v}^{t_{n+1}}e^{-at_{n+1}} \label{eq:7a} \tag{7a}\\ =&\;\hat{v}^{t_{n+1}}e^{-at_{n}}e^{-a\Delta t} \label{eq:7b} \tag{7b}\\ =&\;\left(\hat{v}^{t_n}+\Delta t\; e^{a t_n}\hat{f}(u)\right)e^{-at_{n}}e^{-a\Delta t} \label{eq:7c} \tag{7c}\\ =&\;\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} \label{eq:7d} \tag{7d}\\ =&\;\left( \hat{u}^{t_n}+\Delta t \hat{f}(u) \right)e^{-a\Delta t} \label{eq:7e} \tag{7e}. \end{align*}

There is actually no need to use the variable substitution in \eqref{eq:4a}. We now have an expression for $\hat{u}^{t_{n+1}}$: $$\begin{equation*} \hat{u}^{t_{n+1}}=\left( \hat{u}^{t_n}+\Delta t \hat{f}(u) \right)e^{-a\Delta t}. \label{eq:8} \tag{8} \end{equation*}$$

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}}]. \label{eq:9} \tag{9}$$

### step by step

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

• Calculate the Fourier transform of $u$: $\hat{u}=\mathcal{F}[u]$.
• Have $f(u)$ calculated and then take its Fourier transform: $\hat{f}(u)=\mathcal{F}[f(u)]$.
• 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;} \label{eq:10a} \tag{10a}\\ 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.} \label{eq:10b} \tag{10b} \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 ($N_x\neq N_y$). 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} \label{eq:11} \tag{11}$$

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}. \label{eq:12} \tag{12}$$

Then factor $e^{-a\Delta t}$ equals

$$e^{-a\Delta t}= e^{\left[\epsilon-D(k_x^2+k_y^2)\right]\Delta t}, \label{eq:13} \tag{13}$$

where $k_x^2$ is the element-wise exponentiation of the 2d array $k_{x,2d}$.

• 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}. \label{eq:14} \tag{14}$$
• 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, \label{eq:15} \tag{15}$$ 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. \label{eq:16} \tag{16}$$