A variable nonlinear splitting algorithm for reaction-diffusion systems with self- and cross-diffusion

TL;DR

This paper introduces a novel nonlinear operator splitting algorithm that efficiently models self- and cross-diffusion effects in reaction-diffusion biological systems, ensuring accuracy and stability through rigorous analysis and numerical validation.

Contribution

It presents a new splitting method that directly incorporates both self- and cross-diffusion, improving computational efficiency and stability in biological reaction-diffusion models.

Findings

The method is numerically stable under specified criteria.

The algorithm demonstrates high accuracy in numerical experiments.

It effectively models complex biological diffusion processes.

Abstract

Self- and cross-diffusion are important nonlinear spatial derivative terms that are included into biological models of predator-prey interactions. Self-diffusion models overcrowding effects, while cross-diffusion incorporates the response of one species in light of the concentration of another. In this paper, a novel nonlinear operator splitting method is presented that directly incorporates both self- and cross-diffusion into a computational efficient design. The numerical analysis guarantees the accuracy and demonstrates appropriate criteria for stability. Numerical experiments display its efficiency and accuracy.

Tables

Table 1. Table 1: A absolute maximum error at t = 1 𝑡 1 \displaystyle t=1 as we increase the resolution of the temporal step. Each time that the temporal step is halved we see that the max error is reduced by a fourth, which indicates a second-order convergence of the operator splitting scheme.

$\tau$	$\delta$	Max Error (Dirichlet)	Max Error (Neumman)
$2.5000\times {10}^{-5}$	$1.0000\times {10}^{-2}$	$1.8804\times {10}^{-8}$	$4.4342\times {10}^{-8}$
$1.2500\times {10}^{-5}$	$1.0000\times {10}^{-2}$	$4.6915\times {10}^{-9}$	$1.0839\times {10}^{-9}$
$6.2500\times {10}^{-6}$	$1.0000\times {10}^{-2}$	$1.1694\times {10}^{-9}$	$2.4808\times {10}^{-9}$

Equations

$$u_t-\Delta \left(d_{1}u+s_1{u}^2+c_{12}vu\right)$$

$$v_t-\Delta \left(d_{2}v+s_2{v}^2+c_{21}uv\right)$$

$$u(0,\text{x})=u_0(\text{x}),v(0,\text{x})=v_0(\text{x})\text{x}\in \Omega ,$$

$$=$$

$$\frac{d\text{u}}{dt}$$

$$\frac{d\text{v}}{dt}$$

$$\left(I-\frac{{\tau }_k}{2}\left((P+R)(d_1+c_{12}{D}_{k+1}^{(\text{v})})\right)\right){\text{u}}_{k+1}$$

$$=\left(I+\frac{{\tau }_k}{2}\left((P+R)(d_1+c_{12}{D}_k^{(\text{v})})\right)\right){\text{u}}_k+\frac{{\tau }_k}{2}\left({\text{f}}_{k+1}+{\text{f}}_k\right)$$

$$\left(I-\frac{{\tau }_k}{2}\left((P+R)(d_2+c_{21}{D}_{k+1}^{(\text{u})})\right)\right){\text{v}}_{k+1}$$

$$=\left(I+\frac{{\tau }_k}{2}\left((P+R)(d_2+c_{21}{D}_k^{(\text{u})})\right)\right){\text{v}}_k+\frac{{\tau }_k}{2}\left({\text{g}}_{k+1}+{\text{g}}_k\right)$$

$$\left(I-\frac{d_1{\tau }_k}{2}P\right){\tilde{\text{u}}}^{(1)}$$

$$\left(I-\frac{d_2{\tau }_k}{2}P\right){\tilde{\text{v}}}^{(1)}$$

$$\left(I-\frac{d_1{\tau }_k}{2}R\right){\tilde{\text{u}}}^{(2)}$$

$$\left(I-\frac{d_2{\tau }_k}{2}R\right){\tilde{\text{v}}}^{(2)}$$

$$\left(I-\frac{c_{12}{\tau }_k}{2}P{D}_{k+1}^{(\text{v})}\right){\tilde{\text{u}}}^{(3)}$$

$$\left(I-\frac{c_{21}{\tau }_k}{2}P{D}_{k+1}^{(\text{u})}\right){\tilde{\text{v}}}^{(3)}$$

$$\left(I-\frac{c_{12}{\tau }_k}{2}R{D}_{k+1}^{(\text{v})}\right){\text{u}}_{k+1}$$

$$\left(I-\frac{c_{21}{\tau }_k}{2}R{D}_{k+1}^{(\text{u})}\right){\text{v}}_{k+1}$$

$${\text{u}}_{k+1}$$

$${\text{v}}_{k+1}$$

$$\left(I-\frac{d_1{\tau }_k}{2}P\right)\left(I-\frac{d_1{\tau }_k}{2}R\right)\left(I-\frac{c_{12}{\tau }_k}{2}P{D}_{k+1}^{(\text{v})}\right)\left(I-\frac{c_{12}{\tau }_k}{2}R{D}_{k+1}^{(\text{v})}\right){\text{u}}_{k+1}$$

$$=\left(I+\frac{d_1{\tau }_k}{2}P\right)\left(I+\frac{d_1{\tau }_k}{2}R\right)\left(I+\frac{c_{12}{\tau }_k}{2}P{D}_k^{(\text{v})}\right)\left(I+\frac{c_{12}{\tau }_k}{2}R{D}_k^{(\text{v})}\right){\text{u}}_k$$

$$+\frac{{\tau }_k}{2}({\text{f}}_{k+1}+{\text{f}}_k),$$

$$\left(I-\frac{d_2{\tau }_k}{2}P\right)\left(I-\frac{d_2{\tau }_k}{2}R\right)\left(I-\frac{c_{21}{\tau }_k}{2}P{D}_{k+1}^{(\text{u})}\right)\left(I-\frac{c_{21}{\tau }_k}{2}R{D}_{k+1}^{(\text{u})}\right){\text{v}}_{k+1}$$

$$=\left(I+\frac{d_2{\tau }_k}{2}P\right)\left(I+\frac{d_2{\tau }_k}{2}R\right)\left(I+\frac{c_{21}{\tau }_k}{2}P{D}_k^{(\text{u})}\right)\left(I+\frac{c_{21}{\tau }_k}{2}R{D}_k^{(\text{u})}\right){\text{v}}_k$$

$$+\frac{{\tau }_k}{2}({\text{g}}_{k+1}+{\text{g}}_k),$$

$$I-\frac{\tau }{2}\left(d_1(P+R)+c_{12}(P+R)D_{k+1}^{(\text{v})}\right)+\frac{{\tau }^2}{4}(d_1^{2}PR+(c_{12}^{2}P{D}_{k+1}^{(\text{v})}R$$

$$+d_1{c}_{12}{P}^2+d_1{c}_{12}PR+d_1{c}_{12}RP+d_1{c}_{12}{R}^2)D_{k+1}^{(\text{v})})+\mathcal{O}({\tau }_k^3)$$

$$=I+\frac{\tau }{2}\left(d_1(P+R)+c_{12}(P+R)D_k^{(\text{v})}\right)+\frac{{\tau }^2}{4}(d_1^{2}PR+(c_{12}^{2}P{D}_k^{(\text{v})}R$$

$$+d_1{c}_{12}{P}^2+d_1{c}_{12}PR+d_1{c}_{12}RP+d_1{c}_{12}{R}^2)D_k^{(\text{v})})+\mathcal{O}({\tau }_k^3).$$

$$\left(I-\frac{d_1{\tau }_k}{2}P\right)\left(I-\frac{d_1{\tau }_k}{2}R\right)\left(I-\frac{c_{12}{\tau }_k}{2}P{D}_{k+1}^{(\text{v})}\right).$$

$$\left(I-\frac{d_1{\tau }_k}{2}P\right)\left(I-\frac{d_1{\tau }_k}{2}R\right)\left(I-\frac{c_{12}{\tau }_k}{2}P{D}_{k+1}^{(\text{v})}\right)\left(I-\frac{c_{12}{\tau }_k}{2}R{D}_{k+1}^{(\text{v})}\right){\text{u}}_{k+1}$$

$$=\left(I+\frac{{\tau }_k}{2}\left(d_1(P+R)+c_{12}(P+R)D_k^{(\text{v})}\right)\right){\text{u}}_k+\frac{{\tau }_k}{2}({\text{f}}_{k+1}+{\text{f}}_k)$$

$$+\frac{{\tau }_k^2}{4}\left((RP+R^2+P{D}_k^{(\text{v})}R+P^2+PR)D_k^{(\text{v})}+PR\right){\text{u}}_k+\mathcal{O}({\tau }_k^3),$$

$$\frac{d\text{u}}{dt}$$

$$\frac{d\text{v}}{dt}$$

$$\left(I-\frac{{\tau }_k}{2}\left((P+R)(d_1+s_1{D}_{k+1}^{(\text{u})}+c_{12}{D}_{k+1}^{(\text{v})})\right)\right){\text{u}}_{k+1}$$

$$=\left(I+\frac{{\tau }_k}{2}\left((P+R)(d_1+s_1{D}_k^{(\text{u})}+c_{12}{D}_k^{(\text{v})})\right)\right){\text{u}}_k+\frac{{\tau }_k}{2}\left({\text{f}}_{k+1}+{\text{f}}_k\right),$$

$$\left(I-\frac{{\tau }_k}{2}\left((P+R)(d_2+s_2{D}_{k+1}^{(\text{v})}+c_{21}{D}_{k+1}^{(\text{u})})\right)\right){\text{v}}_{k+1}$$

$$=\left(I+\frac{{\tau }_k}{2}\left((P+R)(d_2+s_2{D}_k^{(\text{v})}+c_{21}{D}_k^{(\text{u})})\right)\right){\text{v}}_k+\frac{{\tau }_k}{2}\left({\text{g}}_{k+1}+{\text{g}}_k\right).$$

$$\left(I-\frac{d_1{\tau }_k}{2}P\right){\tilde{\text{u}}}^{(1)}-{\tau }_k{\text{f}}_k$$

$$\left(I-\frac{d_2{\tau }_k}{2}P\right){\tilde{\text{v}}}^{(1)}-{\tau }_k{\text{g}}_k$$

$$\left(I-\frac{d_1{\tau }_k}{2}R\right){\tilde{\text{u}}}^{(2)}$$

$$\left(I-\frac{d_2{\tau }_k}{2}R\right){\tilde{\text{v}}}^{(2)}$$

$$\left(I-\frac{s_1{\tau }_k}{2}P{D}_{k+1}^{(\text{u})}\right){\tilde{\text{u}}}^{(3)}$$

Full text

A variable nonlinear splitting algorithm for reaction diffusion systems with self- and cross-diffusion

Matthew A. Beauregard111Principal and corresponding author. Email address: beauregama@sfasu.edu. This author was supported in part by internal research grant (No. 150030-26423-150) from Stephen F. Austin State University., Joshua L. Padgett2

1Department of Mathematics $\displaystyle\&$ Statistics, Stephen F. Austin State University, Nacogdoches, TX, 75962

2Department of Mathematics $\displaystyle\&$ Statistics, Texas Tech University, Lubbock, TX, 79409-1042

**Abstract.

**Self- and cross-diffusion are important nonlinear spatial derivative terms that are included into biological models of predator-prey interactions. Self-diffusion models overcrowding effects, while cross-diffusion incorporates the response of one species in light of the concentration of another. In this paper, a novel nonlinear operator splitting method is presented that directly incorporates both self- and cross-diffusion into a computational efficient design. The numerical analysis guarantees the accuracy and demonstrates appropriate criteria for stability. Numerical experiments display its efficiency and accuracy.

Keywords: reaction-diffusion equations, nonlinear operator splitting, self-diffusion, cross-diffusion

1 Introduction

This paper is motivated to develop efficient and accurate numerical approximations to a generalized Shigesada-Kawasaki-Teramoto (SKT) model, which is a two-species predator-prey model first developed in [22] and has been subject to much research ever since. Here, we consider the following self- and cross-diffusion system of equations that models the interaction between the prey ($\displaystyle u$) and predator ($\displaystyle v$),

[TABLE]

defined on $\displaystyle\mathbb{R}^{+}\times\Omega$. Here $\displaystyle\Omega=[0,L]\times[0,L]$, $\displaystyle\Delta$ is the two-dimensional Laplacian operator, and the reactive functions $\displaystyle f$ and $\displaystyle g$ are in $\displaystyle C^{1}(\mathbb{R}^{2})$. We define x to be the spatial coordinate vector in two dimensions. The initial populations are given as

[TABLE]

and are assumed to be nonnegative and in $\displaystyle C^{1}(\Omega)$. The diffusion coefficients are given by $\displaystyle d_{i}>0$. The parameters $\displaystyle c_{12}$ and $\displaystyle c_{21}$ describe the cross-diffusion response, with the parameter $\displaystyle c_{12}$ being nonnegative. This parameter choice encourages a prey to move away from high concentrations of the predator. A similar description is given for $\displaystyle c_{21}>0$. In the event that $\displaystyle c_{21}<0$ the species $\displaystyle v$ move toward higher concentrations of $\displaystyle u$, its primary food source. However, in this case, finite time blow-up of the solution may occur [28], hence, we assume the $\displaystyle c_{21}>0$. The nonnegative self-diffusion coefficients $\displaystyle s_{i}$ model interspecies competition and provides a severe penalty to over-crowding of a particular species.

Since each of the coefficients are assumed nonnegative then it is known that weak solutions exist globally in time when considering homogenous Neumann boundary conditions and reactive functions of Lotka-Volterra type that satisfy $\displaystyle f(0,v)=g(u,0)=0$ [6, 7, 8]. Likewise, classical solutions globally exist under various restrictions on the self and cross diffusion coefficients [20], namely, $\displaystyle s_{1}s_{2}\geq 0,~{}s_{2}>c_{12},~{}s_{1}>c_{21}.$ However, if $\displaystyle d_{1}<d_{2}$ then $\displaystyle s_{1}$ maybe less than $\displaystyle c_{21}$ and solutions will still exist globally in time [18]. In light of these theoretical results, this paper only considers homogeneous Neumann boundary conditions, even though our numerical analysis extends to homogeneous Dirichlet or Neumann-Dirichlet boundary conditions in a straightforward manner.

The inclusion of self- and cross-diffusion terms allow for realistic responses to predator and prey movement and are often incorporated into mathematical models in population biology [4, 11, 26, 27]. Our current efforts consider only a two species model, yet the methods developed herein can by readily extended to more general systems. The considered model is rich in dynamics and stems from the generalizations of the SKT model and the literature is abundant with investigations of variants to this model [1, 7, 8, 9, 17, 18, 19, 23, 25]. However, the development and analysis of accurate and efficient numerical approximations has not been considered.

This paper is organized as follows. In Section 2 we present the nonlinear variable time splitting method where $\displaystyle s_{i}=0$. The algorithm is based on a variant of a splitting method recently developed by the authors for similar equations with self-diffusion in [2]. Section 3 develops the nonlinear splitting algorithm for systems with both self- and cross-diffusion. In both of these sections the methods are overall second-order accurate and their equivalent factorizations are preferable as they allow for parallel implementations and use of high performance computing (HPC) methods. Section 4 details further numerical analysis to provide criteria that guarantee stability. Section 5 provides examples to illustrate the efficiency and accuracy of the algorithm. We conclude and highlight our results in Section 6.

In the following, we define a scheme as computationally efficient if it is second-order accurate in space and time or better and the number of operations per time step is directly proportional to the number of unknowns. In addition, all lowercase bold letters indicate vectors, uppercase letters are used for matrices. Lastly, the $\displaystyle\ell^{2}$-norm is used throughout unless otherwise specified.

2 Splitting algorithm with cross-diffusion only

Given $\displaystyle N\gg 0,$ we may inscribe over $\displaystyle\Omega$ the mesh $\displaystyle\mathcal{D}_{\delta}=\{(x_{i},y_{j})~{}|~{}i,j=0,1,\dots,N+1\},$ where $\displaystyle\delta=L/(N+1)$ and $\displaystyle x_{i}=i\delta$ and $\displaystyle y_{i}=j\delta$ for $\displaystyle i,j=0,1,\dots,N+1.$ Further, define $\displaystyle u_{i,j}(t)$ and $\displaystyle v_{i,j}(t)$ as the approximations to the exact solution $\displaystyle u(x_{i},y_{j},t)$ and $\displaystyle v(x_{i},y_{j},t)$, respectively. Let

[TABLE]

A similar notation is used for the reactive functions $\displaystyle f$ and $\displaystyle g$. Second-order central difference approximations are used to approximate the two-dimensional Laplacian operator in Eqs. (1.1)-(1.2) which yields the semidiscretized equations,

[TABLE]

where $\displaystyle P$, $\displaystyle R$, and $\displaystyle D(\textbf{q})=\mbox{diag}(\textbf{q})$ are $\displaystyle M^{2}\times M^{2}$ matrices, where $\displaystyle M=N+2$. The matrices are defined as $\displaystyle P\mathrel{\mathop{:}}=I_{M}\otimes T$ and $\displaystyle R\mathrel{\mathop{:}}=T\otimes I_{M},$ where $\displaystyle I_{M}$ is the $\displaystyle M\times M$ identity matrix and $\displaystyle T$ is a $\displaystyle M\times M$ tridiagonal matrix with main diagonal $\displaystyle-2/\delta^{2}$ and upper and lower diagonals of $\displaystyle 1/\delta^{2}$ with the exception of $\displaystyle T_{12}=T_{M,M-1}=2$ as a result of the Neumann boundary conditions. The semidiscretized equations are specified at every mesh point, with the Neumann boundary conditions being incorporated via central difference approximations.

A variable time-step second-order Crank-Nicolson method [2] is used to advance the solution in time, that is,

[TABLE]

where $\displaystyle\textbf{u}_{k}$ and $\displaystyle\textbf{v}_{k}$ are approximations to u and v at time $\displaystyle t_{k}=\textstyle\sum_{i=0}^{k-1}\tau_{k}$ and $\displaystyle D_{k}^{(\textbf{q})}=D(\textbf{q}_{k})$ are diagonal matrices whose elements are $\displaystyle\textbf{q}=\textbf{u},\,\textbf{v}.$ Similarly, the reactive functions $\displaystyle\textbf{f}_{k}$ and $\displaystyle\textbf{g}_{k}$ denotes the estimated values of f and g using the approximations of $\displaystyle\textbf{u}_{k}$ and $\displaystyle\textbf{v}_{k}$. To advance the solution in time we propose to solve Eqs. (2.1)-(2.2) through a modified Douglass-Gunn splitting method, similar to that shown in [2], given by

[TABLE]

Note that $\displaystyle(\ref{splitstepa})$-$\displaystyle(\ref{SplitStepFinal})$ involve implicit terms. In order to maintain the desired computational efficiency, the implicit terms are approximated by taking an Euler step, that is,

[TABLE]

By employing the Euler step, we have developed a second-order accurate scheme which remains computationally efficient using standard ADI techinques [3]. The Euler approximation is established in the first two steps (2.3a) and (2.3b). Hence, the overall computational cost is still $\displaystyle\mathcal{O}(N^{2})$. In effect, the splitting method resolves the nonlinearity directly and advances the solution by iterating through one-dimensional subproblems. Careful implementation only requires the use of six $\displaystyle N\times 1$ vectors. The splitting method’s solution is equivalent to the solution of Eqs. (2.1)-(2.2). We demonstrate this equivalence in the following. First, we show that Eqs. (2.1)-(2.2) are equivalent to a useful factorization.

Theorem 2.1**.**

The induced error between the factorizations

[TABLE]

and Eqs. (2.1)-(2.2), respectively, is $\displaystyle\mathcal{O}(\tau_{k}^{3})$.

Proof.

We show the equivalency between Eq. (2.1) and Eq. (2.1). A virtually identical proof can be used to establish the result for v. For ease of exposition, we disregard the reaction function approximations $\displaystyle\textbf{f}_{k}$ and $\displaystyle\textbf{f}_{k+1}.$ Expanding the matrix products on the left- and right-hand sides of Eq. (2.1) yields,

[TABLE]

The desired result is obtained by expanding the $\displaystyle\mathcal{O}(\tau_{k}^{2})$ terms about $\displaystyle\textbf{v}_{k}$ in the left-hand side: that is $\displaystyle D_{k+1}=D_{k}+\mathcal{O}\left(\tau_{k}\right)$. The truncated terms contribute terms of order $\displaystyle\mathcal{O}(\tau^{3})$. Therefore the $\displaystyle\mathcal{O}(\tau^{2})$ terms on both sides cancel. Consequently the factorization is equivalent up to $\displaystyle\mathcal{O}(\tau_{k}^{3})$. ∎

The factorization in Eqs. $\displaystyle(\ref{ufactoredNOADI})$ and $\displaystyle(\ref{ufactoredNNOADI})$ is desirable as it leads to our splitting method given in Eqs. (2.3a)-(2.3h). This is encapsulated in the following theorem.

Theorem 2.2**.**

The induced error between the factorizations given by Eqs. $\displaystyle(\ref{ufactoredNOADI})$ and $\displaystyle(\ref{ufactoredNNOADI})$ and the splitting design of Eqs. (2.3a)-(2.3h) is $\displaystyle\mathcal{O}(\tau_{k}^{3})$.

Proof.

This is a straightforward calculation accomplished by multiplying (2.3g) by the matrix product:

[TABLE]

Expanding and substituting into the result the definitions for $\displaystyle\tilde{\textbf{u}}^{(1)},\tilde{\textbf{u}}^{(2)},$ and $\displaystyle\tilde{\textbf{u}}^{(3)}$ results in

[TABLE]

where we have used the expansions $\displaystyle\textbf{v}_{k+1}=\textbf{v}_{k}+\mathcal{O}\left(\tau_{k}\right)$ and $\displaystyle\textbf{f}_{k+1}=\textbf{f}_{k}+\mathcal{O}(\tau_{k})$. The $\displaystyle\textbf{u}_{k}$ terms on the right-hand side can be factored to yield the desired result. A similar proof can be established for the factorization used to advance v. ∎

The combination of these two theorems establishes that the splitting method is second-order accurate in time. Hence, the overall accuracy of the splitting method combined with the spatial discretization is second-order accurate.

3 Splitting algorithm with self and cross-diffusion

Here, we adopt the same notation as in the previous section and consider the semidiscretization of Eqs. (1.1)-(1.2),

[TABLE]

As before, we advance solution in time using a variable time-step second-order Crank-Nicolson method,

[TABLE]

To advance the solution, just as before, a modified Douglass-Gunn splitting method is utilized to solve Eqs. (3.1)-(3.2). The proposed method is given by

[TABLE]

An Euler step is used to approximate the implicit terms and maintain the computational efficiency. The splitting method’s solution is equivalent to the solution of Eqs. (3.1)-(3.2). The following theorems develop an equivalent factorization of the Eqs. (3.1)-(3.2) and show that this factorization is equivalent to the splitting in Eqs. (3.3a)-(3.3l). The proofs are omitted since they follow a similar pattern to those of Theorems 2.1 and 2.2.

In the following, for brevity, we employ the notations

[TABLE]

and

[TABLE]

Theorem 3.1**.**

The induced error between the factorizations

[TABLE]

and Eqs. (3.1)-(3.2), respectively, is $\displaystyle\mathcal{O}(\tau_{k}^{3})$.

Theorem 3.2**.**

The induced error between the factorizations given by Eqs $\displaystyle(\ref{sc1})$-$\displaystyle(\ref{sc2})$ and the splitting design of Eqs. (3.3a)-(3.3l) is $\displaystyle\mathcal{O}(\tau_{k}^{3})$.

4 Stability and Convergence

We now establish criteria to guarantee the invertibility of specific matrices of the proposed operator splitting method.

Theorem 4.1**.**

Let $\displaystyle\kappa=\max\{d_{1},d_{2},s_{1},s_{2},c_{12},c_{21}\}$. If $\displaystyle\tau_{k}$ is sufficiently small,

[TABLE]

and $\displaystyle\emph{{u}}_{k},\emph{{v}}_{k}\geq 0$ , then, for $\displaystyle i=1$ and $\displaystyle 2$, the matrices

[TABLE]

[TABLE]

are nonsingular.

Proof.

In the case of homogeneous Dirichlet or Neumann boundary conditions we have $\displaystyle\|T\|\leq 4/\delta^{2}.$ Therefore,

[TABLE]

which gives that the matrix $\displaystyle\textstyle I-\frac{d_{1}\tau_{k}}{2}P$ is nonsingular. A similar argument can be established for the matrices $\displaystyle\textstyle I-\frac{d_{2}\tau_{k}}{2}P$ and $\displaystyle\textstyle I-\frac{d_{i}\tau_{k}}{2}R$ for $\displaystyle i=1$ and $\displaystyle 2$.

For the remaining matrices we substitute Taylor expansions of $\displaystyle\textbf{u}_{k+1}$ and $\displaystyle\textbf{v}_{k+1}$ into the matrices and then show that the resulting matrices satisfies the weak-row sum criterion [13] (up to the appropriate order), that is:

The nonzero off-diagonal entries share the same sign and are of opposite sign of the nonzero main diagonal entries and; 2. 2.

The row sums of the matrix are all nonnegative with at least one row sum strictly greater than zero.

Suppose $\displaystyle c_{21}>0$ and consider the matrix $\displaystyle\textstyle I-\frac{c_{21}\tau_{k}}{2}PD_{k+1}^{(\textbf{u})}$. Taylor expanding $\displaystyle\textbf{u}_{k+1}$ shows that

[TABLE]

where

[TABLE]

Since the splitting scheme is second-order accurate we shall neglect the high-order terms in showing that $\displaystyle Q$ satisfies the weak-row sum criterion. Consider the off diagonal entry:

[TABLE]

Utilizing Taylor expansions in the spatial coordinates, we see that for $\displaystyle\tau_{k}$ sufficiently small, $\displaystyle Q_{j,j+1}\leq 0,$ since $\displaystyle\textbf{v}_{k}$ and $\displaystyle\textbf{f}_{k}$ are assumed to be nonnegative. A virtually identical proof can be established for the off-diagonal entries affected by the Neumann boundary conditions. In addition, a similar calculation shows the lower diagonal of $\displaystyle Q$ is not positive. Therefore, $\displaystyle Q_{i,j}\leq 0$ for $\displaystyle i\neq j.$

Next, we consider the row sums and show that these are nonnegative with at least one row sum being strictly positive. Consider the $\displaystyle i^{th}$ row sum of $\displaystyle Q$,

[TABLE]

for $\displaystyle i=1,\dots,M^{2}$. The boundary conditions are $\displaystyle(\textbf{u}_{k})_{-1}=(\textbf{u}_{k})_{1}$, $\displaystyle(\textbf{v}_{k})_{-1}=(\textbf{v}_{k})_{1}$, $\displaystyle(\textbf{u}_{k})_{M^{2}+1}=(\textbf{u}_{k})_{M^{2}-1}$, and $\displaystyle(\textbf{v}_{k})_{M^{2}+1}=(\textbf{v}_{k})_{M^{2}-1}$. Simplifying for $\displaystyle i\neq 1,M^{2}$ yields,

[TABLE]

Clearly, $\displaystyle\tau_{k}$ can be chosen sufficiently small such that the $\displaystyle\mathcal{O}{\left(\tau_{k}^{2}/\delta^{2}\right)}$ term can be made small enough to ensure the result is positive. A similar result can be established for $\displaystyle i=1$ and $\displaystyle M^{2}$. Hence, we have that the weak row sum criterion is satisfied. Likewise, the remaining matrices can be shown to be nonsingular and inverse positive using a virtually identical approach. ∎

We now study the stability of both Eqs. $\displaystyle(\ref{ufactoredNOADI})$-$\displaystyle(\ref{ufactoredNOADI})$ and Eqs. $\displaystyle(\ref{sc1})$-$\displaystyle(\ref{sc2})$. Moreover, we improve upon the previous methodology employed in [2] by significantly reducing the previously imposed regularity conditions on the solution. In order to prove stability of our proposed nonlinear splitting algorithm, we introduce a definition and some lemmas.

Definition 4.1**.**

Let $\displaystyle\|\cdot\|$ be an induced matrix norm. Then the associated logarithmic norm $\displaystyle\mu\,:\,\mathbb{C}^{M\times M}\to\mathbb{R}$ of $\displaystyle A\in\mathbb{C}^{M\times M}$ is defined as

[TABLE]

where $\displaystyle I\in\mathbb{C}^{M\times M}$ is the identity matrix.

Remark 4.1*.*

When the induced matrix norm being considered is the $\displaystyle\ell^{2}-$norm, then $\displaystyle\mu(A)=\max\{\lambda\,:\,\lambda\ \mbox{is an eigenvalue of}\ (A+A^{*})/2\}.$

Lemma 4.1**.**

Let $\displaystyle t\in\mathbb{C}$ and $\displaystyle A\in\mathbb{C}^{M\times M}.$ Then for any induced matrix norm $\displaystyle\|\cdot\|$ we have

[TABLE]

where $\displaystyle E(\cdot)$ is the matrix exponential.

Proof.

See [12]. ∎

Lemma 4.2**.**

Assume that $\displaystyle(\ref{cfl})$ holds and let $\displaystyle L\mathrel{\mathop{:}}=P+R$, $\displaystyle L_{u}\mathrel{\mathop{:}}=d_{1}+s_{1}D_{k}^{(\emph{{u}})}+c_{12}D_{k}^{(\emph{{v}})}$, and $\displaystyle L_{v}\mathrel{\mathop{:}}=d_{2}+s_{2}D_{k}^{(\emph{{v}})}+c_{21}D_{k}^{(\emph{{u}})}$ Then we have

[TABLE]

Proof.

By Lemma 4.1 we have

[TABLE]

thus we need to show that $\displaystyle\mu\left(LL_{u}\right)\leq 0$ to obtain the desired result. To do so, we need to show that the eigenvalues of the matrix $\displaystyle LL_{u}+L_{u}L^{*}$ are nonpositive.

Note that $\displaystyle L_{u}$ is diagonal and by $\displaystyle(\ref{cfl})$ its entries are positive, thus we have $\displaystyle LL_{u}L_{u}^{-1/2}=LL_{u}^{1/2}.$ Then,

$\displaystyle LL_{u}$ is similar to $\displaystyle L_{u}^{1/2}LL_{u}^{1/2}$ with similarity matrix $\displaystyle L_{u}^{-1/2}$ and 2. 2.

$\displaystyle L$ is congruent to $\displaystyle L_{u}^{1/2}LL_{u}^{1/2}$ through the matrix $\displaystyle L_{u}^{-1/2}$.

Therefore by Sylvester’s law of inertia we have that the eigenvalues of $\displaystyle L$ are identical to those of $\displaystyle L_{u}^{1/2}LL_{u}^{1/2}$, which due to similarity, are the same as those of $\displaystyle LL_{u}$ [15]. Clearly, eigenvalues of $\displaystyle L$ are real and nonpositive [24]. Hence, so are the eigenvalues of $\displaystyle LL_{u}$. A similar argument gives that $\displaystyle LL_{v}$ has real and nonpositive eigenvalues. Subsequently, the desired bounds are established. ∎

Remark 4.2*.*

It is worth noting that Lemma 4.2 holds independently of $\displaystyle k,$ and even $\displaystyle\textbf{u}_{k}$ and $\displaystyle\textbf{v}_{k}.$ All that is necessary is that the entries of the diagonal matrix $\displaystyle D$ are positive.

Lemma 4.3**.**

Assume that $\displaystyle(\ref{cfl})$ holds and let $\displaystyle L$, $\displaystyle L_{u}$, and $\displaystyle L_{v}$ be defined as in the previous lemma. Then

[TABLE]

Proof.

For convenience, we employ the notation $\displaystyle X_{k}^{(\textbf{u})}\mathrel{\mathop{:}}=LL_{u}$ and $\displaystyle X_{k}^{(\textbf{v})}\mathrel{\mathop{:}}=LL_{v}.$ By employing the inverse Laplace transform and the second resolvent identity [14], we have

[TABLE]

where $\displaystyle\Gamma$ is a path surrounding the spectrum of the matrices $\displaystyle X_{k}^{(\textbf{u})}$ and $\displaystyle X_{k}^{(\textbf{v})}.$ Taking the norm of the above gives

[TABLE]

We proceed by developing separate bounds for the matrices in the above expression. Note that by the Laplace transform and Lemma 4.2, we have

[TABLE]

We similarly derive the bound

[TABLE]

Further, we have by direct calculation

[TABLE]

Combining the above yields

[TABLE]

which is the desired result. ∎

Lemma 4.4**.**

Assume that $\displaystyle(\ref{cfl})$ holds, $\displaystyle u,v\in H^{1}(\Omega),$ and that

[TABLE]

and

[TABLE]

Then

[TABLE]

for some positive constant $\displaystyle C$.

Proof.

By employing the notation from Lemma 4.3, we have

[TABLE]

and

[TABLE]

Thus, it follows that

[TABLE]

where

[TABLE]

Taking the norm gives

[TABLE]

By our assumptions and [10] we have

[TABLE]

Applying $\displaystyle(\ref{dd12})$ to $\displaystyle(\ref{dd11})$ yields the desired result. ∎

In the following theorem we employ the following product notation for non-commuting matrices $\displaystyle A_{j}\in\mathbb{C}^{M\times M}:$

[TABLE]

Theorem 4.2**.**

Assume that $\displaystyle(\ref{cfl})$ holds, $\displaystyle u,v\in H^{1}(\Omega),$ and let $\displaystyle\tau_{j},\ j=0,\ldots,k,$ be sufficiently small. Then it follows that the schemes $\displaystyle(\ref{ufactoredNOADI})$-$\displaystyle(\ref{ufactoredNNOADI})$ are stable in the sense that there exists $\displaystyle K_{1},K_{2}>0$ such that

[TABLE]

where $\displaystyle\emph{{z}}_{0}^{(\emph{{u}})}=\emph{{u}}_{0}-\tilde{\emph{{u}}}_{0}$ and $\displaystyle\emph{{z}}_{0}^{(\emph{{v}})}=\emph{{v}}_{0}-\tilde{\emph{{v}}}_{0}$ are initial errors, $\displaystyle\emph{{z}}_{k+1}^{(\emph{{u}})}\mathrel{\mathop{:}}=\emph{{u}}_{k+1}-\tilde{\emph{{u}}}_{k+1}$ and $\displaystyle\emph{{z}}_{k+1}^{(\emph{{v}})}\mathrel{\mathop{:}}=\emph{{v}}_{k+1}-\tilde{\emph{{v}}}_{k+1}$ are the $\displaystyle(k+1)$th perturbed error vectors, and $\displaystyle K_{1}$ and $\displaystyle K_{2}$ are independent of $\displaystyle k$ and $\displaystyle\tau_{k}.$

Proof.

We first consider $\displaystyle(\ref{ufactoredNOADI})$ and, for the time, disregard the nonlinear reaction term. Thus, we have

[TABLE]

It will be convenient to appeal to a Magnus-type representation of $\displaystyle(\ref{mag1})$ [5, 16]. That is, by employing the standard nonlinear Magnus-type integration and exponential splitting theory, we have

[TABLE]

where $\displaystyle L\mathrel{\mathop{:}}=P+R,$ $\displaystyle D_{k+1/2}^{(\textbf{v})}=D(\textbf{v}_{k+1/2}),$ and $\displaystyle t_{k+1/2}\mathrel{\mathop{:}}=t_{k}+\tau_{k}/2.$ We now consider perturbations of the numerical solutions $\displaystyle\textbf{u}_{k}$ and $\displaystyle\textbf{v}_{k},$ represented as

[TABLE]

and

[TABLE]

Let $\displaystyle\textbf{z}^{(\textbf{u})}_{k}\mathrel{\mathop{:}}=\textbf{u}_{k}-\tilde{\textbf{u}}_{k}$ and $\displaystyle\textbf{z}^{(\textbf{v})}_{k}\mathrel{\mathop{:}}=\textbf{v}_{k}-\tilde{\textbf{v}}_{k}.$ Iterating $\displaystyle(\ref{mag1})$ by employing $\displaystyle(\ref{mag2})$ yields

[TABLE]

and iterating $\displaystyle(\ref{pert1})$ yields

[TABLE]

Subtracting $\displaystyle(\ref{uu1})$ from $\displaystyle(\ref{u1})$ gives

[TABLE]

where

[TABLE]

with

[TABLE]

We first consider $\displaystyle\Psi_{1}.$ To that end, by Lemma 4.2, we have

[TABLE]

Now we consider $\displaystyle\Psi_{2}.$ By manipulating $\displaystyle\Psi_{2},$ we arrive at the following expression

[TABLE]

Taking the norm of both sides of $\displaystyle(\ref{psi2})$ and employing Lemma 4.3 gives

[TABLE]

where $\displaystyle C$ is a positive constant. By our assumptions and Lemma 4.4, we have

[TABLE]

where $\displaystyle C_{1}$ is a constant independent of $\displaystyle\delta$ and $\displaystyle\tau_{j},\ j=0,\ldots,k.$

Combining the above gives

[TABLE]

where $\displaystyle C_{2}$ and $\displaystyle C_{3}$ are independent of $\displaystyle j,\ \tau_{j},$ $\displaystyle j=0,\ldots,k.$ By employing similar arguments with $\displaystyle(\ref{ufactoredNNOADI})$, it can be shown that $\displaystyle\|\textbf{z}_{k+1}^{(\textbf{v})}\|\leq\|\textbf{z}_{0}^{(\textbf{v})}\|+C_{1}\|\tilde{\textbf{v}}_{0}\|+C_{3}T\equiv K_{2},$ after possibly rescaling $\displaystyle C_{1}$ and $\displaystyle C_{3}.$

We now consider the nonlinear reaction term in the analysis of $\displaystyle(\ref{ufactoredNOADI})$. Recall that $\displaystyle\textbf{f}_{k}=\textbf{f}(\textbf{u}_{k},\textbf{v}_{k}).$ Due to the differentiability of $\displaystyle f,$ we have

[TABLE]

for some $\displaystyle\xi_{k}=s_{1}\textbf{u}_{k}+(1-s_{1})\tilde{\textbf{u}}_{k},\ s_{1}\in[0,1]$ and $\displaystyle\eta_{k}=s_{2}\textbf{v}_{k}+(1-s_{2})\tilde{\textbf{v}}_{k},\ s_{2}\in[0,1].$ Thus, we have

[TABLE]

Since $\displaystyle f$ is bounded, we obtain a similar stability bound as above. The same can be done for $\displaystyle(\ref{ufactoredNNOADI})$ by using the properties of the reaction term $\displaystyle g.$∎

We now present the stability theorem for the scheme $\displaystyle(\ref{sc1})$-$\displaystyle(\ref{sc2})$.

Theorem 4.3**.**

Assume that $\displaystyle(\ref{cfl})$ holds, $\displaystyle u,v\in H^{1}(\Omega),$ and let $\displaystyle\tau_{j},\ j=0,\ldots,k,$ be sufficiently small. Then it follows that the schemes $\displaystyle(\ref{sc1})$-$\displaystyle(\ref{sc2})$ are stable in the sense that there exists $\displaystyle K_{1},K_{2}>0$ such that

[TABLE]

where $\displaystyle\emph{{z}}_{0}^{(\emph{{u}})}=\emph{{u}}_{0}-\tilde{\emph{{u}}}_{0}$ and $\displaystyle\emph{{z}}_{0}^{(\emph{{v}})}=\emph{{v}}_{0}-\tilde{\emph{{v}}}_{0}$ are initial errors, $\displaystyle\emph{{z}}_{k+1}^{(\emph{{u}})}\mathrel{\mathop{:}}=\emph{{u}}_{k+1}-\tilde{\emph{{u}}}_{k+1}$ and $\displaystyle\emph{{z}}_{k+1}^{(\emph{{v}})}\mathrel{\mathop{:}}=\emph{{v}}_{k+1}-\tilde{\emph{{v}}}_{k+1}$ are the $\displaystyle(k+1)$th perturbed error vectors, and $\displaystyle K_{1}$ and $\displaystyle K_{2}$ are independent of $\displaystyle k$ and $\displaystyle\tau_{k}.$

Proof.

This result follows by a proof similar to that of the previous theorem, where we use our developed bounds on $\displaystyle u$ and $\displaystyle v$ and employ the bound $\displaystyle s_{1}D_{k}^{(\textbf{u})}+c_{12}D_{k}^{(\textbf{v})}\leq\max\{s_{1},c_{12}\}D_{k}^{(\textbf{w})},$ where $\displaystyle w_{i}=\max\{\textbf{u}_{i},\textbf{v}_{i}\}$ for $\displaystyle X_{k}^{(\textbf{u})}$. The remaining arguments follow in a similar fashion. ∎

5 Numerical Experiments

In this section we provide illustrative the efficiency and convergence rate of the developed algorithm. All of the computations are carried out on a single HP EliteDesk 800 G1 work station with an Intel®Core(TM) i7-4770 3.40 GHz processor with 16 GB of RAM using Matlab®.

5.1 Example 1.

Consider Eqs. (1.1)-(1.2) with $\displaystyle L=1,$ homogenous Dirichlet boundary conditions, and reaction functions

[TABLE]

It can be shown by direct calculation that the exact solution to this system of partial differential equations is

[TABLE]

Define $\displaystyle u_{i,j}^{(\tau)}$ and $\displaystyle v_{i,j}^{(\tau)}$ as the numerical solutions to $\displaystyle u$ and $\displaystyle v$ at time $\displaystyle T$ with a temporal step size of $\displaystyle\tau$ and location $\displaystyle(i\delta,j\delta)$. Based on the analysis is the previous section, we anticipate that the order of convergence in time is second-order. That is, $\displaystyle\max\{|u_{i,j}^{(\tau)}-u(ih,ij,T)|,|v_{i,j}^{(\tau)}-v(ih,ij,T)|\}\approx C\tau^{p}$ for some arbitrary constant $\displaystyle C$ and $\displaystyle p=2$. We approximate $\displaystyle p$ by,

[TABLE]

Let $\displaystyle\delta=.01$ and $\displaystyle d_{1}=d_{2}=s_{1}=s_{2}=c_{12}=c_{21}=1$. Therefore, based on our CFL condition, we let $\displaystyle\tau=.25\times 10^{-4}<.5\times 10^{-4}.$ We simulate the solution until $\displaystyle T=1$ and estimate the approximate order of convergence to be $\displaystyle p\approx 2.000809.$ This indicates that our nonlinear operator splitting method is converging at the desired second-order rate.

Likewise, if we consider Neumann boundary conditions with the reaction functions

[TABLE]

then the exact solution for $\displaystyle u$ and $\displaystyle v$ is $\displaystyle a+\cos(\pi x)\cos(\pi y)\exp(-\pi^{2}t)$, where $\displaystyle a\in\mathbb{R}$. In the simulations we let $\displaystyle a=1$, which ensures that the initial condition remains nonnegative throughout the computational domain. Using the same resolution as in the Dirichlet case, we determine the order of convergence to be approximately $\displaystyle 1.9945$. In Table 1, we demonstrate the data for the maximum absolute error as we increase the resolution in time. Notice that for a fixed $\displaystyle\delta,$ as we half the temporal step the max error decreases by approximately one-fourth. This is also true if we fix the temporal step size and then halve the spatial step sizes, which is the expected result for a second-order numerical method.

To test the computational efficiency we determine the computational time to complete $\displaystyle 1000$ iterations for a fixed $\displaystyle\tau$ as we increase the spatial resolution, $\displaystyle N$. We determined that the computational time scales as $\displaystyle\mathcal{O}(N^{1.7568}),$ which is an exponent less than the number of unknowns. Figure 1(b) shows the log-log plot of the computational time versus $\displaystyle N$. These results suggests that the procedure is computationally efficient.

5.2 Example 2.

Consider Eqs. (1.1)-(1.2) with Neumann boundary conditions and the following Lotka-Volterra reaction functions,

[TABLE]

In [20] it was shown that for the above Lotka-Volterra reaction function, a global solution exists if the parameters for diffusion satisfy the following conditions:

[TABLE]

Consider a situation that violates the third condition and let $\displaystyle d_{1}=.01,d_{2}=.1,s_{1}=.05,s_{2}=.4,c_{12}=.12,$ and $\displaystyle c_{21}=.06.$ These coefficients satisfy the criteria shown in [18] for global solutions. For the reaction functions’ parameters we choose $\displaystyle a_{1}=1,b_{1}=2,c_{1}=.2,a_{2}=.3,b_{2}=1,$ and $\displaystyle c_{2}=4.$ The initial conditions are chosen to be

[TABLE]

where $\displaystyle a_{i},b_{i},m_{i},$ and $\displaystyle n_{i}$ are positive integers and $\displaystyle L=\pi.$ The amplitudes $\displaystyle\sigma_{i}$ and $\displaystyle\beta_{i}$ were chosen randomly from uniform distribution on the interval $\displaystyle(0,1)$. Regardless of the choice of amplitudes or frequencies the solution converges to spatially homogenous solutions which exists globally in time.

5.3 Example 3.

Consider Eqs. (1.1)-(1.2) with Dirichlet boundary conditions and reaction functions of the form

[TABLE]

In [28] it was shown that finite time blow up can occur if $\displaystyle b_{i}>s_{i}\lambda$ and $\displaystyle a_{i}\geq d_{i}\lambda$ where $\displaystyle\lambda$ is the first eigenvalue of the Laplacian operator and $\displaystyle c_{12}=c_{21}=0$.

Consider the initial conditions

[TABLE]

with parameters $\displaystyle d_{1}=d_{2}=1,s_{1}=s_{2}=.05,a_{1}=a_{2}=3,b_{1}=b_{2}=4$. In such a situation the system of equations are decoupled. The parameters for $\displaystyle u$ and $\displaystyle v$ are chosen to obey the criteria to guarantee blow-up. In order to satisfy Eq. (4.1) the temporal step-size must be reduced as the populations grow in time. We terminate the calculation when $\displaystyle\tau$ is adapted to a minimum of $\displaystyle 10^{-10}$. Indeed, in our simulations the populations are demonstrated to grow without bound. Figure 2 shows a panel of the initial condition coupled with two snapshots of the population during the calculation. Figure 2(b) shows the population at $\displaystyle t=.5$ for which the population has begin to concentrate and grow in the center of the domain. Figure 2(c) shows the population in the iteration prior to termination. At termination of the computation, the population had reached a maximum value of approximately $\displaystyle 1.6370\times 10^{7}$, while the temporal derivative was approximately $\displaystyle 1.0116\times 10^{18}$.

6 Conclusions

In this paper a novel adaptive nonlinear operator splitting method was developed to approximate coupled parabolic partial differential equations that contain self- and cross-diffusion terms. The method extends the algorithm designed in [2] and provides an updated criterion to ensure the method’s stability and the invertibility of the involved matrices. Surprisingly, the inclusion of the cross-diffusion term does not complicate the numerical analysis any further than the inclusion of self-diffusion terms. Hence, convergence and stability follow by employing techniques similar to those developed in our prior work [2]. However, the current research improves upon these techniques and employs techniques that allow for reduced regularity requirements on the solution. Moreover, the analysis does not depend on the boundary conditions utilized and can be readily extended to several problems of particular interest in the mathematical biology community such as in [18] and more recently [21]. Our numerical experiments provide empirical evidence of the anticipated convergence rate for Dirichlet and Neumann boundary conditions, further verifying the efficacy of the proposed method. The computational experiments further provide evidence suggesting that global solutions exist under certain parameter restrictions.

Acknowledgements

The first author would like to explicitly express his gratitude for an internal research grant (No. 150030-26423-150) from Stephen F. Austin State University.

We would like to express our gratitude for the numerous suggestions provided by the reviewers that undoubtedly improved the results and readability of this manuscript.