Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volumeย 2011ย (2011), Article IDย 193781, 27 pages
Research Article

Computing Exponential for Iterative Splitting Methods: Algorithms and Applications

Humboldt-Univeristรคt zu Berlin, 10099 Berlin, Germany

Received 29 November 2010; Accepted 27 January 2011

Academic Editor: Shuyuย Sun

Copyright ยฉ 2011 Jรผrgen Geiser. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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.

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 ๐ฏ 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:๐œ™๐œ•๐‘ก๐‘๐‘–+โˆ‡โ‹…๐…๐‘–๎€ท=๐‘”โˆ’๐‘๐‘–+๐‘๐‘–,im๎€ธ+๐‘˜๐›ผ๎€ทโˆ’๐‘๐‘–+๐‘๐‘–,ad๎€ธโˆ’๐œ†๐‘–,๐‘–๐œ™๐‘๐‘–+๎“๐‘˜=๐‘˜(๐‘–)๐œ†๐‘–,๐‘˜๐œ™๐‘๐‘˜+๎‚๐‘„๐‘–[],๐…,inฮฉร—0,๐‘ก(1.1)๐‘–=๐ฏ๐‘๐‘–โˆ’๐ท๐‘’(๐‘–)โˆ‡๐‘๐‘–+๐œ‚๐„๐‘๐‘–,(1.2)๐œ™๐œ•๐‘ก๐‘๐‘–,im๎€ท๐‘=๐‘”๐‘–โˆ’๐‘๐‘–,im๎€ธ+๐‘˜๐›ผ๎€ท๐‘๐‘–,im,adโˆ’๐‘๐‘–,im๎€ธโˆ’๐œ†๐‘–,๐‘–๐œ™๐‘๐‘–,im+๎“๐‘˜=๐‘˜(๐‘–)๐œ†๐‘–,๐‘˜๐œ™๐‘๐‘˜,im+๎‚๐‘„๐‘–,im[],,inฮฉร—0,๐‘ก(1.3)๐œ™๐œ•๐‘ก๐‘๐‘–,ad=๐‘˜๐›ผ๎€ท๐‘๐‘–โˆ’๐‘๐‘–,ad๎€ธโˆ’๐œ†๐‘–,๐‘–๐œ™๐‘๐‘–,ad+๎“๐‘˜=๐‘˜(๐‘–)๐œ†๐‘–,๐‘˜๐œ™๐‘๐‘˜,ad+๎‚๐‘„๐‘–,ad[],inฮฉร—0,๐‘ก,(1.4)๐œ™๐œ•๐‘ก๐‘๐‘–,im,ad=๐‘˜๐›ผ๎€ท๐‘๐‘–,imโˆ’๐‘๐‘–,im,ad๎€ธโˆ’๐œ†๐‘–,๐‘–๐œ™๐‘๐‘–,im,ad+๎“๐‘˜=๐‘˜(๐‘–)๐œ†๐‘–,๐‘˜๐œ™๐‘๐‘˜,im,ad+๎‚๐‘„๐‘–,im,ad[],๐‘,inฮฉร—0,๐‘ก(1.5)๐‘–(๐ฑ,๐‘ก)=๐‘๐‘–,0(๐ฑ),๐‘๐‘–,ad(๐ฑ,๐‘ก)=0,๐‘๐‘–,im(๐ฑ,๐‘ก)=0,๐‘๐‘–,im,ad๐‘(๐ฑ,๐‘ก)=0,onฮฉ,(1.6)๐‘–(๐ฑ,๐‘ก)=๐‘๐‘–,1(๐ฑ,๐‘ก),๐‘๐‘–,ad(๐ฑ,๐‘ก)=0,๐‘๐‘–,im(๐ฑ,๐‘ก)=0,๐‘๐‘–,im,ad[],(๐ฑ,๐‘ก)=0,on๐œ•ฮฉร—0,๐‘ก(1.7) where the initial value is given as ๐‘๐‘–,0 and we assume Dirichlet boundary conditions with the function ๐‘๐‘–,1(๐‘ฅ,๐‘ก) sufficiently smooth, all other initial and boundary conditions of the other phases zero.๐œ™:effective porosity (โˆ’),๐‘๐‘–:concentration of the ๐‘–th gaseous species in the plasma chamber phase (mol/cm3),๐‘๐‘–,im:concentration of the ๐‘–th gaseous species in the immobile zones of the plasma chamber phase (mol/cm3),๐‘๐‘–,ad:concentration of the ๐‘–th gaseous species in the mobile adsorbed zones of the plasma chamber phase (mol/cm3),๐‘๐‘–,im,ad:concentration of the ๐‘–th gaseous species in the immobile adsorbed zones of the plasma chamber phase (mol/cm3),๐ฏ:velocity through the chamber and porous substrate [4] (cm/nsec),๐ท๐‘’(๐‘–):element-specific diffusions-dispersions tensor (cm2/nsec),๐œ†๐‘–,๐‘–:decay constant of the ๐‘–th species (1/nsec),๎‚๐‘„๐‘–:source term of the ๐‘–th species (mol/(cm3nsec)),๐‘”:exchange rate between the mobile and immobile concentration (1/nsec),๐‘˜๐›ผ:exchange rate between the mobile and adsorbed concentration or immobile and immobile adsorbed concentration (kinetic controlled sorption) (1/nsec),๐„:stationary electric field in the apparatus (V/cm),๐œ‚:the mobility rate in the electric field, see [5] (cm2/nsec).

with ๐‘–=1,โ€ฆ,๐‘€ and ๐‘€ 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 ๐ฏ, that presents the flow-direction and the absolute value of the plasma flux. The velocity field is divergence-free. The decay constant of the ๐‘–th species is denoted by ๐œ†๐‘–. Thereby does ๐‘˜(๐‘–) denote the indices of the other species.

In this paper we concentrate on solving linear evolution equations, such as the differential equation,๐œ•๐‘ก๐‘ข=๐ด๐‘ข,๐‘ข(0)=๐‘ข0,(1.8) where ๐ด 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 Section 2, 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 [๐‘ก๐‘›,๐‘ก๐‘›+1] we solve the following subproblems consecutively for ๐‘–=0,2,โ€ฆ,2๐‘š. (cf. [6, 7]):๐œ•๐‘๐‘–(๐‘ก)๐œ•๐‘ก=๐ด๐‘๐‘–(๐‘ก)+๐ต๐‘๐‘–โˆ’1(๐‘ก),with๐‘๐‘–(๐‘ก๐‘›)=๐‘๐‘›,๐‘0(๐‘ก๐‘›)=๐‘๐‘›,๐‘โˆ’1=0.0,๐œ•๐‘๐‘–+1(๐‘ก)๐œ•๐‘ก=๐ด๐‘๐‘–(๐‘ก)+๐ต๐‘๐‘–+1(๐‘ก),with๐‘๐‘–+1(๐‘ก๐‘›)=๐‘๐‘›,(2.1) where ๐‘๐‘› is the known split approximation at the time level ๐‘ก=๐‘ก๐‘›. The split approximation at the time level ๐‘ก=๐‘ก๐‘›+1 is defined as ๐‘๐‘›+1=๐‘2๐‘š+1(๐‘ก๐‘›+1), (clearly, the function ๐‘๐‘–+1(๐‘ก) depends on the interval [๐‘ก๐‘›,๐‘ก๐‘›+1], too, but, for the sake of simplicity, in our notation we omit the dependence on ๐‘›).

In the following we will analyze the convergence and the rate of convergence of the method (2.1) for ๐‘š tends to infinity for the linear operators ๐ด,๐ตโˆถ๐—โ†’๐—, where we assume that these operators and their sum are generators of the ๐ถ0 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 a ๐ถ0 semigroup based on their underlying graph norms.

Theorem 2.1. Let one consider the abstract Cauchy problem in a Hilbert space ๐—๐œ•๐‘ก๐‘๐‘(๐‘ฅ,๐‘ก)=๐ด๐‘(๐‘ฅ,๐‘ก)+๐ต๐‘(๐‘ฅ,๐‘ก),0<๐‘กโ‰ค๐‘‡,๐‘ฅโˆˆฮฉ,(๐‘ฅ,0)=๐‘0(๐‘ฅ),๐‘ฅโˆˆฮฉ,๐‘(๐‘ฅ,๐‘ก)=๐‘1[],(๐‘ฅ,๐‘ก),๐‘ฅโˆˆ๐œ•ฮฉร—0,๐‘‡(2.2) where ๐ด,๐ตโˆถ๐ท(๐—)โ†’๐— are given linear operators which are generators of the ๐ถ0-semigroup and ๐‘0โˆˆ๐— is a given element. We assume ๐ด, ๐ต are unbounded. Further, we assume the estimations of the unbounded operator ๐ต with sufficient smooth initial conditions (see [8]): โ€–โ€–๐ตexp((๐ด+๐ต)๐œ)๐‘ข0โ€–โ€–โ‰ค๐œ….(2.3)
Further we assume the estimation of the partial integration of the unbounded operator ๐ต (see [8]): โ€–โ€–โ€–๐ต๎€œ๐œ0โ€–โ€–โ€–exp(๐ต๐‘ )๐‘ ๐‘‘๐‘ โ‰ค๐œ๐ถ.(2.4)
Then, we can bound our iterative operator splitting method as โ€–โ€–๎€ท๐‘†๐‘–๎€ธโ€–โ€–โˆ’exp((๐ด+๐ต)๐œ)โ‰ค๐ถ๐œ๐‘–,(2.5) where ๐‘†๐‘– is the approximated solution for the ๐‘–th iterative step and ๐ถ is a constant that can be chosen uniformly on bounded time intervals.

Proof. Let us consider the iteration (2.1) on the subinterval [๐‘ก๐‘›,๐‘ก๐‘›+1].
For the first iterations, we have๐œ•๐‘ก๐‘1(๐‘ก)=๐ด๐‘1๎€ท๐‘ก(๐‘ก),๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป,(2.6) and for the second iteration we have: ๐œ•๐‘ก๐‘2(๐‘ก)=๐ด๐‘1(๐‘ก)+๐ต๐‘2๎€ท๐‘ก(๐‘ก),๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป.(2.7)
In general, we have:
(i)for the odd iterations, ๐‘–=2๐‘š+1 for ๐‘š=0,1,2,โ€ฆ๐œ•๐‘ก๐‘๐‘–(๐‘ก)=๐ด๐‘๐‘–(๐‘ก)+๐ต๐‘๐‘–โˆ’1๎€ท๐‘ก(๐‘ก),๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป,(2.8) where for ๐‘0(๐‘ก)โ‰ก0,(ii)for the even iterations, ๐‘–=2๐‘š for ๐‘š=1,2,โ€ฆ๐œ•๐‘ก๐‘๐‘–(๐‘ก)=๐ด๐‘๐‘–โˆ’1(๐‘ก)+๐ต๐‘๐‘–๎€ท๐‘ก(๐‘ก),๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป.(2.9)
We have the following solutions for the iterative scheme: the solutions for the first two equations are given by the variation of constants:๐‘1๎€ท๐ด๎€ท๐‘ก(๐‘ก)=exp๐‘›+1โˆ’๐‘ก๎€ธ๎€ธ๐‘(๐‘ก๐‘›๎€ท๐‘ก),๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป,๐‘2(๐‘ก)=exp(๐ต๐‘ก)๐‘(๐‘ก๐‘›๎€œ)+๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ต๎€ท๐‘กexp๐‘›+1โˆ’๐‘ ๎€ธ๎€ธ๐ด๐‘1๎€ท๐‘ก(๐‘ )๐‘‘๐‘ ,๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป.(2.10)
For the recurrence relations with even and odd iterations, we have the solutions: for the odd iterations: ๐‘–=2๐‘š+1, for ๐‘š=0,1,2,โ€ฆ,๐‘๐‘–(๐‘ก)=exp(๐ด(๐‘กโˆ’๐‘ก๐‘›))๐‘(๐‘ก๐‘›)+๎€œ๐‘ก๐‘ก๐‘›๐‘กexp๎€ท๎€ท๐‘›+1๎€ธ๐ด๎€ธโˆ’๐‘ ๐ต๐‘๐‘–โˆ’1(๎€ท๐‘ก๐‘ )๐‘‘๐‘ ,๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป.(2.11)
For the even iterations, ๐‘–=2๐‘š, for ๐‘š=1,2,โ€ฆ๐‘๐‘–(๐‘ก)=exp(๐ต(๐‘กโˆ’๐‘ก๐‘›))๐‘(๐‘ก๐‘›)+๎€œ๐‘ก๐‘ก๐‘›๐‘กexp๎€ท๎€ท๐‘›+1๎€ธ๐ต๎€ธโˆ’๐‘ ๐ด๐‘๐‘–โˆ’1(๎€ท๐‘ก๐‘ )๐‘‘๐‘ ,๐‘กโˆˆ๐‘›,๐‘ก๐‘›+1๎€ป.(2.12)
The consistency is given as the following.
For ๐‘’1, we have ๐‘1(๐œ)=exp(๐ด๐œ)๐‘(๐‘ก๐‘›),๐‘(๐œ)=exp((๐ด+๐ต)๐œ)๐‘(๐‘ก๐‘›)=exp(๐ด๐œ)๐‘(๐‘ก๐‘›)+๎€œ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ ๎€ธ๎€ธ๐ตexp(๐‘ (๐ด+๐ต))๐‘(๐‘ก๐‘›)๐‘‘๐‘ .(2.13)
We obtainโ€–โ€–๐‘’1โ€–โ€–=โ€–โ€–๐‘โˆ’๐‘1โ€–โ€–โ‰คโ€–exp((๐ด+๐ต)๐œ)๐‘(๐‘ก๐‘›)โˆ’exp(๐ด๐œ)๐‘(๐‘ก๐‘›)โ€–โ‰ค๐ถ1๐œ๐‘(๐‘ก๐‘›).(2.14)For ๐‘’2 we have (alternating) ๐‘2(๐œ)=exp(๐ต๐œ)๐‘(๐‘ก๐‘›๎€œ)+๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ต๎€ท๐‘กexp๐‘›+1โˆ’๐‘ ๎€ธ๎€ธ๐ดexp(๐‘ ๐ด)๐‘(๐‘ก๐‘›)๐‘‘๐‘ ,๐‘(๐œ)=exp(๐ต๐œ)๐‘(๐‘ก๐‘›๎€œ)+๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ต๎€ท๐‘กexp๐‘›+1โˆ’๐‘ ๎€ธ๎€ธ๐ดexp(๐‘ ๐ด)๐‘(๐‘ก๐‘›+๎€œ)๐‘‘๐‘ ๐‘ก๐‘›+1๐‘ก๐‘›๎€œexp(๐ต๐‘ )๐ด๐‘ก๐‘›+1๐‘กโˆ’๐‘ ๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ โˆ’๐œŒ๎€ธ๎€ธ๐ตexp(๐œŒ(๐ด+๐ต))๐‘(๐‘ก๐‘›)๐‘‘๐œŒ๐‘‘๐‘ .(2.15)
We obtainโ€–โ€–๐‘’2โ€–โ€–โ‰คโ€–โ€–exp((๐ด+๐ต)๐œ)๐‘(๐‘ก๐‘›)โˆ’๐‘2โ€–โ€–โ‰ค๐ถ2๐œ2๐‘(๐‘ก๐‘›).(2.16)
For odd and even iterations, the recursive proof is given in the following. For the odd iterations (only iterations over ๐ด), ๐‘–=2๐‘š+1 for ๐‘š=0,1,2,โ€ฆ, for ๐‘’๐‘–, we have:๐‘๐‘–(๐œ)=exp(๐ด๐œ)๐‘(๐‘ก๐‘›)+๎€œ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ ๎€ธ๎€ธ๐ตexp(๐‘ ๐ด)๐‘(๐‘ก๐‘›+๎€œ)๐‘‘๐‘ ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ 1๐ต๎€œ๎€ธ๎€ธ๐‘ก๐‘›+1โˆ’๐‘ 1๐‘ก๐‘›๐‘กexp๎€ท๎€ท๐‘›+1โˆ’๐‘ 1โˆ’๐‘ 2๎€ธ๐ด๎€ธ๎€ท๐‘ ๐ตexp2๐ด๎€ธ๐‘(๐‘ก๐‘›)๐‘‘๐‘ 2๐‘‘๐‘ 1+๎€œ+โ‹ฏ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ 1๐ต๎€œ๎€ธ๎€ธ๐‘ก๐‘›+1โˆ’๐‘ 1๐‘ก๐‘›๐‘กexp๎€ท๎€ท๐‘›+1โˆ’๐‘ 1โˆ’๐‘ 2๎€ธ๐ด๎€ธ๎€ท๐‘ ๐ตexp2๐ด๎€ธ๐‘(๐‘ก๐‘›)๐‘‘๐‘ 2๐‘‘๐‘ 1+๎€œ+โ‹ฏ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ 1๎€œ๎€ธ๎€ธ๐ตโ‹ฏ๐‘ก๐‘›+1โˆ’โˆ‘๐‘–โˆ’2๐‘—=1๐‘ ๐‘—๐‘ก๐‘›๐‘กexp๎ƒฉ๎ƒฉ๐‘›+1โˆ’๐‘–โˆ’1๎“๐‘—=1๐‘ ๐‘—๎ƒช๐ด๎ƒช๐ต๎€ท๐‘ ร—exp๐‘–๐ด๎€ธ๐‘(๐‘ก๐‘›)๐‘‘๐‘ ๐‘–โˆ’1โ‹ฏ๐‘‘๐‘ 2๐‘‘๐‘ 1,๐‘(๐œ)=exp(๐ด๐œ)๐‘(๐‘ก๐‘›)+๎€œ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ ๎€ธ๎€ธ๐ตexp(๐‘ ๐ด)๐‘(๐‘ก๐‘›)+๎€œ๐‘‘๐‘ ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ 1๐ต๎€œ๎€ธ๎€ธ๐‘ก๐‘›+1โˆ’๐‘ 1๐‘ก๐‘›๐‘กexp๎€ท๎€ท๐‘›+1โˆ’๐‘ 1โˆ’๐‘ 2๎€ธ๐ด๎€ธ๎€ท๐‘ ๐ตexp2๐ด๎€ธ๐‘(๐‘ก๐‘›)๐‘‘๐‘ 2๐‘‘๐‘ 1+๎€œ+โ‹ฏ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ 1๐ต๎€œ๎€ธ๎€ธ๐‘ก๐‘›+1โˆ’๐‘ 1๐‘ก๐‘›๐‘กexp๎€ท๎€ท๐‘›+1โˆ’๐‘ 1โˆ’๐‘ 2๎€ธ๐ด๎€ธ๎€ท๐‘ ๐ตexp2๐ด๎€ธ๐‘(๐‘ก๐‘›)๐‘‘๐‘ 2๐‘‘๐‘ 1+๎€œ+โ‹ฏ๐‘ก๐‘›+1๐‘ก๐‘›๎€ท๐ด๎€ท๐‘กexp๐‘›+1โˆ’๐‘ 1๎€œ๎€ธ๎€ธ๐ตโ‹ฏ๐‘ก๐‘›+1โˆ’โˆ‘๐‘–โˆ’2๐‘—=1๐‘ ๐‘—๐‘ก๐‘›๐‘กexp๎ƒฉ๎ƒฉ๐‘›+1โˆ’๐‘–๎“๐‘—=1๐‘ ๐‘—๎ƒช๐ด๎ƒช๐ต๎€ท๐‘ ร—exp๐‘–๐ด๎€ธ๐‘(๐‘ก๐‘›)๐‘‘๐‘ ๐‘–โ‹ฏ๐‘‘๐‘ 2๐‘‘๐‘ 1.(2.17)
We obtainโ€–โ€–๐‘’๐‘–โ€–โ€–โ‰คโ€–โ€–exp((๐ด+๐ต)๐œ)๐‘(๐‘ก๐‘›)โˆ’๐‘๐‘–โ€–โ€–โ‰ค๐ถ๐œ๐‘–๐‘(๐‘ก๐‘›),(2.18) 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.2. The same idea can be applied to ๐ด=โˆ‡๐ทโˆ‡๐ต=โˆ’๐ฏโ‹…โˆ‡, so that one operator is less unbounded, but we reduce the convergence order โ€–โ€–๐‘’1โ€–โ€–=๐พโ€–๐ตโ€–๐œ๐›ผ1โ€–โ€–๐‘’0โ€–โ€–๎€ท๐œ+๐’ช1+๐›ผ1๎€ธ,โ€–โ€–๐‘’andhence,2โ€–โ€–โ€–โ€–๐‘’=๐พโ€–๐ตโ€–0โ€–โ€–๐œ1+๐›ผ1+๐›ผ2๎€ท๐œ+๐’ช1+๐›ผ1+๐›ผ๎€ธ,(2.19) where 0โ‰ค๐›ผ1, ๐›ผ2<1.

Remark 2.3. If we assume the consistency of ๐’ช(๐œ๐‘š) for the initial value ๐‘’1(๐‘ก๐‘›) and ๐‘’2(๐‘ก๐‘›), we can redo the proof and obtain at least a global error of the splitting methods of ๐’ช(๐œ๐‘šโˆ’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:๐œ•๐‘ก๐œ=๐ด1๐œ+๐ด2๐œ+๐ต1๎€ท๐œโˆ’๐œim๎€ธ+๐ต2๎€ท๐œโˆ’๐œad๎€ธ๐œ•+๐,๐‘ก๐œim=๐ด2๐œim+๐ต1๎€ท๐œim๎€ธโˆ’๐œ+๐ต2๎€ท๐œimโˆ’๐œim,ad๎€ธ+๐im,๐œ•๐‘ก๐œad=๐ด2๐œad+๐ต2๎€ท๐œad๎€ธโˆ’๐œ+๐ad,๐œ•๐‘ก๐œim,ad=๐ด2๐œim,ad+๐ต2๎€ท๐œim,adโˆ’๐œim๎€ธ+๐im,ad,(2.20) where ๐œ=(๐‘1,โ€ฆ,๐‘๐‘š)๐‘‡ is the spatial discretized concentration in the mobile phase, see (1.1), ๐œim=(๐‘1,im,โ€ฆ,๐‘๐‘š,im)๐‘‡ is the concentration in the immobile phase, the some also for the other phase concentrations. ๐ด1 is the stiffness matrix of (1.1), ๐ด2 is the reaction matrix of the right-hand side of (1.1), ๐ต1 and ๐ต2 are diagonal matrices with the exchange of the immobile and kinetic parameters, see (1.4) and (1.5).

Further, ๐,โ€ฆ,๐im,ad are the spatial discretized sources vectors.

Now, we have the following ordinary differential equation:๐œ•๐‘กโŽ›โŽœโŽœโŽœโŽ๐ด๐‚=1+๐ด2+๐ต1+๐ต2โˆ’๐ต1โˆ’๐ต20โˆ’๐ต1๐ด2+๐ต1+๐ต20โˆ’๐ต2โˆ’๐ต20๐ด2+๐ต200โˆ’๐ต20๐ด2+๐ต2โŽžโŽŸโŽŸโŽŸโŽ ๎‚๐‚+๐,(2.21) where ๐‚=(๐œ,๐œim,๐œad,๐œim,ad)๐‘‡ and the right-hand side is given as ๎‚๐=(๐,๐im,๐ad,๐im,ad)๐‘‡.

For such an equation, we apply the decomposition of the matrices:๐œ•๐‘ก๎‚๎‚๐œ•๐‚=๐ด๐‚+๐,๐‘ก๎‚๐ด๐‚=1๎‚๐ด๐‚+2๎‚๐‚+๐,(2.22) where๎‚๐ด1=โŽ›โŽœโŽœโŽœโŽ๐ด1+๐ด20000๐ด20000๐ด20000๐ด2โŽžโŽŸโŽŸโŽŸโŽ ,๎‚๐ด2=โŽ›โŽœโŽœโŽœโŽ๐ต1+๐ต2โˆ’๐ต1โˆ’๐ต20โˆ’๐ต1๐ต1+๐ต20โˆ’๐ต2โˆ’๐ต20๐ต200โˆ’๐ต20๐ต2โŽžโŽŸโŽŸโŽŸโŽ .(2.23)

The equation system is numerically solved by an iterative scheme.

Algorithm 2.4. We divide our time interval [0,๐‘‡] into subintervals [๐‘ก๐‘›,๐‘ก๐‘›+1], where ๐‘›=0,1,โ€ฆ๐‘, ๐‘ก0=0, and ๐‘ก๐‘=๐‘‡.
We start with ๐‘›=0.
(1)The initial conditions are given with ๐‚0(๐‘ก๐‘›+1)=๐‚(๐‘ก๐‘›). We start with ๐‘˜=0.(2)Compute the fix-point iteration scheme given as ๐œ•๐‘ก๐‚๐‘˜=๎‚๐ด1๐‚๐‘˜+๎‚๐ด2๐‚๐‘˜โˆ’1+๎‚๐,(2.24) where ๐‘˜ 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 [๐‘ก๐‘›,๐‘ก๐‘›+1] is given as โ€–โ€–๐‚๐‘˜๎€ท๐‘ก๐‘›+1๎€ธโˆ’๐‚๐‘˜โˆ’1๎€ท๐‘ก๐‘›+1๎€ธโ€–โ€–โ‰คerr,(2.25) 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 ๐‚๎€ท๐‘ก๐‘›+1๎€ธ=๐‚๐‘˜๎€ท๐‘ก๐‘›+1๎€ธ.(2.26)If ๐‘›=๐‘ then we stop and are done.If (2.25) is not fulfilled, we do ๐‘˜=๐‘˜+1 and go to (2).

The error analysis of the schemes are given in the following Theorem.

Theorem 2.5. Let ๐ด,๐ตโˆˆโ„’(๐—) be given linear bounded operators in a Banach space โ„’(๐—). We consider the abstract Cauchy problem: ๐œ•๐‘ก๎‚๎‚๐‚(๐‘ก)=๐ด๐‚(๐‘ก)+๐ต๐‚(๐‘ก),๐‘ก๐‘›โ‰ค๐‘กโ‰ค๐‘ก๐‘›+1,๐‚๎€ท๐‘ก๐‘›๎€ธ=๐‚๐‘›,for๐‘›=1,โ€ฆ,๐‘,(2.27) where ๐‘ก1=0 and the final time is ๐‘ก๐‘=๐‘‡โˆˆโ„+. Then problem (2.27) has a unique solution. For a finite steps with time size ๐œ๐‘›=๐‘ก๐‘›+1โˆ’๐‘ก๐‘›, the iteration (2.24) for ๐‘˜=1,2,โ€ฆ,๐‘ž is consistent with an order of consistency ๐’ช(๐œ๐‘ž๐‘›).

Proof. The outline of the proof is given in [2].

In the next section we describe the computation of the integral formulation with exp-functions.

3. Exponentials to Compute Iterative Schemes

The theoretical ideas can be discussed in the following formulation:๐ท๐ด(๐ต,๐‘ก)[0]๐ท=exp(๐‘ก๐ต),๐ด(๐ต,๐‘ก)[๐‘˜]๎€œ=๐‘˜๐‘ก0exp((๐‘กโˆ’๐‘ )๐ต)๐ด๐ท๐ด(๐ต,๐‘ )[๐‘˜โˆ’1]๐‘‘๐‘ ,(3.1) and the matrix formulation of our two-step scheme is given as๎‚โŽ›โŽœโŽœโŽœโŽœโŽœโŽโŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ ๐ด=๐ด0โ‹ฏโ‹ฏ๐ด๐ต0โ‹ฏ0๐ต๐ดโ‹ฏโ‹ฎโ‹ฑโ‹ฑโ‹ฑ0โ‹ฏ๐ต๐ด,(3.2) the computation of the exp-Matrix can be expressed as:๎‚€๎‚๎‚โŽ›โŽœโŽœโŽœโŽœโŽœโŽœโŽœโŽ๐ทexp๐ด๐‘กโˆถ=exp(๐ด๐‘ก)0โ‹ฏโ‹ฏ๐ดโ‹ฏ๐ท(๐ต)exp(๐ต๐‘ก)0๐ต(๐ด)[2]๐ท2!๐ต(๐ด)[1]๐ท1!exp(๐ด๐‘ก)โ‹ฏโ‹ฎโ‹ฑโ‹ฑโ‹ฑ๐ต(๐ด)[๐‘›]โ‹ฏ๐ท๐‘›!๐ต(๐ด)[1]โŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽŸโŽŸโŽ 1!exp(๐ด๐‘ก).(3.3)

Here, we have to compute the right-hand side as time dependent term, means we evaluate exp(๐ด๐‘ก) and exp(๐ต๐‘ก) 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 assume ๐ท๐ด(๐ต,๐‘ )[๐‘˜โˆ’1]โ‰ˆ๐ท๐ด(๐ต,0)[๐‘˜โˆ’01]).

Later we also follow with more extended schemes.

Computing ๐œ™-Functions:
So the matrix formulation of our scheme is given as๎‚€๎‚๎‚๐‘Œ(๐‘ก)=exp๐ด๐‘ก๐‘Œ(0),(3.5) where ๎‚๐ด is given as, ๎‚โŽ›โŽœโŽœโŽœโŽœโŽœโŽโŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ ๐ด=๐ด0โ‹ฏโ‹ฏ๐ด๐ต0โ‹ฏ0๐ด๐ตโ‹ฏโ‹ฎโ‹ฑโ‹ฑโ‹ฑ0โ‹ฏ๐ต๐ด,(3.6) the computation of the exp-Matrix can be expressed in a first-order scheme and with the assumption of commutations is given as: ๎‚€๎‚๎‚โŽ›โŽœโŽœโŽœโŽœโŽœโŽ๐œ™exp๐ด๐‘กโˆถ=exp(๐ด๐‘ก)0โ‹ฏโ‹ฏโ‹ฏ1๐œ™(๐ต๐‘ก)๐ด๐‘กexp(๐ต๐‘ก)0โ‹ฏโ‹ฏ2(๐ด๐‘ก)๐ต(๐ต+๐ด)๐‘ก2๐œ™1๐œ™(๐ด๐‘ก)๐ต๐‘กexp(๐ด๐‘ก)โ‹ฏโ‹ฏโ‹ฎโ‹ฑโ‹ฑโ‹ฑโ‹ฎ๐‘–(๐ด๐‘ก)๐ต(๐ด+๐ต)๐‘–โˆ’1๐‘ก๐‘–โ‹ฏ๐œ™1โŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ (๐ด๐‘ก)๐ต๐‘กexp(๐ด๐‘ก),(3.7) where we assume ๐‘ข๐‘–(๐‘ )=๐‘ข(๐‘ก๐‘›), mean we approximate the right-hand side term.

For higher-orders we should also include the full derivations of ๐‘1,๐‘2,โ€ฆ.

3.1. Derivation of the Formula

Consider the equationฬ‡๐‘ฅ=๐ด๐‘ฅ+๐‘Ž,๐‘ฅ(0)=0,(3.8) the solution of this equation in terms of the ๐œ™๐‘˜,๐‘˜=0,1, is given as๐‘ฅ=๐‘ก๐œ™1(๐‘ก๐ด)๐‘Ž.(3.9) Similarly, the equationฬ‡๐‘ฅ=๐ด๐‘ฅ+๐‘๐‘ก,๐‘ฅ(0)=0,(3.10) has a solution of the form๐‘ฅ=๐‘ก2๐œ™2(๐‘ก๐ด)๐‘.(3.11) In general, the solution of the equation๐‘กฬ‡๐‘ฅ=๐ด๐‘ฅ+๐‘Ž+๐‘๐‘ก+๐‘22!+โ‹ฏ,๐‘ฅ(0)=0,(3.12) admits the form๐‘ฅ=๐‘ก๐œ™1(๐‘ก๐ด)๐‘Ž+๐‘ก2๐œ™2(๐‘ก๐ด)๐‘+๐‘ก3๐œ™3(๐‘ก๐ด)๐‘+โ‹ฏ.(3.13) In this section, we use the formula (3.21) for iterative schemes given asฬ‡๐‘ข1=๐ด๐‘ข1,ฬ‡๐‘ข2=๐ด๐‘ข1+๐ต๐‘ข2,ฬ‡๐‘ข3=๐ด๐‘ข3+๐ต๐‘ข2,ฬ‡๐‘ข4=โ‹ฏ.(3.14) The solution of the first iteration given by๐‘ข1=๐‘’๐ด๐‘ก๐‘ข0.(3.15) Inserting this into (3.26) and expanding ๐‘’๐ด๐‘ก up to the order 1, we have 2nd order approximation of the exact solution which has a form๐‘ข2=๐‘’๐ต๐‘ก๐‘ข0+๐œ™1(๐‘ก๐ด)๐ด๐‘ก๐‘ข0.(3.16) Similarly, inserting (3.16) into (3.27) and expanding ๐œ™1(๐‘ก๐ด), we have a third order approximation of the exact solution๐‘ข3=๐‘’๐ต๐‘ก๐‘ข0+๐œ™1(๐‘ก๐ต)๐ต๐‘ก๐‘ข0+๐œ™2(๐‘ก๐ต)๐ต(๐ต+๐ด)๐‘ก2๐‘ข0.(3.17) In general for ๐‘–=0,2,4,โ€ฆ, we have the ๐‘=๐‘–+1-th order approximation of the exact solution in terms of the ๐œ™๐‘– function as follows:๐‘ข๐‘–=๐‘’๐ด๐‘ก๐‘ข0+๐œ™1(๐‘ก๐ด)๐ด๐‘ก๐‘ข0+๐œ™2(๐‘ก๐ด)๐ด(๐ต+๐ด)๐‘ก2๐‘ข0+โ‹ฏ+๐œ™๐‘–(๐‘ก๐ด)๐ด(๐ด+๐ต)๐‘–โˆ’1๐‘ข0.(3.18) For ๐‘–=1,3,5,โ€ฆ, we have๐‘ข๐‘–=๐‘’๐ต๐‘ก๐‘ข0+๐œ™1(๐‘ก๐ต)๐ต๐‘ก๐‘ข0+๐œ™2(๐‘ก๐ต)๐ต(๐ต+๐ด)๐‘ก2๐‘ข0+โ‹ฏ+๐œ™๐‘–(๐‘ก๐ต)๐ต(๐ต+๐ด)๐‘–โˆ’1๐‘ก๐‘–๐‘ข0.(3.19)

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 [๐‘ก0,๐‘‡], divide it into ๐‘ subintervals so that time step is โ„Ž=(๐‘‡โˆ’๐‘ก0)/๐‘.

Step 3. On each subinterval, [๐‘ก๐‘›,๐‘ก๐‘›+โ„Ž], ๐‘›=0,1,โ€ฆ,๐‘, use the algorithm by considering initial conditions for each step as ๐‘ข(๐‘ก0)=๐‘ข0, ๐‘ข๐‘–(๐‘ก๐‘›)=๐‘ข๐‘–โˆ’1(๐‘ก๐‘›)=๐‘ข(๐‘ก๐‘›), ๐‘ข2๎€ท๐‘ก๐‘›๎€ธ=๎€ท๐œ™+โ„Ž0(๐ต๐‘ก)+๐œ™1๎€ธ๐‘ข๎€ท๐‘ก(๐ต๐‘ก)๐ด๐‘ก๐‘›๎€ธ๐‘ข3๎€ท๐‘ก๐‘›๎€ธ=๎€ท๐œ™+โ„Ž0(๐ด๐‘ก)+๐œ™1(๐ด๐‘ก)๐ต๐‘ก+๐œ™2(๐ด๐‘ก)๐ต(๐ต+๐ด)๐‘ก2๎€ธ๐‘ข๎€ท๐‘ก๐‘›๎€ธโ‹ฎ๐‘ข2๐‘–๎€ท๐‘ก๐‘›๎€ธ=๎€ท๐œ™+โ„Ž0(๐ต๐‘ก)+๐œ™1(๐ต๐‘ก)๐ด๐‘ก+โ‹ฏ+๐œ™2๐‘–โˆ’1(๐ต๐‘ก)๐ด(๐ด+๐ต)2๐‘–โˆ’1๐‘ก2๐‘–๎€ธ๐‘ข2๐‘–โˆ’1๎€ท๐‘ก๐‘›๎€ธ๐‘ข2๐‘–+1๎€ท๐‘ก๐‘›๎€ธ=๎€ท๐œ™+โ„Ž0(๐ด๐‘ก)+๐œ™1(๐ด๐‘ก)๐ต๐‘ก+โ‹ฏ+๐œ™2๐‘–(๐ด๐‘ก)๐ต(๐ต+๐ด)2๐‘–๐‘ก2๐‘–+1๎€ธ๐‘ข2๐‘–๎€ท๐‘ก๐‘›๎€ธ.(3.20)

Step 4. ๐‘ข๐‘–(๐‘ก๐‘›+โ„Ž)โ†’๐‘ข(๐‘ก๐‘›+โ„Ž).

Step 5. Repeat this procedure for next interval until the desired time ๐‘‡ 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: ๎€บ๎€ปexp(โˆ’๐‘ก๐ด)๐ตexp(๐‘ก๐ด)=๐ต,exp(๐‘ก๐ด),(3.21) where [โ‹…,โ‹…] is the matrix commutator.
The integration is given as๎€œ๐‘ก0๎€บexp(โˆ’๐‘ ๐ด)๐ตexp(๐‘ ๐ด)๐‘‘๐‘ =๐ต,๐œ™1(๎€ป,๐‘ก๐ด)(3.22) the representation is given in [13].
Further, we have the recursive integration๎€บ๐ต,๐œ™๐‘–(๐‘ก๐ด)๐‘ก๐‘–๎€ป=๎€œ๐‘ก0๎€บ๐ต,๐œ™๐‘–โˆ’1(๎€ป๐‘ ๐ด)๐‘‘๐‘ ,(3.23) where ๐œ™๐‘– given defined as ๐œ™0๐œ™(๐ด๐‘ก)=exp(๐ด๐‘ก),๐‘–(๎€œ๐ด๐‘ก)=๐‘ก0๐œ™๐‘–โˆ’1(๐œ™๐ด๐‘ )๐‘‘๐‘ ,๐‘–(๐œ™๐ด๐‘ก)=๐‘–โˆ’1๎€ท๐‘ก(๐ด๐‘ก)โˆ’๐ผ๐‘–โˆ’1๎€ธ/(๐‘–โˆ’1)!๐ด.(3.24)

The iterative scheme with the equations ฬ‡๐‘ข1=๐ด๐‘ข1,(3.25)ฬ‡๐‘ข2=๐ด๐‘ข2+๐ต๐‘ข1,(3.26)ฬ‡๐‘ข3=๐ด๐‘ข3+๐ต๐‘ข2,(3.27)ฬ‡๐‘ข4=โ‹ฏ,(3.28) is solved as๐‘1(๐‘ก)=exp(๐ด๐‘ก)๐‘(๐‘กn๐‘),2(๐‘ก)=๐‘1(๐‘ก)+๐‘1๎€œ(๐‘ก)๐‘ก0๎€บ๎€ป๐‘๐ต,exp(๐‘ ๐ด)๐‘‘๐‘ ,2(๐‘ก)=๐‘1(๐‘ก)+๐‘1๎€บ(๐‘ก)๐ต,๐œ™1๎€ป,๐‘(๐‘ก๐ด)3(๐‘ก)=๐‘2(๐‘ก)+๐‘1๎€œ(๐‘ก)๐‘ก0๎€บ๐ต,exp(๐‘ ๐ด),๐ต,๐œ™1๎€ป๐‘(๐‘ ๐ด)๐‘‘๐‘ ,3(๐‘ก)=๐‘2(๐‘ก)+๐‘1(๐‘ก)๎€ท๎€บ๐ต,exp(๐‘ก๐ด),๐ต,๐œ™2(๎€ป+๎€บ๐‘ก๐ด)๐ต,๐ดexp(๐‘ก๐ด),๐ต,๐œ™3(๎€ท๐‘ก๐‘ก๐ด)๎€ป๎€ธ+๐‘‚3๎€ธ,โ‹ฏ(3.29)

The recursion is given as๐‘๐‘–(๐‘ก)=๐‘๐‘–โˆ’1(๐‘ก)+๐‘1(๎€œ๐‘ก)๐‘ก0๎€บ๎€ป๎€œ๐ต,exp(๐‘ ๐ด)๐‘ 10๎€บ๎€ท๐‘ ๐ต,exp1๐ด๎€œ๎€ธ๎€ป๐‘ 20โ‹ฏ๎€œ๐‘ ๐‘–โˆ’20๎€บ๎€ท๐‘ ๐ต,exp๐‘–โˆ’1๐ด๎€ธ๎€ป๐‘‘๐‘ ๐‘–โˆ’1โ‹ฏ๐‘‘๐‘ 1๐‘‘๐‘ก.(3.30)

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๐‘ข๎…ž๎‚ธ๎‚น(๐‘ก)=1230๐‘ข,๐‘ข(0)=๐‘ข0=๎‚ต01๎‚ถ,(4.1) the exact solution is๐‘ข2๎€ท๐‘’(๐‘ก)=3๐‘กโˆ’๐‘’โˆ’2๐‘ก๎€ธ5.(4.2) We split the matrix as๎‚ธ๎‚น=๎‚ธ๎‚น+๎‚ธ๎‚น123011100120.(4.3)

The Figure 1 presents the numerical errors between the exact and the numerical solution.

Figure 1: Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps.

The Figure 2 presents the CPU time of the standard and the iterative splitting schemes.

Figure 2: CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps.
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 the 10ร—10 ODE system:๐œ•๐‘ก๐‘ข1=โˆ’๐œ†1,1๐‘ข1+๐œ†2,1๐‘ข2+โ‹ฏ+๐œ†10,1๐‘ข10๐œ•,(4.4)๐‘ก๐‘ข2=๐œ†1,2๐‘ข1โˆ’๐œ†2,2(๐‘ก)๐‘ข2+โ‹ฏ+๐œ†10,2๐‘ข10๐œ•,(4.5)โ‹ฎ(4.6)๐‘ก๐‘ข10=๐œ†1,10๐‘ข1+๐œ†2,10(๐‘ก)๐‘ข2+โ‹ฏโˆ’๐œ†10,10๐‘ข10๐‘ข,(4.7)1(0)=๐‘ข1,0,โ€ฆ,๐‘ข10(0)=๐‘ข10,0(initialconditions),(4.8) where ๐œ†1(๐‘ก)โˆˆโ„+ and ๐œ†2(๐‘ก)โˆˆโ„+ are the decay factors and ๐‘ข1,0,โ€ฆ,๐‘ข10,0โˆˆโ„+. We have the time interval ๐‘กโˆˆ[0,๐‘‡].

We rewrite (4.4) in operator notation, we concentrate us to the following equations:๐œ•๐‘ก๐‘ข=๐ด(๐‘ก)๐‘ข+๐ต(๐‘ก)๐‘ข,(4.9) where ๐‘ข1(0)=๐‘ข10=1.0, ๐‘ข2(0)=๐‘ข20=1.0 are the initial conditions, where we have ๐œ†1(๐‘ก)=๐‘ก and ๐œ†2(๐‘ก)=๐‘ก2.

The operators are splitted in the following way, while operator ๐ด, see (4.10), covers the upper diagonal entries and operator ๐ต, 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โŽ›โŽœโŽœโŽœโŽœโŽœโŽœโŽœโŽœโŽโ‹ฎโŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽŸโŽŸโŽŸโŽ ๐ด=โˆ’0.010.010โ‹ฏ0.01โˆ’0.010โ‹ฏ0.010.01โˆ’0.020โ‹ฏ0.010.010.01โˆ’0.030โ‹ฏ0.โˆ’0.0800.โˆ’0.08,(4.10)โŽ›โŽœโŽœโŽœโŽœโŽœโŽœโŽโ‹ฎโŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽŸโŽ ๐ต=โˆ’0.0800.โˆ’โˆ’โˆ’0.010.01000000000.01โˆ’0.01.(4.11)

The Figure 3 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.

Figure 3: Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps.

The Figure 4 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.

Figure 4: CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps.

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 [๐ด,๐ต] and deal with the matrix equation๐‘ข๎…ž๎‚ธ๎‚น(๐‘ก)=101110๐‘ข,๐‘ข(0)=๐‘ข0=๎‚ต11๎‚ถ.(4.12)

In the following, we discuss different splitting ideas.

Version 1
We split the matrix as ๎‚ธ๎‚น=๎‚ธ๎‚น+๎‚ธ๎‚น10111090111109,(4.13) while โ€–[๐ด,๐ต]โ€–โ‰ฅmax{โ€–๐ดโ€–,โ€–๐ตโ€–}.

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 Figure 5 present the numerical errors between the exact and the numerical solution.

Figure 5: Numerical errors of the standard splitting scheme (a) and the iterative schemes (b) with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps.

The Figure 6 present the CPU time of the standard and the iterative splitting schemes.

Figure 6: CPU time of the standard splitting scheme (a) and the iterative schemes (b) with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps.

Further for the one-side,we obtain more improved results for the following splitting.

Version 2
In this version, we have๎‚ธ๎‚น=๎‚ธ๎‚น+๎‚ธ๎‚น10111090191101.(4.14) The Figure 7 presents the numerical errors between the exact and the numerical solution.

Figure 7: Numerical errors of the standard splitting scheme (a) and the iterative schemes based on one-side to operator ๐ด or ๐ต (b) and two-side scheme (c).

A more delicate problem is given for the stiff matrices.

Version 3
In this version, we have๎ƒฌ10411104๎ƒญ=๎ƒฌ104โˆ’101104๎ƒญ+๎‚ธ๎‚นโˆ’11101.(4.15) The Figure 8 present the numerical errors between the exact and the numerical solution.

Figure 8: Numerical errors of the one-side splitting scheme with ๐ด (a), the one-side splitting scheme with ๐ต (b) and the iterative schemes with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps (c).

The Figure 9 present the CPU time of the standard and the iterative splitting schemes.

Figure 9: CPU time of the one-side splitting scheme with ๐ด (a) and the iterative schemes with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps (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 standard ๐ด-๐ต 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:๐œ•๐‘ก๐‘1+โˆ‡โ‹…๐…๐‘1๎€ท=๐‘”โˆ’๐‘1+๐‘1,im๎€ธโˆ’๐œ†1๐‘1[],๐œ•,inฮฉร—0,๐‘ก๐‘ก๐‘2+โˆ‡โ‹…๐…๐‘2๎€ท=๐‘”โˆ’๐‘2+๐‘2,im๎€ธ+๐œ†1๐‘1โˆ’๐œ†2๐‘2[],๐œ•,inฮฉร—0,๐‘ก๐…=๐ฏโˆ’๐ทโˆ‡,๐‘ก๐‘1,im๎€ท๐‘=๐‘”1โˆ’๐‘1,im๎€ธโˆ’๐œ†1๐‘1,im[],๐œ•,inฮฉร—0,๐‘ก๐‘ก๐‘2,im๎€ท๐‘=๐‘”2โˆ’๐‘2,im๎€ธ+๐œ†1๐‘1,imโˆ’๐œ†2๐‘2,im[],๐‘,inฮฉร—0,๐‘ก1(๐ฑ,๐‘ก)=๐‘1,0(๐ฑ),๐‘2(๐ฑ,๐‘ก)=๐‘2,0๐‘(๐ฑ),onฮฉ,1(๐ฑ,๐‘ก)=๐‘1,1(๐ฑ,๐‘ก),๐‘2(๐ฑ,๐‘ก)=๐‘2,1([],๐‘๐ฑ,๐‘ก),on๐œ•ฮฉร—0,๐‘ก1,im(๐ฑ,๐‘ก)=0,๐‘2,im๐‘(๐ฑ,๐‘ก)=0,onฮฉ,1,im(๐ฑ,๐‘ก)=0,๐‘2,im[].(๐ฑ,๐‘ก)=0,on๐œ•ฮฉร—0,๐‘ก(4.16)

In the following, we deal with the semidiscretized equation given with the matrices:๐œ•๐‘กโŽ›โŽœโŽœโŽœโŽ๐‚=๐ดโˆ’ฮ›1ฮ›โˆ’๐บ0๐บ01๐ดโˆ’ฮ›2โˆ’๐บ0๐บ๐บ0โˆ’ฮ›1โˆ’๐บ00๐บฮ›1โˆ’ฮ›2โŽžโŽŸโŽŸโŽŸโŽ โˆ’๐บ๐‚,(4.17) where ๐‚=(๐œ1,๐œ2,๐œ1im,๐œ2im)๐‘‡, while ๐œ1=(๐œ1,1,โ€ฆ,๐œ1,๐ผ) is the solution of the first species in the mobile phase in each spatial discretization point (๐‘–=1,โ€ฆ,๐ผ), the same is also for the other solution vectors.

We have the following two operators for the splitting method๐ท๐ด=ฮ”๐‘ฅ2โ‹…โŽ›โŽœโŽœโŽœโŽœโŽœโŽโŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ +๐‘ฃโˆ’211โˆ’21โ‹ฑโ‹ฑโ‹ฑ1โˆ’211โˆ’2โ‹…โŽ›โŽœโŽœโŽœโŽœโŽœโŽ1โŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ ฮ”๐‘ฅโˆ’11โ‹ฑโ‹ฑโˆ’11โˆ’11โˆˆโ„๐ผร—๐ผ,(4.18) where ๐ผ is the number of spatial points.ฮ›1=โŽ›โŽœโŽœโŽœโŽœโŽœโŽ๐œ†100๐œ†10โ‹ฑโ‹ฑโ‹ฑ0๐œ†100๐œ†1โŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ โˆˆโ„๐ผร—๐ผ,ฮ›2=โŽ›โŽœโŽœโŽœโŽœโŽœโŽ๐œ†200๐œ†20โ‹ฑโ‹ฑโ‹ฑ0๐œ†200๐œ†2โŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ โˆˆโ„๐ผร—๐ผ,โŽ›โŽœโŽœโŽœโŽœโŽœโŽโŽžโŽŸโŽŸโŽŸโŽŸโŽŸโŽ ๐บ=๐‘”00๐‘”0โ‹ฑโ‹ฑโ‹ฑ0๐‘”00๐‘”โˆˆโ„๐ผร—๐ผ.(4.19)

We decouple into the following matrices:๐ด1=โŽ›โŽœโŽœโŽœโŽโŽžโŽŸโŽŸโŽŸโŽ ๐ด0000๐ด0000000000โˆˆโ„4๐ผร—4๐ผ,๎‚๐ด2=โŽ›โŽœโŽœโŽœโŽโˆ’ฮ›1ฮ›0001โˆ’ฮ›20000โˆ’ฮ›1000ฮ›1โˆ’ฮ›2โŽžโŽŸโŽŸโŽŸโŽ โˆˆโ„4๐ผร—4๐ผ,๎‚๐ด3=โŽ›โŽœโŽœโŽœโŽโŽžโŽŸโŽŸโŽŸโŽ โˆ’๐บ0๐บ00โˆ’๐บ0๐บ๐บ0โˆ’๐บ00๐บ0โˆ’๐บโˆˆโ„4๐ผร—4๐ผ.(4.20)

For the operator ๐ด1 and ๐ด2=๎‚๐ด2+๎‚๐ด3, we apply the iterative splitting method.

Based on the decomposition, operator ๐ด1 is only tridiagonal and operator ๐ด2 is block diagonal. Such matrix structure reduce the computation of the exponential operators.

The Figure 10 present the numerical errors between the exact and the numerical solution. Here we obtain optimal results for one-side iterative schemes on operator ๐ต, means we iterate with respect to ๐ต and use ๐ด as right-hand side.

Figure 10: Numerical errors of the one-side splitting scheme with ๐ด (a), the one-side splitting scheme with ๐ต (b) and the iterative schemes with 1,โ€‰โ€ฆโ€‰,โ€‰6 iterative steps (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.


  1. J. Geiser, Discretisation methods for systems of convective-diffusive dispersive-reactive equations and applications, Ph.D. thesis, Universität Heidelberg, Heidelberg, Germany, 2004.
  2. J. Geiser, โ€œDecomposition methods for partial differential equations: theory and applications in multiphysics problems,โ€ in Numerical Analysis and Scientific Computing Series, Magoules and Lai, Eds., CRC Press, Boca Raton, Fla, USA, 2009. View at Google Scholar
  3. M. Hochbruck and A. Ostermann, โ€œExplicit exponential Runge-Kutta methods for semilinear parabolic problems,โ€ SIAM Journal on Numerical Analysis, vol. 43, no. 3, pp. 1069โ€“1090, 2005. View at Publisher ยท View at Google Scholar ยท View at Zentralblatt MATH ยท View at MathSciNet
  4. H. Rouch, โ€œMOCVD research reactor simulation,โ€ in Proceedings of the COMSOL Users Conference, Paris, France, 2006.
  5. L. Rudniak, โ€œNumerical simulation of chemical vapour deposition process in electric field,โ€ Computers and Chemical Engineering, vol. 22, no. 1, supplement, pp. S755โ€“S758, 1998. View at Google Scholar
  6. R. Glowinski, โ€œNumerical methods for fluids,โ€ in Handbook of Numerical Analysis, P. G. Ciarlet and J. Lions, Eds., vol. 9, North-Holland, Amsterdam, The Netherlands, 2003. View at Google Scholar
  7. J. F. Kanney, C. T. Miller, and C. T. Kelley, โ€œConvergence of iterative split-operator approaches for approximating nonlinear reactive problems,โ€ Advances in Water Resources, vol. 26, no. 3, pp. 247โ€“261, 2003. View at Publisher ยท View at Google Scholar
  8. E. Hansen and A. Ostermann, โ€œExponential splitting for unbounded operators,โ€ Mathematics of Computation, vol. 78, no. 267, pp. 1485โ€“1496, 2009. View at Publisher ยท View at Google Scholar ยท View at Zentralblatt MATH
  9. T. Jahnke and C. Lubich, โ€œError bounds for exponential operator splittings,โ€ BIT Numerical Mathematics, vol. 40, no. 4, pp. 735โ€“744, 2000. View at Publisher ยท View at Google Scholar ยท View at Zentralblatt MATH ยท View at MathSciNet
  10. I. Faragó and J. Geiser, โ€œIterative operator-splitting methods for linear problems,โ€ International Journal of Computational Science and Engineering, vol. 3, no. 4, pp. 255โ€“263, 2007. View at Publisher ยท View at Google Scholar
  11. E. Hairer, S. P. Nørsett, and G. Wanner, Solving ordinary differential equations. I, vol. 8 of Springer Series in Computational Mathematics, Springer, Berlin, 1987, Nonstiff problem.
  12. E. Hairer and G. Wanner, Solving Ordinary Differential Equations. II, vol. 14 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 2nd edition, 1996.
  13. W. Magnus, โ€œOn the exponential solution of differential equations for a linear operator,โ€ Communications on Pure and Applied Mathematics, vol. 7, pp. 649โ€“673, 1954. View at Publisher ยท View at Google Scholar ยท View at Zentralblatt MATH