Computing Exponential for Iterative Splitting Methods: Algorithms and Applications

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.

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/cm³),𝑐_𝑖,im:concentration of the 𝑖th gaseous species in the immobile zones of the plasma chamber phase (mol/cm³),𝑐_𝑖,ad:concentration of the 𝑖th gaseous species in the mobile adsorbed zones of the plasma chamber phase (mol/cm³),𝑐_𝑖,im,ad:concentration of the 𝑖th gaseous species in the immobile adsorbed zones of the plasma chamber phase (mol/cm³),𝐯:velocity through the chamber and porous substrate [4] (cm/nsec),𝐷^𝑒(𝑖):element-specific diffusions-dispersions tensor (cm²/nsec),𝜆_𝑖,𝑖:decay constant of the 𝑖th species (1/nsec),𝑄_𝑖:source term of the 𝑖th species (mol/(cm³nsec)),𝑔: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] (cm²/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)=𝑢₀,(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)+𝐵𝑐_𝑖−1(𝑡),with𝑐_𝑖(𝑡^𝑛)=𝑐^𝑛,𝑐₀(𝑡^𝑛)=𝑐^𝑛,𝑐₋₁=0.0,(𝑡)𝜕𝑡=𝐴𝑐^𝑖(𝑡)+𝐵𝑐_𝑖+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 𝐶₀ 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 𝐶₀ semigroup based on their underlying graph norms.

Theorem 2.1. Let one consider the abstract Cauchy problem in a Hilbert space 𝐗𝜕_𝑡𝑐(𝑥,𝑡)=𝐴𝑐𝑐((𝑥𝑥,,𝑡0))+=𝐵𝑐𝑐₀(𝑥,𝑡),0<𝑡≤𝑇,𝑥∈Ω,𝑐(𝑥,𝑡)=𝑐₁(𝑥),𝑥∈Ω,(𝑥,𝑡),𝑥∈𝜕Ω×[0,𝑇],(2.2) where 𝐴,𝐵∶𝐷(𝐗)→𝐗 are given linear operators which are generators of the 𝐶₀-semigroup and 𝑐₀∈𝐗 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((𝐴+𝐵)𝜏)𝑢₀‖‖≤𝜅.(2.3)
Further we assume the estimation of the partial integration of the unbounded operator 𝐵 (see [8]): ‖‖‖𝐵^𝜏0exp(𝐵𝑠)𝑠𝑑𝑠‖‖‖≤𝜏𝐶.(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,(2.6) and for the second iteration we have: 𝜕_𝑡𝑐₂(𝑡)=𝐴𝑐₁(𝑡)+𝐵𝑐₂(𝑡),𝑡∈𝑡^𝑛,𝑡^𝑛+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,(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:𝑐₂𝑐₁(𝑡)=exp𝐴𝑡^𝑛+1−𝑡𝑐(𝑡^𝑛),𝑡∈𝑡^𝑛,𝑡^𝑛+1,(𝑡)=exp(𝐵𝑡)𝑐(𝑡^𝑛)+^𝑡𝑛+1_𝑡𝑛exp𝐵𝑡^𝑛+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 𝑒₁, we have 𝑐(𝜏)=exp((𝐴+𝑐𝐵₁)(𝜏𝜏))=𝑐(e𝑡_𝑛xp(𝐴𝜏)𝑐(𝑡^𝑛),+_{𝑡𝑛+1𝑡𝑛})=exp(𝐴𝜏)𝑐(𝑡^𝑛)exp𝐴𝑡^𝑛+1−𝑠𝐵exp(𝑠(𝐴+𝐵))𝑐(𝑡^𝑛)𝑑𝑠.(2.13)
We obtain‖‖𝑒₁‖‖=‖‖𝑐−𝑐₁‖‖≤‖exp((𝐴+𝐵)𝜏)𝑐(𝑡^𝑛)−exp(𝐴𝜏)𝑐(𝑡^𝑛)‖≤𝐶₁𝜏𝑐(𝑡^𝑛).(2.14)For 𝑒₂ we have (alternating) 𝑐𝑐₂((𝜏𝜏))==eexxpp((𝐵𝐵𝜏𝜏))𝑐𝑐((𝑡𝑡^𝑛_𝑛)+^𝑡𝑛+1_𝑡𝑛exp𝐵𝑡^𝑛+1−𝑠𝐴exp(𝑠𝐴)𝑐(𝑡^𝑛)𝑑𝑠,+)+^𝑡𝑛+1_𝑡𝑛exp𝐵𝑡^𝑛+1−𝑠𝐴exp(𝑠𝐴)𝑐(𝑡^𝑛)𝑑𝑠_{𝑡𝑛+1𝑡𝑛}exp(𝐵𝑠)𝐴^{𝑡𝑛+1𝑡−𝑠}_𝑛exp𝐴𝑡^𝑛+1−𝑠−𝜌𝐵exp(𝜌(𝐴+𝐵))𝑐(𝑡^𝑛)𝑑𝜌𝑑𝑠.(2.15)
We obtain‖‖𝑒₂‖‖≤‖‖exp((𝐴+𝐵)𝜏)𝑐(𝑡^𝑛)−𝑐₂‖‖≤𝐶₂𝜏²𝑐(𝑡^𝑛).(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:𝑐_𝑖(𝜏)=+exp_{𝑡𝑛+1𝑡𝑛}(𝐴𝜏)𝑐(𝑡^𝑛)+exp𝐴𝑡^𝑛+1−𝑠𝐵exp(𝑠𝐴)𝑐(𝑡^𝑛)𝑑𝑠_{𝑡𝑛+1𝑡𝑛}exp𝐴𝑡^𝑛+1−𝑠₁𝐵^{𝑡𝑛+1−𝑠1}_𝑡𝑛exp𝑡^𝑛+1−𝑠₁−𝑠₂𝐴𝐵exp𝑠₂𝐴𝑐(𝑡^𝑛)𝑑𝑠₂𝑑𝑠₁++⋯_{𝑡𝑛+1𝑡𝑛}+exp𝐴𝑡^𝑛+1−𝑠₁𝐵^{𝑡𝑛+1−𝑠1}_𝑡𝑛exp𝑡^𝑛+1−𝑠₁−𝑠₂𝐴𝐵exp𝑠₂𝐴𝑐(𝑡^𝑛)𝑑𝑠₂𝑑𝑠₁+⋯_{𝑡𝑛+1𝑡𝑛}exp𝐴𝑡^𝑛+1−𝑠₁𝐵⋯^{𝑡𝑛+1−∑𝑖−2𝑗=1𝑠𝑗}_𝑡𝑛exp𝑡^𝑛+1−^𝑖−1_𝑗=1𝑠_𝑗×exp𝑠_𝑖𝐴𝐵𝑐(𝜏)=exp(𝐴𝜏)𝑐(𝑡_𝑛𝐴𝑐(𝑡^𝑛)𝑑𝑠_𝑖−1⋯𝑑𝑠₂𝑑𝑠₁,+_{𝑡𝑛+1𝑡𝑛})+exp𝐴𝑡^𝑛+1−𝑠𝐵exp(𝑠𝐴)𝑐(𝑡^𝑛)𝑑𝑠_{𝑡𝑛+1𝑡𝑛}exp𝐴𝑡^𝑛+1−𝑠₁𝐵^{𝑡𝑛+1−𝑠1}_𝑡𝑛exp𝑡^𝑛+1−𝑠₁−𝑠₂𝐴𝐵exp𝑠₂𝐴𝑐(𝑡^𝑛)𝑑𝑠₂𝑑𝑠₁++⋯_{𝑡𝑛+1𝑡𝑛}+exp𝐴𝑡^𝑛+1−𝑠₁𝐵^{𝑡𝑛+1−𝑠1}_𝑡𝑛exp𝑡^𝑛+1−𝑠₁−𝑠₂𝐴𝐵exp𝑠₂𝐴𝑐(𝑡^𝑛)𝑑𝑠₂𝑑𝑠₁+⋯_{𝑡𝑛+1𝑡𝑛}exp𝐴𝑡^𝑛+1−𝑠₁𝐵⋯^{𝑡𝑛+1−∑𝑖−2𝑗=1𝑠𝑗}_𝑡𝑛exp𝑡^𝑛+1−^𝑖_𝑗=1𝑠_𝑗×exp𝑠_𝑖𝐴𝐵𝐴𝑐(𝑡^𝑛)𝑑𝑠_𝑖⋯𝑑𝑠₂𝑑𝑠₁.(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 a‖‖‖‖n𝑒𝑒d₁‖‖he=nc𝐾e,₂‖𝐵‖𝜏^𝛼1‖‖𝑒₀‖‖+𝒪𝜏^1+𝛼1,‖‖=𝐾‖𝐵‖‖‖𝑒₀‖‖𝜏^{1+𝛼1+𝛼2}+𝒪𝜏^1+𝛼1+𝛼,(2.19) where 0≤𝛼₁, 𝛼₂<1.

Remark 2.3. If we assume the consistency of 𝒪(𝜏^𝑚) for the initial value 𝑒₁(𝑡^𝑛) and 𝑒₂(𝑡^𝑛), 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:𝜕𝜕_𝑡𝐜=𝐴₁𝐜+𝐴₂𝐜+𝐵₁𝐜−𝐜_im+𝐵₂𝐜−𝐜_ad+𝐐,_𝑡𝐜_im=𝐴₂𝜕𝐜_im𝑡+𝐵₁𝐜_im−𝐜+𝐵₂𝐜_im−𝐜_im,ad+𝐐_im,𝜕_𝑡𝐜_ad=𝐴₂𝐜_ad+𝐵₂𝐜_ad−𝐜+𝐐_ad,𝐜_im,ad=𝐴₂𝐜_im,ad+𝐵₂𝐜_im,ad−𝐜_im+𝐐_im,ad,(2.20) where 𝐜=(𝑐₁,…,𝑐_𝑚)^𝑇 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. 𝐴₁ is the stiffness matrix of (1.1), 𝐴₂ is the reaction matrix of the right-hand side of (1.1), 𝐵₁ and 𝐵₂ 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:𝜕_𝑡𝐂=⎛⎜⎜⎜⎝𝐴¹+𝐴₂−+𝐵₁𝐵₁+𝐵₂−𝐵₁−𝐵₂0𝐴₂+𝐵₁+𝐵₂0−𝐵₂−0𝐵₂−0𝐵₂𝐴₂+𝐵₂00𝐴₂+𝐵₂⎞⎟⎟⎟⎠𝐂+𝐐,(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:𝜕𝜕_𝑡𝐂=𝐴𝐂+𝐐,_𝑡𝐂=𝐴₁𝐂+𝐴₂𝐂+𝐐,(2.22) where𝐴₁=⎛⎜⎜⎜⎝𝐴¹+0𝐴₂𝐴0₂0000𝐴0₂0000𝐴0₂⎞⎟⎟⎟⎠,𝐴₂=⎛⎜⎜⎜⎝𝐵¹−+𝐵₁𝐵₂−𝐵₁−𝐵₂0𝐵₁+𝐵₂0−𝐵₂−0𝐵₂−0𝐵₂𝐵₂00𝐵₂⎞⎟⎟⎟⎠.(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, and 𝑡^𝑁=𝑇.
We start with 𝑛=0.
(1)The initial conditions are given with 𝐂₀(𝑡^𝑛+1)=𝐂(𝑡^𝑛). We start with 𝑘=0.(2)Compute the fix-point iteration scheme given as 𝜕_𝑡𝐂^𝑘=𝐴₁𝐂^𝑘+𝐴₂𝐂^𝑘−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⁻⁴.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 𝑡₁=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𝐴=⎛⎜⎜⎜⎜⎜⎝⋮𝐴𝐴00⋱⋯𝐵𝐵0⋱⋯𝐵𝐴0⋱⋯⋯⋯𝐴⎞⎟⎟⎟⎟⎟⎠,(3.2) the computation of the exp-Matrix can be expressed as:exp𝐴𝑡∶=⎛⎜⎜⎜⎜⎜⎜⎜⎝e𝐷xp(𝐴𝑡)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:
𝜙_𝑘+1𝜙₀(𝑥)=𝑒^𝑥,(𝑥)=𝜙^𝑘(𝑥)−1/𝑘!𝑥,𝑘≥0.(3.4)
So the matrix formulation of our scheme is given as𝑌(𝑡)=exp𝐴𝑡𝑌(0),(3.5) where 𝐴 is given as, 𝐴=⎛⎜⎜⎜⎜⎜⎝⋮𝐴𝐴00⋱⋯𝐵𝐴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⋯⋯⋯¹𝜙(𝐵𝑡)𝐴𝑡exp(𝐵𝑡)0⋯⋯₂𝜙(𝐴𝑡)𝐵⋮(𝐵+𝐴)𝑡²𝜙₁(⋱𝐴𝑡)𝐵𝑡exp⋱(𝐴𝑡)⋱⋯⋯⋮_𝑖(𝐴𝑡)𝐵(𝐴+𝐵)^𝑖−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 𝑐₁,𝑐₂,….

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𝑥=𝑡𝜙₁(𝑡𝐴)𝑎.(3.9) Similarly, the equation𝑥̇=𝐴𝑥+𝑏𝑡,𝑥(0)=0,(3.10) has a solution of the form𝑥=𝑡²𝜙₂(𝑡𝐴)𝑏.(3.11) In general, the solution of the equation𝑥̇=𝐴𝑥+𝑎+𝑏𝑡+𝑐𝑡²2!+⋯,𝑥(0)=0,(3.12) admits the form𝑥=𝑡𝜙₁(𝑡𝐴)𝑎+𝑡²𝜙₂(𝑡𝐴)𝑏+𝑡³𝜙₃(𝑡𝐴)𝑐+⋯.(3.13) In this section, we use the formula (3.21) for iterative schemes given as𝑢̇₂𝑢̇₁=𝐴𝑢₁,𝑢̇₃=𝐴𝑢₁+𝐵𝑢₂,=𝑢̇₄𝐴𝑢₃+𝐵𝑢₂,=⋯.(3.14) The solution of the first iteration given by𝑢₁=𝑒^𝐴𝑡𝑢₀.(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𝑢₂=𝑒^𝐵𝑡𝑢₀+𝜙₁(𝑡𝐴)𝐴𝑡𝑢₀.(3.16) Similarly, inserting (3.16) into (3.27) and expanding 𝜙₁(𝑡𝐴), we have a third order approximation of the exact solution𝑢₃=𝑒^𝐵𝑡𝑢₀+𝜙₁(𝑡𝐵)𝐵𝑡𝑢₀+𝜙₂(𝑡𝐵)𝐵(𝐵+𝐴)𝑡²𝑢₀.(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:𝑢_𝑖=𝑒^𝐴𝑡𝑢₀+𝜙₁(𝑡𝐴)𝐴𝑡𝑢₀+𝜙₂(𝑡𝐴)𝐴(𝐵+𝐴)𝑡²𝑢₀+⋯+𝜙_𝑖(𝑡𝐴)𝐴(𝐴+𝐵)^𝑖−1𝑢₀.(3.18) For 𝑖=1,3,5,…, we have𝑢_𝑖=𝑒^𝐵𝑡𝑢₀+𝜙₁(𝑡𝐵)𝐵𝑡𝑢₀+𝜙₂(𝑡𝐵)𝐵(𝐵+𝐴)𝑡²𝑢₀+⋯+𝜙_𝑖(𝑡𝐵)𝐵(𝐵+𝐴)^𝑖−1𝑡^𝑖𝑢₀.(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 [𝑡₀,𝑇], divide it into 𝑁 subintervals so that time step is ℎ=(𝑇−𝑡₀)/𝑁.

Step 3. On each subinterval, [𝑡_𝑛,𝑡_𝑛+ℎ], 𝑛=0,1,…,𝑁, use the algorithm by considering initial conditions for each step as 𝑢(𝑡₀)=𝑢₀, 𝑢_𝑖(𝑡_𝑛)=𝑢_𝑖−1(𝑡_𝑛)=𝑢(𝑡_𝑛), 𝑢𝑢₂₃𝑡_𝑛+ℎ=𝜙₀(𝐵𝑡)+𝜙₁(𝐵𝑡)𝐴𝑡𝑢𝑡_𝑛𝑢_2𝑖𝑡_𝑛+ℎ=⋮𝜙₀(𝐴𝑡)+𝜙₁(𝐴𝑡)𝐵𝑡+𝜙₂(𝐴𝑡)𝐵(𝐵+𝐴)𝑡²𝑢𝑡_𝑛𝑢_2𝑖+1𝑡_𝑛+ℎ=𝜙₀(𝐵𝑡)+𝜙₁(𝐵𝑡)𝐴𝑡+⋯+𝜙_2𝑖−1(𝐵𝑡)𝐴(𝐴+𝐵)^2𝑖−1𝑡^2𝑖𝑢_2𝑖−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^𝑡0exp(−𝑠𝐴)𝐵exp(𝑠𝐴)𝑑𝑠=𝐵,𝜙₁(𝑡𝐴),(3.22) the representation is given in [13].
Further, we have the recursive integration𝐵,𝜙_𝑖(𝑡𝐴)𝑡^𝑖=^𝑡0𝐵,𝜙_𝑖−1(𝑠𝐴)𝑑𝑠,(3.23) where 𝜙_𝑖 given defined as 𝜙𝜙₀(𝐴𝑡)=exp(𝐴𝑡),_𝑖𝜙(𝐴𝑡)=^𝑡0𝜙_𝑖−1(𝐴𝑠)𝑑𝑠,_𝑖(𝐴𝑡)=𝜙^𝑖−1(𝐴𝑡)−𝐼𝑡^𝑖−1/(𝑖−1)!𝐴.(3.24)

The iterative scheme with the equations 𝑢̇₂𝑢̇₁=𝐴𝑢₁,(3.25)𝑢̇₃=𝐴𝑢₂+𝐵𝑢₁,(3.26)=𝑢𝐴̇₄𝑢₃+𝐵𝑢₂,(3.27)=⋯,(3.28) is solved as𝑐𝑐₁(𝑡)=exp(𝐴𝑡)𝑐(𝑡ⁿ),₂(𝑡𝑐)=𝑐₁(𝑡)+𝑐₁(𝑡)^𝑡0𝐵,exp(𝑠𝐴)𝑑𝑠,₂𝑐(𝑡)=𝑐₁(𝑡)+𝑐₁(𝑡)𝐵,𝜙₁(𝑡𝐴)₃,𝑐(𝑡)=𝑐₂(𝑡)+𝑐₁(𝑡)^𝑡0𝐵,exp(𝑠𝐴),𝐵,𝜙₁(𝑠𝐴)𝑑𝑠,₃(𝑡)=𝑐₂(𝑡)+𝑐₁(𝑡)𝐵,exp(𝑡𝐴),𝐵,𝜙₂(𝑡⋯𝐴)+𝐵,𝐴exp(𝑡𝐴),𝐵,𝜙₃(𝑡𝐴)+𝑂𝑡³,(3.29)

The recursion is given as𝑐_𝑖(𝑡)=𝑐_𝑖−1(𝑡)+𝑐₁(𝑡)^𝑡0𝐵,exp(𝑠𝐴)^𝑠1₀𝐵,exp𝑠₁𝐴^𝑠2₀⋯^𝑠𝑖−2₀𝐵,exp𝑠_𝑖−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𝑢^(𝑡)=1320𝑢,𝑢(0)=𝑢₀=01,(4.1) the exact solution is𝑢(𝑡)=2𝑒^3𝑡−𝑒^−2𝑡5.(4.2) We split the matrix as1320=1110+0210.(4.3)

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

(a)

(b)

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.

(a)

(b)

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𝑢₁+𝜆_2,1𝑢₂+⋯+𝜆_10,1𝑢₁₀,(4.4)_𝑡𝜕𝑢₂=⋮𝜆_1,2𝑢₁−𝜆_2,2(𝑡)𝑢₂+⋯+𝜆_10,2𝑢₁₀,((44..56))_𝑡𝑢𝑢₁₀=𝜆_1,10𝑢₁+𝜆_2,10(𝑡)𝑢₂+⋯−𝜆_10,10𝑢₁₀,(4.7)₁(0)=𝑢_1,0,…,𝑢₁₀(0)=𝑢_10,0(initialconditions),(4.8) where 𝜆₁(𝑡)∈ℝ⁺ and 𝜆₂(𝑡)∈ℝ⁺ 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 𝑢₁(0)=𝑢₁₀=1.0, 𝑢₂(0)=𝑢₂₀=1.0 are the initial conditions, where we have 𝜆₁(𝑡)=𝑡 and 𝜆₂(𝑡)=𝑡².

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𝐴=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝−000000⋮.....00000.0111111−000000.....00000.0111111−0000...00000.01112−00⋯⋯0..000.011300⋯..0001100⋯..001100..001100..0011−00.08−00.08⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(4.10)𝐵=⎛⎜⎜⎜⎜⎜⎜⎝−0⋮0000.08−00000.0800..000001100..000001100..000001100..000001100..0000011−000..0000.0112−00000....0000.011111−00000....0000.011111⎞⎟⎟⎟⎟⎟⎟⎠.(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.

(a)

(b)

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.

(a)

(b)

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𝑢^(𝑡)=110110𝑢,𝑢(0)=𝑢₀=11.(4.12)

In the following, we discuss different splitting ideas.

Version 1
We split the matrix as 110110=9101+1019,(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.

(a)

(b)

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.

(a)

(b)

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110110=9109+1011.(4.14) The Figure 7 presents the numerical errors between the exact and the numerical solution.

(a)

(b)

(c)

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10⁴1110₄=10⁴1−110₄0−1+1011.(4.15) The Figure 8 present the numerical errors between the exact and the numerical solution.

(a)

(b)

(c)

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.

(a)

(b)

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,im−𝜆₁𝑐₁,inΩ×[0,𝑡_𝑡],𝑐₂+∇𝜕⋅𝐅𝑐₂=𝑔−𝑐₂+𝐅𝑐=𝐯_2,im−+𝐷𝜆∇₁,_𝑡𝑐₁−𝜆₂𝑐₂,inΩ×[0,𝑡],𝜕𝑐_1,im=𝑔𝑐₁−𝑐_1,im−𝜆₁𝑐_1,im,inΩ×[0,𝑡_𝑡],𝑐_2,im𝑐=𝑔𝑐₂−𝑐_2,im+𝜆₁𝑐_1,im−𝜆₂𝑐_2,im,inΩ×[0,𝑡₁],𝑐(𝐱,𝑡)=𝑐_1,0(𝐱),𝑐₂(𝐱,𝑡)=𝑐_2,0(𝐱),onΩ,₁(𝐱,𝑡)=𝑐𝑐_1,1(𝐱,𝑡),𝑐₂(𝐱,𝑡)=𝑐_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:𝜕_𝑡𝐂=⎛⎜⎜⎜⎝𝐴−Λ¹Λ−𝐺0𝐺0₁𝐺𝐴−Λ0₂−𝐺−Λ₁0𝐺0𝐺Λ₁−𝐺0−Λ₂−𝐺⎞⎟⎟⎟⎠𝐂,(4.17) where 𝐂=(𝐜₁,𝐜₂,𝐜_1im,𝐜_2im)^𝑇, while 𝐜₁=(𝐜_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𝐴=𝐷Δ𝑥²+𝑣⋅⎛⎜⎜⎜⎜⎜⎝−12⋱−12⋱11⋱−12−12⎞⎟⎟⎟⎟⎟⎠Δ𝑥⋅⎛⎜⎜⎜⎜⎜⎝−11⋱1⋱−1−111⎞⎟⎟⎟⎟⎟⎠∈ℝ^𝐼×𝐼,(4.18) where 𝐼 is the number of spatial points.Λ₁=⎛⎜⎜⎜⎜⎜⎝𝜆¹0𝜆0₁⋱⋱00⋱𝜆₁0𝜆0₁⎞⎟⎟⎟⎟⎟⎠∈ℝ_𝐼×𝐼Λ₂,=⎛⎜⎜⎜⎜⎜⎝𝜆²0𝜆0₂⋱⋱00⋱𝜆₂0𝜆0₂⎞⎟⎟⎟⎟⎟⎠∈ℝ_𝐼×𝐼𝐺=⎛⎜⎜⎜⎜⎜⎝𝑔0⋱𝑔0⋱00⋱𝑔0𝑔0⎞⎟⎟⎟⎟⎟⎠∈ℝ_𝐼×𝐼,.(4.19)

We decouple into the following matrices:𝐴𝐴₁₂=⎛⎜⎜⎜⎝𝐴000𝐴00000000000⎞⎟⎟⎟⎠∈ℝ^4𝐼×4𝐼,=⎛⎜⎜⎜⎝−Λ¹Λ000₁0−0Λ₂−0Λ₁000Λ₁0−Λ₂⎞⎟⎟⎟⎠∈ℝ_4𝐼×4𝐼𝐴₃,=⎛⎜⎜⎜⎝−𝐺00𝐺−𝐺00𝐺−𝐺00𝐺−𝐺00𝐺⎞⎟⎟⎟⎠∈ℝ^4𝐼×4𝐼.(4.20)

For the operator 𝐴₁ and 𝐴₂=𝐴₂+𝐴₃, we apply the iterative splitting method.

Based on the decomposition, operator 𝐴₁ is only tridiagonal and operator 𝐴₂ 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.

(a)

(b)

(c)

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.