Solving the reaction-diffusion Brusselator system using Generalized Finite Difference Method

1.   Introduction

In this paper we focus on the numerical discretization of the non-linear Brusselator system, which describes oscillatory chemical behavior in autocatalytic reactions. The Brusselator equations read as:

where $x\in \Omega$, time $t > 0$, and U and V are the concentrations of the chemical species (the system is also valid to model biological species). Throughout this paper $\Omega\subset \mathbb{R}^2$ is a bounded domain with a regular boundary.

For its importance in biological or chemical processes we consider solving system (1.1) numerically using the meshless method of Generalized Finite Differences. Recently, several numerical methods have been used for solving reaction-diffusion Brusselator system. For instance, in ^[1] authors used the dual-reciprocity boundary element method (DRBEM) for the numerical solution. Also, a method of lines implemented as a meshless method based on spatial trial spaces spanned by the Newton basis functions have been applied in ^[2]. A non-standard finite difference method has been applied in ^[3], weighted residual methods have been used in ^[4] and the polynomial based differential quadrature method has been considered in ^[5] for solving the Brusselator system.

By the flexibility of our method due to the easy choice of the clouds of points and the steps of time and space, we derive a discretization of the system by means of the Generalized Finite Difference Method (GFDM) and we prove the convergence of the discrete solution to the analytical one. Several numerical examples on the applications of this meshless method over regular and irregular domains are presented in order to illustrate our results. Very recently, a high order accuracy of time discretization technique has been combined with the generalized finite difference method (see ^[6] and ^[7]).

The paper is organized as follows: In Section 2 we recall the analytical foundations of the GFDM and obtain the explicit formulas of the spatial derivatives. Section 3 is devoted to the analytical study of the explicit GFD scheme where we prove the main result of this section, enclosed in Theorem 3.1. In Section 4 we present several numerical examples over regular and irregular domains which show the applicability of the method. Finally, some conclusions are obtained in Section 5.

2.   GFDM: Explicit formulae

The aim of this section is to obtain explicit linear expressions for the approximation of partial derivatives in the points (nodes) of the domain. First of all, an irregular grid or cloud of nodes is generated in the domain. For each one of the nodes of the domain, a star is defined as a set of selected nodes $\textbf{x}_0, \textbf{x}_1, \cdots, \textbf{x}_s$, where $\textbf{x}_0$ is denominated central node of star. In order to select the other nodes of star different criteria can be used, see ^[8,9].

Let $\textbf{x}_0 = (x_0, y_0)$ be the central node and $h_i = x_i-x_0, k_i = y_i-y_0$, where $(x_i, y_i)$ are the coordinates of the $i^{th}$ node of the star. Let us put $U_0 = U(\textbf{x}_0)$ and $U_i = U(\textbf{x}_i),$ then by the Taylor series expansion, for $i = 1, ..., s$, we have

Let us call

and the spatial derivatives

If in (2.1) we don't consider the higher than second-order terms, we may obtain a second-order approximation of $U_i$, which we shall denote by $u_i$. We may then define:

where $\varphi_i = \varphi(h_i, k_i)$ are positive symmetrical weighting functions which decrease in magnitude as the distance to the center increases. Some weighting functions such as potentials or exponentials can be used (see ^[10] for more details).

In order to minimize the error, the norm function given by (2.4) is minimized with respect to the partial derivatives and the following linear system is obtained:

By solving (2.5), the partial derivatives can be obtained as a function of the values $(u_0, u_i, h_i, k_i, \varphi_i)$.

Remark 2.1. Matrix $\boldsymbol{A}$ defined in (2.5) is a positive definite matrix (see ^[10] for a complete proof) and the approximation of the spatial derivatives is of second order $\mathcal{O}(h_i^2, k_i^2)$.

If in accordance with ^[10] we define

then

and (2.7) can be rewritten as

or

where $\boldsymbol{1} = \{1, 1, \cdots, 1\}^T$. Thus, spatial derivatives using GFD as in ^[11] are approximated by

Finally, time derivative is approximated as follows

3.   GFD scheme and convergence

Let us consider a bounded domain $\Omega\subset\mathbb{R}^2$. Then, the GFD scheme for system (1.1) is:

Thus,

In order to prove the main result of the paper concerning the conditional convergence of the GFD scheme for solving system (1.1), we need the following basic results:

Lemma 3.1. Let $A\in\mathfrak{M}_{n\times n}(\mathbb{R})$. If there exists some matrix norm such that $\|A\| < 1,$ then

Lemma 3.2. Assuming $A\in\mathfrak{M}_{n\times n}(\mathbb{R})$, then the following are equivalent:

(i) $\lim\limits_{k\to\infty}A^k = \textbf{0},$

(ii) $\rho(A) < 1$,

where $\rho(\cdot)$ stands for the spectral radius.

Our main result with respect to the proposed numerical scheme is as stated below

Theorem 3.1. Let $(U, V)$ be the exact solution to system (1.1). Then, the GFD explicit scheme given by (3.2) is convergent if the following holds:

Proof.

Notice that, since $U, V$ are the exact solution of system (1.1), they also solve the discrete equation (3.2). Let us denote by $U^n_i$ the exact $U$-solution at time $n$ and node $i$ (respetively $V^n_i$). We take the difference between (3.2) and the same expression for the exact solution. We call $\tilde{u}^n_i = u^n_i-U^n_i$ (similarly for $\tilde{v}^n_i$) and notice the following relations:

Then,

Let us call $\tilde{u} = \max_{i}\{|\tilde{u}^n_i|\}$ (similarly for $\tilde{v}$). Then, we write (3.5) as

Let us rewrite (3.6) as

where

Consider the matrix

Let us consider the matrix norm $N_{1}(\mathfrak{A}) = \max_{i = 1, 2}\{\sum_{j = 1}^2|a_{ij}|\}$. From their definition it is clear that

Then,

inequality (3.10) is equivalent to

and inequalities (3.11) are equivalent to

Applying Lemma 3.1, we have that $\lim_{k\to\infty}\mathfrak{A}^k = \textbf{0}$. Now, by Lemma 3.2, this is equivalent to $\rho(\mathfrak{A}) < 1,$ that is, the greatest absolute value of all eigenvalues of matrix $\mathfrak{A}$ is bounded by 1, which implies the convergence of the explicit scheme under the condition (3.3). □

4.   Results

In this section we illustrate the application of the GFDM for solving the rection-diffusion Brusselator system given by (1.1). We test the method using the regular and irregular clouds of points of Figure 1.

In all the cases considered, the errors has been calculated according to the norms

where $sol(i)$ is the GFD solution in node $i$, $exac(i)$ is the exact value of the solution at node $i$, and $NI$ is the number of nodes of the domain $\Omega$.

For all the numerical results shown in the following sections, we choose $\Delta t$ as stated, which is to say, the condition of convergence is chosen as the minimum value computated using formulaes Eq (3.3) for each star of the domain $\Omega$. Notice that, therefore, each $\Delta t$ may vary from one cloud of points to another since its value depends on the star of nodes. We choose $\Delta t = 0.001$, fulfilling the assumption made in Theorem 3.1. The distance criterion has been used, the number of nodes per star is 8 plus the central node and the weighting function is the inverse of the distance squared, $\varphi(d_i) = \frac{1}{d_i^2}$ with $d_i$ being the distance from node $i$ to the central node of the star.

We present two cases with different boundary conditions: Dirichlet and Neumann.

Case 1: Dirichlet boundary conditions

We consider for this first case the clouds of points 1 and 2 of Figure 1.

By a direct check, the exact solution is (the Dirichlet boundary conditions are chosen so that the equation is fulfilled):

Table 1 shows the error norms in the clouds of points.

Figures 2 and 3 show the plots of the analytical, $(U, V)$, and approximate, $(u, v)$, solutions of case 1 in the clouds of points 1 and 2.

Convergence test

In order to perform a convergence test, we use the clouds of points of Figure 4: 5 (347 nodes), 6 (491 nodes) and 7 (573 nodes). The values of the $l_2(U)$ and $l_2(V)$ errors, and $\Delta t$ are collected in Table 2.

We plot in Figure 5 the $l_2(U)$ and $l_2(V)$ error norms, respectively, for $T = 1\, s$, $\Delta t = 0.001\, s$ and three clouds of points in Figure 4.

Case 2: Neumann boundary conditions

We consider for this first case the clouds of points 3 and 4 of Figure 1.

By a direct check, the exact solution is (the Neumann boundary conditions are chosen so that the equation is fulfilled):

Table 3 shows the error norms in the clouds of points.

Figures 4 and 5 show the plots of the analytical, $(U, V)$, and approximate, $(u, v)$, solutions of case 2 in the clouds of points 3 and 4.

5.   Conclusions

We have studied the convergence and obtained the discretization of system (1.1) using a Generalized Finite Difference Method explicit scheme and the order of convergence of the method is shown. The conditional convergence has been obtained for the numerical scheme. We have presented several numerical cases over different domains. The main advantages of the method are its flexibility in discretizing irregular and complex domains and the simplicity in calculations, saving time and computational resources. As a future line of research, there is the possibility of implementing space-time methods, such as STCM, and comparing them with the current method.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence tools in the creation of this article.

Conflict of interest

Prof. Dr. Francisco Ureña is the Guest Editor of special issue “Applications of Partial Differential Equations to Science and Engineering Problems: Numerical Resolution” for AIMS Mathematics. Prof. Dr. Francisco Ureña was not involved in the editorial review and the decision to publish this article.

The authors declare that they have no conflicts of interest.