Oncolytic viral therapies and the delicate balance between virus-macrophage-tumour interactions: A mathematical approach

1.   Introduction

This paper mainly focuses on constructing and analyzing an efficient energy-preserving finite difference method (EP-FDM) for solving the nonlinear coupled space-fractional Klein-Gordon (KG) equations:

with $(\boldsymbol{x}, t) \in \Omega \times[0, T]$ and the following widely used boundary and initial conditions

Here, $\boldsymbol{x} = (x_{1}, ..., x_{d})^{T}(d = 1, 2, 3) \in \Omega\subset R^{d}$, $\partial \Omega$ is the boundary of $\Omega$, $\bar{\Omega} = \Omega\bigcup\partial \Omega$, $\kappa$ is a constant and $a_{i}, b_{i}, c_{i}$ are all positive constants. $\phi_1, \phi_2, \varphi_1, \varphi_2$ are all known sufficiently smooth functions. $u(\boldsymbol{x}, t)$, $v(\boldsymbol{x}, t)$ are interacting relativistic fields of masses, $\partial_{x_k}^{\alpha_k}u$ and $\partial_{x_k}^{\alpha_k}v$ stand for the Riesz fractional operator with $1 < \alpha_{k}\leq 2, (k = 1, ..., d)$ in $x_{k}$ directions, which are well defined as follows

where $-\infty D_{x_k}^{\alpha_k} u(\boldsymbol{x}, t)$ and ${ }_{x_k} D_{+\infty}^{\alpha_k} u(\boldsymbol{x}, t)$ are the left and right Riemann-Liouville fractional derivative.

Plenty of physical phenomena, such as the long-wave dynamics of two waves, are represented by the system (1.2). For example, these equations are used to study a number of issues in solid state physics, relativistic mechanics, quantum mechanics, and classical mechanics ^[1,2,3,4].

Especially, when $\alpha_{k}$ tends to 2, the fractional derivative $\partial_{x_k}^{\alpha_k}$ would converge to the second-order Laplace operator, and thus Eqs (1.1) and (1.2) reduce to the classical system of multi-dimensional coupled KG equations ^[5,6,7]. The system has the following conserved energy, which is mentioned in detail in ^[11],

where

The coupled KG equations is initially introduced in ^[8] and is applied to model the usual motion of charged mesons within a magnetic field. There have been many works for solving the classical KG equations. Tsutsumi ^[9] considered nonrelativistic approximation of nonlinear KG equations and proved the convergence of solutions rigorously. Joseph ^[10] obtained some exact solutions for these systems. Deng ^[11] developed two kinds of energy-preserving finite difference methods for the systems of coupled sine-Gordon (SG) equations or coupled KG equations in two dimensions. He ^[12] analyzed two kinds of energy-preserving finite element approximation schemes for a class of nonlinear wave equation. Zhu ^[13] developed the finite element method and the mesh-free deep neural network approach in a comparative fashion for solving two types of coupled nonlinear hyperbolic/wave partial differential equations. Deng ^[14] proposed a two-level linearized compact ADI method for solving the nonlinear coupled wave equations. More relevant and significant references can be found in ^[15,16,17].

However, it has been found that fractional derivatives can be used to describe some physical problems with the spatial non-locality of anomalous diffusion. Therefore, more attentions have been paid to fractional KG equations. There are also some related numerical methods for the related models. These methods may be applied to solve the fractional KG systems. For example, Cheng ^[18] constructed a linearized compact ADI scheme for the two-dimensional Riesz space fractional nonlinear reaction-diffusion equations. Wang ^[19] proposed Fourier spectral method to solve space fractional KG equations with periodic boundary condition. Liu ^[20] presented an implicit finite difference scheme for the nonlinear time-space-fractional Schrödinger equation. Cheng ^[21] constructed an energy-conserving and linearly implicit scheme by combining the scalar auxiliary variable approach for the nonlinear space-fractional Schrödinger equations. Similar scalar auxiliary variable approach can also be found in ^[22,23]. Wang et al. ^[24,25] developed some energy-conserving schemes for space-fractional Schrödinger equations. Meanwhile, the equations are also investigated by some analytical techniques, such as the Fourier transform method ^[26], the Mellin transform method ^[27] and so on. Besides, the spatial disccretization of the KG equations usually gives a system of conservative ordinary differential equations. There are also some energy-conserving time discretizations, such as the implicit midpoint method ^[28], some Runge–Kutta methods ^[28,29], relaxation methods ^[30,31,32] and so on ^[33,34]. To the best of our knowledge, there exist few reports on numerical methods for coupled space-fractional KG equations. Most references focus on the KG equations rather than the coupled systems.

The main purpose of this paper is to develop an EP-FDM for the system of nonlinear coupled space-fractional KG equations. Firstly, we transform the coupled systems of KG equations into an equivalent general form and provide energy conservation for the new system. Secondly, we propose a second-order consistent implicit three-level scheme by using the finite difference method to solve problems (1.1) and (1.2). Thirdly, we give the proof of the discrete energy conservation, boundedness of numerical solutions and convergence analysis in discrete $L^2$ norm. More specifically, the results show that the proposed schemes are energy-conserving. And the schemes have second-order accuracy in both the temporal and spatial directions. Finally, numerical experiments are presented to show the performance of our proposed scheme in one and two dimensions. They confirm our obtained theoretical results very well.

The rest of the paper is organized as follows. Some denotations and preliminaries are given in Section 2. An energy-preserving scheme is constructed in Section 3. The discrete conservation law and boundedness of numerical solutions are given in Section 4. The convergence results are given in Section 5. Several numerical tests are offered to validate our theoretical results in Section 6. Finally, some conclusions are given in Section 7.

Throughout the paper, we set $C$ as a general positive constant that is independent of mesh sizes, which may be changed under different circumstances.

2.   Denotations and preliminaries

We first rewrite Eqs (1.1) and (1.2) into an equivalent form

with the widely used boundary and initial conditions

where $G(u, v) = \frac{b_1}{4 c_1} u^4+\frac{b_2}{4 c_2} v^4+\frac{a_1}{2 c_1} u^2+\frac{a_2}{2 c_2} v^2+\frac{1}{2} u^2 v^2$, and $\alpha = 1/c_{1}$, $\beta = \kappa^{2}/c_{1}$, $\gamma = 1/c_{2}$, $\sigma = \kappa^{2}/c_{2}$. A similar treatment is mentioned in ^[11]. The definition of operator $\partial_{x_k}^{\alpha_k}$ is already presented in Eq (1.6), where the left and right Riemann-Liouville fractional derivatives in space of order $\alpha$ are defined as

Theorem 1. Let ${u(\boldsymbol{x}, t), v(\boldsymbol{x}, t)}$ be the solutions of this systems (2.1)–(2.5), the energy conservation law is defined by

Namely, $E(t) = E(0)$, where $\|u(\cdot, t)\|_{L^2}^2 = \int_{\Omega}|u(\boldsymbol{x}, t)|^2 d \boldsymbol{x}$ and $\langle G(u, v), 1 \rangle = \int_{\Omega}G(u, v)d \boldsymbol{x}$.

Proof. Taking inner product of Eqs (2.1) and (2.2) with $u_{t}$ and $v_{t}$, then summing the obtained equations, and finally applying a integration over the time interval $[0, t]$, it yields the required result.

The finite difference method is used to achieve spatial and temporal discretization in this paper. We now denote temporal step size by $\tau$, let $\tau = T/N$, $t_{n} = n\tau$. For a list of functions $\{w^n\}$, we define

Let $\Omega = (a_{1}, b_{1})\times\cdots (a_{d}, b_{d})$, with the given positive integers $M_{1}, \cdots, M_{d}$, for the convenience of subsequent proofs, we have set it uniformly to M, so we get $h_{k} = (b_{k}-a_{k}) / M_{k} = h \ (k = 1, \cdots, d)$ be the spatial stepsizes in $x_{k}$-direction, then the spatial mesh is defined as $\bar{\Omega}_{h} = \{(x_{k_{1}}, x_{k_{2}}, \cdots, x_{k_{d}}) \mid 0 \leq k_{s} \leq M_{s}, s = 1, \cdots, d \}$, where $x_{k_{s}} = a_{s}+k_{s} h_{s}$.

Moreover, we define the space $\mathcal{V}_{h}^{0}$ as follows by using the grid function on $\Omega_{h}$,

where $1\le k_{s}\le M_{s}-1, \ s = 1, \cdots, d, \ 0\leq n\leq N$. Then we write $\delta_{x_{1}} u_{k_{1}\cdots k_{d}} = \frac{u_{k_{1}+1\cdots k_{d}}-u_{k_{1}\cdots k_{d}}}{h}$. Notations $\delta_{x_{s}} u_{k_{1}\cdots k_{d}} \ (s = 2, \cdots, d)$ are defined similarly.

We then introduce the discrete norm, respectively. For $u, v \in \mathcal{V}_{h}^{0}$, denote

Based on the definitions, we give the following lemmas which are important for this paper.

Lemma 1. (^[35]) Suppose $p(x) \in L_1(\mathbb{R})$ and

where $\hat{p}(k)$ is the Fourier transformation of $p(x)$, then for a given $h$, it holds that

where $w_k^{(\alpha)}$ are defined by

where $\lambda_1 = \left(\alpha^2+3 \alpha+2\right) / 12, \quad \lambda_0 = \left(4-\alpha^2\right) / 6, \quad \lambda_{-1} = \left(\alpha^2-3 \alpha+2\right) / 12$ and $g_k^{(\alpha)} = (-1)^k\left(\begin{array}{l}\alpha \\ k\end{array}\right)$.

In addition, we arrange in this section some of the lemmas that are necessary for the demonstration of later theorems in this paper.

Lemma 2. (^[36]) For any two grid functions $u, v\in \mathcal{V}_{h}^{0}$, there exists a linear operator $\Lambda^{\alpha}$ such that $-(\delta_{x}^{\alpha}u, v) = (\Lambda^{\frac{\alpha}{2}}u, \Lambda^{\frac{\alpha}{2}}v)$, where the difference operator $\Lambda^{\frac{\alpha}{2}}$ is defined by $\Lambda^{\frac{\alpha}{2}} u = {\bf L} u$, and matrix ${\bf L}$ satisfying ${\bf C} = {\bf L}^{T}{\bf L}$ is the cholesky factor of matrix ${\bf C} = 1/(2h^{\alpha}\cos(\alpha\pi/2))({\bf P}+{\bf P}^{T})$ with

While for multi-dimensional case, we give a further lemma.

Lemma 3. (^[18]) For any two grid functions $u, v\in \mathcal{V}_{h}^{0}$, there exists a linear operator $\Lambda_{k}^{\frac{\alpha_{k}}{2}}$ such that $-(\delta_{x_{k}}^{\alpha_{k}}u, v) = (\Lambda_{k}^{\frac{\alpha_{k}}{2}}u, \Lambda_{k}^{\frac{\alpha_{k}}{2}}v), \ k = 1, \cdots, d,$ where $\Lambda_{k}^{\frac{\alpha_{k}}{2}}$ is defined by $\Lambda_{k}^{\frac{\alpha_{k}}{2}} u = [2\cos(\alpha_{k}\pi/2)h^{\alpha_{k}}]^{-1/2}{\bf L_{k}} u$, and matrix ${\bf L_{k}}$ is given by $-{\bf I} \otimes \cdots {\bf D}_{\alpha_{{\bf k}}} \otimes {\bf I} =$ $\left[2 \cos \left(\alpha_k \pi / 2\right) h^{\alpha_k}\right]^{-1} {\bf L}_{{\bf k}}^T {\bf L}_{{\bf k}}$. ${\bf I}$ is a unit matrix and matrix ${\bf D}_{\alpha_{{\bf k}}}$ is given by ${\bf D}_{\alpha_{{\bf k}}} = -1 /\left(2 \cos \left(\alpha_k \pi / 2\right) h^{\alpha_k}\right)$ $\left({\bf P}_{{\bf k}}+{\bf P}_{{\bf k}}^T\right)$, ${\bf P}_{{\bf k}}$ is the matrix ${\bf P}$ in the case $\alpha = \alpha_{k}$ as defined in Lemma 2.

Lemma 4. (^[11]) Let $g(x) \in \mathcal{C}^4(I)$, then $\forall x_0 \in I, x_0+\Delta x \in I$, we have

Lemma 5. (^[11]) Let $u(\boldsymbol{x}, t), v(\boldsymbol{x}, t) \in \mathcal{C}^{4, 4}(\Omega \times[0, T])$, and $G(u, v) \in \mathcal{C}^{4, 4}\left(R^1 \times R^1\right)$. Then we have

Lemma 6. For any grid function $u \in \mathcal{V}_h^0$, it holds that

where $C_{p_1}, C_{p_2}, C_{p_3}$ are constants related to p, $l = \min \left\{l_1, \cdots, l_d\right\}$, and d is the dimension of space $\mathcal{V}_h^0$.

Specially, for two-dimensional case, the parameters $C_{p_1} = \frac{2}{p}$, $C_{p_2} = \max \left\{2 \sqrt{2}, \frac{p}{\sqrt{2}}\right\}$ and $C_{p_3} = 1-\frac{2}{p}$ are shown in ^[37,38].

While in the case of three dimensions, $C_{p_1} = \frac{p+6}{4p}$, $C_{p_2} = \max \left\{2 \sqrt{3}, \frac{p}{\sqrt{3}}\right\}$ and $C_{p_3} = \frac{3p-6}{4p}$, the proof is given in Appendix.

Lemma 7. (^[39]) For $M \geq 5, 1 \leq \alpha \leq 2$ and any $v \in \mathcal{V}_h^0$, there exists a positive constant $C_1$, such that

Specially, for multi-dimensional case, it can be written as $\|v\|^2 < C \sum_{k = 1}^d\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}} v\right\|^2$, where $C$ is a positive constant.

Lemma 8. (^[40]) Assume that $\left\{g^n \mid n \geq 0\right\}$ is a nonnegative sequence, $\psi^0 \geq 0$, and the nonnegative sequence $\left\{G^n \mid n \geq 0\right\}$ satisfies

Then it holds that

Lemma 9. For any grid function $u \in V_h^0$, $V_h^0$ is defined in Section 2 for the the three-dimensional case, let $p \leq r \leq q, \alpha \in[0, 1]$ satisfying $\frac{1}{r} = \frac{\alpha}{p}+\frac{1-\alpha}{q}$, then

Proof. By using Hölder inequality, we have

This completes the proof.

3.   The energy-preserving scheme

Now we are ready to construct the fully-discrete numerical scheme for systems (2.1) and (2.2).

With the help of Lemma 1 and for clarity of description, we will denote the space fractional operator under one-dimensional case firstly.

where $w_{k}^{(\alpha)}$ is given in Eq (2.7). In the multi-dimensional case, the definitions of $\delta_{x_{k}}^{\alpha_{k}}$ are similar to it.

For numerically solving systems (2.1)–(2.5), we propose a three-level scheme. We firstly define the following approximations.

Let $u^{n}_{k_{1}\cdots k_{d}} = u(\boldsymbol{x}, t_{n})$ and $v^{n}_{k_{1}\cdots k_{d}} = v(\boldsymbol{x}, t_{n})$, for ease of presentation, we shall henceforth write $u^{n}_{k_{1}\cdots k_{d}}$ for $u^{n}$. Denote numerical solutions of $u^{n}$ and $v^{n}$ by $U^{n}$ and $V^{n}$, respectively.

With the definition of $G(u, v)$ in systems (2.1) and (2.2) and by using Lemma 5, then we have

which is given in ^[11]. Further, using the space fractional operator which is already introduced above and second-order centered finite difference operator to approximate at node $(\boldsymbol{x}, t_{n})$, it holds that

and

where $R_1^n$ and $R_2^n$ are the truncation errors.

Let $u(\boldsymbol{x}, t), v(\boldsymbol{x}, t) \in \mathcal{C}^{4, 4}(\Omega \times[0, T])$. Combining Lemma 4 with Eqs (3.1) and (3.2), the truncation errors can be estimated as follows.

where $C$ is a positive constant and $d$ means the dimension of space.

Omitting the truncation errors in Eqs (3.3)–(3.6), we can get the three-level EP-FDM:

and

where $U^{1}$ and $V^{1}$ are obtained by applying Taylor expansion to expand $u(\boldsymbol{x}, \tau)$ and $v(\boldsymbol{x}, \tau)$ at $(\boldsymbol{x}, 0)$, and by Eq (2.4) we know that $U^0 = \phi_1\left(\boldsymbol{x}\right)$, $V^0 = \phi_2\left(\boldsymbol{x}\right)$.

For contrast, by doing explicit treatment of nonlinear terms $\frac{\partial G}{\partial u}$ and $\frac{\partial G}{\partial v}$, we introduce an explicit scheme as follows

which will be used in Section 6 later.

4.   Boundedness of the numerical solutions and discrete conservation law

In this section, we give the energy conservation of the fully-discrete schemes (3.8)–(3.12) and boundedness of numerical solutions. Here, the lemmas given in Section 2 are applied.

Now, we present the energy conservation of the EP-FDMs (3.8)–(3.12).

Theorem 2. Let $U^n, V^n\in \mathcal{V}_{h}^{0}$ be numerical solutions of the three-level FDMs (3.8)–(3.12). Then, the energy, which is defined by

is conservative. Namely, $E^n = E^0$, for $n = 1, \cdots, N-1,$ where $\Lambda_{k}^{\frac{\alpha_{k}}{2}}$ is already introduced by Lemma 3.

Proof. Multiplying $h^dD_{t}U_{k_{1}\cdots k_{d}}^n$ to both sides of Eq (3.8), summing them over $\Omega_{h}$, by using Lemma 3, we obtain

where the second term can be reduced to

then Eq (4.2) turned into

Similarly, multiplying $h^dD_{t}V_{k_{1}\cdots k_{d}}^n$ to both sides of Eq (3.9), summing them over $\Omega_{h}$, by using Lemma 3, we obtain

Adding up Eqs (4.3) and (4.4) yields that $(E^n-E^{n-1})/\tau = 0$, which infers that $E^n = E^{n-1}$.

By Theorem 2, we present the following estimation.

Theorem 3. Let $U^n, V^n\in \mathcal{V}_{h}^{0}$ be numerical solutions of the EP-FDMs (3.8)–(3.12). Then, the following estimates hold:

where $C$ is a positive constant independent of $\tau$ and $h$ and $1\leq \alpha_{k}\leq 2.$ Specially, when $\alpha_{k} = 2,$ it holds that $|U^n|_{H_{1}}\leq C$, $|V^n|_{H_{1}}\leq C$.

Proof. It follows from Theorem 2, there exists a constant $C$ such that

then, we obtain

By $\|\delta_{t}U^n\|\leq C$, we have $\|U^{n+1}-U^n\|\leq C\tau$, then it is easy to check that

This completes the proof.

5.   Convergence analysis

In this section, the convergence analysis of the proposed scheme is given, which is based on some important lemmas presented in Section 2.

We first give the error equations of the EP-FDMs (3.8) and (3.9). Let $e^n = u^n-U^n$, $\theta^n = v^n-V^n$ and for more readability we denote

By deducting Eqs (3.8) and (3.9) from Eqs (3.3) and (3.4), we have

Before giving a proof of convergence, we provide the following estimates for Eqs (5.1)–(5.2).

Lemma 10. On $\bar{\Omega}_h$, we have

where $C > 0$ is a constant, independent of grid parameters $\tau, h_1, \cdots, h_d$.

Proof. Recalling the definition of $G(u, v)$, we can obtain

Noting that $U^k = u^k-e^k$ and $V^k = v^k-\theta^k\ (k = n-1, n, n+1)$, then we get

When $d = 2$, combining Theorem 3, Lemma 6 with Lemma 7, we can get the estimation of $\|e^m\|_{4}^{4}$, $\|e^m\|_{6}^{6}$, $\|e^m\|_{8}^{8}$, that is

The same reasoning can be used to prove that

Smilarly, when $d = 3$ the results can be found in the same way.

By using Cauchy-Schwarz inequality and the widely used inequality $[(a+b)/2]^s\leq(a^s+b^s)/2 \ (a\geq0, b\geq0, s\geq1)$, multiplying both sides of Eq (5.11) by $h^dD_{t}e^n$, then summing it on whole $\Omega_h$, it follows that

the last inequality is derived by inequalities (5.13) and (5.14), similarly, we can also obtain

combine inequalities (5.15)–(5.17), then we get inequality (5.9) is proved. We can demonstrate that inequality (5.10) is likewise true using techniques similar to inequality (5.9). This completes the proof.

Now we further investigate the accuracy of the proposed scheme with the help of the above lemmas, see Theorem 4.

Theorem 4. Assume that $u(\boldsymbol{x}, t), v(\boldsymbol{x}, t) \in \mathcal{C}^{4, 4}(\Omega \times[0, T])$ are exact solutions of systems (2.1)–(2.5), let $u_{k_{1}\cdots k_{d}}^n = u(\boldsymbol{x}, t)$ and $v_{k_{1}\cdots k_{d}}^n = v(\boldsymbol{x}, t)$, denote numerical solutions by $U_{k_{1}\cdots k_{d}}^n$ and $V_{k_{1}\cdots k_{d}}^n$, define $e^n = u^n-U^n$, $\theta^n = v^n-V^n (1\leq n\leq N)$. Then suppose that $\tau$ is sufficiently small. The error estimates of the EP-FDM are

where $C$ is a positive constant, independent of grid parameters $\tau, h_1, \cdots, h_d$.

Proof. Noting that at every time level, the systems defined in Eqs (3.8) and (3.9) is a linear PDE. Obviously, the existence and uniqueness of the solution can be obtained.

For ease of expression, we write

Apparently, we have that $I^1\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2$.

Multiplying $h^dD_{t}e^n$ and $h^dD_{t}\theta^n$ to both sides of Eqs (5.3) and (5.4), then summing it over the whole $\Omega_h$ respectively. Then adding up the obtained results, it follows that

by using Cauchy-Schwarz inequality, we have

multiplying $2\tau$ to both sides of inequality (5.19), and using Lemma 10, then we get

Thus, $\forall K (2\leq n\leq K\leq N-1)$, summing $n$ from $2$ to $K$, we get

when $C\tau\leq\frac{1}{3}$, inequality (5.21) is turned into

then by using Lemma 8 and inequality (3.7), we obtain

By the definition of $I$, it is easy to conclude that

furthermore, we have

Similarly, $\|\theta^n\|\leq C(\tau^2+h_1^2+\cdots+h_d^2).$ This completes the proof.

6.   Numerical experiments

We carry out several numerical examples to support the theoretical results in this section. All computations are performed with Matlab. Throughout the experiments, the spatial domain is divided into $M$ parts in every direction uniformly, that is, in the 1D case, we set $M_{1} = M$, while in the 2D case, we set $M_{1} = M_{2} = M$, and the time interval $[0, T]$ is also divided uniformly into $N$ parts. Then we use the discrete $L^\infty$-norm to measure the global error of the scheme, namely,

Example 1. Consider the following one-dimensional coupled KG model

with $\Omega = [0, 1]$. The initial and boundary conditions are determined by the exact solutions

as well as the source term $g$. Here, we take $a_1 = a_2 = 1$, $b_1 = -1$, $b_2 = -2$, $c_1 = 1$, $c_2 = 0.5$ and $\kappa = 1$.

The precision of the scheme in spatial direction is first tested by fixing $N = 1000$. We compute the global errors at $T = 1$ with different mesh sizes, and the numerical results with $\alpha = 1.2, 1.5, 1.8$ are listed in Table 1 and Table 2. As can be seen in the table, the proposed scheme can have second order convergence in space, which confirms the results of theoretical analysis in Theorem 4. To track the evolution of the discrete energy, we preserve the initial value condition in this case and set the source term to $g = 0$. Additionally, for the terminal time $T = 50$, we fix $h = 0.05$ and $\tau = 0.05$. The evolutionary trend image for scheme1 (3.8)–(3.12) and explicit scheme2 (3.13)–(3.17) with various $\alpha$ are displayed in Figure 1. Then we further verify that the proposed scheme1 (3.8)–(3.12) preserves the discrete energy very well but scheme2 (3.13)–(3.17) does not.

Example 2. Consider the following two-dimensional coupled KG model

with $\Omega = [0, 2] \times[0, 2]$. The initial and boundary conditions are determined by the exact solutions

as well as the source term $g$. Here, we take $a_1 = a_2 = 1$, $b_1 = -1$, $b_2 = -2$, $c_1 = 1$, $c_2 = 0.5$ and $\kappa = 1$.

Similar to Example 1, we verify the convergence orders of the scheme in spatial direction at $T = 1$. For spatial convergence order, we still set $N = 1000$ and thus the temporal error of the scheme can be negligible. The numerical results are presented in Table 3 and Table 4 with different values of $\alpha_1$ and $\alpha_2$ which are in the x and y directions, respectively. The second-order accuracy of the scheme is achieved. Moreover, for the terminal time $T = 100$, Figure 2 shows the evolution of discrete energy for scheme1 (3.8)–(3.12) and explicit scheme2 (3.13)–(3.17) when $g(x, y, t) = 0$. The figure indicate that the discrete conservation law holds very well if the proposed scheme1 (3.8)–(3.12) are used. In contrast, scheme2 (3.13)–(3.17) cannot preserve the discrete energy. Both tables and figure further confirm the theoretical results.

Example 3. Consider the following two-dimensional coupled KG model

and

with $\Omega = [0, 1] \times[0, 1]$.

Here, we take

and

The scheme1 (3.8)–(3.12) with

are used to Example 3. Figure 3 and Figure 4 show the surfaces of $U_{ij}^n$ and $V_{ij}^n$ at different times, respectively. The significant dynamical evolutionary features of the numerical solutions $U_{ij}^n$ and $V_{ij}^n$, such as radiation and oscillation, can be found in Figure 3 and Figure 4.

7.   Conclusions

In this paper, the three-level energy-preserving scheme is proposed for the space-fractional coupled KG systems. The scheme is derived by using the finite difference method. The discrete conservation law, boundedness of numerical solutions and the global error of the scheme are further discussed. It is shown that the scheme can have second order convergence in both temporal direction and spatial direction. Several numerical examples are performed to support the theoretical results in the paper. Moreover, due to the nonlocal derivative operator and considering that the implicit methods involve Toeplitz matrices, fast methods are fairly meaningful to reduce the computational cost of the proposed scheme; refer to the recent work ^[41,42] for this issue.

Acknowledgments

This work is supported by NSFC (Grant Nos. 11971010, 12001067).

Conflict of interest

The authors declared that they have no conflicts of interest to this work.

AppendixA

In the following, we present the proof of Lemma 6.

Proof of Lemma 6: Obviously, the result holds for $p = 2$. We prove the conclusion for $p > 2$.

For any $m, s = 1, 2, \ldots, M_1-1$, and $m > s$, using mean value theorem, we have

where

Then,

It is easy to verify the above inequality also holds for $m \leq s$. Thus, we have

Multiplying the above inequality by $h_1$ and summing up for $s$ from 1 to $M_1-1$, we have

Dividing the result by $l_1$, and noticing that the above inequality holds for $m = 1, 2, \ldots, M_1-1$, we have

Multiplying the above inequality by $h_2h_3$ and summing over $j, k$, then applying the Cauchy-Schwarz inequality, we obtain

Multiply both sides of inequality (A.1) by $(h_2h_3)^{\frac{1}{2}}$, it follows easily that there exists a constant C such that $(h_2h_3)^{\frac{1}{2}}\leq C$, we obtain

Similarly to the previous analysis, we have

Using the Cauchy-Schwarz inequality, we have

Substituting the above inequality into inequalities (A.2)–(A.4), we have

We now estimate $\|u\|_p^p$,

the last inequality is obtained by inequalities (A.5)–(A.7).

In addition, we set $l = \min \left\{l_1, l_2, l_3\right\}$, by using mean value inequality then we have

Combining inequalities (A.8) and (A.9) yields

We consider the case $p \geq 6$, applying Lemma 9 for $p \geq 6$, it holds

Substituting the above inequality into inequality (A.10), we get

that is

Thus, we have proved the result for $p \geq 6$. Taking $p = 6$ in inequality (A.11) yields

When $2 < p < 6$, using Lemma 9 and inequality (A.12), we have

This completes the proof.