Adaptive Iterative Splitting Methods for Convection-Diffusion-Reaction Equations

: This article proposes adaptive iterative splitting methods to solve Multiphysics problems, which are related to convection–diffusion–reaction equations. The splitting techniques are based on iterative splitting approaches with adaptive ideas. Based on shifting the time-steps with additional adaptive time-ranges, we could embedded the adaptive techniques into the splitting approach. The numerical analysis of the adapted iterative splitting schemes is considered and we develop the underlying error estimates for the application of the adaptive schemes. The performance of the method with respect to the accuracy and the acceleration is evaluated in different numerical experiments. We test the beneﬁts of the adaptive splitting approach on highly nonlinear Burgers’ and Maxwell–Stefan diffusion equations.


Introduction
In this paper, we propose adaptive splitting schemes to solve nonlinear differential equations.We consider spatially discretized convection-diffusion-reaction equations, which we could apply as semi-discretized nonlinear systems of ordinary differential equations.Based on the nonlinearities, it is important to deal with adaptive schemes, whereas we can control the local errors of the underlying schemes, see [1][2][3][4].In general, the splitting methods have local splitting errors, which can be controlled with the time-or spatial steps of the underlying schemes, see [1,2,5,6].

•
Iterative methods, e.g., iterative splitting methods, see [3,10], or Waveform-Relaxation methods, see, e.g., [11,12] For the non-iterative methods, e.g., Lie-Trotter and Strang-splitting schemes, first works exist and discuss the ideas of the adaptive time splitting methods, see [6,8].Although for the iterative methods only some first ideas exist, see [13], based on dealing with different time-step approaches, our new contributions are based on the novel strategy of -shifting, see [14], of the underlying splitting method.We obtain so-called shifted iterative splitting methods, which can be compared with the standard iterative splitting methods, i.e., without the shifting.The local error estimate can be computed with the solution of the shifted and non-shifted iterative splitting approaches.Then, the effective error control is given as a function of the local error estimates and an error-tolerance parameter.
Such a novel adaptive splitting approach is important to reduce the computational time, while we decompose into sub-operators, which can be simulated with faster and more accurate numerical schemes.Further, the novel approach can be optimized with an effective and maximal splitting time step, see also [5,14].We combine these two ideas of adaptivity and splitting into a novel strategy to control and reduce the local splitting error for the iterative splitting methods.Then, we can obtain more accurate results with a maximal time step and reduce the computational cost.
In this paper, we present the novel adaptive splitting techniques as follows.

•
In the first step, we consider the standard splitting techniques with the underlying error analysis, see [3], and • in the second step, we introduce the adaptive techniques, which is based on the -shift technique, see [14], such that we can control the local splitting error.
Then, an error-estimate is computed, such that we are allowed to evaluate a maximum splitting time-step with respect to the shifted and non-shifted iterative splitting approach.The analysis is based on the standard iterative splitting approaches, see [13,15], and the additional adaptive techniques, see [5,6].In the numerical applications, which are based on convection-diffusion-reaction equations, we present the verification and the benefits of the novel adaptive iterative splitting approaches.
The paper is outlined as following.In Section 2, we explain the adaptive splitting approaches.Further, in Section 3, we discuss the error analysis of the adaptive splitting schemes.The applications to different convection-diffusion-reaction equations are done in Section 4. In Section 5, we discuss the theoretical and practical results.

Adaptive Splitting Approaches
We take inspiration for our studies, which are presented below, from real-life simulations of nonlinear convection-diffusion-reaction equations with the help of splitting approaches, see [3,[16][17][18].
We deal with convection-diffusion-reaction equations, which can be written as where f : IR n → IR n is the reaction term; u : IR d × IR → IR n is the solution; and D(u) is the diffusion matrix, which is a tensor of order d × d × n and the velocity vector, which is of order d × n.
In the present paper, we apply spatial discretization methods, such that we consider the spatial discretized partial differential equations with the included boundary conditions, which are given as ordinary differential equations: where u(0) = u 0 is the initial condition.A(u) and B(u) are operators, which are spatially discretized.For example, A(u) is the discretized in space convection operator and B(u) is the discretized in space diffusion operators.For convenience, the nonlinear operators are bounded, e.g., bounded matrices.
In our proposed scheme, we embed the idea of the -shift, which is explained in [14], to the iterative splitting methods, see [3].Then, we obtain a new so called -shifted iterative splitting method, which is a new contribution.Such shifted iterative splitting approach are used to design new adaptive time-splitting methods, which can be applied with maximal time-steps and compute the error-controls of the local splitting-errors.
In the following, we discuss the standard and the shifted splitting approaches.

Standard Splitting Approaches
In this section, we describe the standard splitting methods.

Strang-Marchuk Splitting (SMS)
In the SMS method, in the first step, the operator A is solved in the left half of interval [t n , t n+1 ]; then, in the second step, the operator B is solved in the whole interval [t n , t n+1 ]; and in the third step, the operator A is solved in the right half of the interval [t n , t n+1 ].By the initial conditions, the three subproblems are connected, see where τ = ∆(t n ) = t n+1 − t n is the time step.

Iterative Splitting Methods
The iterative splitting method defined in [13] and extensively studied for ordinary and partial differential equations in [3,10,15,19,20] are alternative operator splitting methods, which are based on iterative techniques.
We apply two versions of the iterative splitting methods: The LIS solves in the first equation the linear part of operator A with the given right hand side of operator B. Then, it solves in the second equation the linear part of operator B with the right hand side of operator A, using the solution of the first equation.The two solver steps are iterated m times before we pass to the next interval.
where i = 1, 2, . . ., m.For the initialisation of the iteration, we start with function u 0 (t), which verifies the initial condition u 0 (0) = u 0 .After, we have performed m iterations of the LIS, we apply the approximated solution u(t n+1 ) = u m (t n+1 ) for the next time-step, till the final step n + 1 = N.
• Quasilinear Iterative Splitting (QIS) The QIS solves in the first equation the nonlinear part of operator A with the given right hand side of operator B. Then, it solves in the second equation the nonlinear part of operator B with the right hand side of operator A, using the solution of the first equation.The two solver steps are iterated m times before we pass to the next interval.
where i = 1, 2, . . ., m.For the initialisation of the iteration, we start with function u 0 (t), which verifies the initial condition u 0 (0) = u 0 .After, we have performed m iterations of the QIS, we apply the approximated solution u(t n+1 ) = u m (t n+1 ) for the next time-step, till the final step n + 1 = N.

Shifted Splitting Approaches for Error Estimations
In this section, we describe the modified splitting methods to apply error-estimates.We propose shifted splitting as a novel method to design error estimates, see [14].

Shifted Strang-Marchuk Splitting (SSMS)
In this method, in interval [t n , t n+1 ], we first solve for operator A with a half time step minus a small , means τ/2 − , then we solve B with the full time-step τ, and again for A with a half time step plus a small , means τ/2 + .The three subproblems are connected by the initial conditions, according to The value is a small fraction of τ, for example = 0.005τ, so that the shifted interval is close to the original one.The error estimate is given as If the error estimate err is higher than a given tolerance η, we redo the computations with a smaller step according to the refinement scheme where we apply ν > 0, near 1 as a security factor.Otherwise, if err ≤ η, we accept the obtained value u Strang (t n+1 ), and proceed with the next time interval.In this case, in order to avoid the usage of unnecessary small time steps, we apply a coarsening scheme where κ is a small positive value depending on the tolerance η.
The Algorithm is given in Algorithm 1: Algorithm 1.

1.
Compute the local time-steps with the Strang and shifted Strang method, means u strang (t n+1 ) and u strang, (t n+1 ).

3.
If err > η, reject the time-step and restart the recent time-interval with the ∆t n = ∆t new obtained from (8).
If err ≤ η, then we are in the error tolerance.Proceed with the next time interval with the increased time step ∆t n+1 = ∆t new given by (9).
Remark 1.The theoretical results of the Algorithm 1, which is based on the -shifted Strang-splitting method, are given in the literature [4,5].Here, the authors applied the shifted Strang-splitting method and could apply a local error estimate of the first and second splitting resolution, see [14].
In the Figure 1, we have the graphically introduction of the shifting ideas.

Shifted Iterative Splitting Methods
In this section, we apply the shifting time-step ideas to the iterative splitting methods.First, we solve for the first iterative step with time-step τ − , then we solve for the second iterative step with time-step τ + .
We modify the two versions of the iterative splitting methods: The SLIS solves in the first equation the linearized part of operator A with the given right hand side of operator B for a minus shift of in the time-step.Then, it solves the second equation the linearized part of operator B with a given right hand side of operator A for a plus shift of in the time-step, using the solution of the first equation.The two solver steps are iterated m times before we pass to the next interval.
where i = 1, 2, . . ., m.For the initialization of the iteration, we start with function u 0 (t), which verifies the initial condition u 0 (0) = u 0 .After we have performed m iterations of the SLIS, we apply the approximated solution u(t n+1 ) = u m (t n+1 ) for the next time-step, until the final step n + 1 = N.
We apply the error-estimates as in Algorithm 2 and then we go to the next time-step.

•
Shifted Quasilinear Iterative Splitting (SQIS) The SQIS solves in the first equation the nonlinear part of operator A with the given right hand side of operator B for a minus shift of in the time-step.Then, it solves the second equation the nonlinear part of operator B with a given right hand side of operator A for a plus shift of in the time-step, using the solution of the first equation.The two solver steps are iterated m times before we pass to the next interval.
where i = 1, 2, . . ., m.For the initialization of the iteration, we start with function u 0 (t), which verifies the initial condition u 0 (0) = u 0 .After we have performed m iterations of the SQIS, we apply the approximated solution u(t n+1 ) = u m (t n+1 ) for the next time-step, till the final step n + 1 = N.
We apply the error-estimates as in Algorithm 2 and then we go to the next time-step.
Here we decided i = 1, but the error estimates also works for i = 1, 2, . . ., m.
The error estimate is given as whereas η is a given error tolerance, e.g., η = 10 −5 .Further the adaptive time-stepping is where we apply ν > 0, near 1 as a security factor.The Algorithm is given in Algorithm 2: Algorithm 2.

1.
We compute the local time-steps with the iterative and shifted iterative method, means u i (t n+1 ) and u i, (t n+1 ) 2.
We compute the error estimation If err ≤ η, then we are in the error tolerance and we accept the time-step means u(t n+1 ) = u i (t n+1 ) and the next time-step is ∆t n+1 = ∆t new .
Otherwise, we reject the time-step and restarted the recent time-interval with ∆t n = ∆t new In the Figure 2, we have the graphically introduction of the shifting ideas.
Unshifted iterative splitting and shifted iterative splitting method.

Error Analysis
The error analysis of the methods are done in the following.We deal with the following assumptions to the nonlinear operators: Assumption 1.
• Estimation of the nonlinear operators: ||B(e i (t where Ã and B are bounded operators, such that the Taylor expansion of the operators can be applied.

•
For the nonlinear operators A and B, we estimated the linearized parts as bounded operators Ã, B : X → X, where X is an appropriate Banach space.Further, we have a Banach-norm for the vector and the matrices, which is given as • .
In Theorem 1, we derive the consistency order of the shifted iterative operator-splitting.
Theorem 1.The operators A, B ∈ L(X) are nonlinear bounded operators with the Assumption 1.The abstract Cauchy problem is given as The abstract Cauchy problem (17) has an existent and unique solution.Then, the shifted iterative splitting method (11) is consistent with the order of the consistency O(τ 2i n ) with i = 1, . . ., m.
Proof.We assume Ã + B ∈ L(X) and we assume that the linear operators are generator of a uniformly continuous semigroup, such that we have a unique solution c(t) = exp(( Ã + B)t)c 0 .
In the following, we consider the local time-interval [t n , t n+1 ]. e i (t) = c(t) − ci (t) and e i+1 (t) = c(t) − c i (t) are the local error functions.
Based on the Assumptions 1, we can assume that the linearized operators Ã and B are generators of the one-parameter C 0 semigroup, which are given as (exp( Ã(t)) t≥0 and (exp( B(t)) t≥0 .
In the following, we can write the abstract Cauchy problem with homogeneous initial conditions as We apply the norms for the vectors and matrices and we can estimate We assume, that ( Ã(t)) t≥0 and ( B(t)) t≥0 are generators of semigroups and we apply the so-called growth estimation.Then, we can estimate exp( Ãt) ≤ K exp(ωt); t ≥ 0, exp( Ãt) ≤ K exp( ωt); t ≥ 0, (22) where the estimations are held for some numbers K ≥ 0 and ω, ω ∈ IR.
In the following, we distinguish between the following two operator-types.
Then, we have the following two estimates of the two groups of operators: and we apply the estimations to ( 21) and obtain the relation • Operators with exponential growth.
Here, we assume that (exp( Ãt)) t≥0 , (exp( Bt)) t≥0 are exponentially growing with some ω > 0, ω > 0. Therefore, we can estimate where Further, we apply The estimations ( 24) and (30) result in and we can apply recursively the error-estimation of Equations ( 31) and (32) and we obtain, Then, we recursively applied the Equation (33) and obtain the proved statement.
Remark 2. Based on the derivation of the error of the shifted iterative method, we obtain the error whereas err = C(∆t 2i−1 ) and also η = C(∆t 2i−1 new ), and we obtain Remark 3. In realistic applications, an optimal relation between the time-step τ n and the number of iterative steps 2i − 1, see Equation (35), is necessary.In practical experiments, we saw that ν > 0, but near to 1 and i ≈ 3, 4, 5 iterations are sufficient, such that we obtain an optimal new time-step ∆t new .To improve the criterion for stopping the iterative processes, we can additionally define an error bound.For example |c i − c i−1 | ≤ err with err = 10 −4 can be used to restrict us to an appropriate low number of iterative steps.
The order of accuracy can be improved by the choice of the initial iteration function, e.g., additional pre-steps with standard splitting approaches, see [13].
Based on our assumption about the initial solutions, we initialize with exact solutions or we apply higher order interpolated split solutions.This assumption allow to derive a theory for the exactness of the iterative methods, see also [13].

Numerical Results
In this section, we present the numerical results based on our novel iterative splitting methods for nonlinear ordinary and partial differential equation.We verified our theoretical results of the error-estimates and applied the shifted iterative splitting methods as a new-solver class.

First Numerical Example: Bernoulli Equation
In the first example, we apply a nonlinear differential equation, which is given as the Bernoulli equation, see For the Bernoulli equation, we can derive analytical solutions as reference solutions, see [15,22].
The analytical solutions are given as .
We apply the following operators for the splitting.
We apply backward Euler method to approximate the derivative in each subinterval [t n , t n+1 ], n = 0, 1, . . ., N, and solve the resulting equation by using the fixed point method and Newton's method with tolerance 10 −12 allowing a maximum of three iterations.The accuracy of the methods is assessed by comparing the numerical result u num with the analytical solution u given by (37).We compute the maximum and mean difference at the nodes t n , according to For the Shifted Strang-Marchuk splitting, we analyze the accuracy (with respect to the analytic solution) and the cost of the algorithm for different tolerances η and coarsening factors 1 + κ, where we have set κ = 4η 1 4 .Taking larger values of κ reduces the number of time intervals, but increases the number of tentative steps, where the ∆t n must be reduced in order to satisfy the error tolerance criterion.The value of κ has been experimentally chosen in order to minimize the total number of steps.
In each splitting step, the differential equation is approximated by back-Euler's method and solved using the fixed point (BEFP) or Newton's (BEN) methods.The cost of the algorithm and the final accuracy depend on the error tolerance η and the coarsening factor 1 + κ.We compare the shifted ABA-operator splitting method with the shifted variants of the iterative splitting methods above considered.
Table 1 shows that the accuracy is roughly proportional to the square root of the tolerance η, and the number of functional evaluations is inversely proportional to the same quantity.Newton's method requires less iterations to fulfill the tolerance; thus, if the number of time steps is similar to that of the fixed point method, it needs less functional evaluations.Nevertheless, Newton's method also evaluates the derivative, which reduces its advantage over the fixed point method.For the Shifted Linear Iterative Splitting method, we take κ = 4 √ η.We obtain similar accuracies to Strang-Marchuk's algorithm working now with higher error tolerances η, as it is shown in Table 2.The accuracy is of the same order as η, whereas the computational cost is slightly higher than the cost of Strang-Marchuk's algorithm.The results for the Shifted Quasilinear Iterative Splitting, see Table 3, are quite alike to the ones of the linear splitting.Increasing the number of iterations, iter = 2, 3, . .., results in a linear increment of the cost without any accuracy improvement.

Second Numerical Example: Mixed Convection-Diffusion and Burgers Equation
In the second numerical example, we apply coupled partial differential equation (PDE).We apply a coupling of a convection-diffusion equation with a Burgers' equation in 2D, which is called mixed convection-diffusion and Burgers equation (MCDB), and given as where the domains are given as Ω = [0, 1] × [0, 1] and T = 1.25.The viscosity is µ.
For such an mixed PDE, we can derive an analytical solution, which is where we can derive the right hand side f (x, y, t).

By considering the following operators
The MCDB Equation ( 38) is splitted into the Burgers' term, A and the convection-diffusion term, B and we obtain the operators: We deal with different viscosities: low viscosity µ = 0.5, high viscosity, µ = 5.The spatial domain is discretized taking a rectangular mesh with n x = n y = 16 intervals and applying standard second order divided difference approximations.The resulting differential system is solved by the same methods as in the previous example.The coarsening strategy applied here when err < η is where err = u i,j is the vector norm of the computed values in the nodes (x i , y j ) at each time step.
For the shifted Strang-Marchuk splitting method we take = 0.05 and different values of η.Table 4 shows the results of solving the equation with low viscosity µ = 0.5 using different tolerances.The solutions of the differential equations are approximated by using back-Euler fixed point method (BEFP) or back-Euler-Newton's method (BEN).Both methods perform similarly in cost and accuracy in this case.
The corresponding results for solving the equation with high viscosity, µ = 5, are shown in Table 5. BEFP requires much more time steps then BEN, but reaches more accuracy.For the linear and quasilinear shifted iterative splitting methods we take = 0.5 and the same coarsening strategy.Tables 6-9 show the cost and the accuracy in the low and high viscosity cases for the shifted linear and quasilinear iterative splitting methods using the back-Euler fixed point method and back-Euler Newton's method as solvers.
The shifted linear and quasilinear iterative splitting methods give similar results in all the considered cases.The behavior of the back-Euler fixed point method is worse in the low viscosity case than in the high viscosity case, as in the shifted Strang-Marchuk splitting method.

Third Numerical Example: Convection-Diffusion-Reaction Equation
In the third numerical example, we deal with a PDE, which is a convection-diffusion-reaction equation in 3D (CDR), see the example in [23]: u(x, y, z, t 0 ) = u 0 (x, y, z), (x, y, z) ∈ Ω, where we have v = (v x , v y , v z ) t , D ∈ IR 3 × IR 3 a diffusion matrix, u ∈ is the velocity field, k is a reaction parameter, and Ω = [0, 4] 3 × [t 0 , T], T = 10.0.We can have a special analytical solution for an instantaneous point source, which is given as:    We have the following parameters.
• instantaneous point source: where we initialise the equation with u 0 (x, y, z) = u ana (x, y, z, t 0 ), • the diffusion parameters are given as D 11 = 0.01, D 22 = 0.02, D 33 = 0.03 all other parameters are 0, the reaction parameter is given as k = 0.1.
By considering the following operators, we decouple into the fast velocity-reaction part and the slow diffusion parts: we split (50) in fast and slow parts The equation is spatially discretized taking a number, n x , n y , n z , of equal subintervals in each direction in Ω, and approximating the spatial derivatives by standard second order divided differences, resulting in a linear differential system.
We first check that the discretization error decreases with the size of the spatial subintervals by solving the differential system using Heun's method and Strang-Marchuk method with different number of spatial subintervals.The numerical results u num in the node points (x i , y j , z k ) are compared with the analytical solution u ana in the same points at the final time T = 10, computing the maximum and the mean absolute differences as before.Table 10 shows that there is no significant difference between both methods.Now we fix the number of spatial subintervals n x = n y = n z = 16, and analyze the performance of the adaptive methods for the CDR example.To estimate the convergence of the methods, we compare their results with the approximation obtained by integrating the differential equation by Heun's method using the same time steps.Table 11 shows the results of the shifted Strang-Marchuk splitting and the shifted linear iterative splitting for different tolerances, η.Lower tolerances produce lower maximum and mean errors but require more time steps.The relationship between the number of time steps and the mean error is depicted in Figure 3.  Remark 4. In the Figure 3, we see the differences in the convergence behaviour between the shifted Strang-Marchuk splitting (SSML) and the shifted linear iterative splitting (SLIS) method.We see in the figure, that the SSML method has only a linear convergence order ≤ 1, the SLIS method has higher order of convergence, here we have at least an order ≥ 2. Therefore the adaptive iterative scheme is much more effective and accurate than the noniterative splitting scheme.The result verified the proposition, that the iterative splitting scheme is a higher order scheme, see [3,15], and that we also conserve the higher order approach in the adaptive version.

Fourth Numerical Example: Nonlinear Diffusion Equation
Our fourth numerical example is a partial differential equation which is nonlinear diffusion equation in 2D, see the example in [24].
An application of such a nonlinear diffusion (NLD) is given by where we have α = 1 Further, we apply with the following parameters in the NLD Equations ( 41) and (42).
The parameters and the initial and boundary conditions are given as: • Uphill example, which is known as semi-degenerated Duncan and Toor experiment, see [25]: The computational domains are given with: Ω = [0, 1] is the spatial domain and [0, T] = [0, 1] is the time domain.

•
The initial conditions are as follows.

1.
Uphill example Asymptotic example • For the boundary conditions, we apply no-flux type conditions: We apply the following splitting of the operators with the one-dimensional spatial derivations: where we have the operators in the following decomposition of the u 1 and u 2 parts with ξ ∈ [0, 1]: where we have ξ = 0.5 a symmetric decomposition.We first check that the non adaptive methods require a very small time step to converge and estimate its convergence by comparing the results doubling successively the number of time steps.The results are shown in Table 12 for the direct integration and for the unshifted Strang-Marchuk method.The errors are computed measuring the difference between the result obtained with a given number of time steps and the result with twice that number at every shared temporal and spatial node.The error estimates for the Strang-Marchuk splitting in the case of 40,000 time steps are not available because the method diverges with 20,000 time steps.
Figure 4 illustrates the uphill phenomenon, where the solutions u 1 and u 2 increase before reaching the stationary state.
The adaptive methods result in an important reduction of the number of time steps obtaining similar error estimations.Tables 13 and 14 show the results for the uphill case and for the asymptotic case of the nonlinear diffusion equation, respectively.Here, the errors are computed by comparing the solution of the shifted methods with the ones obtained by direct integration, using the same time steps as the adaptive method.
The shifted Strang-Marchuk method behaves better for = 0.03, whereas the shifted linear and quasilinear splitting methods work well for = 0.03, except for in Table 15, where the behavior of the considered splitting methods is studied for different splitting weights ξ.

Conclusions and Discussion
We present a novel adaptive iterative splitting approach for partial differential equations of the type convection-diffusion-reaction equation.The numerical analysis shows the convergence of the schemes, while we could apply a shift in time of the methods.In the numerical experiments, we apply different state-of-the-art nonlinear convection-diffusion equations, where we receive benefits in the computational time and also in the accuracy of the methods.The adaptive splitting schemes allow to control the errors of the scheme and reduce the computational time, while we could apply smaller and larger time-steps.

Figure 3 .
Figure 3. Mean error e mean of the shifted Strang-Marchuk splitting (SSMS) and the shifted linear iterative splitting (SLIS) for the CDR equation for different tolerances, η.

Figure 5
depicts the regions in the space-time plane where the uphill phenomenon takes place, that is where N 2 and ∂ x u 2 have the same sign.The equation is solved by the shifted linear iterative splitting with η = 1.0 × 10 −5 .

Figure 4 .
Figure 4. Evolution of the magnitudes u 1 and u 2 in the uphill example.

Figure 5 .
Figure 5. Regions in the space-time domain where N 2 and ∂ x u 2 have the same sign.

Table 2 .
Shifted linear iterative splitting method for Bernoulli's equation.

Table 4 .
Solution of mixed convection-diffusion and Burgers (MCDB) equation for µ = 0.5 using the shifted Strang-Marchuk splitting method.

Table 5 .
Solution of MCDB equation for µ = 5 using the shifted Strang-Marchuk splitting method.

Table 6 .
Results of the shifted linear iterative splitting for the MCDB equation with µ = 0.5.

Table 7 .
Results of the shifted linear iterative splitting for the MCDB equation with µ = 5.

Table 8 .
Results of the shifted quasi linear iterative splitting for the MCDB equation with µ = 0.5.

Table 9 .
Results of the shifted quasi linear iterative splitting for the MCDB equation with µ = 5.

Table 10 .
Differences between the analytical solution of convection-diffusion-reaction (CDR) equation and the numerical solutions obtained by direct (Heun) integration and by Strang-Marchuk method with different number of spatial subintervals.

Table 11 .
Results of the shifted Strang-Marchuk's splitting, (SSMS), and the shifted linear iterative splitting, SLIS, for the CDR equation for different tolerances.

Table 12 .
Cost and error estimates of the Heun integration (HI) and the Strang-Marchuk (SM) splitting for the nonlinear diffusion equation.