# Computing Exponential for Iterative Splitting Methods: Algorithms and Applications

Authors: Jürgen Geiser
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

Iterative splitting methods have a huge amount to compute matrix exponential. Here, the
acceleration and recovering of higher-order schemes can be achieved. From a theoretical point of
view, iterative splitting methods are at least alternating Picards fix-point iteration schemes.
For practical applications, it is important to compute very fast matrix exponentials. In this
paper, we concentrate on developing fast algorithms to solve the iterative splitting scheme.
First, we reformulate the iterative splitting scheme into an integral notation of matrix exponential.
In this notation, we consider fast approximation schemes to the integral formulations,
also known asϕ-functions. Second, the error analysis is explained and applied to the integral
formulations. The novelty is to compute cheaply the decoupled exp-matrices and apply only
cheap matrix-vector multiplications for the higher-order terms. In general, we discuss an elegant way of embedding recently survey on methods for computing
matrix exponential with respect to iterative splitting schemes.
We present numerical benchmark examples, that compared standard splitting schemes
with the higher-order iterative schemes. A real-life application in contaminant transport as a
two phase model is discussed and the fast computations of the operator splitting method is
explained.


---

## Body

## 1. Introduction

We are motivated to solving multiple phase problems that arose of transport problem in porous media. In the last years, the interest in numerical simulations with multiphase problems, that can be used to model potential damage events has significantly increased in the area of final repositories for chemical or radioactive waste.Here, the modeling of the underlying porous media, which is the geosphere, spooled with water, is important, and we apply mobile and immobile pores as a two-phase problem in the media. With such a model, we achieve more realistic scenarios of the transported species, see [1]. The transport is structured in a convection and a diffusion portion. The convection portion is determined by the velocity vector v and indicates the spatial direction of the concentrations with the velocity. The diffusion portion is the spatial diffusion process of the concentrations, see [1].The sorption allegorizes the exchange between the solute (mobile) pollutant and at the surface sorbed (immobile) pollutants and appears in the temporal term as well as in the reaction term. The  reaction is reversible. The equilibrium-sorption therefore can be declared as a coefficient in the specific terms.We concentrate on simplified models and taken into account all real conditions for achieving statements of such realistic simulations of potential damage event. Here we have the delicate problem of coupled partial and ordinary differential equations and their underlying multiscale problems. While we deal with splitting methods and decomposing such scale problems, we can overcome such scaling problems, see [2].Moreover the computational part is important, while we dealing with large matrices and standard splitting schemes are expensive with respect to compute the exponential matrices. We solve the computational problem with a novel iterative splitting scheme, that concentrate on developing fast algorithms to solve an integral formulation of matrix exponentials. In such a scheme, we consider fast approximation schemes to the integral formulations, also known asϕ-functions, see [3]. The novelty is to compute cheaply the decoupled exp matrices and apply only cheap matrix-vector multiplications for the higher-order terms. In general, we discuss an elegant way of embedding recent survey on methods for computing matrix exponential with respect to iterative splitting schemes.In the following, we describe our model problem. The model equation for the multiple phase equations are given as coupled partial and ordinary differential equations:(1.1)ϕ∂tci+∇⋅Fi=g(-ci+ci,im)+kα(-ci+ci,ad)-λi,iϕci+∑k=k(i)λi,kϕck+Q̃i,inΩ×[0,t],(1.2)Fi=vci-De(i)∇ci+ηEci,(1.3)ϕ∂tci,im=g(ci-ci,im)+kα(ci,im,ad-ci,im)-λi,iϕci,im+∑k=k(i)λi,kϕck,im+Q̃i,im,inΩ×[0,t],(1.4)ϕ∂tci,ad=kα(ci-ci,ad)-λi,iϕci,ad+∑k=k(i)λi,kϕck,ad+Q̃i,ad,inΩ×[0,t],(1.5)ϕ∂tci,im,ad=kα(ci,im-ci,im,ad)-λi,iϕci,im,ad+∑k=k(i)λi,kϕck,im,ad+Q̃i,im,ad,inΩ×[0,t],(1.6)ci(x,t)=ci,0(x),ci,ad(x,t)=0,ci,im(x,t)=0,ci,im,ad(x,t)=0,on⁡Ω,(1.7)ci(x,t)=ci,1(x,t),ci,ad(x,t)=0,ci,im(x,t)=0,ci,im,ad(x,t)=0,on⁡∂Ω×[0,t],
where the initial value is given as ci,0 and we assume Dirichlet boundary conditions with the function ci,1(x,t) sufficiently smooth, all other initial and boundary conditions of the other phases zero.ϕ:
effective porosity(-),ci:
concentration of theith gaseous species in the plasma chamber phase (mol/cm3),ci,im:
concentration of theith gaseous species in the immobile zones of the plasma chamber phase (mol/cm3),ci,ad:
concentration of theith gaseous species in the mobile adsorbed zones of the plasma chamber phase (mol/cm3),ci,im,ad:
concentration of theith gaseous species in the immobile adsorbed zones of the plasma chamber phase (mol/cm3),v:
velocity through the chamber and porous substrate [4] (cm/nsec),De(i):
element-specific diffusions-dispersions tensor (cm2/nsec),λi,i:
decay constant of theith species (1/nsec),Q̃i:
source term of theith species (mol/(cm3nsec)),g:
exchange rate between the mobile and immobile concentration (1/nsec),kα:
exchange rate between the mobile and adsorbed concentration or immobile and immobile adsorbed concentration (kinetic controlled sorption) (1/nsec),E:
stationary electric field in the apparatus (V/cm),η:
the mobility rate in the electric field, see [5] (cm2/nsec).withi=1,…,M and M denotes the number of components.The parameters in (1.1) are further described, see also [1].The effective porosity is denoted byϕ and declares the portion of the porosities of the aquifer, that is filled with plasma, and we assume a nearly fluid phase. The transport term is indicated by the Darcy velocity v, that presents the flow-direction and the absolute value of the plasma flux. The velocity field is divergence-free. The decay constant of the ith species is denoted by λi. Thereby does k(i) denote the indices of the other species.In this paper we concentrate on solving linear evolution equations, such as the differential equation,(1.8)∂tu=Au,u(0)=u0,
where A can be an unbounded operator. We assume to achieve large-scale differential equation, which are delicate to solve with standard solvers.Our main focus will be to consider and contrast higher-order algorithms derived from standard schemes as for example Strang-splitting schemes.We propose iterative splitting schemes as a solver scheme which is simple to implement and parallelisible.A rewriting to integral formulation, allows to reduce the computation to numerical approximation schemes. While standard schemes has the disadvantage to compute commutators to achieve higher-order schemes, we could speed up our schemes by recursive algorithms. Iterative schemes can be seen as Successive approximations, which are based on recursive integral formulations in which an iterative method is enforce the time dependency.The paper is outlined as follows. In Section2, we present the underlying splitting methods. The algorithms for the exponentials are given in Section 3. Numerical verifications are given in Section 4. In Section 5, we briefly summarize our results.


## 2. Iterative Splitting Method

The following algorithm is based on the iteration with fixed-splitting discretization step-sizeτ, namely, on the time interval [tn,tn+1] we solve the following subproblems consecutively for i=0,2,…,2m. (cf. [6, 7]):(2.1)∂ci(t)∂t=Aci(t)+Bci-1(t),withci(tn)=cn,c0(tn)=cn,c-1=0.0,∂ci+1(t)∂t=Aci(t)+Bci+1(t),withci+1(tn)=cn,
where cn is the known split approximation at the time level t=tn. The split approximation at the time level t=tn+1 is defined as cn+1=c2m+1(tn+1), (clearly, the function ci+1(t) depends on the interval [tn,tn+1], too, but, for the sake of simplicity, in our notation we omit the dependence on n).In the following we will analyze the convergence and the rate of convergence of the method (2.1) for m tends to infinity for the linear operators A,B:X→X, where we assume that these operators and their sum are generators of the C0 semigroups. We emphasize that these operators are not necessarily bounded, so the convergence is examined in a general Banach space setting.The novelty of the convergence results are the reformulation in integral-notation. Based on this, we can assume to have bounded integral operators which can be estimated and given in a recursive form. Such formulations are known in the work of [8, 9] and estimations of the kernel part with the exponential operators are sufficient to estimate the recursive formulations.


### 2.1. Error Analysis

We present the results of the consistency of our iterative method. We assume for the system of operator the generator of aC0 semigroup based on their underlying graph norms.Theorem 2.1.
Let one consider the abstract Cauchy problem in a Hilbert spaceX(2.2)∂tc(x,t)=Ac(x,t)+Bc(x,t),0<t≤T,x∈Ω,c(x,0)=c0(x),x∈Ω,c(x,t)=c1(x,t),x∈∂Ω×[0,T],
where A,B:D(X)→X are given linear operators which are generators of the C0-semigroup and c0∈X is a given element. We assume A, B are unbounded. Further, we assume the estimations of the unbounded operator B with sufficient smooth initial conditions (see [8]):
(2.3)‖Bexp⁡((A+B)τ)u0‖≤κ.
Further we assume the estimation of the partial integration of the unbounded operatorB (see [8]):(2.4)‖B∫0τexp⁡(Bs)sds‖≤τC.
Then, we can bound our iterative operator splitting method as(2.5)‖(Si-exp⁡((A+B)τ))‖≤Cτi,
where Si is the approximated solution for the ith iterative step and C is a constant that can be chosen uniformly on bounded time intervals.Proof.
Let us consider the iteration (2.1) on the subinterval [tn,tn+1].
For the first iterations, we have(2.6)∂tc1(t)=Ac1(t),t∈(tn,tn+1],
and for the second iteration we have:
(2.7)∂tc2(t)=Ac1(t)+Bc2(t),t∈(tn,tn+1].
In general, we have:
(i)
for the odd iterations,i=2m+1 for m=0,1,2,…(2.8)∂tci(t)=Aci(t)+Bci-1(t),t∈(tn,tn+1],
where for c0(t)≡0,(ii)
for the even iterations,i=2m for m=1,2,…(2.9)∂tci(t)=Aci-1(t)+Bci(t),t∈(tn,tn+1].
We have the following solutions for the iterative scheme: the solutions for the first two equations are given by the variation of constants:(2.10)c1(t)=exp⁡(A(tn+1-t))c(tn),t∈(tn,tn+1],c2(t)=exp⁡(Bt)c(tn)+∫tntn+1exp⁡(B(tn+1-s))Ac1(s)ds,t∈(tn,tn+1].
For the recurrence relations with even and odd iterations, we have the solutions: for the odd iterations:i=2m+1, for m=0,1,2,…,(2.11)ci(t)=exp⁡(A(t-tn))c(tn)+∫tntexp⁡((tn+1-s)A)Bci-1(s)ds,t∈(tn,tn+1].
For the even iterations,i=2m, for m=1,2,…(2.12)ci(t)=exp⁡(B(t-tn))c(tn)+∫tntexp⁡((tn+1-s)B)Aci-1(s)ds,t∈(tn,tn+1].
The consistency is given as the following.
Fore1, we have
(2.13)c1(τ)=exp⁡(Aτ)c(tn),c(τ)=exp⁡((A+B)τ)c(tn)=exp⁡(Aτ)c(tn)+∫tntn+1exp⁡(A(tn+1-s))Bexp⁡(s(A+B))c(tn)ds.
We obtain(2.14)‖e1‖=‖c-c1‖≤‖exp⁡((A+B)τ)c(tn)-exp⁡(Aτ)c(tn)‖≤C1τc(tn).Fore2 we have (alternating)
(2.15)c2(τ)=exp⁡(Bτ)c(tn)+∫tntn+1exp⁡(B(tn+1-s))Aexp⁡(sA)c(tn)ds,c(τ)=exp⁡(Bτ)c(tn)+∫tntn+1exp⁡(B(tn+1-s))Aexp⁡(sA)c(tn)ds+∫tntn+1exp⁡(Bs)A∫tntn+1-sexp⁡(A(tn+1-s-ρ))Bexp⁡(ρ(A+B))c(tn)dρds.
We obtain(2.16)‖e2‖≤‖exp⁡((A+B)τ)c(tn)-c2‖≤C2τ2c(tn).
For odd and even iterations, the recursive proof is given in the following. For  the odd iterations (only iterations overA), i=2m+1 for m=0,1,2,…, for ei, we have:(2.17)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-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)B×exp⁡(siA)c(tn)dsi-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-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=1isj)A)B×exp⁡(siA)c(tn)dsi⋯ds2ds1.
We obtain(2.18)‖ei‖≤‖exp⁡((A+B)τ)c(tn)-ci‖≤Cτic(tn),
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.2.
The same idea can be applied toA=∇D∇B=-v·∇, so that one operator is less unbounded, but we reduce the convergence order
(2.19)‖e1‖=K‖B‖τα1‖e0‖+O(τ1+α1),andhence,‖e2‖=K‖B‖‖e0‖τ1+α1+α2+O(τ1+α1+α),
where 0≤α1, α2<1.Remark 2.3.
If we assume the consistency of𝒪(τm) for the initial value e1(tn) and e2(tn), we can redo the proof and obtain at least a global error of the splitting methods of 𝒪(τm-1).


### 2.2. Splitting Method to Couple Mobile, Immobile, and Adsorbed Parts

The motivation of the splitting method are based on the following observations.(i)
The mobile phase is semidiscretized with fast finite volume methods and can be stored into a stiffness-matrix. We achieve large time steps, if we consider implicit Runge-Kutta methods of lower order (e.g., implicit Euler) as a time discretization method.(ii)
The immobile, adsorbed and immobile-adsorbed phases are purely ordinary differential equations and each of them is cheap to solve with explicit Runge-Kutta schemes.(iii)
The ODEs can be seen as perturbations and can be solved all explicit in a fast iterative scheme.For the full equation we consider the following matrix notation:(2.20)∂tc=A1c+A2c+B1(c-cim)+B2(c-cad)+Q,∂tcim=A2cim+B1(cim-c)+B2(cim-cim,ad)+Qim,∂tcad=A2cad+B2(cad-c)+Qad,∂tcim,ad=A2cim,ad+B2(cim,ad-cim)+Qim,ad,
where c=(c1,…,cm)T is the spatial discretized concentration in the mobile phase, see (1.1), cim=(c1,im,…,cm,im)T is the concentration in the immobile phase, the some also for the other phase concentrations. A1 is the stiffness matrix of (1.1), A2 is the reaction matrix of the right-hand side of (1.1), B1 and B2 are diagonal matrices with the exchange of the immobile and kinetic parameters, see (1.4) and (1.5).Further,Q,…,Qim,ad are the spatial discretized sources vectors.Now, we have the following ordinary differential equation:(2.21)∂tC=(A1+A2+B1+B2-B1-B20-B1A2+B1+B20-B2-B20A2+B200-B20A2+B2)C+Q̃,
where C=(c,cim,cad,cim,ad)T and the right-hand side is given as Q̃=(Q,Qim,Qad,Qim,ad)T.For such an equation, we apply the decomposition of the matrices:(2.22)∂tC=ÃC+Q̃,∂tC=Ã1C+Ã2C+Q̃,
where(2.23)Ã1=(A1+A20000A20000A20000A2),Ã2=(B1+B2-B1-B20-B1B1+B20-B2-B20B200-B20B2).The equation system is numerically solved by an iterative scheme.Algorithm 2.4.
We divide our time interval[0,T] into subintervals [tn,tn+1], where n=0,1,…N, t0=0, and tN=T.
We start withn=0.
(1)
The initial conditions are given withC0(tn+1)=C(tn). We start with k=0.(2)
Compute the fix-point iteration scheme given as(2.24)∂tCk=Ã1Ck+Ã2Ck-1+Q̃,
where k is the iteration index, see [10]. For the time integration, we apply Runge-Kutta methods as ODE solvers, see [11, 12].(3)
The stop criterion for the time interval[tn,tn+1] is given as
(2.25)‖Ck(tn+1)-Ck-1(tn+1)‖≤err⁡,
where ∥·∥ is the maximum norm over all components of the solution vector. err is a given error bound, for example, err=10-4.If (2.25) is fulfilled, we have the result
(2.26)C(tn+1)=Ck(tn+1).Ifn=N then we stop and are done.If (2.25) is not fulfilled, we do k=k+1 and go to (2).The error analysis of the schemes are given in the following Theorem.Theorem 2.5.
LetA,B∈ℒ(X) be given linear bounded operators in a Banach space ℒ(X). We consider the abstract Cauchy problem:
(2.27)∂tC(t)=ÃC(t)+B̃C(t),tn≤t≤tn+1,C(tn)=Cn,forn=1,…,N,
where t1=0 and the final time is tN=T∈ℝ+. Then problem (2.27) has a unique solution. For a finite steps with time size τn=tn+1-tn, the iteration (2.24) for k=1,2,…,q is consistent with an order of consistency 𝒪(τnq).Proof.
The outline of the proof is given in [2].In the next section we describe the computation of the integral formulation withexp-functions.


## 2.1. Error Analysis

We present the results of the consistency of our iterative method. We assume for the system of operator the generator of aC0 semigroup based on their underlying graph norms.Theorem 2.1.
Let one consider the abstract Cauchy problem in a Hilbert spaceX(2.2)∂tc(x,t)=Ac(x,t)+Bc(x,t),0<t≤T,x∈Ω,c(x,0)=c0(x),x∈Ω,c(x,t)=c1(x,t),x∈∂Ω×[0,T],
where A,B:D(X)→X are given linear operators which are generators of the C0-semigroup and c0∈X is a given element. We assume A, B are unbounded. Further, we assume the estimations of the unbounded operator B with sufficient smooth initial conditions (see [8]):
(2.3)‖Bexp⁡((A+B)τ)u0‖≤κ.
Further we assume the estimation of the partial integration of the unbounded operatorB (see [8]):(2.4)‖B∫0τexp⁡(Bs)sds‖≤τC.
Then, we can bound our iterative operator splitting method as(2.5)‖(Si-exp⁡((A+B)τ))‖≤Cτi,
where Si is the approximated solution for the ith iterative step and C is a constant that can be chosen uniformly on bounded time intervals.Proof.
Let us consider the iteration (2.1) on the subinterval [tn,tn+1].
For the first iterations, we have(2.6)∂tc1(t)=Ac1(t),t∈(tn,tn+1],
and for the second iteration we have:
(2.7)∂tc2(t)=Ac1(t)+Bc2(t),t∈(tn,tn+1].
In general, we have:
(i)
for the odd iterations,i=2m+1 for m=0,1,2,…(2.8)∂tci(t)=Aci(t)+Bci-1(t),t∈(tn,tn+1],
where for c0(t)≡0,(ii)
for the even iterations,i=2m for m=1,2,…(2.9)∂tci(t)=Aci-1(t)+Bci(t),t∈(tn,tn+1].
We have the following solutions for the iterative scheme: the solutions for the first two equations are given by the variation of constants:(2.10)c1(t)=exp⁡(A(tn+1-t))c(tn),t∈(tn,tn+1],c2(t)=exp⁡(Bt)c(tn)+∫tntn+1exp⁡(B(tn+1-s))Ac1(s)ds,t∈(tn,tn+1].
For the recurrence relations with even and odd iterations, we have the solutions: for the odd iterations:i=2m+1, for m=0,1,2,…,(2.11)ci(t)=exp⁡(A(t-tn))c(tn)+∫tntexp⁡((tn+1-s)A)Bci-1(s)ds,t∈(tn,tn+1].
For the even iterations,i=2m, for m=1,2,…(2.12)ci(t)=exp⁡(B(t-tn))c(tn)+∫tntexp⁡((tn+1-s)B)Aci-1(s)ds,t∈(tn,tn+1].
The consistency is given as the following.
Fore1, we have
(2.13)c1(τ)=exp⁡(Aτ)c(tn),c(τ)=exp⁡((A+B)τ)c(tn)=exp⁡(Aτ)c(tn)+∫tntn+1exp⁡(A(tn+1-s))Bexp⁡(s(A+B))c(tn)ds.
We obtain(2.14)‖e1‖=‖c-c1‖≤‖exp⁡((A+B)τ)c(tn)-exp⁡(Aτ)c(tn)‖≤C1τc(tn).Fore2 we have (alternating)
(2.15)c2(τ)=exp⁡(Bτ)c(tn)+∫tntn+1exp⁡(B(tn+1-s))Aexp⁡(sA)c(tn)ds,c(τ)=exp⁡(Bτ)c(tn)+∫tntn+1exp⁡(B(tn+1-s))Aexp⁡(sA)c(tn)ds+∫tntn+1exp⁡(Bs)A∫tntn+1-sexp⁡(A(tn+1-s-ρ))Bexp⁡(ρ(A+B))c(tn)dρds.
We obtain(2.16)‖e2‖≤‖exp⁡((A+B)τ)c(tn)-c2‖≤C2τ2c(tn).
For odd and even iterations, the recursive proof is given in the following. For  the odd iterations (only iterations overA), i=2m+1 for m=0,1,2,…, for ei, we have:(2.17)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-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)B×exp⁡(siA)c(tn)dsi-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-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=1isj)A)B×exp⁡(siA)c(tn)dsi⋯ds2ds1.
We obtain(2.18)‖ei‖≤‖exp⁡((A+B)τ)c(tn)-ci‖≤Cτic(tn),
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.2.
The same idea can be applied toA=∇D∇B=-v·∇, so that one operator is less unbounded, but we reduce the convergence order
(2.19)‖e1‖=K‖B‖τα1‖e0‖+O(τ1+α1),andhence,‖e2‖=K‖B‖‖e0‖τ1+α1+α2+O(τ1+α1+α),
where 0≤α1, α2<1.Remark 2.3.
If we assume the consistency of𝒪(τm) for the initial value e1(tn) and e2(tn), we can redo the proof and obtain at least a global error of the splitting methods of 𝒪(τm-1).


## 2.2. Splitting Method to Couple Mobile, Immobile, and Adsorbed Parts

The motivation of the splitting method are based on the following observations.(i)
The mobile phase is semidiscretized with fast finite volume methods and can be stored into a stiffness-matrix. We achieve large time steps, if we consider implicit Runge-Kutta methods of lower order (e.g., implicit Euler) as a time discretization method.(ii)
The immobile, adsorbed and immobile-adsorbed phases are purely ordinary differential equations and each of them is cheap to solve with explicit Runge-Kutta schemes.(iii)
The ODEs can be seen as perturbations and can be solved all explicit in a fast iterative scheme.For the full equation we consider the following matrix notation:(2.20)∂tc=A1c+A2c+B1(c-cim)+B2(c-cad)+Q,∂tcim=A2cim+B1(cim-c)+B2(cim-cim,ad)+Qim,∂tcad=A2cad+B2(cad-c)+Qad,∂tcim,ad=A2cim,ad+B2(cim,ad-cim)+Qim,ad,
where c=(c1,…,cm)T is the spatial discretized concentration in the mobile phase, see (1.1), cim=(c1,im,…,cm,im)T is the concentration in the immobile phase, the some also for the other phase concentrations. A1 is the stiffness matrix of (1.1), A2 is the reaction matrix of the right-hand side of (1.1), B1 and B2 are diagonal matrices with the exchange of the immobile and kinetic parameters, see (1.4) and (1.5).Further,Q,…,Qim,ad are the spatial discretized sources vectors.Now, we have the following ordinary differential equation:(2.21)∂tC=(A1+A2+B1+B2-B1-B20-B1A2+B1+B20-B2-B20A2+B200-B20A2+B2)C+Q̃,
where C=(c,cim,cad,cim,ad)T and the right-hand side is given as Q̃=(Q,Qim,Qad,Qim,ad)T.For such an equation, we apply the decomposition of the matrices:(2.22)∂tC=ÃC+Q̃,∂tC=Ã1C+Ã2C+Q̃,
where(2.23)Ã1=(A1+A20000A20000A20000A2),Ã2=(B1+B2-B1-B20-B1B1+B20-B2-B20B200-B20B2).The equation system is numerically solved by an iterative scheme.Algorithm 2.4.
We divide our time interval[0,T] into subintervals [tn,tn+1], where n=0,1,…N, t0=0, and tN=T.
We start withn=0.
(1)
The initial conditions are given withC0(tn+1)=C(tn). We start with k=0.(2)
Compute the fix-point iteration scheme given as(2.24)∂tCk=Ã1Ck+Ã2Ck-1+Q̃,
where k is the iteration index, see [10]. For the time integration, we apply Runge-Kutta methods as ODE solvers, see [11, 12].(3)
The stop criterion for the time interval[tn,tn+1] is given as
(2.25)‖Ck(tn+1)-Ck-1(tn+1)‖≤err⁡,
where ∥·∥ is the maximum norm over all components of the solution vector. err is a given error bound, for example, err=10-4.If (2.25) is fulfilled, we have the result
(2.26)C(tn+1)=Ck(tn+1).Ifn=N then we stop and are done.If (2.25) is not fulfilled, we do k=k+1 and go to (2).The error analysis of the schemes are given in the following Theorem.Theorem 2.5.
LetA,B∈ℒ(X) be given linear bounded operators in a Banach space ℒ(X). We consider the abstract Cauchy problem:
(2.27)∂tC(t)=ÃC(t)+B̃C(t),tn≤t≤tn+1,C(tn)=Cn,forn=1,…,N,
where t1=0 and the final time is tN=T∈ℝ+. Then problem (2.27) has a unique solution. For a finite steps with time size τn=tn+1-tn, the iteration (2.24) for k=1,2,…,q is consistent with an order of consistency 𝒪(τnq).Proof.
The outline of the proof is given in [2].In the next section we describe the computation of the integral formulation withexp-functions.


## 3. Exponentials to Compute Iterative Schemes

The theoretical ideas can be discussed in the following formulation:(3.1)DA(B,t)[0]=exp⁡(tB),DA(B,t)[k]=k∫0texp⁡((t-s)B)ADA(B,s)[k-1]ds,
and the matrix formulation of our two-step scheme is given as(3.2)Ã=(A0⋯⋯AB0⋯0BA⋯⋮⋱⋱⋱0⋯BA),
the computation of the exp-Matrix can be expressed as:(3.3)exp⁡(Ãt):=(exp⁡(At)0⋯⋯DA(B)exp⁡(Bt)0⋯DB(A)[2]2!DB(A)[1]1!exp⁡(At)⋯⋮⋱⋱⋱DB(A)[n]n!⋯DB(A)[1]1!exp⁡(At)).Here, we have to compute the right-hand side as time dependent term, means we evaluateexp(At) and exp(Bt) as a Taylor expansion. Then the integral formulation can be done with polynomials and we derive the expansion of such integral operators with the ϕ-function.Theϕ-functions are defined as the integration of the exp-functions. We can analytically derive ϕ-functions with respect to the integral-formulation of a matrix exp-function.In the following we reduce to a approximation of the fixed right-hand side (means we assumeDA(B,s)[k-1]≈DA(B,0)[k-01]).Later we also follow with more extended schemes.Computingϕ-Functions:
(3.4)ϕ0(x)=ex,ϕk+1(x)=ϕk(x)-1/k!x,k≥0.
So the matrix formulation of our scheme is given as(3.5)Y(t)=exp⁡(Ãt)Y(0),
where Ã is given as,
(3.6)Ã=(A0⋯⋯AB0⋯0AB⋯⋮⋱⋱⋱0⋯BA),
the computation of the exp-Matrix can be expressed in a first-order scheme and with the assumption of commutations is given as:
(3.7)exp⁡(Ãt):=(exp⁡(At)0⋯⋯⋯ϕ1(Bt)Atexp⁡(Bt)0⋯⋯ϕ2(At)B(B+A)t2ϕ1(At)Btexp⁡(At)⋯⋯⋮⋱⋱⋱⋮ϕi(At)B(A+B)i-1ti⋯ϕ1(At)Btexp⁡(At)),
where we assume ui(s)=u(tn), mean we approximate the right-hand side term.For higher-orders we should also include the full derivations ofc1,c2,….


### 3.1. Derivation of the Formula

Consider the equation(3.8)ẋ=Ax+a,x(0)=0,
the solution of this equation in terms of the ϕk,k=0,1, is given as(3.9)x=tϕ1(tA)a.
Similarly, the equation(3.10)ẋ=Ax+bt,x(0)=0,
has a solution of the form(3.11)x=t2ϕ2(tA)b.
In general, the solution of the equation(3.12)ẋ=Ax+a+bt+ct22!+⋯,x(0)=0,
admits the form(3.13)x=tϕ1(tA)a+t2ϕ2(tA)b+t3ϕ3(tA)c+⋯.
In this section, we use the formula (3.21) for iterative schemes given as(3.14)u̇1=Au1,u̇2=Au1+Bu2,u̇3=Au3+Bu2,u̇4=⋯.
The solution of the first iteration given by(3.15)u1=eAtu0.
Inserting this into (3.26) and expanding eAt up to the order 1, we have 2nd order approximation of the exact solution which has a form(3.16)u2=eBtu0+ϕ1(tA)Atu0.
Similarly, inserting (3.16) into (3.27) and expanding ϕ1(tA), we have a third order approximation of the exact solution(3.17)u3=eBtu0+ϕ1(tB)Btu0+ϕ2(tB)B(B+A)t2u0.
In general for i=0,2,4,…, we have the p=i+1-th order approximation of the exact solution in terms of the ϕi function as follows:(3.18)ui=eAtu0+ϕ1(tA)Atu0+ϕ2(tA)A(B+A)t2u0+⋯+ϕi(tA)A(A+B)i-1u0.
For i=1,3,5,…, we have(3.19)ui=eBtu0+ϕ1(tB)Btu0+ϕ2(tB)B(B+A)t2u0+⋯+ϕi(tB)B(B+A)i-1tiu0.


### 3.2. Algorithms

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


### 3.3. Two-Side Scheme

For the implementation of the integral formulation, we have to deal with parallel ideas, which means we select independent parts of the formulation.Step 1.
Determine the order of the method by fixed iteration number.Step 2.
Consider the time interval[t0,T], divide it into N subintervals so that time step is h=(T-t0)/N.Step 3.
On each subinterval,[tn,tn+h], n=0,1,…,N, use the algorithm by considering initial conditions for each step as u(t0)=u0, ui(tn)=ui-1(tn)=u(tn),
(3.20)u2(tn+h)=(ϕ0(Bt)+ϕ1(Bt)At)u(tn)u3(tn+h)=(ϕ0(At)+ϕ1(At)Bt+ϕ2(At)B(B+A)t2)u(tn)⋮u2i(tn+h)=(ϕ0(Bt)+ϕ1(Bt)At+⋯+ϕ2i-1(Bt)A(A+B)2i-1t2i)u2i-1(tn)u2i+1(tn+h)=(ϕ0(At)+ϕ1(At)Bt+⋯+ϕ2i(At)B(B+A)2it2i+1)u2i(tn).Step 4.
ui(tn+h)→u(tn+h).Step 5.
Repeat this procedure for next interval until the desired timeT is reached.


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

For the one-side scheme, we taken into account of the following commutator relation.Definition 3.1.
The following relation is given:(3.21)exp⁡(-tA)Bexp⁡(tA)=[B,exp⁡(tA)],
where [·,·] is the matrix commutator.
The integration is given as(3.22)∫0texp⁡(-sA)Bexp⁡(sA)ds=[B,ϕ1(tA)],
the representation is given in [13].
Further, we have the recursive integration(3.23)[B,ϕi(tA)ti]=∫0t[B,ϕi-1(sA)]ds,
where ϕi given defined as
(3.24)ϕ0(At)=exp⁡(At),ϕi(At)=∫0tϕi-1(As)ds,ϕi(At)=ϕi-1(At)-I(ti-1/(i-1)!)A.The iterative scheme with the equations(3.25)u̇1=Au1,(3.26)u̇2=Au2+Bu1,(3.27)u̇3=Au3+Bu2,(3.28)u̇4=⋯,
is solved as(3.29)c1(t)=exp⁡(At)c(tn),c2(t)=c1(t)+c1(t)∫0t[B,exp⁡(sA)]ds,c2(t)=c1(t)+c1(t)[B,ϕ1(tA)],c3(t)=c2(t)+c1(t)∫0t[B,exp⁡(sA),B,ϕ1(sA)]ds,c3(t)=c2(t)+c1(t)([B,exp⁡(tA),B,ϕ2(tA)]+[B,Aexp⁡(tA),B,ϕ3(tA)])+O(t3),⋯The recursion is given as(3.30)ci(t)=ci-1(t)+c1(t)∫0t[B,exp⁡(sA)]∫0s1[B,exp⁡(s1A)]∫0s2⋯∫0si-2[B,exp⁡(si-1A)]dsi-1⋯ds1dt.Remark 3.2.
For the novel notation, we have embedded the commutator to the computational scheme. For such a scheme we could save to compute additional the commutators.


## 3.1. Derivation of the Formula

Consider the equation(3.8)ẋ=Ax+a,x(0)=0,
the solution of this equation in terms of the ϕk,k=0,1, is given as(3.9)x=tϕ1(tA)a.
Similarly, the equation(3.10)ẋ=Ax+bt,x(0)=0,
has a solution of the form(3.11)x=t2ϕ2(tA)b.
In general, the solution of the equation(3.12)ẋ=Ax+a+bt+ct22!+⋯,x(0)=0,
admits the form(3.13)x=tϕ1(tA)a+t2ϕ2(tA)b+t3ϕ3(tA)c+⋯.
In this section, we use the formula (3.21) for iterative schemes given as(3.14)u̇1=Au1,u̇2=Au1+Bu2,u̇3=Au3+Bu2,u̇4=⋯.
The solution of the first iteration given by(3.15)u1=eAtu0.
Inserting this into (3.26) and expanding eAt up to the order 1, we have 2nd order approximation of the exact solution which has a form(3.16)u2=eBtu0+ϕ1(tA)Atu0.
Similarly, inserting (3.16) into (3.27) and expanding ϕ1(tA), we have a third order approximation of the exact solution(3.17)u3=eBtu0+ϕ1(tB)Btu0+ϕ2(tB)B(B+A)t2u0.
In general for i=0,2,4,…, we have the p=i+1-th order approximation of the exact solution in terms of the ϕi function as follows:(3.18)ui=eAtu0+ϕ1(tA)Atu0+ϕ2(tA)A(B+A)t2u0+⋯+ϕi(tA)A(A+B)i-1u0.
For i=1,3,5,…, we have(3.19)ui=eBtu0+ϕ1(tB)Btu0+ϕ2(tB)B(B+A)t2u0+⋯+ϕi(tB)B(B+A)i-1tiu0.


## 3.2. Algorithms

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


## 3.3. Two-Side Scheme

For the implementation of the integral formulation, we have to deal with parallel ideas, which means we select independent parts of the formulation.Step 1.
Determine the order of the method by fixed iteration number.Step 2.
Consider the time interval[t0,T], divide it into N subintervals so that time step is h=(T-t0)/N.Step 3.
On each subinterval,[tn,tn+h], n=0,1,…,N, use the algorithm by considering initial conditions for each step as u(t0)=u0, ui(tn)=ui-1(tn)=u(tn),
(3.20)u2(tn+h)=(ϕ0(Bt)+ϕ1(Bt)At)u(tn)u3(tn+h)=(ϕ0(At)+ϕ1(At)Bt+ϕ2(At)B(B+A)t2)u(tn)⋮u2i(tn+h)=(ϕ0(Bt)+ϕ1(Bt)At+⋯+ϕ2i-1(Bt)A(A+B)2i-1t2i)u2i-1(tn)u2i+1(tn+h)=(ϕ0(At)+ϕ1(At)Bt+⋯+ϕ2i(At)B(B+A)2it2i+1)u2i(tn).Step 4.
ui(tn+h)→u(tn+h).Step 5.
Repeat this procedure for next interval until the desired timeT is reached.


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

For the one-side scheme, we taken into account of the following commutator relation.Definition 3.1.
The following relation is given:(3.21)exp⁡(-tA)Bexp⁡(tA)=[B,exp⁡(tA)],
where [·,·] is the matrix commutator.
The integration is given as(3.22)∫0texp⁡(-sA)Bexp⁡(sA)ds=[B,ϕ1(tA)],
the representation is given in [13].
Further, we have the recursive integration(3.23)[B,ϕi(tA)ti]=∫0t[B,ϕi-1(sA)]ds,
where ϕi given defined as
(3.24)ϕ0(At)=exp⁡(At),ϕi(At)=∫0tϕi-1(As)ds,ϕi(At)=ϕi-1(At)-I(ti-1/(i-1)!)A.The iterative scheme with the equations(3.25)u̇1=Au1,(3.26)u̇2=Au2+Bu1,(3.27)u̇3=Au3+Bu2,(3.28)u̇4=⋯,
is solved as(3.29)c1(t)=exp⁡(At)c(tn),c2(t)=c1(t)+c1(t)∫0t[B,exp⁡(sA)]ds,c2(t)=c1(t)+c1(t)[B,ϕ1(tA)],c3(t)=c2(t)+c1(t)∫0t[B,exp⁡(sA),B,ϕ1(sA)]ds,c3(t)=c2(t)+c1(t)([B,exp⁡(tA),B,ϕ2(tA)]+[B,Aexp⁡(tA),B,ϕ3(tA)])+O(t3),⋯The recursion is given as(3.30)ci(t)=ci-1(t)+c1(t)∫0t[B,exp⁡(sA)]∫0s1[B,exp⁡(s1A)]∫0s2⋯∫0si-2[B,exp⁡(si-1A)]dsi-1⋯ds1dt.Remark 3.2.
For the novel notation, we have embedded the commutator to the computational scheme. For such a scheme we could save to compute additional the commutators.


## 4. Numerical Examples

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


### 4.1. First Example: Matrix Problem

For another example, consider the matrix equation(4.1)u′(t)=[1230]u,u(0)=u0=(01),
the exact solution is(4.2)u(t)=2(e3t-e-2t)5.
We split the matrix as(4.3)[1230]=[1110]+[0120].The Figure1 presents the numerical errors between the exact and the numerical solution.Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)The Figure2 presents the CPU time of the standard and the iterative splitting schemes.CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)


### 4.2. Second Experiment:10×10 Matrix

We deal in the first with an ODE and separate the complex operator in two simpler operators.We deal with the10×10 ODE system:(4.4)∂tu1=-λ1,1u1+λ2,1u2+⋯+λ10,1u10,(4.5)∂tu2=λ1,2u1-λ2,2(t)u2+⋯+λ10,2u10,(4.6)⋮(4.7)∂tu10=λ1,10u1+λ2,10(t)u2+⋯-λ10,10u10,(4.8)u1(0)=u1,0,…,u10(0)=u10,0(initialconditions),
where λ1(t)∈ℝ+ and λ2(t)∈ℝ+ are the decay factors and u1,0,…,u10,0∈ℝ+. We have the time interval t∈[0,T].We rewrite (4.4) in operator notation, we concentrate us to the following equations:(4.9)∂tu=A(t)u+B(t)u,
where u1(0)=u10=1.0, u2(0)=u20=1.0 are the initial conditions, where we have λ1(t)=t and λ2(t)=t2.The operators are splitted in the following way, while operatorA, see (4.10), covers the upper diagonal entries and operator B, see (4.11), covers the lower diagonal entries of the ODE system (4.4)–(4.7). Such a decomposition allows to reduce the computation of the matrix exponential to a upper or lower matrix and speedup the solver process(4.10)A=(-0.010.010⋯0.01-0.010⋯0.010.01-0.020⋯0.010.010.01-0.030⋯⋮0.010.010.010.010.010.010.010.01-0.0800.010.010.010.010.010.010.010.010-0.08),(4.11)B=(-0.0800.010.010.010.010.010.010.010.010-0.080.010.010.010.010.010.010.010.01⋮0000000-0.020.010.0100000000-0.010.01000000000.01-0.01).The Figure3 present the numerical errors between the exact and the numerical solution. Here we have the best approximation to the exact solution with a sixth-order iterative scheme (6 iterative steps), but no additionally more computational time for finer time step.Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)The Figure4 present the CPU time of the standard and the iterative splitting schemes. Compared to standard splitting schemes, the iterative schemes are cheaper to compute and higher in accuracy. We can also choose smaller time steps to obtain more accurate results and use nearly the same computational resource.CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)Remark 4.1.
The computational results show the benefit of the iterative schemes. We save computational time and achieve higher-order accuracy. The one-side and two-side schemes have the same results.


### 4.3. Third Example: Commutator Problem

We assume to have a large norm of the commutator[A,B] and deal with the matrix equation(4.12)u′(t)=[101110]u,u(0)=u0=(11).In the following, we discuss different splitting ideas.Version 1
We split the matrix as(4.13)[101110]=[9011]+[1109],
while ∥[A,B]∥≥max{∥A∥,∥B∥}.Here we have to deal with highly noncommutative operators and the computational speedup is given in the iterative schemes, while the commutator is not needed to obtain more accurate results, while the standard splitting schemes deal with such commutators for the error reduction.The Figure5 present the numerical errors between the exact and the numerical solution.Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)The Figure6 present the CPU time of the standard and the iterative splitting schemes.CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)Further for the one-side,we obtain more improved results for the following splitting.Version 2
In this version, we have(4.14)[101110]=[9019]+[1101].
The Figure 7 presents the numerical errors between the exact and the numerical solution.Numerical errors of the standard splitting scheme (a) and the iterative schemes based on one-side to operatorA or B (b) and two-side scheme (c).
(a)(b)(c)A more delicate problem is given for the stiff matrices.Version 3
In this version, we have(4.15)[10411104]=[104-101104-1]+[1101].
The Figure 8 present the numerical errors between the exact and the numerical solution.Numerical errors of the one-side splitting scheme withA (a), the one-side splitting scheme with B (b) and the iterative schemes with 1, … , 6 iterative steps (c).
(a)(b)(c)The Figure9 present the CPU time of the standard and the iterative splitting schemes.CPU time of the one-side splitting scheme withA (a) and the iterative schemes with 1, … , 6 iterative steps (b).
(a)(b)Remark 4.2.
The iterative schemes with fast computations of the exponential matrices have a speedup. The constant CPU time of the iterative schemes shows that it benefit instead of the expensive standard schemes. Also for stiff problems with multi iterative steps, we reach the same results of the standardA-B or Strang-Splitting schemes, while we are independent of the commutator estimation.


### 4.4. Two-Phase Example

The next example is a simplified real-life problem for a multiphase transport-reaction equation. We deal with mobile and immobile pores in the porous media, such simulations are given for waste scenarios.We concentrate on the computational benefits of a fast computation of the iterative scheme, given with matrix exponentials.The equation is given as:(4.16)∂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].In the following, we deal with the semidiscretized equation given with the matrices:(4.17)∂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 discretization point (i=1,…,I), the same is also for the other solution vectors.We have the following two operators for the splitting method(4.18)A=DΔx2⋅(-211-21⋱⋱⋱1-211-2)+vΔx⋅(1-11⋱⋱-11-11)∈RI×I,
where I is the number of spatial points.(4.19)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈RI×I,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈RI×I,G=(g00g0⋱⋱⋱0g00g)∈RI×I.We decouple into the following matrices:(4.20)A1=(A0000A0000000000)∈R4I×4I,Ã2=(-Λ1000Λ1-Λ20000-Λ1000Λ1-Λ2)∈R4I×4I,Ã3=(-G0G00-G0GG0-G00G0-G)∈R4I×4I.For the operatorA1 and A2=Ã2+Ã3, we apply the iterative splitting method.Based on the decomposition, operatorA1 is only tridiagonal and operator A2 is block diagonal. Such matrix structure reduce the computation of the exponential operators.The Figure10 present the numerical errors between the exact and the numerical solution. Here we obtain optimal results for one-side iterative schemes on operator B, means we iterate with respect to B and use A as right-hand side.Numerical errors of the one-side splitting scheme withA (a), the one-side splitting scheme with B (b) and the iterative schemes with 1, … , 6 iterative steps (c).
(a)(b)(c)Remark 4.3.
For all iterative schemes, we can reach faster results as for the The iterative schemes with fast computations of the exponential matrices standard schemes. With 4-5 iterative steps we obtain more accurate results as we did for the expensive standard schemes. With one-side iterative schemes we reach the best convergence results.


## 4.1. First Example: Matrix Problem

For another example, consider the matrix equation(4.1)u′(t)=[1230]u,u(0)=u0=(01),
the exact solution is(4.2)u(t)=2(e3t-e-2t)5.
We split the matrix as(4.3)[1230]=[1110]+[0120].The Figure1 presents the numerical errors between the exact and the numerical solution.Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)The Figure2 presents the CPU time of the standard and the iterative splitting schemes.CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)


## 4.2. Second Experiment:10×10 Matrix

We deal in the first with an ODE and separate the complex operator in two simpler operators.We deal with the10×10 ODE system:(4.4)∂tu1=-λ1,1u1+λ2,1u2+⋯+λ10,1u10,(4.5)∂tu2=λ1,2u1-λ2,2(t)u2+⋯+λ10,2u10,(4.6)⋮(4.7)∂tu10=λ1,10u1+λ2,10(t)u2+⋯-λ10,10u10,(4.8)u1(0)=u1,0,…,u10(0)=u10,0(initialconditions),
where λ1(t)∈ℝ+ and λ2(t)∈ℝ+ are the decay factors and u1,0,…,u10,0∈ℝ+. We have the time interval t∈[0,T].We rewrite (4.4) in operator notation, we concentrate us to the following equations:(4.9)∂tu=A(t)u+B(t)u,
where u1(0)=u10=1.0, u2(0)=u20=1.0 are the initial conditions, where we have λ1(t)=t and λ2(t)=t2.The operators are splitted in the following way, while operatorA, see (4.10), covers the upper diagonal entries and operator B, see (4.11), covers the lower diagonal entries of the ODE system (4.4)–(4.7). Such a decomposition allows to reduce the computation of the matrix exponential to a upper or lower matrix and speedup the solver process(4.10)A=(-0.010.010⋯0.01-0.010⋯0.010.01-0.020⋯0.010.010.01-0.030⋯⋮0.010.010.010.010.010.010.010.01-0.0800.010.010.010.010.010.010.010.010-0.08),(4.11)B=(-0.0800.010.010.010.010.010.010.010.010-0.080.010.010.010.010.010.010.010.01⋮0000000-0.020.010.0100000000-0.010.01000000000.01-0.01).The Figure3 present the numerical errors between the exact and the numerical solution. Here we have the best approximation to the exact solution with a sixth-order iterative scheme (6 iterative steps), but no additionally more computational time for finer time step.Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)The Figure4 present the CPU time of the standard and the iterative splitting schemes. Compared to standard splitting schemes, the iterative schemes are cheaper to compute and higher in accuracy. We can also choose smaller time steps to obtain more accurate results and use nearly the same computational resource.CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)Remark 4.1.
The computational results show the benefit of the iterative schemes. We save computational time and achieve higher-order accuracy. The one-side and two-side schemes have the same results.


## 4.3. Third Example: Commutator Problem

We assume to have a large norm of the commutator[A,B] and deal with the matrix equation(4.12)u′(t)=[101110]u,u(0)=u0=(11).In the following, we discuss different splitting ideas.Version 1
We split the matrix as(4.13)[101110]=[9011]+[1109],
while ∥[A,B]∥≥max{∥A∥,∥B∥}.Here we have to deal with highly noncommutative operators and the computational speedup is given in the iterative schemes, while the commutator is not needed to obtain more accurate results, while the standard splitting schemes deal with such commutators for the error reduction.The Figure5 present the numerical errors between the exact and the numerical solution.Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)The Figure6 present the CPU time of the standard and the iterative splitting schemes.CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1, … , 6 iterative steps.
(a)(b)Further for the one-side,we obtain more improved results for the following splitting.Version 2
In this version, we have(4.14)[101110]=[9019]+[1101].
The Figure 7 presents the numerical errors between the exact and the numerical solution.Numerical errors of the standard splitting scheme (a) and the iterative schemes based on one-side to operatorA or B (b) and two-side scheme (c).
(a)(b)(c)A more delicate problem is given for the stiff matrices.Version 3
In this version, we have(4.15)[10411104]=[104-101104-1]+[1101].
The Figure 8 present the numerical errors between the exact and the numerical solution.Numerical errors of the one-side splitting scheme withA (a), the one-side splitting scheme with B (b) and the iterative schemes with 1, … , 6 iterative steps (c).
(a)(b)(c)The Figure9 present the CPU time of the standard and the iterative splitting schemes.CPU time of the one-side splitting scheme withA (a) and the iterative schemes with 1, … , 6 iterative steps (b).
(a)(b)Remark 4.2.
The iterative schemes with fast computations of the exponential matrices have a speedup. The constant CPU time of the iterative schemes shows that it benefit instead of the expensive standard schemes. Also for stiff problems with multi iterative steps, we reach the same results of the standardA-B or Strang-Splitting schemes, while we are independent of the commutator estimation.


## 4.4. Two-Phase Example

The next example is a simplified real-life problem for a multiphase transport-reaction equation. We deal with mobile and immobile pores in the porous media, such simulations are given for waste scenarios.We concentrate on the computational benefits of a fast computation of the iterative scheme, given with matrix exponentials.The equation is given as:(4.16)∂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].In the following, we deal with the semidiscretized equation given with the matrices:(4.17)∂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 discretization point (i=1,…,I), the same is also for the other solution vectors.We have the following two operators for the splitting method(4.18)A=DΔx2⋅(-211-21⋱⋱⋱1-211-2)+vΔx⋅(1-11⋱⋱-11-11)∈RI×I,
where I is the number of spatial points.(4.19)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈RI×I,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈RI×I,G=(g00g0⋱⋱⋱0g00g)∈RI×I.We decouple into the following matrices:(4.20)A1=(A0000A0000000000)∈R4I×4I,Ã2=(-Λ1000Λ1-Λ20000-Λ1000Λ1-Λ2)∈R4I×4I,Ã3=(-G0G00-G0GG0-G00G0-G)∈R4I×4I.For the operatorA1 and A2=Ã2+Ã3, we apply the iterative splitting method.Based on the decomposition, operatorA1 is only tridiagonal and operator A2 is block diagonal. Such matrix structure reduce the computation of the exponential operators.The Figure10 present the numerical errors between the exact and the numerical solution. Here we obtain optimal results for one-side iterative schemes on operator B, means we iterate with respect to B and use A as right-hand side.Numerical errors of the one-side splitting scheme withA (a), the one-side splitting scheme with B (b) and the iterative schemes with 1, … , 6 iterative steps (c).
(a)(b)(c)Remark 4.3.
For all iterative schemes, we can reach faster results as for the The iterative schemes with fast computations of the exponential matrices standard schemes. With 4-5 iterative steps we obtain more accurate results as we did for the expensive standard schemes. With one-side iterative schemes we reach the best convergence results.


## 5. Conclusions and Discussions

In this work, we have presented a very economical and practical method by successive approximations. Here the idea to decouple the expensive computation of matrix exponential via iterative splitting schemes has the benefit of less computational time. While the error analysis present stable methods and higher-order schemes, the applications show the speedup with the iterative schemes. In, future, we concentrate on matrix dependent scheme, based on iterative splitting algorithms.



---

*Source: 193781-2011-03-29.xml*