Embedded Zassenhaus Expansion to Splitting Schemes: Theory and Multiphysics Applications

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.


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

International Journal of Differential Equations
We focus our attention on the case of two linear operators (i.e., we consider the Cauchy problem) as follows: whereby the initial function 0 is given, and and are assumed to be bounded linear operators in the Banach space X with , : X → X. 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 (0) is the initial condition.
Here one of the main ideas to compute such delicate solutions is to separate or decouple the term exp( + ) into simpler and parallel computational expression, for example, exp( ) exp( ); 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.

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.

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.
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 init ( ) 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 1 ( ) of the iterative splitting schemes can be improved by a cheap computable initial process 1 ( ) ≈ exp(( + ) ). We will next discuss the fast convergent Zassenhaus formula for the initial process.

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 L( , ) is the free Lie algebra generated by and . The kernel function of solution (1), exp(( + ) ), can be uniquely decomposed as where ( , ) ∈ L( , ) 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 exp (− ) exp (( + ) ) = exp ( + ( )) , where ( ) is a Lie polynomial of degree >1, and Further we have exp (− ) exp ( + ( )) = exp ( 2 2 +̃( )) , wherẽ( ) 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 exp( ) of the scheme is much less than that of the main operators exp( ) exp( ), 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 ( , ) ∈ L( , ) is a homogeneous Lie polynomial in and of degree ; see [21]. The BCH formula is given as 4 International Journal of Differential Equations 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 2 of the Zassenhaus formula where = + , and D is the domain of the operators.
(ii) Condition for the second term 3 of the Zassenhaus formula where = + , and D 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 Ω = (0, 1) 2 . Then we need the following regularity result: (i) condition for the first term 2 of the Zassenhaus formula (ii) condition for the second term 3 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].

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 O( ), 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 O( + ).
Combining both splitting schemes, we have a local error of International Journal of Differential Equations 5 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 Δ = +1 − with = 0, . . . , − 1. Further ( 0 ) is the initial condition of the Cauchy problem (1), and we compute ( 1 ), . . . , ( ) successively. We have the following two steps; first we apply the Zassenhaus formula as initialization of 0 ( +1 ) for the iterative scheme; second we apply the iterative splitting scheme to the initial value problem (1) and compute the final solution ( +1 ) as follows.
(1) We start with = 2 and compute the Zassenhaus formula.
(2) We compute the matrix norm for the Zassenhauss exponents .
(3) If the error of the iterative scheme err iter ≤ err 2 , then we are finished; else = +1, and we compute the next iterative step, till err iter ≤ err 2 .

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.

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 Iterative splitting with one iterative step and one Zassenhaus step as prestep 3 Iterative splitting with one iterative step and two Zassenhaus steps as presteps 4 Iterative splitting with two iterative steps and two Zassenhaus steps as presteps 5 Iterative splitting with three iterative steps and two Zassenhaus steps as presteps 6 Iterative splitting with four iterative steps and two Zassenhaus steps as presteps 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 [ , ] ̸ = 0. 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 timesteps 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 = 3, . . . , 6, as expected. see [9]. The equations are coupled with the reaction terms and are presented as follows:  For the initialisation of the kinetic part, we set 0 = 0. 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 = 0 at the inflow boundary. The transport part is given by the velocity v ∈ R and is piecewise constant; see [2]. The exchange between the mobile and immobile part is given by .

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

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.
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 1 and 2 given in Section 2.1.
The submatrices are diff, = Δ 2 1 = Δ 2 ⋅ ( where is the number of spatial points and Δ is the spatial step size. Consider We have the commutator where we assume that [ 1 , 2 ] ≈ 0. 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 Δ = 0.1 ≫ 0, so that we scale the singular entry ( 2 − 1 )/Δ 2 ≈ 1 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 [sec], such that we need at least for such benchmark problems about 100-1000 [sec]. 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. 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 1 , 2 , and Δ with ( 2 − 1 )/Δ 2 ≈ 1, we can control the delicate blowup of the computed commutators. By the way, if we neglect the restriction and Δ with ( 2 − 1 )/Δ 2 → ∞, the commutator blows up and we lose the benefits of such an algorithm. This is a limitation of applying such an embedded Zassenhaus formula.

Two Phase Example.
The next example is a more delicate real life problem involving a multiphase transportreaction 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.
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̃= 1 + 2 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 Δ = 0.1 ≫ 0, so that we scale the singular entry /Δ 2 ≈ 1 and V/Δ ≈ 1. 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 onesided iterative schemes on the operator̃.

Remark 16.
We obtain the same results as in the onephase 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̃.

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 matrixdependent schemes, that switch between noniterative and iterative schemes based on Zassenhaus products.