1 Fitzhugh-Nagumo Equation
1.1 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.
1.2 The code
1.2.1 import packages
1.2.2 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[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 a1.2.3 run simulation, save snapshots
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)1.2.4 make movie, delete snapshots
1.3 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.
1.3.1 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, \tag{1}
where f(u) is a nonlinear function. First, we compute the Fourier transform of : \hat{u}_t=\epsilon \hat{u} + \hat{f}(u)-k^2D\hat{u}, \tag{2} where the hat denotes the Fourier transform.
We rearrange in the following way: \hat{u}_t+a\hat{u}=\hat{f}(u), \tag{3} where a=-\epsilon +k^2D, and now we make a variable substitution
\begin{align} \hat{v}(k,t)&=\;\hat{u}(k,t)\,e^{at} \tag{4a}\\ \hat{v}_t&=\;\hat{u}_te^{at}+a\hat{u}\,e^{at}. \tag{4b} \end{align}
We multiply by e^{at} and we finally get \hat{v}_t=e^{at}\hat{f}(u). \tag{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). \tag{6}
What we really want is \hat{u}, which, according to , is given by
\begin{align} \displaystyle\hat{u}^{t_{n+1}}=&\;\hat{v}^{t_{n+1}}e^{-at_{n+1}} \tag{7a}\\ =&\;\hat{v}^{t_{n+1}}e^{-at_{n}}e^{-a\Delta t} \tag{7b}\\ =&\;\left(\hat{v}^{t_n}+\Delta t\; e^{a t_n}\hat{f}(u)\right)e^{-at_{n}}e^{-a\Delta t} \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} \tag{7d}\\ =&\;\left( \hat{u}^{t_n}+\Delta t \hat{f}(u) \right)e^{-a\Delta t} \tag{7e}. \end{align}
There is actually no need to use the variable substitution in . 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}. \tag{8}
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}}]. \tag{9}
1.3.2 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
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), \quad \text{if N is even;} \tag{10a}\\ k&=2\pi\cdot\left[ 0,1,\cdots,\tfrac{N-1}{2},-\tfrac{N-1}{2},\cdots,-1 \right]/(N\, \delta x), \quad \text{if N is odd.} \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} \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}. \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}, \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}. \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}}].
1.3.3 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, \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. \tag{16}