International Journal of Differential Equations

Volume 2013, Article ID 314290, 11 pages

http://dx.doi.org/10.1155/2013/314290

## Embedded Zassenhaus Expansion to Splitting Schemes: Theory and Multiphysics Applications

Institute of Physics, Felix-Hausdorff-Street 6, 17489 Greifswald, Germany

Received 25 March 2013; Revised 5 August 2013; Accepted 19 August 2013

Academic Editor: Shuyu Sun

Copyright © 2013 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

We present some operator splitting methods improved by the use of the Zassenhaus product and designed for applications to multiphysics problems. We treat iterative splitting methods that can be improved by means of the Zassenhaus product formula, which is a sequential splitting scheme. The main idea for reducing the computation time needed by the iterative scheme is to embed fast and cheap Zassenhaus product schemes, since the computation of the commutators involved is very cheap, since we are dealing with nilpotent matrices. We discuss the coupling ideas of iterative and sequential splitting techniques and their convergence. While the iterative splitting schemes converge slowly in their first iterative steps, we improve the initial convergence rates by embedding the Zassenhaus product formula. The applications are to multiphysics problems in fluid dynamics. We consider phase models in computational fluid dynamics and analyse how to obtain higher order operator splitting methods based on the Zassenhaus product. The computational benefits derive from the use of sparse matrices, which arise from the spatial discretisation of the underlying partial differential equations. Since the Zassenhaus formula requires nearly constant CPU time due to its sparse commutators, we have accelerated the iterative splitting schemes.

#### 1. Introduction

Our motivation to study the operator splitting methods comes from models in fluid dynamics, for example problems in bioremediation [1] or radioactive contaminants [2]. Such multiphysics problems are delicate, and the solver methods can be accelerated by decoupling the different physical behaviours; see [3, 4]. Based on the splitting error of all the splitting schemes, we have to take into account higher order ideas; see [5, 6]. While standard splitting methods deal with lower order convergence, see Lie splitting and Strang splitting [7, 8], we propose a combination of iterative splitting methods with an embedded Zassenhaus product formula. Such a combination allows reducing the splitting error and also reducing the CPU time needed; see [9]. Theoretically, we combine fixed point schemes (iterative splitting methods) with sequential splitting schemes (Zassenhaus products), which are connected with the theory of Lie groups and Lie algebras. Based on that relation, we can construct higher order splitting schemes for an underlying Lie algebra and improve the convergence results with cheap iterative schemes.

Historically, the efficiency of decoupling different physical processes into simpler processes, for example, convection and reaction, has helped to accelerate the solution process and is discussed in [10].

We propose the following ideas.(i)Iterative splitting schemes are based on fixed point schemes, for example, waveform relaxation, which linearly improve the convergence order. Based on reducing the integral operators to cheaply computable matrices, they can be seen as solver methods; see [9]. (ii)Zassenhaus formula is based on nested commutators, which are the main keys to deriving higher order standard splitting schemes (e.g., Lie-Trotter and Strang splitting). They are simple to compute because of their nilpotent structure; see [11].In this paper we study the following ordinary differential equations, that can arose of semidiscretized partial differential equations, for example, multiphysics problem, see [2, 3], or theoretical physics; see [12, 13].

We focus our attention on the case of two linear operators (i.e., we consider the Cauchy problem) as follows: whereby the initial function is given, and and are assumed to be bounded linear operators in the Banach space with . In realistic applications, the operators correspond to discretised matrices, for example, such as convection and diffusion matrices; see [2]. Further, they can be given matrices, for example, such as in perturbation theory [14].

For all such problems, the solution of the Cauchy problem (1) is given as where is the initial condition.

Here one of the main ideas to compute such delicate solutions is to separate or decouple the term into simpler and parallel computational expression, for example, ; see [10, 15]. To reduce the splitting error for such decompositions, iterative splitting method and products of exponential of noncommuting operators, see [16, 17], are important. We concentrate on iterative splitting methods and improve their convergence order with embedded Zassenhaus expansion; see [18].

While iterative splitting methods are cheap to implement and can be computed fast, one of their drawbacks is due to their initial process. Here, Zassenhaus expansion is an attractive tool to overcome such a drawback and improve the initial solutions to higher accuracy. While on the other hand, we deal only with lower nested commutators of the Zassenhaus formula, the drawbacks of stringent regularity conditions on the operators are less; see [19, 20].

The outline of the present paper is as follows. The splitting methods are discussed in Section 2. In Section 3, we present the numerical experiments and the benefits of the higher order splitting methods. Finally, we discuss future work in the area of iterative and noniterative splitting methods.

#### 2. Splitting Methods

We consider the following splitting schemes: (i)iterative splitting schemes, (ii)Zassenhaus formula.Their ideas and their underlying convergence analysis, based on the Cauchy problem (1), are discussed in the next subsections.

##### 2.1. Iterative Operator Splitting

In the following, we develop two ideas of the iterative splitting schemes.

Iterative splitting with respect to one operator (also called a one-sided scheme) as follows: Iterative splitting with respect to two operators (also called a two-sided scheme) as follows:

Theorem 1. *Let one considers the abstract Cauchy problem given in (1). One obtains for the one-sided iterative operator splitting method (4) the following accuracy:
**
where is the approximate solution for the th iterative step and is a constant that can be chosen uniformly on bounded time intervals.*

*Proof. *The proof is given for and with the consistency error , 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. *The accuracy of the initialisation is important for conserving or improving the underlying iterative splitting scheme.

Here we have the following initialisation schemes:
So the slow convergence in the initial process of the iterative splitting schemes can be improved by a cheap computable initial process . We will next discuss the fast convergent Zassenhaus formula for the initial process.

##### 2.2. Zassenhaus Formula (Sequential Splitting)

The Zassenhaus formula is an extension to the exponential splitting schemes and can be given in Theorem 3; see also [20].

Theorem 3. *One assumes is the free Lie algebra generated by and . The kernel function of solution (1), , can be uniquely decomposed as
**
where is a homogeneous Lie polynomial of and of degree .*

*Proof. *The sketch of the proof is given in [20]. The idea is based on the existence of the BCH (Baker-Campbell-Hausdorff) formula, and we obtain
where is a Lie polynomial of degree >1, and
Further we have
where is a Lie polynomial of degree >2, and
We can repeat the recursive process and obtain by complete induction the results.

The first Zassenhaus polynomials are given as Here we see the benefit of such schemes in the computation of the commutators. If we assume that we have quickly computable commutators, for example, nilpotent matrices or sparse matrices, the computation time needed for the adjacent (or perturbed) operators of the scheme is much less than that of the main operators , and so the scheme is very effective and can be used as an initial solution of the iterative schemes.

*Remark 4. *We can also generalize the application of the Zassenhaus formula to decomposition methods with existing BCH formulas, for example,(i)- splitting or Lie-Trotter splitting, (ii)Strang splitting.

*Example 5. *(1)* A-B or Lie-Trotter Splitting*. For the Lie-Trotter splitting, there exists a BCH formula, and we have derived the polynomials in Theorem 3.

(2)* Strang Splitting*. There exists also a BCH formula, which is given as
where the is a homogeneous Lie polynomial in and of degree ; see [21].

The BCH formula is given as
and so there exists a new product based on the Zassenhaus polynomials as follows:
and we obtain a higher order scheme; see also [22].

*Remark 6. *The application of the Zassenhaus formula to partial differential equations, for example, second-order parabolic evolution equations, needs additional stringent regularity conditions on the operators. Due to the nested commutators, which are needed to apply the Zassenhaus formula, one needs at least higher regularity conditions.

We need some functional analytical results on the domains of the operators. The domains of the operators have to satisfy the following compatibility conditions for the first Zassenhaus exponents. (i) Condition for the first term of the Zassenhaus formula
where , and is the domain of the operators. (ii) Condition for the second term of the Zassenhaus formula
where , and is the domain of the operators.

*Example 7. *If we assume to apply the Zassenhaus formula to a second-order diffusion equation and we assume to split each spatial direction, see the functional analytical framework [23], where the domains are given as
where .

Then we need the following regularity result: (i) condition for the first term of the Zassenhaus formula
(ii) condition for the second term of the Zassenhaus formula
In the example, we need a very high regular domain to apply a first or second Zassenhaus exponent. This is one of the limitations and drawbacks of applying the Zassenhaus formula. In our applications, we restrict on applying the first or at least the second Zassenhaus exponent. We apply the third and fourth Zassenhaus exponents to the benchmark problems given as matrix equations without such restrictions to the domains; see [20].

##### 2.3. Embedding the Zassenhaus Expansion into an Iterative Splitting Scheme

We now discuss the embedding of the Zassenhaus formula into the iterative operator splitting schemes.

*Definition 8. *We apply the Lie-Trotter splitting with the embedded Zassenhaus formula, which is given as
where the sequential Zassenhaus operators are defined as
and we obtain an accuracy of , and is the number of Zassenhaus components.

In the following theorem, we present the improvement of the iterative splitting scheme with the Zassenhaus expansion.

Theorem 9. *One solves the initial value problem (1). We assume bounded and constant operators , .**The initialisation process is done with the Zassenhaus formula
**
where is given in (25).**Furthermore, the improved solutions are embedded in the iterative splitting scheme (4) and one has after iterative steps the following result:
**
This has improved the error of the iterative scheme to .*

*Proof. *The solution of the iterative splitting scheme (4) is
where , given in Theorem 1.

The initialisation is given by the Zassenhaus formula
Combining both splitting schemes, we have a local error of

*Remark 10. *We obtain an acceleration of the iterative scheme by applying cheap higher order Zassenhaus formula. Such a benefit allows reducing to 2-3 iterative steps, while the initial solution is of a higher order. Taking into account the sparsity or nilpotency of the commutators, we now have a fast iterative splitting scheme.

*Remark 11. *For our numerical experiment, we apply the Zassenhaus formulaat the beginning of each integration step, as the so-called preprocessor. Therefore, we obtain more accurate starting conditions for the integration steps of the iterative splitting scheme. Based on the regularity conditions of the nested commutator, we apply only 1-2 Zassenhaus exponents.

In the following, we explain the final algorithm of the embedded Zassenhaus formula to the iterative splitting scheme.

*Algorithm 12. *We compute time steps, where each time step is given as with . Further is the initial condition of the Cauchy problem (1), and we compute successively. We have the following two steps; first we apply the Zassenhaus formula as initialization of for the iterative scheme; second we apply the iterative splitting scheme to the initial value problem (1) and compute the final solution as follows. (1) We start with and compute the Zassenhaus formula. (2) We compute the matrix norm for the Zassenhauss exponents . If , then we start with the iterative splitting scheme with , and go to (3), else and we go to (1). (3) If the error of the iterative scheme , then we are finished; else , and we compute the next iterative step, till .

#### 3. Numerical Examples

In this section, we discuss examples of the use of embedding Zassenhaus product methods in iterative splitting schemes. In the first example, we demonstrate somewhat artificially how the proposed Zassenhaus splitting method avoids the splitting error up to a certain order. The next example shows the benefit of the combined splitting schemes for applications to multiphysics problems.

In the following, we consider a benchmark example to validate the theoretical results and real life problems to apply our schemes.

We applied the following splitting schemes, given in Table 1.

##### 3.1. Benchmark Problem: Matrix Examples

In the following, the idea of the matrix example is to present the decomposition into two simpler matrices. In context of the multiphysics problems, even such simple problems can be seen as test examples. We assume that each simple matrix is presenting a basic single physics problems, for example, reaction process between two species, while the combination of the two matrices results in a more complicate physics problem, for example, coupled reaction process between two species; see [2].

We consider the matrix equation the exact solution of which is We split the matrix as where . For this simple example, we did not obtain nilpotent matrices for the commutators.

We apply 10–100 time steps for the computations. The final integration time of each single time step was less than <1 [sec], such that we need at least for such benchmark problems about 1–10 [sec]. Figure 1 presents the numerical errors, that is, the difference between the exact and the numerical solution.

Figure 2 presents the CPU times needed by the standard and the iterative splitting schemes. The application of the Zassenhaus formula did not increase, while the next Zassenhaus coefficients are at least nearly nilpotent matrices, and we can neglect the computational time. Moreover different time-steps did not influence the computations of the Zassenhaus exponents, which is a preprocess. First we compute and save all the Zassenhaus exponents, and later we multiply with the necessary time steps.

*Remark 13. *We see that the errors decrease significantly with increasing order of approximation for the initialisation process from the Zassenhaus splitting (we apply 1–3 Zassenhaus coefficients). Here we apply commutators, which did not reduce their matrix rank; therefore we have the same computation time as for the pure iterative schemes. Overall, we have their benefits in the application of the Zassenhaus product schemes to the standard Lie-Trotter or Strang-Marchuk splitting. Similar results are given with more iterative steps , as expected.

##### 3.2. Real Life Problems: Multiphase Examples

In the next examples, we simulate coupled transport and reaction equations with different phases, so-called multiphase problems; see [9]. The equations are coupled with the reaction terms and are presented as follows: where is the number of equations and is the index of each component. The unknown mobile concentrations are taken in , where is the spatial dimension. The unknown immobile concentrations are taken in , where is the spatial dimension. The retardation factors are constant, and . The kinetic part is given by the factors . They are constant and . For the initialisation of the kinetic part, we set . The kinetic part is linear and irreversible, so the successors have only one predecessor. The initial conditions are given for each component as constants or linear impulses. For the boundary conditions, we have trivial inflow and outflow conditions with at the inflow boundary. The transport part is given by the velocity and is piecewise constant; see [2]. The exchange between the mobile and immobile part is given by .

In the following, we discuss the one-phase and two-phase examples.

###### 3.2.1. One-Phase Example

The next example is a simplified real life problem with a multiphase transport-reaction equation. We treat the mobile and immobile pores in a porous media. Such simulations are made about waste scenarios; see [2, 24].

We concentrate on the computational benefits of a fast computation of the mixed iterative scheme with the Zassenhaus formula.

The one phase equation is Here the parameters take the values , , , , and .

We now turn to the finite difference schemes and semidiscretise the equation, with its convection and diffusion operators as follows: We obtain the two matrices and decouple the diffusion and convection parts as follows: We apply the splitting method to the operators and given in Section 2.1.

The submatrices are where is the number of spatial points and is the spatial step size. Consider We have the commutator where we assume that . This is a nilpotent commutator, and the computation time needed for the commutators of such operators is less than that for the full operators. Furthermore, we assume that we have a bounded spatial step size, which is given by , so that we scale the singular entry and evade a blowup in such matrices.

We have obtained the optimal results of a basically constant CPU time for decreasing time steps, which increase in the standard scheme.

We apply 10–100 timesteps for the computations with 10–100 spatial points. The final integration time of each single time step was about 1–10 , such that we need at least for such benchmark problems about 100–1000 .

Figure 3 presents the numerical error, that is, the difference between the exact and the numerical solution.

Figure 4 presents the CPU time of the standard and the iterative splitting schemes. Based on the preprocess of computing the Zassenhaus exponents at the beginning and store the matrices, we could save additional computational time with larger time steps.

*Remark 14. *For the iterative schemes with embedded Zassenhaus products, we can reach improved results more quickly. The benefits come from the constant CPU time needed for calculating the nilpotent commutators in the Zassenhaus formula. So it makes sense to have 1–3 Zassenhaus steps before starting the iterative scheme. At that point, only 2-3 iterative steps are needed to obtain more accurate results than those from the expensive standard schemes. Here, too, we obtain the best results with the one-sided iterative schemes.

*Remark 15. *Based on our restriction to , , and with , we can control the delicate blowup of the computed commutators. By the way, if we neglect the restriction and with , the commutator blows up and we lose the benefits of such an algorithm. This is a limitation of applying such an embedded Zassenhaus formula.

###### 3.2.2. Two Phase Example

The next example is a more delicate real life problem involving a multiphase transport-reaction equation. We treat mobile and immobile pores in a porous media. Such simulations are made for waste scenarios; see [2].

We concentrate on the computational benefits of a fast computation of the commutators and their use for the initialisation of the iterative splitting scheme.

The equation is Here the parameters are , , , and .

We next treat the semidiscretised equation given by the matrices where , while is the solution of the first species in the mobile phase in each spatial discretisation point , and the same for the other solution vectors.

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

Consider

We decouple this into the following matrices: For the operators and and we apply the iterative splitting method given in Section 2.1.

Based on the decomposition, is block diagonal and is tridiagonal.

The commutator is where we have a sparser operator than the main operator , and therefore we save computation time in computing such operators. Further we assume that we have a bounded spatial step size, which is given by , so that we scale the singular entry and .

Figure 5 presents the numerical error, that is, the difference between the exact and the numerical solution and the CPU time required. The optimal results are from the one-sided iterative schemes on the operator .

*Remark 16. *We obtain the same results as in the one-phase example, because of the basically constant CPU time for decreasing time steps of the iterative scheme; we obtain optimal results when applying 1–3 Zassenhaus steps. Therefore, with the embedded Zassenhaus formula, we can reach improved results more quickly than with the standard schemes. With 2-3 more iterative steps, we obtain more accurate results. With the one-sided iterative schemes, we reach the best convergence results when using the more delicate operator .

#### 4. Conclusions and Discussion

In this paper, we have presented a novel splitting scheme combining the ideas of iterative and sequential schemes based on the Zassenhaus formula. Here the idea is to decouple the expensive iterative computation of the schemes based only on the matrix exponential, obtaining simpler embedded Zassenhaus schemes based on nilpotent commutators. Such combinations have the benefit of reducing the computation time because the commutators can be computed very cheaply. On the other hand, simple linear iterative steps can be done very cheaply if they are initialised with a sufficiently accurate starting solution, and this accelerates the solvers. The error analysis showed that these are stable methods for higher order schemes. In the applications, we were able to demonstrate the speedup from the Zassenhaus enhanced method and were able to apply it to multiphysics applications. In the future, we will concentrate on linear and nonlinear matrix-dependent schemes, that switch between noniterative and iterative schemes based on Zassenhaus products.

#### References

- R. E. Ewing, “Up-scaling of biological processes and multiphase ow in porous media,” in
*IIMA Volumes in Mathematics and Its Applications*, vol. 295, pp. 195–215, Springer, New York, NY, USA, 2002. View at Google Scholar - J. Geiser,
*Discretization and Simulation of Systems for Convection-Diffusion-Dispersion Reactions with Applications in Groundwater Contamination*, Groundwater Modelling, Management and Contamination, Nova Science; Monograph, New York, NY, USA, 2008. - N. Antonic, C. J. van Duijn, W. Jäger, and A. Mikelic,
*Multiscale Problems in Science and Technology: Challenges to Mathematical Analysis and Perspectives*, Springer, Berlin, Germany, 2002. - E. Weinan,
*Principle of Multiscale Modelling*, Cambridge University Press, Cambridge, UK, 2011. - S. Blanes and F. Casas, “On the necessity of negative coefficients for operator splitting schemes of order higher than two,”
*Applied Numerical Mathematics*, vol. 54, no. 1, pp. 23–37, 2005. View at Publisher · View at Google Scholar · View at MathSciNet - E. Hansen and A. Ostermann, “High order splitting methods for analytic semigroups exist,”
*BIT. Numerical Mathematics*, vol. 49, no. 3, pp. 527–542, 2009. View at Publisher · View at Google Scholar · View at MathSciNet - G. Strang, “On the construction and comparison of difference schemes,”
*SIAM Journal on Numerical Analysis*, vol. 5, pp. 506–517, 1968. View at Publisher · View at Google Scholar · View at MathSciNet - G. I. Marchuk, “Some applicatons of splitting-up methods to the solution of problems in mathematical physics,”
*Aplikace Matematiky*, vol. 1, pp. 103–132, 1968. View at Google Scholar - J. Geiser, “Computing exponential for iterative splitting methods: algorithms and applications,”
*Journal of Applied Mathematics*, vol. 2011, Article ID 193781, 27 pages, 2011. View at Publisher · View at Google Scholar - R. I. McLachlan and G. R. W. Quispel, “Splitting methods,”
*Acta Numerica*, vol. 11, pp. 341–434, 2002. View at Publisher · View at Google Scholar · View at MathSciNet - J. Geiser, G. Tanoğlu, and N. Gücüyenen, “Higher order operator splitting methods via Zassenhaus product formula: theory and applications,”
*Computers & Mathematics with Applications*, vol. 62, no. 4, pp. 1994–2015, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - G. H. Weiss and A. A. Maradudin, “The Baker-Hausdorff formula and a problem in crystal physics,”
*Journal of Mathematical Physics*, vol. 3, pp. 771–777, 1962. View at Publisher · View at Google Scholar · View at MathSciNet - R. M. Wilcox, “Exponential operators and parameter differentiation in quantum physics,”
*Journal of Mathematical Physics*, vol. 8, pp. 962–982, 1967. View at Publisher · View at Google Scholar · View at MathSciNet - J. Geiser, “An iterative splitting approach for linear integro-differential equations,”
*Applied Mathematics Letters*, vol. 26, no. 11, pp. 1048–1052, 2013. View at Google Scholar - F. Casas and A. Murua, “An efficient algorithm for computing the Baker-Campbell-Hausdorff series and some of its applications,”
*Journal of Mathematical Physics*, vol. 50, no. 3, Article ID 033513, 2009. View at Publisher · View at Google Scholar · View at MathSciNet - J. Geiser,
*Iterative Splitting Methods for Differential Equations*, Chapman & Hall/CRC Numerical Analysis and Scientific Computing, CRC Press, Boca Raton, Fla, USA, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - S. Vandewalle,
*Parallel Multigrid Waveform Relaxation for Parabolic Problems*, Teubner, Stuttgart, Germany, 1993. - J. Geiser and G. Tanoglu, “Operator-splitting methods via the Zassenhaus product formula,”
*Applied Mathematics and Computation*, vol. 217, no. 9, pp. 4557–4575, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - F. Bayen, “On the convergence of the Zassenhaus formula,”
*Letters in Mathematical Physics*, vol. 3, no. 3, pp. 161–167, 1979. View at Publisher · View at Google Scholar · View at MathSciNet - F. Casas, A. Murua, and M. Nadinic, “Efficient computation of the Zassenhaus formula,”
*Computer Physics Communications*, vol. 183, no. 11, pp. 2386–2391, 2012. View at Publisher · View at Google Scholar · View at MathSciNet - E. Hairer, C. Lubich, and G. Wanner,
*Geometric Numerical Integration*, vol. 31 of*Springer Series in Computational Mathematics*, Springer, Berlin, Germany, 2002. View at MathSciNet - H. Yoshida, “Construction of higher order symplectic integrators,”
*Physics Letters A*, vol. 150, no. 5–7, pp. 262–268, 1990. View at Publisher · View at Google Scholar · View at MathSciNet - E. Faou, A. Ostermann, and K. Schratz, “Analysis of exponential splitting methods for inhomogeneous parabolic equations,” http://arxiv.org/abs/1212.5827.
- J. Geiser, “Discretization methods with embedded analytical solutions for convection-diffusion dispersion-reaction equations and applications,”
*Journal of Engineering Mathematics*, vol. 57, no. 1, pp. 79–98, 2007. View at Publisher · View at Google Scholar · View at MathSciNet