Skip to main content Skip to navigation

2D Simulation

Rescaling the domain

In order to run the Barkley model such that spirals are produced on a general square domain, some rescaling is necessary. For convenience, define \( \Omega_c \) for \( c \in \mathbb R ^+ \) to be the interval \( [0,c ]^2 \).
We know that the Barkley equation with parameters

\begin{equation*} \Gamma =\Omega_{150} , ~~ \epsilon = {1 \over {50}}, ~~ a=1, ~~b= {1 \over {100} } ~~c={3 \over 4 } \end{equation*}

and initial conditions

\begin{equation*} u_0 = \mathbb I _{ {y \over {150}} > {{51} \over {100}}} , ~~ v_0 = \mathbb I _{ {y \over {150}} > {{51} \over {100}}} \end{equation*}

gives a spiral \( u: \Omega_{150} \to [0,1] \) when evaluated using the finite difference method. Now, we want to change the domain of \( u \). We can rescale to \( \Omega _L \) so that we have an equation for \( \tilde u , \tilde v : \Omega_L \to [0,1] \) where

\begin{align*} \tilde u (t,x) &= u \left( t , {{150} \over {L}} x \right) \\ \tilde v (t,x) &= v \left( t , {{150} \over {L}} x \right) \end{align*}

which are solutions to the rescaled Barkley equation

\begin{align*} \partial_t \tilde u- { {L^2} \over {150^2}} \Delta \tilde u &=f \left( \tilde u, \tilde v \right) \\ \partial_t \tilde v &= g \left( \tilde u,\tilde v \right) \end{align*}

where we have chosen the appropriate time derivative and Laplace-Beltrami operator for our problem and recalled that \( \mathbf{v} \equiv 0\) as the surface is stationary.

Stability of the 2D scheme

Ultimately, we wish to simulate spiral waves on a unit sphere. The surface area of a unit sphere is \( 4 \pi \). We want our spherical domain and square domain to have the same area. The area of a square is \( L^2 \) so we require

\begin{equation*} L^2= 4 \pi \end{equation*}

which implies that our square domain is \( \Omega_{3.5} \) (since \(L = \sqrt{2 \pi} \approx 3.5\) ). From the calculations above, a reasonable value of the diffusion coefficient \(a \) is

\begin{align*} a &= { {4 \pi } \over { 150^2 }} \\ & \approx {1 \over {1790.49}} \end{align*}

Unfortunately this does not give us stable results, so we change the coefficient to be

\begin{align*} a &= {1 \over {179.049}}, \end{align*}

which does provide stability.

We also tested

\begin{equation*} a= {2 \over {1790.49}} , {1 \over {500}} , {1 \over {350}} , {1 \over {250}} \end{equation*}

all of which produced instability.

Finite element versus finite difference

We know that we can obtain results using both a semi-implicit method and an explicit method. We are also able to compute results using the finite differences method. In this case the numerical scheme can be again given as a \(\theta \) scheme with explicit left hand-side

\begin{equation*} u^{n+1} -a \tau \theta {\Delta_ \Gamma }^ h u^{n+1} =\left( 1-\theta \right) a \tau{\Delta_ \Gamma }^ h u^n + u^n + \tau f \left(u^n , v^n \right) \end{equation*}

where \( u^n \) is defined on a finite spatial grid \( u^n =u^n (i,j,k) \) and thus

\begin{equation*} {\Delta_ \Gamma }^ h u^n (i,j,k ) = \partial_x^2 u^n (i ,j,k) + \partial_y^2 u^n (i ,j,k) \end{equation*}


\begin{equation*} \partial_x^2 u^n (i ,j,k) ={{u^n (i+1,j,k ) -4 u^n(i,j,k) + u^n(i-1,j,k ) }\over {h^2}} \end{equation*}

using the forward diference and backward difference

\begin{align*} \partial_{x+} u^n (i,j,k ) & = {{ u^n(i+1,j,k) - u^n (i,j,k) } \over {h_x}}\\ \partial_{x-} u^n (i,j,k ) & = {{ u^n(i,j,k) - u^n (i-1,j,k) } \over {h}} \end{align*}

and that

\begin{equation*} \partial_x^2 u^n (i,j,k) = \partial_{x+} \partial_{x-}u^n (i,j,k) \end{equation*}

and \( \partial_y u^n \) is defined similarly

The numerical scheme can also be given as in a semi-implicit form. If we define

\begin{equation*} u_{th} = {{v^n +b } \over c} \end{equation*}

Then our equation is:

\begin{cases} \left( 1+ \tau \epsilon \left( 1- u^n \right) \left( u_{th} -u^n \right) \right) u^{n+1}- a \tau {\Delta_ \Gamma }^ h u^{n+1} =u^n ~~~~~~~~~~~~~~~~~~~~~~~~ u_{th} \ge u^n \\ \left( 1 -\tau \epsilon \left( u_{th} -u^n \right) \right) u^{n+1} - a \tau {\Delta_ \Gamma }^ h u^{n+1} =u^n - \tau \left( u_{th} -u^n \right) u^n ~~~~~~~~~ u_{th} \le u^n \end{cases}

FD explicit methodFem explicit method

Fem explicit methodFEM semi-implicit

Tests on \(\Omega_{3.5}\) with \(a=\frac{1}{1790.49}\), time=10

Top left: using FD explicit RHS method.

Top right: using FEM explicit RHS method.

Bottom left: using FD semi-implicit RHS method

Bottom right: using FEM semi-implicit RHS method.

Spiral on a \(\Omega_{3.5}\), \(a=1/1790.49\).

Spiral on a \(\Omega_{3.5}\), \(a=1/179.049\).