New efficient substepping methods for exponential timestepping

Exponential integrators are time stepping schemes which exactly solve the linear part of a semilinear ODE system. This class of schemes requires the approxima- tion of a matrix exponential in every step, and one successful modern method is the Krylov subspace projection method. We investigate the effect of breaking down a single timestep into arbitrary multiple substeps, recycling the Krylov subspace to minimise costs. For these recyling based schemes we analyse the lo- cal error, investigate them numerically and show they can be applied to a large system with 106 unknowns. We also propose a new second order integrator that is found using the extra information from the substeps to form a corrector to increase the overall order of the scheme. This scheme is seen to compare favourably with other order two integrators.


Introduction
We consider the numerical integration of a large system of semilinear ODEs of the form with u, F (t, u(t)) ∈ R N and L ∈ R N ×N a matrix.Equation (1) arises, for example, from the spatial discretisation of reaction-diffusion-advection equations.An increasingly popular method for approximating the solution of semlinar ODE systems such as (1) are exponential integrators.These are a class of schemes which approximate (1) by exactly solving the linear part and are characterised by requiring the evaluation or approximation of a matrix exponential function of L at each timestep.A major class of exponential integrators are the multistep Exponential Time Differencing (ETD) schemes, first developed in [1], other classes include the Exponential Euler Midpoint method [2] and Rosenbrock type methods [3,4].For an overview of exponential integrators see [5,6] and other useful references can be found in [7].
Our schemes are based on the standard exponential integrator ETD1, which can be written as where ϕ 1 (∆tL) is defined shortly; u etd n ≈ u(t n ) at discrete times t n = n∆t for fixed ∆t > 0, F etd n ≡ F (t n , u etd n )0 and n ∈ N. ETD1 is globally first order, and is derived from (1) by variation of constants and approximating F (t, u(t)) by the constant F etd n over one timestep.See for example [5,6,7,21] for more detail.It is useful to introduce the additional notation g(t) ≡ Lu(t) + F (t, u(t)) and g etd n ≡ Lu etd n + F (t n , u etd n ).
The function ϕ 1 is part of a family of matrix exponential functions defined by ϕ 0 (z) = e z , ϕ 1 (z) = z −1 (e z − I) , and in general where I is the identity matrix.These ϕ−functions appear in all exponential integrator schemes; see [19].In particular we use ϕ 1 , and for brevity we introduce the following notation p τ ≡ τ ϕ 1 (τ L).
We can then re-write (2) as We consider the Krylov projection method for approximating terms like p ∆t g etd n in (2).In the Krylov method, this term is approximated on a Krylov subspace defined by the vector g n and the matrix L. Typically the subspace is recomputed, in the form of a matrix of basis vectors V m , every time the solution vector, u n in (2), is updated (and thus also g n ).This is done using a call to the Arnoldi algorithm [16], and is often the most expensive part of each step.It is possible to 'recycle' this matrix at least once, as demonstrated in [22] for the exponential Euler method (EEM) (see [5] and (A.1)).In this paper we investigate this possibility further and use it to construct new methods based on ETD1 and in Appendix A we show how to construct the general recycling method for EEM.
We examine the effect of splitting the single step of (2) of length ∆t in to S substeps of length δt = ∆t S , through which the Krylov subspace and matrices are recycled.By deriving expressions for the local error, we show that the scheme remains locally second order for any number S of substeps, and that the leading term of the local error decreases.This gives a method based on recyling the Krylov subspace for S substeps.We then obtain a second method using the extra information from the substeps to form a corrector to increase the overall order of the scheme.
The paper is arranged as follows.In Section 2 we describe the Krylov subspace projection method for approximating the action of ϕ−functions on vectors.In Section 3 we describe the concept of recycling the Krylov subspace across substeps in order to increase the accuracy of the ETD1 based scheme, and show that the leading term of the local error of the scheme decreases as the number of substeps uses increases.We then prove a lemma to express the local error expression at arbitrary order.With this information about the local error expansion, and the extra information from the substeps taken, it is possible to construct correctors for the scheme the increase the accuracy and local order of the scheme.We demonstrate one simple such corrector in Section 4. Numerical examples demonstrating the effectiveness of this scheme are presented in Section 5.

The Krylov Subspace Projection Method and ETD1
We describe the Krylov subspace projection method for approximating ϕ 1 (∆tL) in (2).We motivate this by showing how the leading powers of ∆tL in L are captured by the subspace.The series definition of ϕ 1 (∆tL) is, The challenge in applying the scheme (2) is to efficiently compute, or approximate, the action of ϕ 1 on the vector g etd n .The sum in (6) is useful in motivating a polynomial Krylov subspace approximation.The m-dimensional Krylov subspace for the matrix L and vector g ∈ R N is defined by: Approximating sum in ( 6) by the first m terms is equivalent to approximation in the subspace K m (L, g etd n ) in (7).We now review some simple results about the general subspace K m (L, g), with arbitrary vector g, before using the results with g = g etd n to demonstrate how they are used in the evaluation of (2).The Arnoldi algorithm (see e.g.[19], [16]) is used to produce an orthonormal basis {v 1 , . . ., v n } for the space K m (L, g) such that It produces two matrices V m ∈ R N ×m , whose columns are the v k , and an upper Hessenburg matrix H m ∈ R m×m .The matrices L,H m and V m are related by see, for example, equation ( 2) in [16], of which ( 9) is a consequence.From (9) it follows that, For any x represents the orthogonal projection into the space K m (L, g).Therefore, since L k g ∈ K m (L, g n ), we also have that We now consider the relationship between L k g and Proof.We have that m , the projector into K(L, g), so that (10) can be more briefly written V m H m V T m = πLπ.Then by (10) we find Now consider using the vector g = g etd n , to generate the subspace K m (L, g etd n ), and the corresponding matrices H m , V m , by the Arnoldi algorithm.By Lemma 2.1 we have that, up to k = m, Thus, inserting the approximation the first m terms are correctly approximated.The Krylov approximation is then Let us introduce a shorthand notation for the Krylov approximation of the ϕ−function.Analogous to (4), for Using ( 13) and ( 14) we then approximate (5) by u etd n+1 = u etd n + p∆t g etd n .
The key here is that the ϕ 1 (∆tH m ) now needs to be evaluated instead of ϕ 1 (∆tL).m is chosen such that m ≪ N , and a classical method such as a rational Padé is used for ϕ 1 (∆tH m ), which would be prohibitively expensive for ϕ 1 (∆tL) for large N .One step of the ETD1 scheme (2), under the approximation ϕ 1 (∆tL) ≈ V m ϕ 1 (∆tH m )V t m , becomes where e 1 is the first basis vector in R m .

Recycling the Krylov subspace
In the Krylov subspace projection method described in §2, the subspace K m (L, g n ) and thus the matrices H m and V m depend on g n .At each step it is understood that a new subspace must be formed, and H m , V m be regenerated by the Arnoldi method, since g n changes.In [22] it is demonstrated that splitting the timestep into two substeps, and recycling H m and V m , i.e. recycling the Krylov subspace, can be viable (in that it does not decrease the local order of the scheme, and apparently decreases the error).We expand on this concept with a more detailed analysis of the effect of this kind of recycled substepping applied to the locally second order ETD1 scheme (2) (EEM is considered in Appendix A).We replace a single step of length ∆t of (15) with S substeps of length δt, such that ∆t = Sδt.We denote the approximations used in this scheme analogously to the notation for ETD1 earlier, without the etd superscript, and introduce the following notation to keep track of substeps and ).
For j = 1 we calculate H m , V m , from g n , and for the remaining S − 1 steps, where the matrices H m and V m are not re-calculated for any substep, j > 1.
We call substeps of the form (17) 'recycled steps' and substeps of the form (16) 'initial steps'.The approximation to u(t n + ∆t) at the end of the step of length ∆t is then given by The recycling steps ( 16), ( 17) can be succinctly expressed using the definition of pτ ;

The local error of the recycling scheme
We now derive an expression for the local error of the scheme defined by ( 16), (17) and prove that the leading term decreases with the number of substeps S. We use the local error assumption that u n = u(t n ) and make an assumption about the accuracy of the initial Krylov approximation with respect to the error of the scheme.Let the error in the polynomial Krylov approximation over a single step (including the error from the approximation of ϕ 1 (τ H m ) using, e.g.Padé approximation), with subspace dimension of size m, be given by E m n+1 , so that, pτ Assumption 3.1.The Krylov approximation error E m n+1 is much less than the error of ETD1, and thus does not affect the leading term of the local error of ETD1.
Bounds on E m n+1 can be found in for example [17,16].Practically, we can always reduce ∆t or increase m until Assumption 3.1 is satisfied.For the local error of the recycling scheme, the following result will be used.Lemma 3.2.For any τ 1 , τ 2 ∈ R, and any vector v ∈ R N , and the same relation holds for the Krylov approximations, that is, Proof.We prove the second equation.The first can be proved using an almost identical argument, replacing pτ by p τ where appropriate.By the definitions of pτ , and ϕ By (9) this becomes pτ2 (Lp which is pτ1+τ2 v as desired. Without recycling substeps, a single ETD1 step (2) of length ∆t, using the polynomial Krylov approximation, would be: To examine the local error we compare u etd n+1 with the u n+1 obtained after some number S of recycled substeps.We can write , where R S n+1 represents the deviation from ( 22) over one step.Then we have: Lemma 3.3.The approximation u n+ j S produced by j substeps of the recycling scheme ( 19), (20), satisfies with Proof.By induction.For j = 1, u n+ 1 S is given by ( 19) and R S n+ 1 (23) holds for some j ≥ 1.Then u n+ j+1 S is obtained by a step of (20).Using (23) we find, and since Lu n = g n − F (t n , u(t n )), by the induction hypothesis, Thus by Lemma 3.2 we have that, To complete the proof we need to show: By the induction hypothesis that (24) holds for j, Hence the lemma is proved.
Using (23) we now express the leading order term of the local error in terms of S. First we examine the leading order term of R S n+1 .
Lemma 3.4.The term R S n+ j S in Lemma 3.3, when expanded in powers of ∆t, satisfies Proof.By induction.Since ( 27) is true for j = 1.Now assume the result holds for some j.Then we can express the term F n+ j S follows: We thus have that We then insert (28) into the inductive expression (27) for R S n+ j S and use the Using the induction assumption (27), Noting that ∆t = Sδt we can write O(∆t 3 ) as O(δt 3 ).Collecting terms we have, The lemma follows since j(j−1) The leading local error term of the ETD1 scheme without substeps is well known to be ∆t 2 2 dF (t) dt (see [23]), so that we can finally recover the leading term from Lemma 3.3.
Corollary 3.5.The leading term of the recycling scheme after j steps is Corollary 3.6.The local error u(t n + ∆t) − u n+1 of an ETD1 Krylov recycling scheme is second order for any number S of recycled substeps.Moreover, the local error after j recycled steps is In particular or in terms of ∆t It is interesting to compare (31) with the leading term of the local error of regular ETD1, ∆t 2 part in (31) is the projection of the ETD1 error into K, multiplied by a factor S−1 S ≤ 1.Thus, in the leading term, according to (31), the recycling scheme reduces the error of ETD1 by effectively eliminating the part of the error which lives in K.In the limit S → ∞, the entirety of the error in K will be eliminated.The effectiveness of the recycling scheme therefore depends on how much of df (t)  dt can be found in K. Corollary 3.6 shows that using S > 1 recycled substeps is advantageous over the basic ETD1 scheme, in the sense of reducing the magnitude of the leading local error term, whenever where ||•|| is a given vector norm.We show in Lemma 3.8 that increasing S will decrease the Euclidean norm || • || 2 of the leading term of the local error.First we require a result on V m V T m , the projector into the Krylov Supspace K.
Remark 3.7.Let x = 0 be a vector such that Proof.An elementary result for orthogonal projectors (see, e.g.[24]) is that ) and the definition of the Euclidean norm.Equation (33) is a generalisation of (34) as can be shown as follows.
x, and then, noting that Using (34) to substitute for ( dt = 0. Let E S1 be the local error using the recycling scheme over a timestep of length ∆t with S 1 ≥ 1 substeps, and E S2 the local error with S 2 substeps with S 2 > S 1 .Then, Proof.The local errors E S k , k = 1, 2 are given in Corollary 3.6.
Note that we have that V m V T m x = 0 from the assumptions.This is because, dt is already entirely within K.Then, We have that 1 To prove the lemma we apply (33) to x, with γ in place of α.If we have that From Lemma (3.8) we see that S recycled Krylov substeps not only maintains the local error order of the ETD1 scheme, but also decreases the 2-norm of the leading term with increasing S. Note that the leading term does not tend towards zero as S → ∞, but towards a constant.We thus expect diminishing returns in the increase in accuracy with increasing S, and the existence of an optimal S for efficiency.

Using the additional substeps for correctors
We now establish a new second order scheme based on a finite difference approximation to the derivative of the nonlinear term F and the recycling scheme given in ( 16), (17).
The first step is to expand the local error for the standard ETD1 scheme.Using variation of constants and a Taylor series expansion of F (t, u(t)), the exact solution of (1) can be expressed as a power series (see for example [5,23]) with Under the local error assumption u n = u(t n ), the local error of the ETD1 step given in ( 2) is Since the approximation from a substepping scheme is related to the approximation from the ETD1 scheme (over one step) by u n+1 = u etd n+1 + R S n+1 , we have the local error for the recycling scheme: The terms of error expression (38) at arbitrary order can be found using (36), Lemma 3.3, and the information on Krylov projection methods in §2.We see that the expansion consists of terms involving the value of F or derivatives thereof at various substeps.These terms can be approximated by finite differences of the values for F at the different substeps, and used as a corrector to eliminate terms for the error.We consider extrapolation in the leading error in the case of two substeps, that is S ≡ 2. Assume that the error from the Krylov approximation, E m n is negligible compared to E n and R n , so that it does not introduce any terms at the first and second and third order expansion of E n and R n .Then we can express exactly the leading second and third order error terms.First we have the leading terms of E etd n from (37), We also have the leading terms of R 2 n+1 (from two substeps, recall ( 24)) Note that the terms in (40) are an order higher than written since ).We then have that The idea now is as follows.Define a corrected approximation: In (42), C is a corrector intended to cancel out some of the leading terms in (41).The term ∆t 2 V m V T m (F n+ 1 2 − F n ) is the only leading term in (41) to involve the matrix V m , and so is added directly to the corrected approximation (42) to allow C to be free of dependence on the matrix V m .Indeed, C will be a linear combination of the the three function values of F , F n , F n+ 1  2 and F n+1 , available at the end of the full step.The approximation to u produced by substeps of the scheme, and thus also to F , is locally second order.We define the C term as follows, with coefficients α, β, γ to be chosen later.
where we have used that F n = F (t n , u(t n )) (under the local error assumptions), ) and so on.The new term ∆t 3 E c is introduced to represent the O(∆t 3 ) error in writing ∆tF n+ 1  2 as ∆tF t n + ∆t 2 , and so on.From (43), we must choose the coefficients to satisfy the two conditions α + β + γ = 0, and With these values of the parameters, the local error of the corrected approximation is We have three coefficients to determine, and two constraints.We are therefore in a position to pick another constraint to reduce the new leading error in (44).
It would be helpful to know the form of the error term E c , introduced by the approximation of F in (43).We have: using Corollary 3.6.We also have Substituting into (44) , − ∆t 3 ∂F ∂u (47) We have the option here to use the final constraint to eliminate the coefficient of F (2) in the leading term: Note that E c cannot be eliminated without taking the inverse of V m V T m , so this is not an efficient option.It can be seen that the values that satisfy the three constraints are: Of course E c also depends on the values of α, β, γ, so the magnitude of the third order term will be affected by the choice of these values also through E c .
With the choices given above, we have the numerical scheme and the E c term in (47) becomes: Here we have used all the extra information from the two substeps to completely eliminate the lowest order from the local error, and a part of the new leading order term for the scheme.A more thorough use of the error expressions in the lemmas here may give rise to recycling schemes that use more substeps and are able to completely eliminate higher order terms from the error, leading to a kind of new exponential Runge-Kutta framework involving recycled Krylov subspaces.Below we demonstrate the efficacy of our two-step corrected recycling scheme with numerical examples.In Appendix A we show how to apply the analysis of the substepping method to the locally third order exponential integrator scheme EEM.

Numerical Results
Here we examine the performance of the recyling scheme (19), (20) and the corrector scheme (48) (for the first two examples).The PDEs investigated in these experiments are all advection-diffusion-reaction equations, which are converted into semilinear ODE systems (1) by spatial disretisation before our timestepping schemes are applied, see for example, [21] for more details.We use the v(x, t) ∈ R to represent the solution to the PDE, while u(t) ∈ R N represents the solution of the corresponding ODE system.The spatial discretisation is a simple finite volume method in all the examples.In examples 2 and 3, the grid was using code from MRST [25].We compare the second order corrector scheme (48) to both the standard second order exponential integrator (ETD2; refer to for example Equation ( 6) in [1]) and standard second order Rosenbrock scheme (ROS2; the same as the simplest Rosenbrock scheme described in [3,4]).For the first two experiments, the error is estimated by comparison with a low ∆t comparison solve u comp with ETD2.Our ETD2 implementation uses phipm.m[19] for each timestep; which requires the following parameters: an initial Krylov subspace dimension m, and an error tolerance.For our comparison solve runs, we used m = 30 and 10 −7 respectively for these parameters.The first two experiments are also found in [21]; see this for more details.For the third experiment a comparison solve was prohibitively time consuming, so error was instead estimated by differencing successive results.We estimate the error in a disrecte aproximation of the L 2 (Ω) norm, where D is the computational domain.

Allen-Cahn type Reaction Diffusion
We approximate the solution to the PDE, The (1D) spatial domain for this experiment was Ω = [0, 100] This was discretised into a grid of N = 100 cells.We imposed no flow boundary conditions, i.e., ∂u ∂x = 0 where x = 0 or x = 100.There was a uniform diffusivity field of D(x) = 1.0.The initial condition was u(x, 0) = cos 2πx N and we solved to a final time T = 1.0.
In Figure A.1 a) and c) we show the estimated error against ∆t, for the recycling scheme with varying number of substeps S, (S = 1, 2, 5, 10, 50, 100).Note that S = 1 is the standard ETD1 integtrator.The behaviour is as expected; increasing S decreases the error and the scheme is first order.The diminishing returns of increasing S (see (31)) can also be observed; for example compare the significant increase in accuracy in increasing S from 1 to 5, with the lesser increase in accuracy in increasing S from 5 to 10. Figure A.1 shows this more emphatically -the increase in accuracy in increasing S from 10 to 50 is significant, but the effect of increasing S from 50 to 100 is very small.The limiting value of the error with respect to S discussed above is clearly close to being reached here.In Figure A.1 b) and d) we plot estimated error against cputime to demonstrate the efficiency of the scheme with varying S. In this case increasing S appears to increase efficiency until an optimal S is reached, after which it decreases, as predicted.Figure A.1 d) shows that the optimal S lies between 50 and 100 for this system.In Figure A.2 we examine the 2-step corrector (48).Plot a) shows estimated error against ∆t.The corrector scheme is second order as intended, and has quite high accuracy compared to the other two schemes, possibly due to the heuristic attempt to decrease the error in the leading term (see discussion in §4).In plot b) we see that the 2-step corrector is of comparable efficiency to ROS2.
In ) indicates that the second order, 2-step corrector method can produce error more than one order of magnitude smaller than the first order recycling scheme with S = 10.

Fracture system with Langmuir-type reaction
We approximate the solution to the PDE, where D(x) is the diffusivity and V (x) is the velocity.In this example a single layer of cells is used, making the problem effectively two dimensional.The domain is Ω = 10 × 10 × 10 metres, divided into 100 × 100 × 1 cells of equal size.We impose no-flow boundary conditions on every edge.The initial condition imposed is initial v(x) = 0 everywhere except at x = (4.95,9.95) T where v(x) = 1.The diffusivity D in the grid varies with x, in a way intended to model a fracture in the medium.A subset of the cells in the 2D grid were chosen to be cells in the fracture.These cells were chosen by a weighted random walk through the grid (weighted to favour moving in the positive y-direction so that the fracture would bisect the domain).This process started on an initial cell which was marked as being in the fracture, then randomly chose a neighbour of the cell and repeated the process.We set the diffusivity to be D = 100 on the fracture and D = 0.1 elsewhere.There is also a constant velocity V field in the system, uniformly one in the x-direction and zero in the other directions in the domain, i.e., v(x) = (1, 0, 0) T , to the right in For sufficiently low ∆t we have the predicted results, with the error being first order with respect to ∆t, and decreasing as S increases.For ∆t too large, this is not the case.Here the Krylov subspace dimension m is most likely the limiting factor as Assumption 3.1 becomes invalid.In Figure A.4 b) and d) we show the efficiency by plotting the estimated error against cputime.For ∆t low enough that the substepping schemes are effective, the scheme with 10 substeps is the most efficient.We can see the existence of an optimal S for efficiency, as predicted, in Figure A.4 d), where the scheme using S = 50 is more efficient than the scheme using S = 100.Any increase in accuracy by increasing S from 50 to 100 is extremely small (indeed, it is unnoticeable in Figure A.4 c), and not enough to offset the increase in cputime.In fact, Figure A.4 d) shows that for this experiment the scheme using S = 10 is more efficient than both the S = 50 and S = 100 schemes.Figure A.4 c) shows that the S = 10 scheme is also slightly more accurate than both.This is likely because at S 10 the improvement in accuracy is already close to the limiting value, and greatly increasing S to 50 or 100 only accumulates rounding errors without further benefit.Figure A.4 a) shows that the improvement from S = 1 to S = 10 is quite significant on its own.
In shows that the second order corrector scheme can be almost three orders of magnitude more accurate for a fixed cputime than the first order recycling scheme with S = 10.

Large 2D Example with Random Fields
In this example the 2D grid models a domain with physical dimensions 100× 100 × 10; the grid is split into a 1000 × 1000 × 1 cells.The model equation is the same as the previous example (50).The diffusivity is kept constant at D = 0.01, while a random velocity V field is used.For this we generated a random permativity field K, which was then used to generate a corresponding pressure field and then a velocity field in a standard way, using Darcy's Law, see [21,25].The pressure p field was determined by the permativity field and the Dirichlet boundary conditions p = 1 where y = 0 and p = 0 where y = 100.The initial conditions for v were zero everywhere, and the boundary conditions were the same as for p, v = 1 where y = 0 and v = 0 where y = 100.The final time was T = 500.
Due to the large size of system (10 6 unknowns) we only examine the recycling scheme and for a system of this size, it was necessary to increase m to 100 to prevent the Krylov error from being dominant.The results are shown in Figure A.7.We see that, for ∆t sufficiently low, increasing S decreases the error and increases the efficiency of the scheme.The improvement in efficiency between S = 5 and S = 10 is marginal; the optimal S for this example would not be much greater than 10.

Conclusions
We have extended the notion of recycling a Krylov subspace for increased accuracy in the sense of [22].We have applied this new method to the first order ETD1 scheme and examined the effect of taking an arbitrary number S of substeps.The local error has been expressed in terms of S, and the expression shows that the local error will decrease with S down to a finite limit.The discussion in Appendix A examines construction for EEM.Results suggest that there maybe an optimal S for a maximal efficiency increase and some preliminary analysis in this direction may be fuond in [21].Convergence and existence of an optimal S > 1 has been demonstrated with numerical experiments.Additional information from the substeps was used to form a corrector and a second order scheme.This was shown to be comparable to, or slightly better than, ETD2 and ROS2 in our tests.
The schemes currently rely on Assumption 3.1, essentially requiring that ∆t be sufficiently small and m be sufficiently large, to be effective.Numerical experiments have shown how having ∆t too large can cause the schemes to become inaccurate as the error of the initial Krylov approximation becomes significant.It is already well established how the Krylov approximation error can be controlled by adapting m and the use of non-recycling substeps.Applying these techniques to the schemes presented here in future work would allow them to be effective over wider ∆t ranges.
We now combine the leading term of the remainder R and the known local error of EEM 1 6 ∆t 3 g T Ĵg.
(see, for example, [21]) to find the local error of the new recycling scheme.
Corollary Appendix A.3.The leading term of the local error of the S step recycling scheme for EEM at the end of a timestep is From this we can predict similar properties to the ETD1 recycling scheme.This extends the work of [22], where the recycling substepping EEM scheme was used for a single substep.
Figure A.1 b) we see that for the same cputime, increasing S from 1 to 10 decreases the estimated error by roughly an order of magnitude.We can see in Figure A.1 d) that increasing S from 10 to 50 can further decrease error for a fixed cputime, though less significantly.Comparing a fixed cputime in Figure A.1 b) and Figure A.2 b

Figure A. 3 .
The initial condition was c(x) = 0 everywhere except at x = (4.95,9.95) T where c(x) = 1.In Figure A.3 we show the final state of the system at time T = 2.4.The result in plot a) was produced with the 2-step recycling scheme with a timestep ∆t = 2.4 × 10 −4 .Plot b) shows the high accuracy comparison ETD2 solve, produced with ∆t = 2.4 × 10 −5 .In Figure A.4 we demonstrate the effect of increasing the number of substeps S on the error.Figure A.4 a) shows estimated error against timestep ∆t, for schemes using S = 1, 2, 5, 10 substeps, while Figure A.4 c) shows the same for schemes using S = 10, 50, 100.Recall that S = 1 is the standard ETD1 integtrator.

Figure A. 5
we compare the 2-step corrector scheme against the two other second order exponential integrators, ETD2 and ROS2.Figure A.5 a)  shows estimated error against ∆t, and we see that, like Figure A.4 a), the Krylov recycling scheme does not function as intended above a certain ∆t threshold; again this is due to the timestep being too large with respect to m.The standard exponential integrators do not have this problem, as their timesteps are driven by phipm.m,which takes extra (non-recycled, linear) substeps to achieve a desired error.Below the ∆t threshold, the 2-step corrector scheme functions exactly as intended, exhibiting second order convergence and high accuracy.In Figure A.5 b) we can observe that the 2-step corrector scheme is more efficient than the other two schemes for lower ∆t, and of comparable efficiency for larger ∆t.It is interesting to compare Figure A.4 a) and Figure A.5 a) and note the threshold ∆t for the corrector scheme seems to be lower than for the substepping schemes.In Figure A.4 b) we can again see that for a fixed cputime, increasing S from 1 to 10 decreases error by roughly one order of magnitude; however Figure A.4 d) shows no improvement in increasing S from 10 to 50.Comparing Figure A.4 b) and Figure A.5 b)

Figure A. 6 :
Figure A.6: Result for the example in §5.3, in which solute enters through the lower boundary and flows according to a random velocity field.Produced by the 10-step recycling scheme with ∆t = 0.2441, i.e. 2048 steps.a) Shows the system at the final time T = 500; the axes indicates the physical dimensions (i.e., the domain is 100 × 100 metres).b) Is a corresponding contour plot for clarity; the axes indicate the cells in the finite volume grid (i.e., the grid has 1000 cells along each side).

Figure A. 7 :
Figure A.7: Results for the substepping scheme applied to the large Langmuir type reaction system in §5.3.(a) Estimated errors against timestep ∆t and (b) displays estimated error against cputime, showing efficiency.Note that for the largest timestep the error is dominated by the Krylov error as m is too small for the given ∆t (c.f.Assumption 3.1 )