Numerical investigation of two models of nonlinear fractional reaction subdiffusion equations

We consider new numerical schemes to solve two different systems of nonlinear fractional reaction subdiffusion equations. These systems of equations model the reversible reaction A+B⇌C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A+B \rightleftharpoons C$$\end{document} in the presence of anomalous subdiffusion. The first model is based on the Henry & Wearne [1] model where the reaction term is added to the subdiffusion equation. The second model is based on the model by Angstmann, Donnelly & Henry [2] which involves a modified fractional differential operator. For both models the Keller Box method [3] along with a modified L1 scheme (ML1), adapted from the Oldham and Spanier L1 scheme [4], are used to approximate the spatial and fractional derivatives respectively. Numerical prediction of both models were compared for a number of examples given the same initial and boundary conditions and the same anomalous exponents. From the results, we see similar short time behaviour for both models predicted. However for long times the solution of the second model remains positive whilst the Henry & Wearne based–model predictions may become negative.


Introduction
Fractional reaction subdiffusion equations have been derived from Continuous Time Random Walk models that take into account the effect of anomalous subdiffusion on the reaction process by using long-tailed waiting time densities [1,2,[5][6][7]. For a variety of cases, master equations that control the temporal development of particle density have been derived. Henry and Wearne introduced this with mass action reaction terms added to the sub-diffusion term [1]. The study of the solution of a fractional reaction-subdiffusion equations has become more prominent and important since there is growing estimation that anomalous diffusion is in fact ubiquitous in nature [8]. But analytic solutions of such equations are seldom available and so numerical techniques are needed. Even when they available they often involve special functions, such as the Fox (H-function) function [9] and the Mittag-Leffler function [10], which are difficult to routinely evaluate. Hence numerical solution schemes are required.
There are a number of numerical schemes that have been used to solve fractional reaction-diffusion equations, [11][12][13][14][15][16][17][18][19][20][21][22]. Recently Angstmann, Donnelly, Henry & Nichols [23] introduced a novel numerical scheme for solving a fractional Fokker-Planck equation where the numerical scheme was constructed by a discrete time and space stochastic process. They found that the probability density can be evaluated and used to approximate the solution of the fractional Fokker-Planck equation. An explicit numerical method is also developed by Angstmann et al. [24] is based on a stochastic process and has been used to solve a fractional reaction subdiffusion equation. They showed, in the diffusion limit, the master equation recovers the fractional partial differential equation. In [25] a similar approach was also used to solve sets of fractional ordinary differential equation for a modified SIR epidemic model. In this article we consider two models of reversible reactions in the presence of anomalous subdiffusion, which we solve numerically and compare their predictions. In these models we let A, B and C be three chemical species undergoing a reversible reaction, A + B C. In the absence of diffusion, the concentration of A, B, and C can be modelled by the following system of the reaction kinetic equations where k 1 and k −1 are the forward and reverse reaction rates respectively. These equations correspond to the reaction A + B → C, if k −1 = 0, [26][27][28][29]. Reversible reactions, in the presence of subdiffusion has been modelled by the system of fractional reaction-diffusion equations by using the CTRW model as in the work of Henry & Wearne [1], which we will denote as Model 1, (1.6) where the reaction term R(x, t, A, B, C) is defined by (1.7) These equations involve the fractional Riemann-Liouville derivative operator [10,30], where 0 < γ < 1. The solution of the Model 1 in the case C → A + B was found by Langlands, Henry & Wearne [31] in the infinite domain. In their results a negative value was predicted which is physically unrealistic. We also consider a second model, which we will denote as Model 2, which is based upon a more recent model from Angstmann, Donnelly & Henry [2] which has a modified fractional operator. For the reversible reaction A + B C we have the following model (1.10) This system includes a non-standard fractional derivative operator L which takes into account species may react and therefore be removed before diffusion can take place. The main difference between the two models is Model 2 contain a non-standard fractional derivative operator (1.11) but it is not include in Model 1. Anomalous nodes in a network is an example for both models. Fedotov and Stage demonstrated that anomalous cumulative inertia overpowers highly connected nodes in attracting network individuals. This contradicts the classical result that people tend to gather at high-order nodes [32]. The Keller Box method is an implicit numerical scheme which is second order accurate in both space and time when applied to the standard diffusion equation [33]. In the Keller Box method the order of spatial derivative in the equation is rewritten as the first derivative of an introduced additional variable. This method can be used for more general equations where we cannot rewrite the fractional partial differential equation with a Caputo derivative on the left. In this work, the current methods for approximating fractional derivatives will need to be modified to approximate the operator in equation (1.11). We consider numerical solutions for equations (1.4) -(1.6) and equations (1.8) -(1.10), by applying the Keller Box method with the modification of the L1 scheme as in [34,35]. The main contribution is the discretization of the modified fractional operator, there may be a similar method was used in [7].

Derivation of the numerical method for both models
In this section a numerical scheme for solving both models is developed based upon the Keller Box method. As in [34,35], we denote the spatial grid points by x i , for We also use the equally-spaced temporal points as t j = j t, for j = 0, 1, ..., M with the constant time step of t = T /M. In the following sections we discuss in detail the derivation of the scheme for Model 2. The scheme for Model 1 can be developed in a similar manner. We approximate the Riemann-Liouville fractional derivative by using the ML1 scheme developed in [34,35], given by where σ γ and the weights β j (γ ) and μ j−k (γ ) are given by The ML1 scheme has been shown to be convergent of order [35]. The similar scheme to ML1 scheme was used in [36] but with different weights and include the estimation at the midpoints t k+1/2 , where k = 0, 1, 2, . . . , j. For the scheme given in [36] the convergence order was similar to that in [35].

The scheme for model 2
In this section we consider the second model given by equations (1.8) -(1.10). These equations include the non-standard fractional derivative operator in equation (1.11). The current methods for approximating fractional derivatives will need to be modified to approximate the operator in equation (1.11). To do this we also define the auxiliary variables y 1 and y 2 as Taking the derivative with respect to t of equations (2.3) we then find the following governing differential equations and These equations are supplemented by the initial conditions y 1 (x, 0) = 1 and y 2 (x, 0) = 1. In the following we denote (2.9) Equation (2.7) is supplemented by (2.10) Approximating the first order spatial and time derivatives in equations (2.8) and (2.10), by using the centred finite difference scheme, we then obtain the equations (2.12) We then approximate the term [R] by the corresponding spatial averages at i − 1 and i, and so equations (2.11) and (2.12) then become and Solving equation (2.14) for v j i−1 and then combining with equation (2.13) gives Using a similar process to above, except now replacing i with i + 1 in equations (2.8) and (2.10) and eliminating v j i+1 , we obtain the equation Combining these two equations then gives After using the ML1 approximation of the fractional derivative, given by equation (2.1), in equation (2.17), replacing the terms at t = t j+ 1 2 by their corresponding temporal averages and using the notation in equation (2.6), we then have the following equation, given in the case of constant grid spacing Here, at point p = i, i − 1, and i + 1, the terms H j p (y 1 , A) and φ and where σ γ is defined in equation (2.2), and d is given by In the above we have used the condition y 1 (x, 0) = 1. The corresponding equation for a uniform mesh for species B in equation (1.9) is Finally we find the approximations for the equations for the auxiliary variables y 1 (x, t) and y 2 (x, t) given in equations (2.4) and (2.5). Approximating the first order time derivatives by a centred finite difference and approximating the values at t = t j+ 1 2 by their temporal average, we then obtain the equations

The scheme for Model 1
Using a similar procedure that was used for Model 2 as given in Section 2.1, the corresponding equations for Model 1, equations (1.4), (1.5) and (1.6), are

Accuracy of the numerical method for both models
In this section we consider the order of accuracy of the numerical schemes in Sections 2.1 and 2.2. We suppose that the continuous problem for both models, equations (1.4) -(1.6) and equations (1.8) -(1.10), has a smooth solution such that A(

Accuracy of the numerical method for Model 2
We now determine the truncation error of KBML1 scheme for Model 2. Equation (2.18) can be rewritten as Note in the above, the following identity was used After adding and subtracting the term 2D t γ −1 with p = i, i − 1, i + 1, equation (3.2) then becomes Note the second term on the right-hand side of equation ( with j i (y 1 , A) = y 1 and Taking the Taylor series expansion around the point x i , t j+ 1 2 , we then obtain ∂ A ∂t (3.10) Evaluating the spatial derivatives and then combining the common terms we find ∂ A ∂t (3.12) Each of the terms of the form are of order t 1+γ accurate in time given the ML1 scheme is O( t 1+γ ) [35]. We then obtain ∂ A ∂t (3.14) Using a similar argument to that in [34], we see that the truncation error, τ

Accuracy of the numerical scheme for Model 1
Similar to Section 3.1, we find the truncation error for equations (2.27), (2.28) and (2.29) are of order 1 + γ in time and second order in space.

Consistency of the numerical schemes
The numerical schemes for solving Model 1 and Model 2 are consistent, as the truncation error approaches zero as t → 0 and x → 0. Hence the Keller Box method with ML1 scheme is consistent with the original systems of fractional reaction diffusion equations.

Numerical examples and results
In this section we investigate the solution of both Models 1 and 2, for Examples 1 and 2 on the domain 0 ≤ x ≤ 1 and 0 ≤ t ≤ 1 and for Example 3 on the domain 0 ≤ x ≤ 1 and 0 ≤ t ≤ 5. We estimate the order of convergence numerically for the KBML1 method by computing the maximum norm of the error between the numerical estimate and the approximate "exact" solution at the time t = 1 for Examples 1 and 2, and at t = 5 for Example 3. The approximate "exact" solution was found by using a large number of grid points and time steps. The approximate order of convergence in x, R 1 , was estimated by computing and the approximate order of convergence in t, R 2 , was estimated by computing and where L = 1. For Model 2 we also include the initial conditions y 1 (x, 0) = 1 and y 2 (x, 0) = 1. For both models we set the fractional exponent as γ = 1/2, the forward reaction rate as k 1 = 1, the backward reaction rate as k −1 = 1, and the diffusion coefficient as D = 1. Here the numerical solutions of Model 1 and Model 2, where found using the numerical schemes, given in Sections 2.2 and 2.1 respectively, with t = 0.001 and x = 0.02. Figure 1 shows the results for A and C for Model 1. We see the reverse reaction C → A + B dominates with A and B being produced from C at a faster rate than the forward reaction A + B → C. Note the predicted values of B (not shown) are the same as those as A in this example. We also see that C decays to a homogeneous steady state, whilst A and B increase to a homogeneous steady state. The numerical solution of Model 2 under the same initial and boundary conditions is shown in Figure 2. Here we see similar behaviour to the results of Model 1. It should be noted, given the initial conditions are the same A(x, 0) = B(x, 0), then A(x, t) = B(x, t) for both models.
In Figure 3 we give a comparison of the predicted values of A(x, t) for both models at the points x = 0.5 and x = 0.9. From this figure we see similar asymptotic behaviour from both models with Model 2 (solid lines) predicting a slightly faster decay to the homogeneous state compared to Model 1 (dashed lines). We also see in Figure 4 the difference = A 1 (0.5, t) − A 2 (0.5, t) at x = 0.5. From this we see Model 1 predicts a slightly higher value than Model 2 at the midpoint.

Example 2
In this next example we consider the solution of the two fractional partial differential equation model where there is no reverse reaction, so we have the reaction  Figures 5 and 6 we show the numerical solution of Model 1 and Model 2 by using the KBML1 scheme subject to the boundary and initial conditions mentioned above. We again see similar asymptotic behaviour with the solution of A and B decaying to a homogeneous state for both models. A comparison is given in Figure 7 of the predicted values of A(x, t) and B(x, t) for both models at the points x = 0.3 and x = 0.9. There we see both models have a similar asymptotic behaviour, but Model 2 predicts a slightly faster decay to the homogeneous state. The differences between the predicted values for A and for B found from both models at x = 0.5 are shown in Figure 8. From this figure we again see the results from Model 1 predicts a slightly higher value at x = 0.5 than Model 2 for both species A and B.  Figure 9 we show the numerical predictions of both models. The predicted values of C from both models quickly decay to a homogeneous state before decaying to zero. A comparison of the values of C for Model 1 and Model 2 with γ = 0.5 at the points x = 0.2 and x = 0.9 is shown in Figure 10 with k −1 = 2, and in Figure 11 with k −1 = 5. We see in Figures 10 and 11, that the predicted solution of Model 1 becomes negative whilst Model 2's predictions remain positive. The solution becomes negative earlier as the value of k −1 is increased from 2 to 5. These results replicate the negative predicted results from [31] in the infinite domain case for Model 1. The negative prediction is physically unrealistic which suggests Model 2 with the modified operator is the better model to use.

Estimation of the order of convergence
In this section we estimate the order of convergence for Example 1. Before that we first demonstrate the predictions converge by using numerical experiments. To do this we calculate the difference between predictions in Figure 12 at x = 0.5, when different time steps, t, is used for Model 1 and Model 2 respectively over the solution interval 0 ≤ t ≤ 1. The value of 1 is the difference between the estimates if t = 10 −2 (100 time steps) and when the time step is t = 10 −5 is used (10 5 time steps). The value of 2 is the difference between when t = 10 −3 and t = 10 −5 time steps are used, and the value of 3 is the difference between when t = 10 −4 and t = 10 −5 time steps are used. We see the difference between the numerical predictions for species A decreases as the time step is decreased and appear to converge to zero as indicated by the arrows in each figure. We also found similar behaviour for Examples 2 and 3 as given in Figures 13 -15. The error and order of convergence estimates are next found from applying the numerical scheme on Model 1 and Model 2 for species A and C for Example 1, species A and B for Example 2, and species C for Example 3. The error is approximated by comparing these numerical results with an approximate "exact solution" found using a long run with a large number of time steps, with t = 1.25×10 −4 , and a large number of grid points, with x = 5 × 10 −4 . The results are given in Tables 1 and 2 for species A in Example 1 for γ = 0.1, 0.5, 0.9 and time t = 1.0. We kept t fixed at 10 −3 to estimate the convergence in space and varying x. To estimate the convergence in time we kept x fixed at 10 −3 whilst varying t. We see both numerical schemes, for Model 1 and Model 2, appear to be order O( x 2 ) which compares well with the accuracy analysis. Using a similar process, we found a similar convergence order in the space for species C for Example 1 and species A and B for Example 2 (not shown).
The obtained numerical convergence order estimates however do not appear to match up with the expected order of 1 + γ in time but the errors do decrease as t is decreased showing convergence for both models. We note that the convergence was slightly better for Model 2 than that for Model 1. Using a similar procedure, we estimate the convergence order using Example 3 for γ = 0.1, 0.5, 0.9 and time t = 5.0 as shown in Tables 3 and 4. Here we kept t = 0.1 fixed and varied x to estimate the convergence order in space. Likewise, we kept x = 10 −3 fixed and varied t to estimate the convergence order in time. From the results in Tables 3 and 4, we see similar convergence order in space and time as was shown in Example 1.
Stability and convergence are very important requirements for a robust numerical scheme. These models involve a very complicated set of equations and so it is not possible to obtain exact analytical conditions for the stability and convergence. However, we have tested our scheme under different time and grid steps and found the solutions do converge.

Conclusions
In this work, we extended the KBML1 scheme given in [34] to the case of systems of nonlinear fractional partial differential equations. We considered two models of a     Table 4 Numerical convergence order in x for Model 1 and Model 2, for Example 3 for species C(x, t), and R2 is the order of convergence  [2]. The accuracy of the KBML1 method for both models, by using the Taylor series, was found to be order 1 + γ in time and second order in space. We note that Model 2 takes longer to run computationally when compared to Model 1. This is to be expected since we need to solve a system of five differential equations in Model 2 rather than three in Model 1. The convergence of the KBML1 scheme for both models has been demonstrated numerically noting the error in the numerical predictions decreased as the time step was decreased. However, the expected 1 + γ order of convergence was not observed but was greater than 1. The spatial convergence order of two was confirmed for both model schemes. The two models were compared for a number of examples given the same initial and boundary conditions and with the same anomalous exponent. From the results found, similar behaviour for both models was predicted for short times. However, we see for longer times the solution of Model 2 remains positive whilst the Model 1 predictions may become negative which is physically unrealistic as shown for the case of the reaction C → A + B. This verifies numerically the analytic results in [31]. We conclude that using Model 2 with the modified operator is more realistic than using Model 1 if we require the solution to remain positive.