# Fitzhugh-Nagumo Equation
> Simulating the Fitzhugh-Nagumo equations using a spectral method

- toc: true 
- badges: true
- comments: false
- categories: [jupyter]

> youtube: https://youtu.be/5au-G5FuI_A

[https://youtu.be/5au-G5FuI_A](https://youtu.be/5au-G5FuI_A)

## 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](http://in.bgu.ac.il/en/bidr/SIDEER/DSEEP/Ehud_Meron/Documents/Papers/HaMe94prl.pdf).

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](yairmau.com).

## The code

### import packages


```python
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


```python
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
```

### run simulation, save snapshots


```python
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)
```

### make movie, delete snapshots


```python
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"](http://aranne5.bgu.ac.il/others/MauYair3.pdf).

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"](https://pdfs.semanticscholar.org/6bb0/b794a07b0c9b19e7c7ab03a9a523eac2141d.pdf), 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
```python3
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}
$$




```python

```

be92f1c872179bd6e9259148612fc2f2a024fdd8
ipynb
Jupyter Notebook
_notebooks/2020-01-01-fhn.ipynb
yairmau/website
74d8296d1536a2b65746e089951c5a11c643c4ed
Apache-2.0
2020-05-30T13:17:26.000Z
2021-12-23T14:45:11.000Z
_notebooks/2020-01-01-fhn.ipynb
yairmau/website
74d8296d1536a2b65746e089951c5a11c643c4ed
Apache-2.0
2020-09-26T07:34:07.000Z
2021-08-05T12:54:43.000Z
_notebooks/2020-01-01-fhn.ipynb
yairmau/website
74d8296d1536a2b65746e089951c5a11c643c4ed
Apache-2.0
2020-07-14T11:30:02.000Z
2020-07-14T11:30:02.000Z
Qwen/Qwen-72B
1. YES
2. YES