Journal of Applied Mathematics

Volume 2011 (2011), Article ID 193781, 27 pages

http://dx.doi.org/10.1155/2011/193781

## Computing Exponential for Iterative Splitting Methods: Algorithms and Applications

Humboldt-Univeristät zu Berlin, 10099 Berlin, Germany

Received 29 November 2010; Accepted 27 January 2011

Academic Editor: Shuyu Sun

Copyright © 2011 Jürgen Geiser. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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 fix-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.

#### 1. 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 and indicates the spatial direction of the concentrations with the velocity. The diffusion portion is the spatial diffusion process of the concentrations, see [1].

The sorption allegorizes the exchange between the solute (mobile) pollutant and at the surface sorbed (immobile) pollutants and appears in the temporal term as well as in the reaction term. The reaction is reversible. The equilibrium-sorption therefore can be declared as a coefficient in the specific terms.

We concentrate on simplified models and taken into account all real conditions for achieving statements of such realistic simulations of potential damage event. Here we have the delicate problem of coupled partial and ordinary differential equations and their underlying multiscale problems. While we deal with splitting methods and decomposing such scale problems, we can overcome such scaling problems, see [2].

Moreover the computational part is important, while we dealing with large matrices and standard splitting schemes are expensive with respect to compute the exponential matrices. We solve the computational problem with a novel iterative splitting scheme, that concentrate on developing fast algorithms to solve an integral formulation of matrix exponentials. In such a scheme, we consider fast approximation schemes to the integral formulations, also known as -functions, see [3]. The novelty is to compute cheaply the decoupled matrices and apply only cheap matrix-vector multiplications for the higher-order terms. In general, we discuss an elegant way of embedding recent survey on methods for computing matrix exponential with respect to iterative splitting schemes.

In the following, we describe our model problem. The model equation for the multiple phase equations are given as coupled partial and ordinary differential equations:
where the initial value is given as and we assume Dirichlet boundary conditions with the function sufficiently smooth, all other initial and boundary conditions of the other phases zero.:effective porosity ,:concentration of the th gaseous species in the plasma chamber phase (mol/cm^{3}),:concentration of the th gaseous species in the immobile zones of the plasma chamber phase (mol/cm^{3}),:concentration of the th gaseous species in the mobile adsorbed zones of the plasma chamber phase (mol/cm^{3}),:concentration of the th gaseous species in the immobile adsorbed zones of the plasma chamber phase (mol/cm^{3}),:velocity through the chamber and porous substrate [4] (cm/nsec),:element-specific diffusions-dispersions tensor (cm^{2}/nsec),:decay constant of the th species (1/nsec),:source term of the th species (mol/(cm^{3}nsec)),:exchange rate between the mobile and immobile concentration (1/nsec),:exchange rate between the mobile and adsorbed concentration or immobile and immobile adsorbed concentration (kinetic controlled sorption) (1/nsec),:stationary electric field in the apparatus (V/cm),:the mobility rate in the electric field, see [5] (cm^{2}/nsec).

with and 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 , that presents the flow-direction and the absolute value of the plasma flux. The velocity field is divergence-free. The decay constant of the th species is denoted by . Thereby does denote the indices of the other species.

In this paper we concentrate on solving linear evolution equations, such as the differential equation, where 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.

#### 2. Iterative Splitting Method

The following algorithm is based on the iteration with fixed-splitting discretization step-size , namely, on the time interval we solve the following subproblems consecutively for . (cf. [6, 7]): where is the known split approximation at the time level . The split approximation at the time level is defined as , (clearly, the function depends on the interval , too, but, for the sake of simplicity, in our notation we omit the dependence on ).

In the following we will analyze the convergence and the rate of convergence of the method (2.1) for tends to infinity for the linear operators , where we assume that these operators and their sum are generators of the 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.

##### 2.1. Error Analysis

We present the results of the consistency of our iterative method. We assume for the system of operator the generator of a semigroup based on their underlying graph norms.

Theorem 2.1. *Let one consider the abstract Cauchy problem in a Hilbert space **
where are given linear operators which are generators of the -semigroup and is a given element. We assume , are unbounded. Further, we assume the estimations of the unbounded operator with sufficient smooth initial conditions (see [8]):
**Further we assume the estimation of the partial integration of the unbounded operator (see [8]): **Then, we can bound our iterative operator splitting method as
**
where is the approximated solution for the th iterative step and is a constant that can be chosen uniformly on bounded time intervals.*

* Proof. *Let us consider the iteration (2.1) on the subinterval .

For the first iterations, we have
and for the second iteration we have:

In general, we have:

(i)for the odd iterations, for
where for ,(ii)for the even iterations, for

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

For the recurrence relations with even and odd iterations, we have the solutions: for the odd iterations: , for ,

For the even iterations, , for * The consistency is given as the following. *

For , we have

We obtainFor we have (alternating)

We obtain

For odd and even iterations, the recursive proof is given in the following. For the odd iterations (only iterations over ), for , for , we have:

We obtain
where is the number of iterative steps.

The same idea can be applied to the even iterative scheme and also for alternating and .

*Remark 2.2. *The same idea can be applied to , so that one operator is less unbounded, but we reduce the convergence order
where , .

*Remark 2.3. *If we assume the consistency of for the initial value and , we can redo the proof and obtain at least a global error of the splitting methods of .

##### 2.2. 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: where is the spatial discretized concentration in the mobile phase, see (1.1), is the concentration in the immobile phase, the some also for the other phase concentrations. is the stiffness matrix of (1.1), is the reaction matrix of the right-hand side of (1.1), and are diagonal matrices with the exchange of the immobile and kinetic parameters, see (1.4) and (1.5).

Further, are the spatial discretized sources vectors.

Now, we have the following ordinary differential equation: where and the right-hand side is given as .

For such an equation, we apply the decomposition of the matrices: where

The equation system is numerically solved by an iterative scheme.

*Algorithm 2.4. *We divide our time interval into subintervals , where , , and .

We start with .

(1)The initial conditions are given with . We start with .(2)Compute the fix-point iteration scheme given as
where 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 is given as
where is the maximum norm over all components of the solution vector. is a given error bound, for example, .If (2.25) is fulfilled, we have the result
If then we stop and are done.If (2.25) is not fulfilled, we do and go to (2).

The error analysis of the schemes are given in the following Theorem.

Theorem 2.5. *Let be given linear bounded operators in a Banach space . We consider the abstract Cauchy problem:
**
where and the final time is . Then problem (2.27) has a unique solution. For a finite steps with time size , the iteration (2.24) for is consistent with an order of consistency .*

* Proof. *The outline of the proof is given in [2].

In the next section we describe the computation of the integral formulation with -functions.

#### 3. 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 and 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 -functions. We can analytically derive -functions with respect to the integral-formulation of a matrix -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.

*Computing -Functions:*

So the matrix formulation of our scheme is given as
where 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:
where we assume , mean we approximate the right-hand side term.

For higher-orders we should also include the full derivations of .

##### 3.1. Derivation of the Formula

Consider the equation the solution of this equation in terms of the , is given as Similarly, the equation has a solution of the form In general, the solution of the equation admits the form In this section, we use the formula (3.21) for iterative schemes given as The solution of the first iteration given by Inserting this into (3.26) and expanding up to the order 1, we have 2nd order approximation of the exact solution which has a form Similarly, inserting (3.16) into (3.27) and expanding , we have a third order approximation of the exact solution In general for , we have the -th order approximation of the exact solution in terms of the function as follows: For , we have

##### 3.2. Algorithms

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

##### 3.3. 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 , divide it into subintervals so that time step is .

*Step 3. *On each subinterval, , , use the algorithm by considering initial conditions for each step as , ,

*Step 4. *.

*Step 5. *Repeat this procedure for next interval until the desired time is reached.

##### 3.4. One-Side Scheme (Alternative Notation with Commutators)

For the one-side scheme, we taken into account of the following commutator relation.

*Definition 3.1. *The following relation is given:
where is the matrix commutator.

The integration is given as
the representation is given in [13].

Further, we have the recursive integration
where given defined as

The iterative scheme with the equations is solved as

The recursion is given as

*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.

#### 4. Numerical Examples

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

##### 4.1. First Example: Matrix Problem

For another example, consider the matrix equation the exact solution is 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.

##### 4.2. Second Experiment: Matrix

We deal in the first with an ODE and separate the complex operator in two simpler operators.

We deal with the ODE system: where and are the decay factors and . We have the time interval .

We rewrite (4.4) in operator notation, we concentrate us to the following equations: where , are the initial conditions, where we have and .

The operators are splitted in the following way, while operator , see (4.10), covers the upper diagonal entries and operator , 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.

##### 4.3. Third Example: Commutator Problem

We assume to have a large norm of the commutator and deal with the matrix equation

In the following, we discuss different splitting ideas.

*Version 1*

We split the matrix as
while .

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
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
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 or Strang-Splitting schemes, while we are independent of the commutator estimation.

##### 4.4. 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: where , while is the solution of the first species in the mobile phase in each spatial discretization point , the same is also for the other solution vectors.

We have the following two operators for the splitting method where is the number of spatial points.

We decouple into the following matrices:

For the operator and , we apply the iterative splitting method.

Based on the decomposition, operator is only tridiagonal and operator 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 , means we iterate with respect to and use 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.

#### 5. 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.

#### References

- J. Geiser,
*Discretisation methods for systems of convective-diffusive dispersive-reactive equations and applications*, Ph.D. thesis, Universität Heidelberg, Heidelberg, Germany, 2004. - J. Geiser, “Decomposition methods for partial differential equations: theory and applications in multiphysics problems,” in
*Numerical Analysis and Scientific Computing Series*, Magoules and Lai, Eds., CRC Press, Boca Raton, Fla, USA, 2009. View at Google Scholar - M. Hochbruck and A. Ostermann, “Explicit exponential Runge-Kutta methods for semilinear parabolic problems,”
*SIAM Journal on Numerical Analysis*, vol. 43, no. 3, pp. 1069–1090, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - H. Rouch, “MOCVD research reactor simulation,” in
*Proceedings of the COMSOL Users Conference*, Paris, France, 2006. - L. Rudniak, “Numerical simulation of chemical vapour deposition process in electric field,”
*Computers and Chemical Engineering*, vol. 22, no. 1, supplement, pp. S755–S758, 1998. View at Google Scholar - R. Glowinski, “Numerical methods for fluids,” in
*Handbook of Numerical Analysis*, P. G. Ciarlet and J. Lions, Eds., vol. 9, North-Holland, Amsterdam, The Netherlands, 2003. View at Google Scholar - J. F. Kanney, C. T. Miller, and C. T. Kelley, “Convergence of iterative split-operator approaches for approximating nonlinear reactive problems,”
*Advances in Water Resources*, vol. 26, no. 3, pp. 247–261, 2003. View at Publisher · View at Google Scholar - E. Hansen and A. Ostermann, “Exponential splitting for unbounded operators,”
*Mathematics of Computation*, vol. 78, no. 267, pp. 1485–1496, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - T. Jahnke and C. Lubich, “Error bounds for exponential operator splittings,”
*BIT Numerical Mathematics*, vol. 40, no. 4, pp. 735–744, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - I. Faragó and J. Geiser, “Iterative operator-splitting methods for linear problems,”
*International Journal of Computational Science and Engineering*, vol. 3, no. 4, pp. 255–263, 2007. View at Publisher · View at Google Scholar - E. Hairer, S. P. Nørsett, and G. Wanner,
*Solving ordinary differential equations. I*, vol. 8 of*Springer Series in Computational Mathematics*, Springer, Berlin, 1987, Nonstiff problem. - E. Hairer and G. Wanner,
*Solving Ordinary Differential Equations. II*, vol. 14 of*Springer Series in Computational Mathematics*, Springer, Berlin, Germany, 2nd edition, 1996. - W. Magnus, “On the exponential solution of differential equations for a linear operator,”
*Communications on Pure and Applied Mathematics*, vol. 7, pp. 649–673, 1954. View at Publisher · View at Google Scholar · View at Zentralblatt MATH