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

Authors: Jürgen Geiser
Journal: International Journal of Differential Equations
 (2013)

---

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


---

## Body

## 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:(1)∂c(t)∂t=Ac(t)+Bc(t),witht∈[0,T],c(0)=c0,
whereby the initial function c0 is given, and A and B are assumed to be bounded linear operators in the Banach space X with A,B: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
(2)c(t)=exp(A+B)c(0),witht∈[0,T],
where c(0) is the initial condition.Here one of the main ideas to compute such delicate solutions is to separate or decouple the termexp(A+B) into simpler and parallel computational expression, for example, exp(A)exp(B); 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 Section2. 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:(3)∂ci(t)∂t=Aci(t)+Bci-1(t),withci(tn)=cn,i=1,2,…,m.
Iterative splitting with respect to two operators (also called a two-sided scheme) as follows:
(4)∂ci(t)∂t=Aci(t)+Bci-1(t),withci(tn)=cn,(5)∂ci+1(t)∂t=Aci(t)+Bci+1(t),withci+1(tn)=cn,i=1,3,…,2m+1.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:
(6)∥Si-exp((A+B)τ)∥≤Cτi,
where Si is the approximate solution for the ith iterative step and C is a constant that can be chosen uniformly on bounded time intervals.Proof.
The proof is given fori=1,2,… and with the consistency error ei(τ)=c(τ)-ci(τ), we have
(7)ci(τ)=exp(Aτ)c(tn)+∫tntn+1exp(A(tn+1-s))Bexp(sA)c(tn)ds+∫tntn+1exp(A(tn+1-s1))B×∫tntn+1-s1exp((tn+1-s1-s2)A)×Bexp(s2A)c(tn)ds2ds1+⋯+∫tntn+1exp(A(tn+1-s1))B⋯∫tntn+1-∑j=1i-2sjexp((tn+1-∑j=1i-1sj)A)×Bcinit(si)dsidsi-1,…,ds2ds1,c(τ)=exp(Aτ)c(tn)+∫tntn+1exp(A(tn+1-s))Bexp(sA)c(tn)ds+∫tntn+1exp(A(tn+1-s1))B×∫tntn+1-s1exp((tn+1-s1-s2)A)×Bexp(s2A)c(tn)ds2ds1+⋯+∫tntn+1exp(A(tn+1-s1))B⋯∫tntn+1-∑j=1i-1sjexp((tn+1-∑j=1isj)A)×Bexp((si(A+B))c(tn)dsi,…,ds2ds1.
We obtain
(8)∥ei∥≤∥exp((A+B)τ)c(tn)-Si(τ)cinit(τ)∥≤Cτimaxsi∈[0,τ]∥exp(si(A+B))c(tn)-cinit(si)∥≤Cτi∥exp(τ(A+B))c(tn)-cinit(τ)∥,
where i is the number of iterative steps.
The same idea can be applied to the even iterative scheme and also for alternatingA and B.Remark 2.
The accuracy of the initialisationcinit(τ) is important for conserving or improving the underlying iterative splitting scheme.
Here we have the following initialisation schemes:(9)c1(τ)=exp(Aτ)c(tn)⟶∥ei∥≤Cτic(tn),c1(τ)=exp(Aτ)exp(Bτ)c(tn)⟶∥ei∥≤Cτi+1c(tn).
So the slow convergence in the initial process c1(τ) of the iterative splitting schemes can be improved by a cheap computable initial process c1(τ)≈exp((A+B)τ). 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 Theorem3; see also [20].Theorem 3.
One assumesℒ(A,B) is the free Lie algebra generated by A and B. The kernel function of solution (1), exp((A+B)t), can be uniquely decomposed as
(10)exp((A+B)t)=exp(At)exp(Bt)πi=1jπk=2m×exp(Ck(A,B)tk)+𝒪(tm+1),
where Cj(A,B)∈ℒ(A,B) is a homogeneous Lie polynomial of A and B of degree j.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
(11)exp(-At)exp((A+B)t)=exp(Bt+D(t)),
where D(t) is a Lie polynomial of degree >1, and
(12)exp((A+B)t)=exp(At)exp(Bt)+𝒪(t2).
Further we have
(13)exp(-Bt)exp(Bt+D(t))=exp(C2t2+D~(t)),
where D~(t) is a Lie polynomial of degree >2, and
(14)exp((A+B)t)=exp(At)exp(Bt)exp(C2t2)+𝒪(t3).
We can repeat the recursive process and obtain by complete induction the results.The first Zassenhaus polynomials are given as(15)C2(A,B)=12[B,A],C3(A,B)=13[[B,A],B]+16[[B,A],A],C4(A,B)=124[[[B,A],A],A]+18[[[B,A],B],A]+18[[[B,A],B],B].
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(Cktk) of the scheme is much less than that of the main operators exp(At)exp(Bt), 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)
A-B 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
(16)exp(At2)exp(Bt)exp(At2)=exp(tS1+t3S3+t5S5+⋯),
where the Si(A,B)∈ℒ(A,B) is a homogeneous Lie polynomial in A and B of degree i; see [21].
The BCH formula is given as(17)exp((A2+B2)t)=Πk=2∞exp(C~ktk)exp(A2t)exp(B2t),exp((B2+A2)t)=exp(B2t)exp(A2t)Πk=2∞exp(Cktk),
and so there exists a new product based on the Zassenhaus polynomials as follows:
(18)Πk=3∞exp(Dktk)=Πk=2∞exp(C~ktk)Πk=2∞exp(Cktk),
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 termC2 of the Zassenhaus formula
(19)𝒟(L2)⊂𝒟(A2)∩𝒟(AB)∩𝒟(BA)∩𝒟(B2),whereL=A+B, and 𝒟 is the domain of the operators.(ii)
Condition for the second termC3 of the Zassenhaus formula
(20)𝒟(L3)⊂𝒟(A3)∩𝒟(A2B)∩𝒟(ABA)∩𝒟(AB2)∩𝒟(B3),whereL=A+B, 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
(21)𝒟(L)=H2(Ω)∩H01(Ω),𝒟(A)={c∈L2(Ω);∂xxc,∂xc∈L2Ω,u(0,y)=u(1,y)=0,∀y∈(0,1)c∈L2(Ω);∂xxc,∂xc∈L2Ω,},𝒟(B)={c∈L2(Ω);∂yyc,∂yc∈L2Ω,u(x,0)=u(x,1)=0,∀x∈(0,1)c∈L2(Ω);∂yyc,∂yc∈L2Ω,},
where Ω=(0,1)2.
Then we need the following regularity result:(i)
condition for the first termC2 of the Zassenhaus formula
(22)𝒟(L2)⊂H4(Ω),(ii)
condition for the second termC3 of the Zassenhaus formula
(23)𝒟(L3)⊂H6(Ω).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(24)EZassen,Comp,1(t)=exp(At)exp(Bt),EZassen,Comp,j(t)=exp(Cjtj),forj∈2,…,i,
where the sequential Zassenhaus operators are defined as
(25)EZassen,i(t)=Πj=1iEZassen,Comp,j(t),
and we obtain an accuracy of 𝒪(ti), and i 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 A, B.
The initialisation process is done with the Zassenhaus formula(26)ci(t)=EZassen,i(t)c0,
where EZassen,i(t) is given in (25).
Furthermore, the improved solutions are embedded in the iterative splitting scheme (4) and one has after j iterative steps the following result:
(27)ci+j(t)=Eiter,j(t)EZassen,i(t)c0.
This has improved the error of the iterative scheme to 𝒪(ti+j).Proof.
The solution of the iterative splitting scheme (4) is
(28)ci+j(t)=Eiter,i(t)cinit,j,
where Si(t)=Eiter,i(t), given in Theorem 1.
The initialisation is given by the Zassenhaus formula(29)cinit,j(t)=EZassen,j(t)c0.
Combining both splitting schemes, we have a local error of
(30)ei+j(t)=∥c(t)-ci+j(t)∥=∥c(t)-Eiter,i(t)EZassen,j(t)c0∥≤𝒪(ti+j)c0.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 computeN time steps, where each time step Δtn is given as Δtn=tn+1-tn with n=0,…,N-1. Further c(t0) is the initial condition of the Cauchy problem (1), and we compute c(t1),…,c(tN) successively. We have the following two steps; first we apply the Zassenhaus formula as initialization of c0(tn+1) for the iterative scheme; second we apply the iterative splitting scheme to the initial value problem (1) and compute the final solution c(tn+1) as follows.(1)
We start withj=2 and compute the Zassenhaus formula.(2)
We compute the matrix norm for the Zassenhauss exponentsCj.IferrZassj≤err1, then we start with the iterative splitting scheme with c0(t)=cZassj, i=1 and go to (3), else j=j+1 and we go to (1).(3)
If the error of the iterative schemeerriter≤err2, then we are finished; else i=i+1, and we compute the next iterative step, till erriter≤err2.


## 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:(3)∂ci(t)∂t=Aci(t)+Bci-1(t),withci(tn)=cn,i=1,2,…,m.
Iterative splitting with respect to two operators (also called a two-sided scheme) as follows:
(4)∂ci(t)∂t=Aci(t)+Bci-1(t),withci(tn)=cn,(5)∂ci+1(t)∂t=Aci(t)+Bci+1(t),withci+1(tn)=cn,i=1,3,…,2m+1.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:
(6)∥Si-exp((A+B)τ)∥≤Cτi,
where Si is the approximate solution for the ith iterative step and C is a constant that can be chosen uniformly on bounded time intervals.Proof.
The proof is given fori=1,2,… and with the consistency error ei(τ)=c(τ)-ci(τ), we have
(7)ci(τ)=exp(Aτ)c(tn)+∫tntn+1exp(A(tn+1-s))Bexp(sA)c(tn)ds+∫tntn+1exp(A(tn+1-s1))B×∫tntn+1-s1exp((tn+1-s1-s2)A)×Bexp(s2A)c(tn)ds2ds1+⋯+∫tntn+1exp(A(tn+1-s1))B⋯∫tntn+1-∑j=1i-2sjexp((tn+1-∑j=1i-1sj)A)×Bcinit(si)dsidsi-1,…,ds2ds1,c(τ)=exp(Aτ)c(tn)+∫tntn+1exp(A(tn+1-s))Bexp(sA)c(tn)ds+∫tntn+1exp(A(tn+1-s1))B×∫tntn+1-s1exp((tn+1-s1-s2)A)×Bexp(s2A)c(tn)ds2ds1+⋯+∫tntn+1exp(A(tn+1-s1))B⋯∫tntn+1-∑j=1i-1sjexp((tn+1-∑j=1isj)A)×Bexp((si(A+B))c(tn)dsi,…,ds2ds1.
We obtain
(8)∥ei∥≤∥exp((A+B)τ)c(tn)-Si(τ)cinit(τ)∥≤Cτimaxsi∈[0,τ]∥exp(si(A+B))c(tn)-cinit(si)∥≤Cτi∥exp(τ(A+B))c(tn)-cinit(τ)∥,
where i is the number of iterative steps.
The same idea can be applied to the even iterative scheme and also for alternatingA and B.Remark 2.
The accuracy of the initialisationcinit(τ) is important for conserving or improving the underlying iterative splitting scheme.
Here we have the following initialisation schemes:(9)c1(τ)=exp(Aτ)c(tn)⟶∥ei∥≤Cτic(tn),c1(τ)=exp(Aτ)exp(Bτ)c(tn)⟶∥ei∥≤Cτi+1c(tn).
So the slow convergence in the initial process c1(τ) of the iterative splitting schemes can be improved by a cheap computable initial process c1(τ)≈exp((A+B)τ). 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 Theorem3; see also [20].Theorem 3.
One assumesℒ(A,B) is the free Lie algebra generated by A and B. The kernel function of solution (1), exp((A+B)t), can be uniquely decomposed as
(10)exp((A+B)t)=exp(At)exp(Bt)πi=1jπk=2m×exp(Ck(A,B)tk)+𝒪(tm+1),
where Cj(A,B)∈ℒ(A,B) is a homogeneous Lie polynomial of A and B of degree j.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
(11)exp(-At)exp((A+B)t)=exp(Bt+D(t)),
where D(t) is a Lie polynomial of degree >1, and
(12)exp((A+B)t)=exp(At)exp(Bt)+𝒪(t2).
Further we have
(13)exp(-Bt)exp(Bt+D(t))=exp(C2t2+D~(t)),
where D~(t) is a Lie polynomial of degree >2, and
(14)exp((A+B)t)=exp(At)exp(Bt)exp(C2t2)+𝒪(t3).
We can repeat the recursive process and obtain by complete induction the results.The first Zassenhaus polynomials are given as(15)C2(A,B)=12[B,A],C3(A,B)=13[[B,A],B]+16[[B,A],A],C4(A,B)=124[[[B,A],A],A]+18[[[B,A],B],A]+18[[[B,A],B],B].
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(Cktk) of the scheme is much less than that of the main operators exp(At)exp(Bt), 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)
A-B 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
(16)exp(At2)exp(Bt)exp(At2)=exp(tS1+t3S3+t5S5+⋯),
where the Si(A,B)∈ℒ(A,B) is a homogeneous Lie polynomial in A and B of degree i; see [21].
The BCH formula is given as(17)exp((A2+B2)t)=Πk=2∞exp(C~ktk)exp(A2t)exp(B2t),exp((B2+A2)t)=exp(B2t)exp(A2t)Πk=2∞exp(Cktk),
and so there exists a new product based on the Zassenhaus polynomials as follows:
(18)Πk=3∞exp(Dktk)=Πk=2∞exp(C~ktk)Πk=2∞exp(Cktk),
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 termC2 of the Zassenhaus formula
(19)𝒟(L2)⊂𝒟(A2)∩𝒟(AB)∩𝒟(BA)∩𝒟(B2),whereL=A+B, and 𝒟 is the domain of the operators.(ii)
Condition for the second termC3 of the Zassenhaus formula
(20)𝒟(L3)⊂𝒟(A3)∩𝒟(A2B)∩𝒟(ABA)∩𝒟(AB2)∩𝒟(B3),whereL=A+B, 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
(21)𝒟(L)=H2(Ω)∩H01(Ω),𝒟(A)={c∈L2(Ω);∂xxc,∂xc∈L2Ω,u(0,y)=u(1,y)=0,∀y∈(0,1)c∈L2(Ω);∂xxc,∂xc∈L2Ω,},𝒟(B)={c∈L2(Ω);∂yyc,∂yc∈L2Ω,u(x,0)=u(x,1)=0,∀x∈(0,1)c∈L2(Ω);∂yyc,∂yc∈L2Ω,},
where Ω=(0,1)2.
Then we need the following regularity result:(i)
condition for the first termC2 of the Zassenhaus formula
(22)𝒟(L2)⊂H4(Ω),(ii)
condition for the second termC3 of the Zassenhaus formula
(23)𝒟(L3)⊂H6(Ω).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(24)EZassen,Comp,1(t)=exp(At)exp(Bt),EZassen,Comp,j(t)=exp(Cjtj),forj∈2,…,i,
where the sequential Zassenhaus operators are defined as
(25)EZassen,i(t)=Πj=1iEZassen,Comp,j(t),
and we obtain an accuracy of 𝒪(ti), and i 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 A, B.
The initialisation process is done with the Zassenhaus formula(26)ci(t)=EZassen,i(t)c0,
where EZassen,i(t) is given in (25).
Furthermore, the improved solutions are embedded in the iterative splitting scheme (4) and one has after j iterative steps the following result:
(27)ci+j(t)=Eiter,j(t)EZassen,i(t)c0.
This has improved the error of the iterative scheme to 𝒪(ti+j).Proof.
The solution of the iterative splitting scheme (4) is
(28)ci+j(t)=Eiter,i(t)cinit,j,
where Si(t)=Eiter,i(t), given in Theorem 1.
The initialisation is given by the Zassenhaus formula(29)cinit,j(t)=EZassen,j(t)c0.
Combining both splitting schemes, we have a local error of
(30)ei+j(t)=∥c(t)-ci+j(t)∥=∥c(t)-Eiter,i(t)EZassen,j(t)c0∥≤𝒪(ti+j)c0.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 computeN time steps, where each time step Δtn is given as Δtn=tn+1-tn with n=0,…,N-1. Further c(t0) is the initial condition of the Cauchy problem (1), and we compute c(t1),…,c(tN) successively. We have the following two steps; first we apply the Zassenhaus formula as initialization of c0(tn+1) for the iterative scheme; second we apply the iterative splitting scheme to the initial value problem (1) and compute the final solution c(tn+1) as follows.(1)
We start withj=2 and compute the Zassenhaus formula.(2)
We compute the matrix norm for the Zassenhauss exponentsCj.IferrZassj≤err1, then we start with the iterative splitting scheme with c0(t)=cZassj, i=1 and go to (3), else j=j+1 and we go to (1).(3)
If the error of the iterative schemeerriter≤err2, then we are finished; else i=i+1, and we compute the next iterative step, till erriter≤err2.


## 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 Table1.Table 1
Applied splitting schemes.
Abbreviation in the examples
Splitting methods
A-B
A-B splitting
Strang
Strang splitting
c
1
Iterative splitting with one iterative step
c
2
Iterative splitting with one iterative step and one Zassenhaus step as prestep
c
3
Iterative splitting with one iterative step and two Zassenhaus steps as presteps
c
4
Iterative splitting with two iterative steps and two Zassenhaus steps as presteps
c
5
Iterative splitting with three iterative steps and two Zassenhaus steps as presteps
c
6
Iterative splitting with four iterative steps and two Zassenhaus steps as presteps


### 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(31)u′(t)=[1230]u,u(0)=u0=(01),t∈[0,1],
the exact solution of which is
(32)u(t)=2(e3t-e-2t)5.
We split the matrix as
(33)[1230]=[1110]+[0120]=A+B,
where [A,B]≠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]. Figure1 presents the numerical errors, that is, the difference between the exact and the numerical solution.Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)Figure2 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.CPU times of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)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 stepsi=3,…,6, 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:
(34)∂tRiui+∇·vui=-λiRiui+λi-1Ri-1ui-1+β(-ui+gi),inΩ×(0,T),ui,0(x)=ui(x,0),onΩ,∂tRigi=-λiRigi+λi-1Ri-1gi-1+β(-gi+ui),inΩ×(0,T),gi,0(x)=gi(x,0),onΩ,i=1,…,m,
where m is the number of equations and i is the index of each component. The unknown mobile concentrations ui=ui(x,t) are taken in Ω×(0,T)⊂ℝn×ℝ+, where n is the spatial dimension. The unknown immobile concentrations gi=gi(x,t) are taken in Ω×(0,T)⊂ℝn×ℝ+, where n is the spatial dimension. The retardation factors Ri are constant, and Ri≥0. The kinetic part is given by the factors λi. They are constant and λi≥0. 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 i as constants or linear impulses. For the boundary conditions, we have trivial inflow and outflow conditions with ui=0 at the inflow boundary. The transport part is given by the velocity v∈ℝn 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(35)∂tc1+∇·F1c1=-λ1c1,inΩ×[0,t],∂tc2+∇·F2c2=λ1c1-λ2c2,inΩ×[0,t],Fi=vi-Di∇,i=1,2,c1(x,t)=c1,0(x),c2(x,t)=c2,0(x),onΩ,c1(x,t)=c1,1(x,t),c2(x,t)=c2,1(x,t),on∂Ω×[0,t].
Here the parameters take the values v1=0.1, v2=0.05, D1=0.01, D2=0.005, and λ1=λ2=0.1.We now turn to the finite difference schemes and semidiscretise the equation, with its convection and diffusion operators as follows:(36)∂tc=(A1+A2)c.
We obtain the two matrices and decouple the diffusion and convection parts as follows:
(37)A1=(Adiff00Adiff)∈ℝ2I×2I,A2=(AConv00AConv)+(-Λ10Λ1-Λ2)∈ℝ2I×2I.
We apply the splitting method to the operators A1 and A2 given in Section 2.1.The submatrices are(38)Adiff,i=DiΔx2R1=DiΔx2·(-211-21⋱⋱⋱1-211-2)∈ℝI×I,Aconv,i=-viΔxR2=-viΔx·(1-11⋱⋱-11-11)∈ℝI×I,
where I is the number of spatial points and Δx is the spatial step size. Consider
(39)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈ℝI×I,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈ℝI×I.
We have the commutator
(40)[A1,A2]=(00D2-D1Δx2R1Λ10)∈ℝ2I×2I,
where we assume that [R1,R2]≈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 Δx=0.1≫0, so that we scale the singular entry (D2-D1)/Δx2≈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].Figure3 presents the numerical error, that is, the difference between the exact and the numerical solution.Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)
(c)Figure4 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.CPU time of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)
(c)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 toD1, D2, and Δx with (D2-D1)/Δx2≈1, we can control the delicate blowup of the computed commutators. By the way, if we neglect the restriction and Δx with (D2-D1)/Δx2→∞, 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(41)∂tc1+∇·Fc1=g(-c1+c1,im)-λ1c1,inΩ×[0,t],∂tc2+∇·Fc2=g(-c2+c2,im)+λ1c1-λ2c2,inΩ×[0,t],F=v-D∇,∂tc1,im=g(c1-c1,im)-λ1c1,im,inΩ×[0,t],∂tc2,im=g(c2-c2,im)+λ1c1,im-λ2c2,im,inΩ×[0,t],c1(x,t)=c1,0(x),c2(x,t)=c2,0(x),onΩ,c1(x,t)=c1,1(x,t),c2(x,t)=c2,1(x,t),on∂Ω×[0,t],c1,im(x,t)=0,c2,im(x,t)=0,onΩ,c1,im(x,t)=0,c2,im(x,t)=0,on∂Ω×[0,t].
Here the parameters are v=0.1, D=0.01, λ1=λ2=0.1, and g=0.01.We next treat the semidiscretised equation given by the matrices(42)∂tC=(A-Λ1-G0G0Λ1A-Λ2-G0GG0-Λ1-G00GΛ1-Λ2-G)C,
where C=(c1,c2,c1im,c2im)T, while c1=(c1,1,…,c1,I) is the solution of the first species in the mobile phase in each spatial discretisation point (i=1,…,I), and the same for the other solution vectors.We have the following two operators for the splitting method:(43)A=DΔx2·(-211-21⋱⋱⋱1-211-2),+vΔx·(1-11⋱⋱-11-11)∈ℝI×I,
where I is the number of spatial points.Consider(44)Λ1=(λ100λ10⋱⋱⋱0λ100λ1),Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈ℝI×I,G=(g00g0⋱⋱⋱0g00g)∈ℝI×I.We decouple this into the following matrices:(45)A~=(A0000A0000000000)∈ℝ4I×4I,B1=(-Λ1000Λ1-Λ20000-Λ1000Λ1-Λ2),B2=(-G0G00-G0GG0-G00G0-G)∈ℝ4I×4I.
For the operators A~ and B~=B1+B2 and B we apply the iterative splitting method given in Section 2.1.Based on the decomposition,A~ is block diagonal and B~ is tridiagonal.The commutator is(46)[A~,B~]=(00AG0000AG-AG0000-AG00)∈ℝ4I×4I,
where we have a sparser operator than the main operator B~, 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 Δx=0.1≫0, so that we scale the singular entry D/Δx2≈1 and v/Δx≈1.Figure5 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 B~.Numerical errors and CPU time of the standard splitting scheme and the iterative schemes with Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)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 operatorB~.


## 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(31)u′(t)=[1230]u,u(0)=u0=(01),t∈[0,1],
the exact solution of which is
(32)u(t)=2(e3t-e-2t)5.
We split the matrix as
(33)[1230]=[1110]+[0120]=A+B,
where [A,B]≠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]. Figure1 presents the numerical errors, that is, the difference between the exact and the numerical solution.Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)Figure2 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.CPU times of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)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 stepsi=3,…,6, 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:
(34)∂tRiui+∇·vui=-λiRiui+λi-1Ri-1ui-1+β(-ui+gi),inΩ×(0,T),ui,0(x)=ui(x,0),onΩ,∂tRigi=-λiRigi+λi-1Ri-1gi-1+β(-gi+ui),inΩ×(0,T),gi,0(x)=gi(x,0),onΩ,i=1,…,m,
where m is the number of equations and i is the index of each component. The unknown mobile concentrations ui=ui(x,t) are taken in Ω×(0,T)⊂ℝn×ℝ+, where n is the spatial dimension. The unknown immobile concentrations gi=gi(x,t) are taken in Ω×(0,T)⊂ℝn×ℝ+, where n is the spatial dimension. The retardation factors Ri are constant, and Ri≥0. The kinetic part is given by the factors λi. They are constant and λi≥0. 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 i as constants or linear impulses. For the boundary conditions, we have trivial inflow and outflow conditions with ui=0 at the inflow boundary. The transport part is given by the velocity v∈ℝn 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(35)∂tc1+∇·F1c1=-λ1c1,inΩ×[0,t],∂tc2+∇·F2c2=λ1c1-λ2c2,inΩ×[0,t],Fi=vi-Di∇,i=1,2,c1(x,t)=c1,0(x),c2(x,t)=c2,0(x),onΩ,c1(x,t)=c1,1(x,t),c2(x,t)=c2,1(x,t),on∂Ω×[0,t].
Here the parameters take the values v1=0.1, v2=0.05, D1=0.01, D2=0.005, and λ1=λ2=0.1.We now turn to the finite difference schemes and semidiscretise the equation, with its convection and diffusion operators as follows:(36)∂tc=(A1+A2)c.
We obtain the two matrices and decouple the diffusion and convection parts as follows:
(37)A1=(Adiff00Adiff)∈ℝ2I×2I,A2=(AConv00AConv)+(-Λ10Λ1-Λ2)∈ℝ2I×2I.
We apply the splitting method to the operators A1 and A2 given in Section 2.1.The submatrices are(38)Adiff,i=DiΔx2R1=DiΔx2·(-211-21⋱⋱⋱1-211-2)∈ℝI×I,Aconv,i=-viΔxR2=-viΔx·(1-11⋱⋱-11-11)∈ℝI×I,
where I is the number of spatial points and Δx is the spatial step size. Consider
(39)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈ℝI×I,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈ℝI×I.
We have the commutator
(40)[A1,A2]=(00D2-D1Δx2R1Λ10)∈ℝ2I×2I,
where we assume that [R1,R2]≈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 Δx=0.1≫0, so that we scale the singular entry (D2-D1)/Δx2≈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].Figure3 presents the numerical error, that is, the difference between the exact and the numerical solution.Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)
(c)Figure4 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.CPU time of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)
(c)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 toD1, D2, and Δx with (D2-D1)/Δx2≈1, we can control the delicate blowup of the computed commutators. By the way, if we neglect the restriction and Δx with (D2-D1)/Δx2→∞, 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(41)∂tc1+∇·Fc1=g(-c1+c1,im)-λ1c1,inΩ×[0,t],∂tc2+∇·Fc2=g(-c2+c2,im)+λ1c1-λ2c2,inΩ×[0,t],F=v-D∇,∂tc1,im=g(c1-c1,im)-λ1c1,im,inΩ×[0,t],∂tc2,im=g(c2-c2,im)+λ1c1,im-λ2c2,im,inΩ×[0,t],c1(x,t)=c1,0(x),c2(x,t)=c2,0(x),onΩ,c1(x,t)=c1,1(x,t),c2(x,t)=c2,1(x,t),on∂Ω×[0,t],c1,im(x,t)=0,c2,im(x,t)=0,onΩ,c1,im(x,t)=0,c2,im(x,t)=0,on∂Ω×[0,t].
Here the parameters are v=0.1, D=0.01, λ1=λ2=0.1, and g=0.01.We next treat the semidiscretised equation given by the matrices(42)∂tC=(A-Λ1-G0G0Λ1A-Λ2-G0GG0-Λ1-G00GΛ1-Λ2-G)C,
where C=(c1,c2,c1im,c2im)T, while c1=(c1,1,…,c1,I) is the solution of the first species in the mobile phase in each spatial discretisation point (i=1,…,I), and the same for the other solution vectors.We have the following two operators for the splitting method:(43)A=DΔx2·(-211-21⋱⋱⋱1-211-2),+vΔx·(1-11⋱⋱-11-11)∈ℝI×I,
where I is the number of spatial points.Consider(44)Λ1=(λ100λ10⋱⋱⋱0λ100λ1),Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈ℝI×I,G=(g00g0⋱⋱⋱0g00g)∈ℝI×I.We decouple this into the following matrices:(45)A~=(A0000A0000000000)∈ℝ4I×4I,B1=(-Λ1000Λ1-Λ20000-Λ1000Λ1-Λ2),B2=(-G0G00-G0GG0-G00G0-G)∈ℝ4I×4I.
For the operators A~ and B~=B1+B2 and B we apply the iterative splitting method given in Section 2.1.Based on the decomposition,A~ is block diagonal and B~ is tridiagonal.The commutator is(46)[A~,B~]=(00AG0000AG-AG0000-AG00)∈ℝ4I×4I,
where we have a sparser operator than the main operator B~, 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 Δx=0.1≫0, so that we scale the singular entry D/Δx2≈1 and v/Δx≈1.Figure5 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 B~.Numerical errors and CPU time of the standard splitting scheme and the iterative schemes with Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)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 operatorB~.


## 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(35)∂tc1+∇·F1c1=-λ1c1,inΩ×[0,t],∂tc2+∇·F2c2=λ1c1-λ2c2,inΩ×[0,t],Fi=vi-Di∇,i=1,2,c1(x,t)=c1,0(x),c2(x,t)=c2,0(x),onΩ,c1(x,t)=c1,1(x,t),c2(x,t)=c2,1(x,t),on∂Ω×[0,t].
Here the parameters take the values v1=0.1, v2=0.05, D1=0.01, D2=0.005, and λ1=λ2=0.1.We now turn to the finite difference schemes and semidiscretise the equation, with its convection and diffusion operators as follows:(36)∂tc=(A1+A2)c.
We obtain the two matrices and decouple the diffusion and convection parts as follows:
(37)A1=(Adiff00Adiff)∈ℝ2I×2I,A2=(AConv00AConv)+(-Λ10Λ1-Λ2)∈ℝ2I×2I.
We apply the splitting method to the operators A1 and A2 given in Section 2.1.The submatrices are(38)Adiff,i=DiΔx2R1=DiΔx2·(-211-21⋱⋱⋱1-211-2)∈ℝI×I,Aconv,i=-viΔxR2=-viΔx·(1-11⋱⋱-11-11)∈ℝI×I,
where I is the number of spatial points and Δx is the spatial step size. Consider
(39)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈ℝI×I,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈ℝI×I.
We have the commutator
(40)[A1,A2]=(00D2-D1Δx2R1Λ10)∈ℝ2I×2I,
where we assume that [R1,R2]≈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 Δx=0.1≫0, so that we scale the singular entry (D2-D1)/Δx2≈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].Figure3 presents the numerical error, that is, the difference between the exact and the numerical solution.Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)
(c)Figure4 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.CPU time of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)
(c)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 toD1, D2, and Δx with (D2-D1)/Δx2≈1, we can control the delicate blowup of the computed commutators. By the way, if we neglect the restriction and Δx with (D2-D1)/Δx2→∞, 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(41)∂tc1+∇·Fc1=g(-c1+c1,im)-λ1c1,inΩ×[0,t],∂tc2+∇·Fc2=g(-c2+c2,im)+λ1c1-λ2c2,inΩ×[0,t],F=v-D∇,∂tc1,im=g(c1-c1,im)-λ1c1,im,inΩ×[0,t],∂tc2,im=g(c2-c2,im)+λ1c1,im-λ2c2,im,inΩ×[0,t],c1(x,t)=c1,0(x),c2(x,t)=c2,0(x),onΩ,c1(x,t)=c1,1(x,t),c2(x,t)=c2,1(x,t),on∂Ω×[0,t],c1,im(x,t)=0,c2,im(x,t)=0,onΩ,c1,im(x,t)=0,c2,im(x,t)=0,on∂Ω×[0,t].
Here the parameters are v=0.1, D=0.01, λ1=λ2=0.1, and g=0.01.We next treat the semidiscretised equation given by the matrices(42)∂tC=(A-Λ1-G0G0Λ1A-Λ2-G0GG0-Λ1-G00GΛ1-Λ2-G)C,
where C=(c1,c2,c1im,c2im)T, while c1=(c1,1,…,c1,I) is the solution of the first species in the mobile phase in each spatial discretisation point (i=1,…,I), and the same for the other solution vectors.We have the following two operators for the splitting method:(43)A=DΔx2·(-211-21⋱⋱⋱1-211-2),+vΔx·(1-11⋱⋱-11-11)∈ℝI×I,
where I is the number of spatial points.Consider(44)Λ1=(λ100λ10⋱⋱⋱0λ100λ1),Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈ℝI×I,G=(g00g0⋱⋱⋱0g00g)∈ℝI×I.We decouple this into the following matrices:(45)A~=(A0000A0000000000)∈ℝ4I×4I,B1=(-Λ1000Λ1-Λ20000-Λ1000Λ1-Λ2),B2=(-G0G00-G0GG0-G00G0-G)∈ℝ4I×4I.
For the operators A~ and B~=B1+B2 and B we apply the iterative splitting method given in Section 2.1.Based on the decomposition,A~ is block diagonal and B~ is tridiagonal.The commutator is(46)[A~,B~]=(00AG0000AG-AG0000-AG00)∈ℝ4I×4I,
where we have a sparser operator than the main operator B~, 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 Δx=0.1≫0, so that we scale the singular entry D/Δx2≈1 and v/Δx≈1.Figure5 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 B~.Numerical errors and CPU time of the standard splitting scheme and the iterative schemes with Zassenhaus formula, where we have at least1,…,4 iterative steps.
(a)
(b)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 operatorB~.


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



---

*Source: 314290-2013-09-24.xml*