Fitzhugh-Nagumo Equation
Simulating the Fitzhugh-Nagumo equations using a spectral method
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.
import numpy as np
import matplotlib.pyplot as plt
import os
from matplotlib import rcParams
rcParams['font.family'] = 'monospace'
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
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)
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 beThe frequency bin vector $k$ looks like:
numpy.fft.fftfreq(N, dx).
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:
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} $$