Iterative Operator-splitting Methods with Embedded Discretization Schemes

In this paper we describe a computation of iterative operator-splitting method, which are known as competitive splitting methods, see [10] and [11]. We derived a closed form, based on commutators for the iterative method. The time discretization schemes apply extrapolation schemes and Pade approximations to the exp-functions. Spatial discretization schemes considered Lax-Wendroff methods and are combined with the iterative schemes. The error analysis describe the approximation errors. Numerical examples of ordinary and partial differential equations support the fast computation ideas.


Introduction
In this paper we concentrate on approximation to the solution of the linear evolution equation where L, A and B are unbounded operators and T ∈ R + .We solve our equation with the following numerical scheme: and c 0 (t n ) = c n , c −1 = 0.0, with c i+1 (t n ) = c n , where c(t n ) is the approximation at t = t n .
The numerical method (2)-( 3) can be written algebraic in a 2-stage iterative splitting scheme: c i (t) = exp(At)c(0) + t 0 exp(As)Bc i−1 ds, (4) c i+1 (t) = exp(B t)c(0) + t 0 exp(Bs)Ac i ds, (5) where i = 1, 3, 5, . . .and the initial or start solution is given as c 0 (t) = 0 or constant.Further we have the conditions, that c n is the known split approximation at the time-level t = t n .The split approximation at the time-level t = t n+1 is defined as c n+1 = c 2m+1 (t n+1 ).(Clearly, the function c i+1 (t) depends on the interval [t n , t n+1 ], too, but, for the sake of simplicity, in our notation we omit the dependence on n.) Based on our motivation to design effective algorithms for large equation systems.The problem arose in the field of optimizing the computation of the iteration steps of very large systems of differential equations fixed on time-scale and on one discretization method.Here novel ideas in computing e x p-functions can be used, so we apply the multi-product expansion as an extrapolation method, see [8].The iterative splitting method can be formulated as a generalization of a waveform relaxation approach, see [5].So that also numerical schemes of sparse solvers can be taken into account to such a method.For partial differential equations, the spatial discretization methods are very important and we concentrate on efficient embedded discretization schemes, that overcome the problem of order restriction, while using lower order spatial discretization schemes, see [12].
Historically, effective computational methods can derived by considering the local character of each equation part.So in the last years the ideas of splitting into simpler equations are established, see [16], [7] and [14].We concentrate on choosing extrapolation to obtain higher order schemes without loosing efficiency in computing the operators.
The outline of the paper is as follows.The operator-splitting-method is introduced and the error-analysis of the operator-splitting method is presented in Section 2. A closed form is discussed in Section 4, where we discuss an efficient computation of the iterative splitting method with based on extrapolation methods.In Section 6 we present the numerical results for the methods.Finally we discuss future works in the area of iterative methods.

Error analysis
The following algorithm is based on the iteration with fixed-splitting discretization step-size τ, namely, on the time-interval [t n , t n+1 ] we solve the following sub-problems consecutively for i = 0, 2, . . ., 2m. (cf.[16]): and c 0 (t n ) = c n , c −1 = 0.0, with where c n is the known split approximation at the time-level t = t n .The split approximation at the time-level t = t n+1 is defined as c n+1 = c 2m+1 (t n+1 ).(Clearly, the function c i+1 (t) depends on the interval [t n , t n+1 ], too, but, for the sake of simplicity, in our notation we omit the dependence on n.)

Two-step iterative schemes for operators
In the following we discuss the consistency of the 2 stage iterative method, taken into account to iterate over both operators.

Theorem 1. Let us consider the abstract Cauchy problem in a Banach space
where A, B : X → X are given linear operators which are generators of the analytical semigroups and c 0 ∈ X is a given element.We assume dom(B) ⊂ dom(A), so we are restricted to balance the operators.Further, we assume the estimations of an unbounded operator, see [15]: where α ∈ (0, 1) and we assume B 1−α is the infinitesimal generator of an analytical semigroup for all α ∈ (0, 1), see [4].
Further we assume: where α, β , p, q ∈ (0, 1) and The error of the first time-step is of accuracy (τ i A α n ), where τ n = t n+1 − t n and we have equidistant time-steps, with n = 1, . . ., N .Further i A are the iterative steps with operator A.
For the first iterations we have: and for the second iteration we have: In general we have: For the odd iterations: i = 2m + 1 for m = 0, 1, 2, . . .
We have the following solutions for the iterative scheme: The solutions for the first two equations are given by the variation of constants: For the recursive even and odd iterations we have the solutions: For the odd iterations: i = 2m + 1 for m = 0, 1, 2, . . .
For e 2 we have: We obtain: For odd and even iterations, the recursive proof is given in the following.In the next steps, we shift t n → 0 and t n+1 → τ n for simpler calculations, see [15].The initial conditions are given with c(0) = c(t n ).
For the odd iterations means the iteration over operator A: i = 2m + 1, with m = 0, 1, 2, . .., we obtain for c i and c: By shifting 0 → t n and τ n → t n+1 , we obtain our result: where α = min i j=1 {α i } and 0 ≤ α i < 1 and i A is the number of odd iteration steps means over the operator A.
The same proof idea can be applied to the even iterative scheme.
Where for the even scheme we could not obtain a higher order, see: and we have t n exp(B(t n+1 − s 1 ))A. . .
So we could not improve the order in the weaker iterations, so we need at least some strong iterative steps.

Remark 1. An example can be assumed with
We apply two iterative steps with operator A and have the following local errors: and hence where 0 ≤ α < 1 and C, C are constants independent of τ n .

Convergence results
To derive convergent results, we apply stability and consistency done in the previous sections.
The convergence result for the assumption B = A 1−α and apply the iteration to operator A (one-side iteration).
We obtain the following results: where m are number of the local steps, C is a constant independent of τ n .
For i A = 2: where m are number of the local steps, C is a constant independent of τ n .
The same recursive argument can be done for arbitrary i A ∈ N + and we obtain where m are number of the local steps, C is a constant independent of τ n and i A are the number of iterative steps.

Remark 2. Local we have an error O(t i
This means we need sufficient iterative steps i A .At least for α ∈ (0, 1) and an assumed convergence (t c ) with order c > 0, we need In numerical experiments we obtain much more better results and achieve second or third order methods with two and three iterative steps.
In the next section we describe the computation of the integral formulation with exp-functions.

Computation of the iterative splitting schemes: Closed formulation
In the last years, the computational effort to compute integral with exp-function has increased, we present a closed form, and re-substitute the integral with closed functions.Such benefits accelerate the computation and made the ideas to parallelize, see [2] and [8].
Recursion.We study the stability of the linear system ( 6) and ( 7), based on different closed formulations.We consider the suitable vector norm • on R M , together with its induced operator norm.The matrix exponential of Z ∈ R M ×M is denoted by exp(Z).We assume that: where K A , K B ∈ R + are given as the growth estimation of the exponential functions, see [5].
It can be shown that the system (1) implies exp(τ n (A + B))c n ≤ K c n and is itself stable.
For more transparency of the splitting scheme ( 6) and ( 7), we consider a wellconditioned system of eigenvectors whereby we can consider the eigenvalues λ 1 of A and λ 2 of B instead of the operators A, B themselves.
We assume that all initial values c i (t n ) = c approx (t n ) with i = 0, 1, 2, . . ., are as where m is the order, see [5].
Further we assume λ 1 = λ 2 , otherwise we do not consider the iterative splitting method, while the time-scales are equal, see [7].
A(α)-stability.We define z k = τλ k , k = 1, 2. We start with c 0 (t) = u n and we obtain: where S m is the stability function of the scheme with m-iterations.
Let us consider the A(α)-stability given by the following eigenvalues in a wedge: For the A-stability we have The stability of the splitting schemes are given in the following theorems with respect to A and A(α)-stability.

Computation of the iterative splitting methods: Closed formulation with integral computations
A further computation of the iterative schemes are given by the variation of constants, see for exponential splitting schemes [15].
To obtain analytical solutions of the differential equations: . . .
where c(t n ) is the initial condition and A, B are bounded operators, the initialization is with c 1,iter (t) = exp(B t) exp(At)c(t n ) is a first order splitting scheme.
The application of the variation of constants is given as: We apply the numerical integration of the integral with Trapezoidal rule for the first integral and Simpson's rule for the second integral and obtain: where we compute c 1,iter (t n+1 ), c 2,iter (t n+1 ), . . ., and s The forth order method can also be computed with Bode's or Romberg's rules: Example.
where we have to compute the subinterval results: We we have to compute t ∈ [0, T ], with t 0 , t 1 , . . ., t N and N number of time steps, where the time steps are equidistant of τ = t j − t j−1 , j = 1, . . ., N .
The generalization is given with Romberg's extrapolation scheme is given in the following algorithm.
(i) We start with n = 0 and the initial condition c(0), and starting solution We compute the time interval t n , t n+1 and the solution c(t n+1 ) is obtained by: (a) We start with i = 2 where and We compute the integrals of the functions f A,i−1 , f B,i by: where ) we increase i = i + 1, till i = I and we go to (iii) (iii) The result is given as c(t n+1 ) = c I (t n+1 ), we increase n = n + 1 and goto (ii), if n = N we are finished.
Remark 3. The same recurrent argument can be applied to the next iterative scheme.A higher numerical integration method is necessary.Here we have only to apply matrix multiplications and can skip the time-consuming integral computations.Only two evaluations for the exponential function for A and B are necessary.The main disadvantage of computing the iterative scheme exactly are the time-consuming inverse matrices.These can be skipped with numerical methods.
We have the following assumptions for the stability formulations: The stability of the methods are given in the following Theorem 2 Theorem 2. We have the following stability for the integral formulated iterative schemes: For the stability function S i,iter of iterative splitting schemes we have with ω ∈ [0, 1] and the initial conditions are c(t n ) = c n and i is the iteration index.
Proof.We proof the stability of S 1,iter .
For the extrapolation schemes we have the stability function: For both possibilities z 1 → −∞ and z 2 → −∞ we have S 1,iter → 0. For the higher iteration steps we taken into account the assumptions (59)-(60).
Based on this assumptions, we write for i = 2 For both possibilities z 1 , z 2 → −∞ we have S 2,iter → 0. Same proof idea is used for the higher iterative steps.

Spatial discretization schemes
To apply our method to partial differential equations, we have to consider higher order spatial discretization schemes.
In the following we discuss an efficient scheme based on Lax-Wendroff.

The Lax-Wendroff scheme in one dimension
The 2nd order Lax-Wendroff scheme for the 1D advection diffusion equation is: where ν = v∆t ∆x is the Courant-Friedrich-Levy number and Since the Lax-Wendroff scheme is made for hyperbolic PDEs only we have to examine the stability closely and carry out a von-Neumann stability analysis.We encounter the following stability condition for the Fourier coefficients of u j (t n ) Obviously this statement is always true for right-hand sides ≥ 1 so it is sufficient to estimate the solutions of In Figure 5.1 we see the pairs (ν, z) for which the statement holds.Of course only the region {(ν, z), ν ≥ 0, z ≥ 1} is of practical interest.The figure suggests that for all ν ≥ 0 there are z ≥ 1 leading to a stable scheme.For the small D(z, v, ∆t) for which this is true we may say that the problem is convection dominated.

Generalization to two dimensions
The advection-diffusion equation in two dimensions: Application of the Lax-Wendroff scheme yields the second order accurate formula: We use central differences to achieve a second order scheme.This is similar to the 1D scheme (77) except for the new cross term incorporating ∂ x y u: introducing the Courant-Friedrich-Levy numbers and the constants z x , z y This finite difference scheme is now brought into conservation form * : At this point we introduced the flux limiter φ(θ i ) which is used to handle steep gradients of u(x, y, t) where the Lax-Wendroff scheme adds spurious oscillations (see also [6]).θ i = is a measure for the slope of u and is estimated in the x and in the y direction respectively.For our purposes we chose the van Leer limiter:

Dimensional splitting
During the numerical experiments we will compare the discussed Lax-Wendroff scheme to the dimensional splitting which is also translated into a second order finite difference scheme, filtering out the numerical viscosity . In this case we use an implicit BDF2 method to achieve the second order in time yielding: with the spatial discretization We suppressed the dependencies u i j (t n ) on the right hand side

Advection-diffusion splitting
Another way of splitting the advection-diffusion equation is the following: This splitting will also be used during the experiments.The implementation of the operators is very similar to the above finite differencing schemes.However numerical aspects demand the use of an explicit method namely the Adams-Bashforth method.For further discussions on the practicability of finite differencing schemes in conjunction with the iterative splitting see [10].

Numerical experiments
In the following we present numerical experiments with the closed computable iterative splitting methods and their benefits to standard schemes.

First experiment
We deal with the 2-dimensional advection-diffusion equation and periodic boundary conditions with the parameters The given advection-diffusion problem has an analytical solution which we will use as a convenient initial function: We apply dimensional splitting to our problem where We use a 1st order upwind scheme for ∂ ∂ x and a 2nd order central difference scheme for 2 we achieve a 2nd order finite difference scheme because the new diffusion constant eliminates the first order error (i.e. the numerical viscosity) of the Taylor expansion of the upwind scheme.L y u is derived in the same way.
We apply a BDF5 method to gain 5th order accuracy in time: Our aim is to compare the iterative splitting method with AB-splitting.Since [A x , A y ] = 0 there is no splitting-error for the AB-splitting and therefore we cannot expect to achieve better results with the iterative splitting in terms of general numerical accuracy.Instead we will show that the iterative splitting out competes AB-splitting regarding the computational effort and round-off-errors.But first there are some remarks which have to be made concerning the special behavior of both methods when combined with high-order Runge-Kutta and BDF methods.

Splitting and schemes of high order in time.
Concerning AB-Splitting.The principle of AB-splitting is well known and simple.The equation du d t = Au + Bu is broken up into which are connected via u n+1 (t) = u n+1/2 (t + ∆t).This is pointed out in figure (6.2).AB-splitting works very well for any given one-step method like the Crank-Nicholson-Scheme.Not taking into account the splitting-error (which is an error in time) it is also compatible with high order schemes such as explicit/implicit Runge-Kutta-schemes.
Things look different if one tries to use a multi-step method like the implicit BDF or the explicit Adams method with AB-splitting, these cannot be properly applied as is shown by the following example: Choose for instance a BDF2 method which, in case of du/d t = f (u), has the scheme So the first step of the AB-splitting looks like: Clearly u n+1/2 (t) = u n (t) but what is u n+1/2 (t − ∆t)?This is also shown in figure (6.2) and it is obvious that we wont have knowledge about u n+1/2 (t − ∆t) unless we compute it separately which means additional computational effort.This overhead even increases dramatically when we move to a multi-step method of higher order.The mentioned problems with the AB-splitting will not occur with a higher order Runge-Kutta method since only knowledge of u n (t) is needed.

Remarks about the iterative splitting.
The BDF methods apply very well to the iterative splitting.Let us recall at this point that this method, although being a real splitting scheme, always remains a combination of the operators A and B so no steps have to be done into one direction only ‡ .
In particular we do a subdivision of our given time-discretization t j = t 0 + j∆t into I parts.So we have subintervals t j,i = t j + i∆t/I , 0 ≤ i ≤ I on which we solve ‡ As we will see there is an exception to this. the following equations iteratively: u −1/I is either 0 or a reasonable approximation § while u 0 = u(t j ) and u 1 = u(t j +∆t).The crucial point here is that we only know our approximations at given times which don't happen to be the times at which a Runge-Kutta method needs to know them.Therefore, in case of a RK method, the values of the approximations have to be interpolated with at least the accuracy one wishes to attain with the splitting and this means a lot of additional computational effort.We may summarize our results now in Table 6.1 that shows which methods are practicable for each kind of splitting scheme.¶ Numerical results.After resolving the technical aspects of this issue we can now proceed to the actual computations.The question which arises is which of the splitting methods has the least computational effort since we can expect them to solve the problem with more or less the same accuracy if we use practicable methods with equal order because [A x , B x ] = 0. We tested the dimensional splitting of the 2d-advection-diffusion equation with the AB-splitting combined with a 5th order RK method after Dormand and Prince and with the iterative splitting in conjunction with a BDF5 scheme.We used 40 × 40-and 80 × 80-grids and completed n t time-steps with each of which subdivided into 10 smaller steps until we reached time t end = 0.6 which is sufficient to see the main effects.The iterative splitting was done with 2 iterations which was already enough to attain the desired order.In Tables 7.2 and 6.3 the errors at time t end and the computation times are shown.§ In fact the order of the approximation is not of much importance if we fulfill a sufficient number of iterations.In case of u −1/I = 0 we have the exception that a step in A-direction is done while B is left out.The error of this step vanishes after a few but mostly only one iteration ¶ In favor of the iterative splitting scheme take also into the account that AB-splitting may be used along with the mentioned high order methods but cannot maintain the order if [A, B] = 0 while the iterative splitting re-establishes the maximum order of the scheme when a sufficient number of iterations is done.As we can see, the error of the iterative splitting reaches the AB-splitting error after a certain number of time-steps and stays below it for all additional steps we accomplish.Of course the error cannot sink under a certain amount which is governed by the spatial discretization.It is to be noticed that while the computation time used for the iterative splitting is always about 20%-40% less than that of the AB-splitting the accuracy is, with a sufficient number of timesteps, slightly better than that of the AB-splitting.This is due to the roundoff error which is higher for the Runge-Kutta method because of the greater amount of basic operations needed to compute the RK steps.
A future task will be to introduce non-commuting operators in order to show the superiority of the iterative splitting over the AB-splitting when the order in time is reduced due to the splitting error.

Second experiment
In the second experiment we consider the linear and nonlinear advectiondiffusion equation in two dimensions.

Two dimensional advection-diffusion equation. We deal with the two dimensional advection-diffusion equation with the parameters
The code for both methods is kept in the simplest possible form.

One and two step iterative schemes for advection-diffusion equations.
We use the splitting (96) and apply the one and two step iterative splitting respectively.
Results are shown in Tables 7.3 and 7.4.We see in case of the one step method that the choice of A and B makes a difference for 1 Iteration.This is due to the dominance of the advection part.In general there is no improvement with higher iteration numbers.In this case the operators are defined as follows: B = −u i−1 ∂ y + D∂ y y (105) u i−1 denotes the (i − 1)th iteration.This takes care of the coupling between the iteration steps.We use a finite differencing scheme of third order in space (central differences) and time (explicit Adams-Bashforth) and compare the results with a scheme of only second order.We also compare the iterative scheme with a standard Strang-Marchuk splitting of second order.Numerical results are shown in Table 7.6.Remark 4. In the second experiment, we can improve iterative splitting methods with embedded spatial discretization methods.Higher order schemes as Lax-Wendroff schemes are embedded as dimensional splitting schemes to our iterative splitting method.Also nonlinear equations can be embedded.To compute the iterative scheme, we have only to deal with spatial discretization schemes on each dimensions.At least less iterative steps are needed to achieve higher order results.

Conclusions and Discussions
We have presented an iterative operator-splitting method computed with extrapolation schemes.We have analyzed the splitting error for the operators.
Under weak assumptions we could proof the higher order error bounds.Closed formulations allow to compute the delicate iterative scheme efficient.We embedded spatial discretization schemes to the iterative splitting method and achieved higher order results.Numerical examples confirm the applications to differential equations and achieve the theoretical results.In the future we will focus us on the development of improved operator-splitting methods with respect to their application in timedependent and nonlinear differential equations.

Figure 5 . 1 .
Figure 5.1.(ν, z) inside the blue area satisfy the stability condition.The region of interest is also shown in detail.

Table 6 . 2 .
Errors and computation times of AB-splitting and iterative splitting for a 40 × 40-grid.Number of steps Error AB Error It. spl.AB computation time It.spl.computation time

Table 6 . 3 .
Errors and computation times of AB-splitting and iterative splitting for a 80 × 80-grid.

Table 7 . 3 .
Errors for different diffusion constants on a 40×40-grid with the one step iterative splitting with A = −v∇ and B = D∆.

Table 7 . 4 .
Errors for different diffusion constants on a 40×40-grid with the one step iterative splitting with A = D∆ and B = −v∇.

One and two step iterative schemes for 2D Burgers equation.
In our last experiment we apply the iterative splitting to the 2D Burgers equation:

Table 7 . 5 .
Errors for different diffusion constants on a 40×40-grid with the two step iterative splitting.

Table 7 . 6 .
Burgers equation on a 40×40-grid with the two step iterative splitting and schemes of second and third order (t 0 = 0.25, T = 30).