Computing Exponential for Iterative Splitting Methods: Algorithms and Applications

. Iterative splitting methods have a huge amount to compute matrix exponential. Here, the acceleration and recovering of higher-order schemes can be achieved. From a theoretical point of view, iterative splitting methods are at least alternating Picards ﬁx-point iteration schemes. For practical applications, it is important to compute very fast matrix exponentials. In this paper, we concentrate on developing fast algorithms to solve the iterative splitting scheme. First, we reformulate the iterative splitting scheme into an integral notation of matrix exponential. In this notation, we consider fast approximation schemes to the integral formulations, also known as φ - functions. Second, the error analysis is explained and applied to the integral formulations. The novelty is to compute cheaply the decoupled exp-matrices and apply only cheap matrix-vector multiplications for the higher-order terms. In general, we discuss an elegant way of embedding recently survey on methods for computing matrix exponential with respect to iterative splitting schemes. We present numerical benchmark examples, that compared standard splitting schemes with the higher-order iterative schemes. A real-life application in contaminant transport as a two phase model is discussed and the fast computations of the operator splitting method is explained.


Introduction
We are motivated to solving multiple phase problems that arose of transport problem in porous media.In the last years, the interest in numerical simulations with multiphase problems, that can be used to model potential damage events has significantly increased in the area of final repositories for chemical or radioactive waste.
Here, the modeling of the underlying porous media, which is the geosphere, spooled with water, is important, and we apply mobile and immobile pores as a two-phase problem in the media.With such a model, we achieve more realistic scenarios of the transported species, see 1 .The transport is structured in a convection and a diffusion portion.The convection portion is determined by the velocity vector v and indicates the spatial direction where the initial value is given as c i,0 and we assume Dirichlet boundary conditions with the function c i,1 x, t sufficiently smooth, all other initial and boundary conditions of the other phases zero.
φ: effective porosity − , c i : concentration of the ith gaseous species in the plasma chamber phase mol/cm 3 , c i,im : concentration of the ith gaseous species in the immobile zones of the plasma chamber phase mol/cm 3 , c i,ad : concentration of the ith gaseous species in the mobile adsorbed zones of the plasma chamber phase mol/cm 3 , c i,im,ad : concentration of the ith gaseous species in the immobile adsorbed zones of the plasma chamber phase mol/cm 3 , v: velocity through the chamber and porous substrate 4 cm/nsec , D e i : element-specific diffusions-dispersions tensor cm 2 /nsec , λ i,i : decay constant of the ith species 1/nsec , Q i : source term of the ith species mol/ cm 3 nsec , g: exchange rate between the mobile and immobile concentration 1/nsec , k α : exchange rate between the mobile and adsorbed concentration or immobile and immobile adsorbed concentration kinetic controlled sorption 1/nsec , E: stationary electric field in the apparatus V/cm , η: the mobility rate in the electric field, see 5 cm 2 /nsec .
with i 1, . . ., M and M denotes the number of components.
The parameters in 1.1 are further described, see also 1 .
The effective porosity is denoted by φ and declares the portion of the porosities of the aquifer, that is filled with plasma, and we assume a nearly fluid phase.The transport term is indicated by the Darcy velocity v, that presents the flow-direction and the absolute value of the plasma flux.The velocity field is divergence-free.The decay constant of the ith species is denoted by λ i .Thereby does k i denote the indices of the other species.
In this paper we concentrate on solving linear evolution equations, such as the differential equation, where A can be an unbounded operator.We assume to achieve large-scale differential equation, which are delicate to solve with standard solvers.
Our main focus will be to consider and contrast higher-order algorithms derived from standard schemes as for example Strang-splitting schemes.
We propose iterative splitting schemes as a solver scheme which is simple to implement and parallelisible.
A rewriting to integral formulation, allows to reduce the computation to numerical approximation schemes.While standard schemes has the disadvantage to compute commutators to achieve higher-order schemes, we could speed up our schemes by recursive algorithms.Iterative schemes can be seen as Successive approximations, which are based on recursive integral formulations in which an iterative method is enforce the time dependency.
The paper is outlined as follows.In Section 2, we present the underlying splitting methods.The algorithms for the exponentials are given in Section 3. Numerical verifications are given in Section 4. In Section 5, we briefly summarize our results.

Iterative Splitting Method
The following algorithm is based on the iteration with fixed-splitting discretization step-size τ, namely, on the time interval t n , t n 1 we solve the following subproblems consecutively for i 0, 2, . . ., 2m. cf.6, 7 : where c n is the known split approximation at the time level t t n .The split approximation at the time level t t n 1 is defined as c n 1 c 2m 1 t n 1 , clearly, the function c i 1 t depends on the interval t n , t n 1 , too, but, for the sake of simplicity, in our notation we omit the dependence on n .
In the following we will analyze the convergence and the rate of convergence of the method 2.1 for m tends to infinity for the linear operators A, B : X → X, where we assume that these operators and their sum are generators of the C 0 semigroups.We emphasize that these operators are not necessarily bounded, so the convergence is examined in a general Banach space setting.
The novelty of the convergence results are the reformulation in integral-notation.Based on this, we can assume to have bounded integral operators which can be estimated and given in a recursive form.Such formulations are known in the work of 8, 9 and estimations of the kernel part with the exponential operators are sufficient to estimate the recursive formulations.

Error Analysis
We present the results of the consistency of our iterative method.We assume for the system of operator the generator of a C 0 semigroup based on their underlying graph norms.
where A, B : D X → X are given linear operators which are generators of the C 0 -semigroup and c 0 ∈ X is a given element.We assume A, B are unbounded.Further, we assume the estimations of the unbounded operator B with sufficient smooth initial conditions (see [8]):

2.3
Further we assume the estimation of the partial integration of the unbounded operator B (see [8]): Then, we can bound our iterative operator splitting method as where S i is the approximated solution for the ith iterative step and C is a constant that can be chosen uniformly on bounded time intervals.
Proof.Let us consider the iteration 2.1 on the subinterval t n , t n 1 .
For the first iterations, we have

2.9
We have the following solutions for the iterative scheme: the solutions for the first two equations are given by the variation of constants:

2.10
For the recurrence relations with even and odd iterations, we have the solutions: for the odd iterations: i 2m 1, for m 0, 1, 2, . ..,

2.12
The consistency is given as the following.
For e 1 , we have

2.13
We obtain

2.14
For e 2 we have alternating

2.15
We obtain

2.16
For odd and even iterations, the recursive proof is given in the following.For the odd iterations only iterations over A , i 2m 1 for m 0, 1, 2, . .., for e i , we have:

2.17
We obtain where i is the number of iterative steps.
The same idea can be applied to the even iterative scheme and also for alternating A and B.
Remark 2.2.The same idea can be applied to A ∇D∇ B −v • ∇, so that one operator is less unbounded, but we reduce the convergence order and hence, Remark 2.3.If we assume the consistency of O τ m for the initial value e 1 t n and e 2 t n , we can redo the proof and obtain at least a global error of the splitting methods of O τ m−1 .

Splitting Method to Couple Mobile, Immobile, and Adsorbed Parts
The motivation of the splitting method are based on the following observations.
i The mobile phase is semidiscretized with fast finite volume methods and can be stored into a stiffness-matrix.We achieve large time steps, if we consider implicit Runge-Kutta methods of lower order e.g., implicit Euler as a time discretization method.
ii The immobile, adsorbed and immobile-adsorbed phases are purely ordinary differential equations and each of them is cheap to solve with explicit Runge-Kutta schemes.
iii The ODEs can be seen as perturbations and can be solved all explicit in a fast iterative scheme.
For the full equation we consider the following matrix notation:

2.20
where c c 1 , . . ., c m T is the spatial discretized concentration in the mobile phase, see 1.1 , T is the concentration in the immobile phase, the some also for the other phase concentrations.A 1 is the stiffness matrix of 1.1 , A 2 is the reaction matrix of the righthand side of 1.1 , B 1 and B 2 are diagonal matrices with the exchange of the immobile and kinetic parameters, see 1.4 and 1.5 .Further, Q, . . ., Q im,ad are the spatial discretized sources vectors.Now, we have the following ordinary differential equation: where C c, c im , c ad , c im,ad T and the right-hand side is given as For such an equation, we apply the decomposition of the matrices:

2.23
The equation system is numerically solved by an iterative scheme.
We start with n 0.
1 The initial conditions are given with C 0 t n 1 C t n .We start with k 0.
2 Compute the fix-point iteration scheme given as where k is the iteration index, see 10 .For the time integration, we apply Runge-Kutta methods as ODE solvers, see 11, 12 .
3 The stop criterion for the time interval t n , t n 1 is given as where • is the maximum norm over all components of the solution vector.err is a given error bound, for example, err 10 −4 .
If 2.25 is fulfilled, we have the result

2.26
If n N then we stop and are done.
If 2.25 is not fulfilled, we do k k 1 and go to 2 .
The error analysis of the schemes are given in the following Theorem.
Theorem 2.5.Let A, B ∈ L X be given linear bounded operators in a Banach space L X .We consider the abstract Cauchy problem:

2.27
where t 1 0 and the final time is t N T ∈ Ê .Then problem 2.27 has a unique solution.For a finite steps with time size τ n t n 1 − t n , the iteration 2.24 for k 1, 2, . . ., q is consistent with an order of consistency O τ q n .
Proof.The outline of the proof is given in 2 .
In the next section we describe the computation of the integral formulation with expfunctions.

Exponentials to Compute Iterative Schemes
The theoretical ideas can be discussed in the following formulation: and the matrix formulation of our two-step scheme is given as the computation of the exp-Matrix can be expressed as: Here, we have to compute the right-hand side as time dependent term, means we evaluate exp At and exp Bt as a Taylor expansion.Then the integral formulation can be done with polynomials and we derive the expansion of such integral operators with the φfunction.
The φ-functions are defined as the integration of the exp-functions.We can analytically derive φ-functions with respect to the integral-formulation of a matrix exp-function.
In the following we reduce to a approximation of the fixed right-hand side means we assume Later we also follow with more extended schemes.

3.4
So the matrix formulation of our scheme is given as where A is given as, the computation of the exp-Matrix can be expressed in a first-order scheme and with the assumption of commutations is given as: exp At : where we assume u i s u t n , mean we approximate the right-hand side term.For higher-orders we should also include the full derivations of c 1 , c 2 , . ...

Derivation of the Formula
Consider the equation ẋ Ax a, x 0 0, 3.8 the solution of this equation in terms of the φ k , k 0, 1, is given as x tφ 1 tA a.

3.9
Similarly, the equation ẋ Ax bt, x 0 0, 3.10 has a solution of the form x t 2 φ 2 tA b.

3.11
In general, the solution of the equation ẋ Ax a bt c t 2  2! • • • , x 0 0, 3.12 admits the form In this section, we use the formula 3.21 for iterative schemes given as

3.14
The solution of the first iteration given by u 1 e At u 0 .

3.15
Inserting this into 3.26 and expanding e At up to the order 1, we have 2nd order approximation of the exact solution which has a form u 2 e Bt u 0 φ 1 tA Atu 0 .

3.16
Similarly, inserting 3.16 into 3.27 and expanding φ 1 tA , we have a third order approximation of the exact solution

3.17
In general for i 0, 2, 4, . .., we have the p i 1-th order approximation of the exact solution in terms of the φ i function as follows:

Algorithms
In the following, we discuss the one-and two-side algorithms.

Two-Side Scheme
For the implementation of the integral formulation, we have to deal with parallel ideas, which means we select independent parts of the formulation.
Step 1. Determine the order of the method by fixed iteration number.
Step 2. Consider the time interval t 0 , T , divide it into N subintervals so that time step is h T − t 0 /N.
Step 3. On each subinterval, t n , t n h , n 0, 1, . . ., N, use the algorithm by considering initial conditions for each step as u t 0 u 0 , u i t n u i−1 t n u t n ,

3.20
Step 4. u i t n h → u t n h .
Step 5. Repeat this procedure for next interval until the desired time T is reached.

One-Side Scheme (Alternative Notation with Commutators)
For the one-side scheme, we taken into account of the following commutator relation.

3.29
The recursion is given as

3.30
Remark 3.2.For the novel notation, we have embedded the commutator to the computational scheme.For such a scheme we could save to compute additional the commutators.

Numerical Examples
In the following, we deal with numerical example to verify the theoretical results.

First Example: Matrix Problem
For another example, consider the matrix equation

4.2
We split the matrix as The Figure 1 presents the numerical errors between the exact and the numerical solution.
The Figure 2 presents the CPU time of the standard and the iterative splitting schemes.

Second Experiment: 10 × 10 Matrix
We deal in the first with an ODE and separate the complex operator in two simpler operators.We deal with the 10 × 10 ODE system:  where λ 1 t ∈ Ê and λ 2 t ∈ Ê are the decay factors and u 1,0 , . . ., u 10,0 ∈ Ê .We have the time interval t ∈ 0, T .We rewrite 4.4 in operator notation, we concentrate us to the following equations: where u 1 0 u 10 1.0, u 2 0 u 20 1.0 are the initial conditions, where we have λ 1 t t and λ 2 t t 2 .The operators are splitted in the following way, while operator A, see 4.10 , covers the upper diagonal entries and operator B, see 4.11 , covers the lower diagonal entries of the ODE system 4.4 -4.7 .Such a decomposition allows to reduce the computation of the matrix exponential to a upper or lower matrix and speedup the solver process The Figure 3 present the numerical errors between the exact and the numerical solution.Here we have the best approximation to the exact solution with a sixth-order iterative scheme 6 iterative steps , but no additionally more computational time for finer time step.
The Figure 4 present the CPU time of the standard and the iterative splitting schemes.Compared to standard splitting schemes, the iterative schemes are cheaper to compute and higher in accuracy.We can also choose smaller time steps to obtain more accurate results and use nearly the same computational resource.
Remark 4.1.The computational results show the benefit of the iterative schemes.We save computational time and achieve higher-order accuracy.The one-side and two-side schemes have the same results.

Third Example: Commutator Problem
We assume to have a large norm of the commutator A, B and deal with the matrix equation In the following, we discuss different splitting ideas.

Version 1
We split the matrix as  Here we have to deal with highly noncommutative operators and the computational speedup is given in the iterative schemes, while the commutator is not needed to obtain more accurate results, while the standard splitting schemes deal with such commutators for the error reduction.
The Figure 5 present the numerical errors between the exact and the numerical solution.
The Figure 6 present the CPU time of the standard and the iterative splitting schemes.Further for the one-side,we obtain more improved results for the following splitting.

Version 2
In this version, we have

4.14
The Figure 7 presents the numerical errors between the exact and the numerical solution.
A more delicate problem is given for the stiff matrices.

Version 3
In this version, we have

4.15
The Figure 8 present the numerical errors between the exact and the numerical solution.
The Figure 9 present the CPU time of the standard and the iterative splitting schemes.
Remark 4.2.The iterative schemes with fast computations of the exponential matrices have a speedup.The constant CPU time of the iterative schemes shows that it benefit instead of the expensive standard schemes.Also for stiff problems with multi iterative steps, we reach the same results of the standard A-B or Strang-Splitting schemes, while we are independent of the commutator estimation.

Two-Phase Example
The next example is a simplified real-life problem for a multiphase transport-reaction equation.We deal with mobile and immobile pores in the porous media, such simulations are given for waste scenarios.
We concentrate on the computational benefits of a fast computation of the iterative scheme, given with matrix exponentials.
The equation is given as: In the following, we deal with the semidiscretized equation given with the matrices: .17 where C c 1 , c 2 , c 1im , c 2im T , while c 1 c 1,1 , . . ., c 1,I is the solution of the first species in the mobile phase in each spatial discretization point i 1, . . ., I , the same is also for the other solution vectors.
We have the following two operators for the splitting method

4.18
where I is the number of spatial points.

4.20
For the operator A 1 and A 2 A 2 A 3 , we apply the iterative splitting method.Based on the decomposition, operator A 1 is only tridiagonal and operator A 2 is block diagonal.Such matrix structure reduce the computation of the exponential operators.
The Figure 10 present the numerical errors between the exact and the numerical solution.Here we obtain optimal results for one-side iterative schemes on operator B, means we iterate with respect to B and use A as right-hand side.Remark 4.3.For all iterative schemes, we can reach faster results as for the The iterative schemes with fast computations of the exponential matrices standard schemes.With 4-5 iterative steps we obtain more accurate results as we did for the expensive standard schemes.With one-side iterative schemes we reach the best convergence results.

Conclusions and Discussions
In this work, we have presented a very economical and practical method by successive approximations.Here the idea to decouple the expensive computation of matrix exponential via iterative splitting schemes has the benefit of less computational time.While the error analysis present stable methods and higher-order schemes, the applications show the speedup with the iterative schemes.In, future, we concentrate on matrix dependent scheme, based on iterative splitting algorithms.

Figure 1 :
Figure 1: Numerical errors of the standard splitting scheme a and the iterative schemes b with 1, . . ., 6 iterative steps.

Figure 2 :
Figure 2: CPU time of the standard splitting scheme a and the iterative schemes b with 1, . . ., 6 iterative steps.

Strang c 1 c 2 c 3 c 4 c 5 c 6 ABStrang c 1 c 2 c 3 c 4 c 5 c 6 AB,Figure 3 :Figure 4 :
Figure 3: Numerical errors of the standard splitting scheme a and the iterative schemes b with 1, . . ., 6 iterative steps.
while A, B ≥ max{ A , B }.

Figure 5 :Figure 6 :
Figure 5: Numerical errors of the standard splitting scheme a and the iterative schemes b with 1, . . ., 6 iterative steps.

Figure 7 :
Figure 7: Numerical errors of the standard splitting scheme a and the iterative schemes based on one-side to operator A or B b and two-side scheme c .

Figure 8 :
Figure 8: Numerical errors of the one-side splitting scheme with A a , the one-side splitting scheme with B b and the iterative schemes with 1, . . ., 6 iterative steps c .

Figure 9 :
Figure 9: CPU time of the one-side splitting scheme with A a and the iterative schemes with 1, . . ., 6 iterative steps b .

Figure 10 :
Figure 10: Numerical errors of the one-side splitting scheme with A a , the one-side splitting scheme with B b and the iterative schemes with 1, . . ., 6 iterative steps c . 2