Embedded Zassenhaus Expansion to Splitting Schemes: Theory and Multiphysics Applications

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: 𝜕𝑐(𝑡)(𝑡)(𝑡)[](0)𝜕𝑡=𝐴𝑐+𝐵𝑐,with𝑡∈0,𝑇,𝑐=𝑐₀, (1) 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 𝑐(𝑡)=exp(𝐴+𝐵)𝑐(0),with𝑡∈[0,𝑇], (2) 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.

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: 𝜕𝑐_𝑖(𝑡)𝜕𝑡=𝐴𝑐_𝑖(𝑡)+𝐵𝑐_𝑖−1(𝑡),with𝑐_𝑖𝑡^𝑛=𝑐^𝑛,𝑖=1,2,…,𝑚. (3) Iterative splitting with respect to two operators (also called a two-sided scheme) as follows: 𝜕𝑐_𝑖(𝑡)𝜕𝑡=𝐴𝑐_𝑖(𝑡)+𝐵𝑐_𝑖−1(𝑡),with𝑐_𝑖𝑡^𝑛=𝑐^𝑛,𝜕𝑐_𝑖+1(𝑡)𝜕𝑡=𝐴𝑐_𝑖(𝑡)+𝐵𝑐_𝑖+1(𝑡),with𝑐_𝑖+1𝑡^𝑛=𝑐^𝑛,𝑖=1,3,…,2𝑚+1. (4)(5)

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: 𝑆_𝑖−exp((𝐴+𝐵)𝜏)≤𝐶𝜏^𝑖, (6) 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 𝑖=1,2,… and with the consistency error 𝑒_𝑖(𝜏)=𝑐(𝜏)−𝑐_𝑖(𝜏), we have 𝑐_𝑖(𝜏)=exp(𝐴𝜏)𝑐𝑡^𝑛+^𝑡𝑛+1_𝑡𝑛exp𝐴𝑡^𝑛+1−𝑠𝐵exp(𝑠𝐴)𝑐𝑡^𝑛+𝑑𝑠_{𝑡𝑛+1𝑡𝑛}exp𝐴𝑡^𝑛+1−𝑠₁𝐵×^{𝑡𝑛+1−𝑠1}_𝑡𝑛exp𝑡^𝑛+1−𝑠₁−𝑠₂𝐴×𝐵exp𝑠₂𝐴𝑐𝑡^𝑛𝑑𝑠₂𝑑𝑠₁+⋯+^𝑡𝑛+1_𝑡𝑛exp𝐴𝑡^𝑛+1−𝑠₁𝐵⋯^{𝑡𝑛+1−∑𝑖−2𝑗=1𝑠𝑗}_𝑡𝑛exp𝑡^𝑛+1−^𝑖−1_𝑗=1𝑠_𝑗𝐴×𝐵𝑐_init𝑠_𝑖𝑑𝑠_𝑖𝑑𝑠_𝑖−1,…,𝑑𝑠₂𝑑𝑠₁,𝑐(𝜏)=exp(𝐴𝜏)𝑐𝑡^𝑛+^𝑡𝑛+1_𝑡𝑛exp𝐴𝑡^𝑛+1−𝑠𝐵exp(𝑠𝐴)𝑐𝑡^𝑛+𝑑𝑠_{𝑡𝑛+1𝑡𝑛}exp𝐴𝑡^𝑛+1−𝑠₁𝐵×^{𝑡𝑛+1−𝑠1}_𝑡𝑛exp𝑡^𝑛+1−𝑠₁−𝑠₂𝐴×𝐵exp𝑠₂𝐴𝑐𝑡^𝑛𝑑𝑠₂𝑑𝑠₁+⋯+^𝑡𝑛+1_𝑡𝑛exp𝐴𝑡^𝑛+1−𝑠₁𝐵⋯^{𝑡𝑛+1−∑𝑖−1𝑗=1𝑠𝑗}_𝑡𝑛exp𝑡^𝑛+1−^𝑖_𝑗=1𝑠_𝑗𝐴×𝐵exp𝑠_𝑖(𝐴+𝐵)𝑐𝑡^𝑛𝑑𝑠_𝑖,…,𝑑𝑠₂𝑑𝑠₁. (7) We obtain 𝑒_𝑖≤exp((𝐴+𝐵)𝜏)𝑐𝑡^𝑛−𝑆_𝑖(𝜏)𝑐_init(𝜏)≤𝐶𝜏^𝑖max_{𝑠𝑖∈[]0,𝜏}exp𝑠_𝑖(𝐴+𝐵)𝑐𝑡^𝑛−𝑐_init𝑠_𝑖≤𝐶𝜏^𝑖exp(𝜏(𝐴+𝐵))𝑐𝑡^𝑛−𝑐_init(𝜏), (8) 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: 𝑐₁(𝜏)=exp(𝐴𝜏)𝑐𝑡^𝑛⟶𝑒_𝑖≤𝐶𝜏^𝑖𝑐𝑡^𝑛,𝑐₁(𝜏)=exp(𝐴𝜏)exp(𝐵𝜏)𝑐𝑡^𝑛⟶𝑒_𝑖≤𝐶𝜏^𝑖+1𝑐𝑡^𝑛. (9) So the slow convergence in the initial process 𝑐₁(𝜏) of the iterative splitting schemes can be improved by a cheap computable initial process 𝑐₁(𝜏)≈exp((𝐴+𝐵)𝜏). 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), exp((𝐴+𝐵)𝑡), can be uniquely decomposed as exp((𝐴+𝐵)𝑡)=exp(𝐴𝑡)exp(𝐵𝑡)𝜋^𝑗_𝑖=1𝜋^𝑚_𝑘=2×exp𝐶_𝑘(𝐴,𝐵)𝑡^𝑘+𝒪𝑡^𝑚+1, (10) 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 exp(−𝐴𝑡)exp((𝐴+𝐵)𝑡)=exp(𝐵𝑡+𝐷(𝑡)), (11) where 𝐷(𝑡) is a Lie polynomial of degree >1, and exp((𝐴+𝐵)𝑡)=exp(𝐴𝑡)exp(𝐵𝑡)+𝒪𝑡². (12) Further we have exp(−𝐵𝑡)exp(𝐵𝑡+𝐷(𝑡))=exp𝐶₂𝑡²+𝐷(𝑡), (13) where 𝐷(𝑡) is a Lie polynomial of degree >2, and exp((𝐴+𝐵)𝑡)=exp(𝐴𝑡)exp(𝐵𝑡)exp𝐶₂𝑡²+𝒪𝑡³. (14) We can repeat the recursive process and obtain by complete induction the results.

The first Zassenhaus polynomials are given as 𝐶₂()=1𝐴,𝐵2[],𝐶𝐵,𝐴₃()=1𝐴,𝐵3]]+1[[𝐵,𝐴,𝐵6]],𝐶[[𝐵,𝐴,𝐴₄()=1𝐴,𝐵]]]+124[[[𝐵,𝐴,𝐴,𝐴8]]]+1[[[𝐵,𝐴,𝐵,𝐴8[[[𝐵,𝐴],𝐵],𝐵]. (15) 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 exp𝐴𝑡2exp(𝐵𝑡)exp𝐴𝑡2=exp𝑡𝑆₁+𝑡³𝑆₃+𝑡⁵𝑆₅+⋯, (16) where the 𝑆_𝑖(𝐴,𝐵)∈ℒ(𝐴,𝐵) is a homogeneous Lie polynomial in 𝐴 and 𝐵 of degree 𝑖; see [21].
The BCH formula is given as 𝐴exp2+𝐵2𝑡=Π^∞_𝑘=2exp𝐶_𝑘𝑡^𝑘𝐴exp𝐵2𝑡exp,𝐵2𝑡exp2+𝐴2𝑡𝐵=exp𝐴2𝑡expΠ2𝑡^∞𝑘=2exp𝐶_𝑘𝑡^𝑘, (17) and so there exists a new product based on the Zassenhaus polynomials as follows: Π^∞_𝑘=3exp𝐷_𝑘𝑡^𝑘=Π^∞_𝑘=2exp𝐶_𝑘𝑡^𝑘Π^∞_𝑘=2exp𝐶_𝑘𝑡^𝑘, (18) 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 𝒟𝐿²⊂𝒟𝐴²∩𝒟(𝐴𝐵)∩𝒟(𝐵𝐴)∩𝒟𝐵², (19) where 𝐿=𝐴+𝐵, and 𝒟 is the domain of the operators.  (ii) Condition for the second term 𝐶₃ of the Zassenhaus formula 𝒟𝐿³⊂𝒟𝐴³∩𝒟𝐴²𝐵∩𝒟(𝐴𝐵𝐴)∩𝒟𝐴𝐵²∩𝒟𝐵³, (20) 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 𝒟(𝐿)=𝐻²(Ω)∩𝐻¹₀(Ω),𝒟(𝐴)=𝑐∈𝐿²(Ω);𝜕_𝑥𝑥𝑐,𝜕_𝑥𝑐∈𝐿²𝑢(,𝒟(𝐵)=Ω,0,𝑦=𝑢1,𝑦=0,∀𝑦∈0,1)𝑐∈𝐿²(Ω);𝜕_𝑦𝑦𝑐,𝜕_𝑦𝑐∈𝐿²𝑢()()(,Ω,𝑥,0=𝑢𝑥,1=0,∀𝑥∈0,1) (21) where Ω=(0,1)².
Then we need the following regularity result:  (i) condition for the first term 𝐶₂ of the Zassenhaus formula 𝒟𝐿²⊂𝐻⁴(Ω), (22) (ii) condition for the second term 𝐶₃ of the Zassenhaus formula 𝒟𝐿³⊂𝐻⁶(Ω). (23)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 𝐸_{Zassen,Comp,1}(𝑡)()(),𝐸=exp𝐴𝑡exp𝐵𝑡_{Zassen,Comp,𝑗}(𝑡)=exp𝐶_𝑗𝑡^𝑗,for𝑗∈2,…,𝑖, (24) where the sequential Zassenhaus operators are defined as 𝐸_Zassen,𝑖(𝑡)=Π^𝑖_𝑗=1𝐸_{Zassen,Comp,𝑗}(𝑡), (25) 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 𝑐_𝑖(𝑡)=𝐸_{𝑍𝑎𝑠𝑠𝑒𝑛,𝑖}(𝑡)𝑐₀, (26) 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: 𝑐_𝑖+𝑗(𝑡)=𝐸_{𝑖𝑡𝑒𝑟,𝑗}(𝑡)𝐸_{𝑍𝑎𝑠𝑠𝑒𝑛,𝑖}(𝑡)𝑐₀. (27) This has improved the error of the iterative scheme to 𝒪(𝑡^𝑖+𝑗).

Proof. The solution of the iterative splitting scheme (4) is 𝑐_𝑖+𝑗(𝑡)=𝐸_iter,𝑖(𝑡)𝑐_init,𝑗, (28) where 𝑆_𝑖(𝑡)=𝐸_iter,𝑖(𝑡), given in Theorem 1.
The initialisation is given by the Zassenhaus formula 𝑐_init,𝑗(𝑡)=𝐸_Zassen,𝑗(𝑡)𝑐₀. (29) Combining both splitting schemes, we have a local error of 𝑒_𝑖+𝑗(𝑡)=𝑐(𝑡)−𝑐_𝑖+𝑗(𝑡=𝑐(𝑡))−𝐸_iter,𝑖(𝑡)𝐸_Zassen,𝑗(𝑡)𝑐₀≤𝒪𝑡^𝑖+𝑗𝑐₀. (30)

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 𝑐(𝑡⁰) 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 𝑐₀(𝑡^𝑛+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 𝐶_𝑗. If err_Zass𝑗≤err₁, then we start with the iterative splitting scheme with 𝑐₀(𝑡)=𝑐_Zass𝑗, 𝑖=1 and go to (3), else 𝑗=𝑗+1 and we go to (1). (3) If the error of the iterative scheme err_iter≤err₂, then we are finished; else 𝑖=𝑖+1, and we compute the next iterative step, till err_iter≤err₂.

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.

Table 1 

Applied splitting schemes.

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 𝑢^′(𝑡)=(0)1230𝑢,𝑢=𝑢₀=01,𝑡∈[0,1], (31) the exact solution of which is 𝑢(𝑡)=2𝑒^3𝑡−𝑒^−2𝑡5. (32) We split the matrix as =[phone]=𝐴+𝐵, (33) 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.

(a)

(b)

Figure 1 

Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least 1,…,4 iterative steps.

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.

(a)

(b)

Figure 2 

CPU times of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least 1,…,4 iterative 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.

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: 𝜕_𝑡𝑅_𝑖𝑢_𝑖+∇⋅𝐯𝑢_𝑖=−𝜆_𝑖𝑅_𝑖𝑢_𝑖+𝜆_𝑖−1𝑅_𝑖−1𝑢_𝑖−1+𝛽−𝑢_𝑖+𝑔_𝑖(),𝑢,inΩ×0,𝑇_𝑖,0(𝑥)=𝑢_𝑖()𝜕𝑥,0,onΩ,_𝑡𝑅_𝑖𝑔_𝑖=−𝜆_𝑖𝑅_𝑖𝑔_𝑖+𝜆_𝑖−1𝑅_𝑖−1𝑔_𝑖−1+𝛽−𝑔_𝑖+𝑢_𝑖(),𝑔,inΩ×0,𝑇_𝑖,0(𝑥)=𝑔_𝑖(𝑥,0),onΩ,𝑖=1,…,𝑚, (34) where 𝑚 is the number of equations and 𝑖 is the index of each component. The unknown mobile concentrations 𝑢_𝑖=𝑢_𝑖(𝑥,𝑡) are taken in Ω×(0,𝑇)⊂ℝ^𝑛×ℝ⁺, where 𝑛 is the spatial dimension. The unknown immobile concentrations 𝑔_𝑖=𝑔_𝑖(𝑥,𝑡) are taken in Ω×(0,𝑇)⊂ℝ^𝑛×ℝ⁺, where 𝑛 is the spatial dimension. The retardation factors 𝑅_𝑖 are constant, and 𝑅_𝑖≥0. The kinetic part is given by the factors 𝜆_𝑖. They are constant and 𝜆_𝑖≥0. For the initialisation of the kinetic part, we set 𝜆₀=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 𝐯∈ℝ^𝑛 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 𝜕_𝑡𝑐₁+∇⋅𝐅₁𝑐₁=−𝜆₁𝑐₁[],𝜕,inΩ×0,𝑡_𝑡𝑐₂+∇⋅𝐅₂𝑐₂=𝜆₁𝑐₁−𝜆₂𝑐₂[],𝐅,inΩ×0,𝑡_𝑖=𝐯_𝑖−𝐷_𝑖𝑐∇,𝑖=1,2,₁(𝐱,𝑡)=𝑐_1,0(𝐱),𝑐₂(𝐱,𝑡)=𝑐_2,0(𝐱)𝑐,onΩ,₁(𝐱,𝑡)=𝑐_1,1(𝐱,𝑡),𝑐₂(𝐱,𝑡)=𝑐_2,1(),[].𝐱,𝑡on𝜕Ω×0,𝑡 (35) Here the parameters take the values 𝑣₁=0.1, 𝑣₂=0.05, 𝐷₁=0.01, 𝐷₂=0.005, and 𝜆₁=𝜆₂=0.1.

We now turn to the finite difference schemes and semidiscretise the equation, with its convection and diffusion operators as follows: 𝜕_𝑡𝐜=𝐴₁+𝐴₂𝐜. (36) We obtain the two matrices and decouple the diffusion and convection parts as follows: 𝐴₁=𝐴_diff00𝐴_diff∈ℝ^2𝐼×2𝐼,𝐴₂=𝐴_Conv00𝐴_Conv+−Λ₁0Λ₁−Λ₂∈ℝ^2𝐼×2𝐼. (37) We apply the splitting method to the operators 𝐴₁ and 𝐴₂ given in Section 2.1.

The submatrices are 𝐴_diff,𝑖=𝐷_𝑖Δ𝑥²𝑅₁=𝐷_𝑖Δ𝑥²⋅−211−21⋱⋱⋱1−211−2∈ℝ^𝐼×𝐼,𝐴_conv,𝑖𝑣=−^𝑖𝑅Δ𝑥²𝑣=−^𝑖⋅1Δ𝑥−11⋱⋱−11−11∈ℝ^𝐼×𝐼, (38) where 𝐼 is the number of spatial points and Δ𝑥 is the spatial step size. Consider Λ₁=𝜆₁00𝜆₁0⋱⋱⋱0𝜆₁00𝜆₁∈ℝ^𝐼×𝐼,Λ₂=𝜆₂00𝜆₂0⋱⋱⋱0𝜆₂00𝜆₂∈ℝ^𝐼×𝐼. (39) We have the commutator 𝐴₁,𝐴₂=𝐷00₂−𝐷₁Δ𝑥²𝑅₁Λ₁0∈ℝ^2𝐼×2𝐼, (40) where we assume that [𝑅₁,𝑅₂]≈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 (𝐷₂−𝐷₁)/Δ𝑥²≈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.

(a)

(b)

(c)

Figure 3 

Numerical errors of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least 1,…,4 iterative steps.

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.

(a)

(b)

(c)

Figure 4 

CPU time of the standard splitting scheme and the iterative schemes with the Zassenhaus formula, where we have at least 1,…,4 iterative 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 (𝐷₂−𝐷₁)/Δ𝑥²≈1, 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 𝜕_𝑡𝑐₁+∇⋅𝐅𝑐₁=𝑔−𝑐₁+𝑐_1,𝑖𝑚−𝜆₁𝑐₁[],𝜕,inΩ×0,𝑡_𝑡𝑐₂+∇⋅𝐅𝑐₂=𝑔−𝑐₂+𝑐_2,𝑖𝑚+𝜆₁𝑐₁−𝜆₂𝑐₂[],𝜕,inΩ×0,𝑡𝐅=𝐯−𝐷∇,_𝑡𝑐_1,𝑖𝑚=𝑔𝑐₁−𝑐_1,𝑖𝑚−𝜆₁𝑐_1,𝑖𝑚[],𝜕,inΩ×0,𝑡_𝑡𝑐_2,𝑖𝑚=𝑔𝑐₂−𝑐_2,𝑖𝑚+𝜆₁𝑐_1,𝑖𝑚−𝜆₂𝑐_2,𝑖𝑚[],𝑐,inΩ×0,𝑡₁(𝐱,𝑡)=𝑐_1,0(𝐱),𝑐₂(𝐱,𝑡)=𝑐_2,0(𝐱)𝑐,onΩ,₁(𝐱,𝑡)=𝑐_1,1(𝐱,𝑡),𝑐₂(𝐱,𝑡)=𝑐_2,1(),[],𝑐𝐱,𝑡on𝜕Ω×0,𝑡_1,𝑖𝑚(𝐱,𝑡)=0,𝑐_2,𝑖𝑚()𝑐𝐱,𝑡=0,onΩ,_1,𝑖𝑚(𝐱,𝑡)=0,𝑐_2,𝑖𝑚(𝐱,𝑡)=0,on𝜕Ω×[0,𝑡]. (41) Here the parameters are 𝑣=0.1, 𝐷=0.01, 𝜆₁=𝜆₂=0.1, and 𝑔=0.01.

We next treat the semidiscretised equation given by the matrices 𝜕_𝑡𝐂=𝐴−Λ₁Λ−𝐺0𝐺0₁𝐴−Λ₂−𝐺0𝐺𝐺0−Λ₁−𝐺00𝐺Λ₁−Λ₂−𝐺𝐂, (42) where 𝐂=(𝐜_𝟏,𝐜_𝟐,𝐜_𝟏𝑖𝑚,𝐜_𝟐𝑖𝑚)^𝑇, while 𝐜_𝟏=(𝑐_1,1,…,𝑐_1,𝐼) is the solution of the first species in the mobile phase in each spatial discretisation point (𝑖=1,…,𝐼), and the same for the other solution vectors.

We have the following two operators for the splitting method: 𝐷𝐴=Δ𝑥²⋅,+𝑣−211−21⋱⋱⋱1−211−2⋅1Δ𝑥−11⋱⋱−11−11∈ℝ^𝐼×𝐼, (43) where 𝐼 is the number of spatial points.

Consider Λ₁=𝜆₁00𝜆₁0⋱⋱⋱0𝜆₁00𝜆₁,Λ₂=𝜆₂00𝜆₂0⋱⋱⋱0𝜆₂00𝜆₂∈ℝ^𝐼×𝐼,𝐺=𝑔00𝑔0⋱⋱⋱0𝑔00𝑔∈ℝ^𝐼×𝐼. (44)

We decouple this into the following matrices: 𝐴=𝐴0000𝐴0000000000∈ℝ^4𝐼×4𝐼,𝐵₁=−Λ₁Λ000₁−Λ₂0000−Λ₁000Λ₁−Λ₂,𝐵₂=−𝐺0𝐺00−𝐺0𝐺𝐺0−𝐺00𝐺0−𝐺∈ℝ^4𝐼×4𝐼. (45) 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 𝐴,𝐵=00𝐴𝐺0000𝐴𝐺−𝐴𝐺0000−𝐴𝐺00∈ℝ^4𝐼×4𝐼, (46) 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 𝐷/Δ𝑥²≈1 and 𝑣/Δ𝑥≈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 one-sided iterative schemes on the operator 𝐵.

(a)

(b)

Figure 5 

Numerical errors and CPU time of the standard splitting scheme and the iterative schemes with Zassenhaus formula, where we have at least 1,…,4 iterative steps.

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.