Numerical Methods
Numerical Methods
To solve the Barkley model numerically we need some numerical methods to discretise in space and time. Two schemes for the numerical approximation are considered: the \(\theta\)-scheme and the semi-implicit method.
Finite Difference method
Time discretisation
Starting with the Barkley equation, we discretise in time with the \(\theta\)-scheme with timestep \(\tau\) to give
\begin{align} T_\theta u^{n+1} &= T_{\theta-1}u^n + \tau \tilde f^n \\ v^{n+1} &= v^n + \tau \tilde g^n,\end{align}
where \(T_\theta\) is the linear operator defined by
\begin{equation*} (T_\theta u)(x)=u(x)-\theta\tau a\Delta_\Gamma u(x), \end{equation*}
\(\tilde f^n=\tilde f^n(u^n,u^{n+1},v^n,v^{n+1})\) is the term involving \(f\) and \(\tilde g^n=\tilde g^n(u^n,u^{n+1},v^n,v^{n+1})\) is the one involving \(g\) (which will be discussed later). We also need to discretise the space, which the following section addresses.
Space discretisation
We mainly focus on the excitation variable \(u\) because the equation for \(v\) is much more straightforward. When we are only simulating the equations on squares with side length \(L\), FD is a simple and fast way to produce results.
By dividing each coordinate axis into \(N+1\) evenly-spaced points which we denote by \((x_i)_0^N\) and \((y_i)_0^N\), the solution is approximated by interpolations from the vectors \((u^n(x_i,y_j))_{i,j=0}^N\). Both the Laplacian are approximated using the central difference.
Then \(T_\theta\) is approximated by matrices \(A_{\theta,h}=I-\theta\tau a\Delta_h\) operating on vectors \(u^n\), hence if (and only if) \(\tilde f^n\) is linear with respect to \(u^{n+1}\), then the discretisation of \(u\) is reduced to a system of linear equations: \begin{equation}\label{eq:tmp} A_{\theta,h}u^{n+1}=A_{\theta-1,h}u^{n}+\tau \tilde f^n. \end{equation}
Surface Finite Element method
We then discretise in space using surface finite elment method (FE). Multiplying by a test function \(\phi\colon\mathcal{G}_T\to\mathbb{R}\) \(\left(\text{where }\mathcal{G}_T = \bigcup_{t\in[0,T]}(\{t\}\times\Gamma_t)\right)\) such that \(\phi(t,\cdot) \in H^1(\Gamma_t)\) and integrating, we obtain the weak formulation of \(u\):
\begin{align*} \int_{\Gamma_t}\dot u\phi+\int_{\Gamma_t}u\phi\nabla_{\Gamma_t}\cdot \mathbf v+\int_{\Gamma_t}a\nabla_{\Gamma_t}u\cdot\nabla_{\Gamma_t}\phi=\int_{\Gamma_t}f(u,v)\phi. \end{align*}
By using the Leibniz formula, this becomes
\begin{align*} \frac{d}{dt}\int u\phi+\int a\nabla_{\Gamma}u\cdot\nabla_{\Gamma}\phi=\int f(u,v)\phi+\int u\dot\phi. \end{align*}
The work done by Elliott et. al. shows that by considering a set of \(m\)-dimensional spaces \(V_t\subset H^1(\Gamma_t)\), a set of basis functions \((Z_i(t,\cdot))_{i\le m}\) can be chosen such that \(\dot Z_i=0\) for all \(i\). Therefore plugging \(\phi=Z_i\) into the previous equation, we have
\begin{align*} \frac{d}{dt}\int uZ_i+\int a\nabla_{\Gamma}u\cdot\nabla_{\Gamma}Z_i=\int f(u,v)Z_i. \end{align*}
Therefore discretising in time and applying the \(\theta\)-scheme, we have
\begin{align*} \int_{\Gamma_{(n+1)\tau}}u^{n+1}Z_i^{n+1}+\tau \theta a\nabla_\Gamma u^{n+1}\cdot\nabla_\Gamma Z^{n+1}_i=\int_{\Gamma_{n\tau}}u^nZ_i^n+\tau (\theta-1)a\nabla_\Gamma u^n\cdot\nabla_\Gamma Z_i^n+\tilde{F}^n, \end{align*}
where \(u^n=u(n\tau,\cdot)\), \(Z^n=Z(n\tau,\cdot)\) and \(\tilde{F}^n\) is the integral containing \(f\).
This yields
\begin{equation}\label{eq:u_fem} a_\theta^{n+1}(u^{n+1},Z_i^{n+1})=a_{\theta-1}^n(u^n,Z_i^n)+\tilde F^n, \end{equation}
where \(a_\theta^n\colon H^1(\Gamma_{n\tau})\times H^1(\Gamma_{n\tau})\to\mathbb{R}\) is a bilinear form defined by
\begin{align*} a_\theta^n(\xi,\eta)=\int_{\Gamma_{n\tau}}(\xi\eta+\tau \theta a\nabla_\Gamma \xi\cdot\nabla_\Gamma\eta). \end{align*}
FEM aims at approximating \(u^n\) by its projection on \(V\). This allows us to write \(u^n\) as
\begin{align*} u^n &\approx P_V(u^n)= \sum_{i=1} ^{N} \alpha_i^n Z_i^n \\ \end{align*}
with coefficients \(\alpha_i^n\in \mathbb{R}\), which are to be solved for. For simplicity we identify \(u^n\) with \(P_V(u^n)\). Then the bilinear form equation can be rewritten as:
\begin{equation*} \sum_{j=1}^{m}\alpha_j^{n+1}a_\theta^{n+1}(Z_j^{n+1},Z_i^{n+1})=\sum_{j=1}^m\alpha_j^na_{\theta-1}^n(Z_j^n,Z_i^n)+\tilde F^n \end{equation*}
for \(i = 1, ..., m.\) Hence if \(\tilde F^n\) is linear with respect to \(\alpha_j^{n+1}\), \(j\le m\), then the above equation is again reduced to a system of linear equations.
Approximating the right hand side
Now let us discuss how to deal with \(\tilde f^n\) and \(\tilde F^n\) terms in the time and space discretisation. If we work strictly in line with \(\theta\)-method, then the respectively \(\tilde f^n\) in the time-discretisation and \(\tilde F^n\) in the spatial discretisation are
\begin{align}\label{eq:rhs_dis} \tilde f^n&= \theta f(u^{n+1},v^{n+1})+ (1-\theta)f(u^n,v^n).\\ \nonumber \tilde F^n&=\tau\theta\int_{\Gamma_{(n+1)\tau}}f(u^{n+1},v^{n+1})Z^{n+1}_i+\tau(1-\theta)\int_{\Gamma_{n\tau}}f(u^{n},v^{n})Z^{n}_i\\ \nonumber &=\tau\theta\int_{\Gamma_{(n+1)\tau}}f\left(\sum_j\alpha_jZ^{n+1}_j,v^{n+1}\right)Z^{n+1}_i+\tau(1-\theta)\int_{\Gamma_{n\tau}}f(u^{n},v^{n})Z^{n}_i. \end{align}
When \(\theta>0\), \(f\) is nonlinear with respect to \(u^{n+1}\), which would inflict intolerable pains if one tries to implement (as we no longer have a simple matrix equation), especially when we set \(\theta=1\) for sake of stability. We consider two ways of solving this problem.
Explicit method
One solution is to consider \(\tilde f^n\) and \(\tilde F^n\) fully explicitly regardless of \(\theta\) in time discretisation: \begin{align*} \tilde f^n&=f(u^n,v^n)\\ \tilde F^n&=\tau\int_{\Gamma_{n\tau}}f(u^{n},v^{n})Z^{n}_i. \end{align*}
Then the the FD method at the end of the section above on space discretisation becomes a system of linear equations:
\begin{equation*} A_{\theta,h}u^{n+1}=A_{\theta-1,h}u^n+f(u^n,v^n). \end{equation*}
And the FE method:
\begin{equation*} \sum_{j=1}^m\alpha_j^{n+1}a_\theta^{n+1}(Z_j^{n+1},Z_i^{n+1})=\sum_{j=1}^m\alpha_j^na_{\theta-1}^n(Z_j^n,Z_i^n)+\int f(u^n,v^n)Z_i^n;\quad\forall i=1,\dots,m. \end{equation*}
We call this the explict RHS method. It requires a small time step to ensure stability.
Semi-implicit method
An alternative method to the explicit method is to test with a semi-implicit RHS, which is half-way between fully explicit (as described in the previous solution) and fully implicit (which is what we would obtain in the \(\theta-\)scheme with \(\theta=1\)). It has the advantage of allowing for a larger timestep without succumbing to numerical instability unlike the explicit RHS. The idea is as follows. First, consider the FD method. Since our \(f\) is a factorized cubic polynomial, take
\begin{equation*} \tilde f^n=\frac{1}{\epsilon}u^{n_1}(1-u^{n_2})\left(u^{n_3}-\frac{v^{n_3}+b}{c}\right), \end{equation*}
where \(n_1,n_2,n_3\) are either \(n\) or \(n+1\). When \(n_1=n_2=n_3=n\), \(\tilde f^n\) is the fully explicit RHS described in previous section. When \(n_1=n_2=n_3=n+1\), \(\tilde f^n\) is the fully implicit RHS, i.e., take \(\theta=1\) in the definition given in the section on approximating the right-hand side. The semi-implicit method takes two of the \(n_i\) to be \(n\) and the remaining one to be \(n+1\) to make \(\tilde f^n\) linear in \(u^{n+1}\). Exactly which \(n_i\) is to be taken as \(n+1\) depends on the excitability threshold
\begin{equation*} u^n_{th} = {{v^n+b} \over c}. \end{equation*}
When \(u^n_{th}\ge u^n\), take \(n_1=n+1\); when \(u^n_{th}< u^n\), we take \(n_2=n+1\). Then the right-hand side function splits and our equationvbecomes:
\begin{align*} \left(A_{\theta,h}+\frac{\tau}{\epsilon}u^n(u^n-u^n_{th})I\right)u^{n+1}&=A_{\theta-1,h}u^{n}+\frac{\tau}{\epsilon}u^n(u^n-u^n_{th})\;\;\text{ when }u^n_{th}< u^n,\\ \left(A_{\theta,h}-\frac{\tau}{\epsilon}(1-u^n)(u^n-u^n_{th})I\right)u^{n+1}&=A_{\theta-1,h}u^{n}\;\;\;\; \ \;\;\;\; \;\;\;\; \;\;\;\; \;\;\;\; \text{ when }u^n_{th}\ge u^n, \end{align*}
where \(I\) is the identity matrix. For the FE method a corresponding semi-implicit \(\tilde F^n\) is
\begin{align*} \tilde F^n=\begin{cases} \frac{\tau}{\epsilon}\int_{\Gamma_{(n+1)\tau}}u^{n+1}(1-\bar u^n)(\bar u^n-\bar u^n_{th}),&\text{ if }u^n_{th}\ge u^n;\\ \frac{\tau}{\epsilon}\int_{\Gamma_{n\tau}}u^n(u^n-u^n_{th})-\frac{\tau}{\epsilon}\int_{\Gamma_{(n+1)\tau}}\bar u^n(\bar u^n-\bar u^n_{th})u^{n+1},&\text{ if }u^n_{th}< u^n, \end{cases} \end{align*}
where the \(\bar u^n=u^n\circ\Phi_{(n+1)\tau}^{n\tau}\) and \(\bar u^n_{th}=u^n_{th}\circ\Phi_{(n+1)\tau}^{n\tau}\) are \(u^n\) and \(u^n_{th}\) taking values at mapped points. Then our variational form becomes a system of linear equations:
\begin{align*} \sum_{j=1}^m\alpha_j^{n+1}\left(a^{n+1}_\theta(Z_j^{n+1},Z_i^{n+1})-\frac{\tau}{\epsilon}\int Z^{n+1}_j(1-\bar u^n)(\bar u^n-\bar u^n_{th})\right)&=\sum_{j=1}^m\alpha_j^na_{\theta-1}^n(Z_j^n,Z_i^n)\\\text{ when }u^n_{th}\ge u^n, \text{ and}\\ \sum_{j=1}^m\alpha_j^{n+1}\left(a^{n+1}_\theta(Z_j^{n+1},Z_i^{n+1})+\frac{\tau}{\epsilon}\int Z^{n+1}_j\bar u^n(\bar u^n-\bar u^n_{th})\right)&=\sum_{j=1}^m\alpha_j^na_{\theta-1}^n(Z_j^n,Z_i^n)\\ &\;\;\;\;+\frac{\tau}{\epsilon}\int u^n(u^n-u^n_{th})\\\text{ when }u^n_{th}< u^n \end{align*}
for all \(i=1,\dots,m.\)
The recovery variable \(v\)
Solving \(v\) is a much simpler task than \(u\) due to its ODE-like equation. When using FD for space discretisation, \(g^n\) is treated fully explicitly:
\begin{equation*} \tilde g^n(u^n,v^n,u^{n+1},v^{n+1})=g(u^n,v^n). \end{equation*} Hence the scheme for \(v\) is explicit: \begin{equation*} v^{n+1}=v^{n}+\tau g(u^n,v^n). \end{equation*}
When using the FE method, we also treat \(g\) explicitly and we treat \(v\) in the same way as \(u\), i.e.,
\begin{equation*} \sum_{j=1}^m\beta_j^{n+1}\int Z^{n+1}_jZ^{n+1}_i=\sum_{j=1}^m\beta_j^n\int Z^n_jZ^n_i+\tau\int g(u^n,v^n)Z^n_i,\quad\forall i=1,\dots,m, \end{equation*}
where \(\beta_j^n\)'s are the coordinates of \(v^n\) in \(V_{n\tau}\):
\begin{equation*} v^n=\sum_{j=1}^m\beta_j^nZ^n_j. \end{equation*}
For the finite element simulations that we create in the following sections, we make use of the Distributed and Unified Numerics Environment
(DUNE) software. Material on the finite element method can be found here, and for information about the iterative solvers used, see here. The theory of the DUNE grid is described here and here.
