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.
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>')
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 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 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}}]. $$
To implement this technique, one just has to follow the steps below:
numpy.fft.fftfreq(N, dx)
. The frequency bin vector $k$ looks like: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 $k_x^2$ is the element-wise exponentiation of the 2d array $k_{x,2d}$.
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. $$