2.1 Initial value problem

(1)

(2)

1 Introduction

Ordinary and partial differential equations appear in many areas of science and engineering as mathematical models of a large variety of natural phenomena. For example in fluid dynamics, air pollution transport or reaction kinetics the main objective of mathematical research is to find thesolution of the given system of equations. In many cases the evolution of a quantity in time is sought, where at a given point of time the quantity can be a real number, a set of numbers or—in the most abstract form—a function of other variables, such as a distribution in space. Even in the simplest case to find the symbolic solution of the given differential equation is impossible, therefore methods that can provide an approximate solution are in the focus of many scientific work. A large variety of numerical methods have been developed that perform satisfactorily for different types of problems, however no universal numerical scheme exists that would generate an approximation of the required accuracy in reasonable time. Thus to handle a complex system containing terms of different characteristics one needs additional tools.

Operator splitting is a family of methods that are widely used to solve large systems of differential equations that describe the evolution of the given variables in time. The splitting procedure replaces the original complex problem with a set of subproblemsthat are simpler from the aspect of numerical treatment. For the subproblems different, suitable numerical methods can be used, thus approximation of prescribed accuracy can be produced in reasonable time for each subprob- lem. Then a solution for the original problem can be derived from the ones of the subproblems.

Operator splitting was proposed the first time in 1957: Bagrinovskii and Godunov [3] applied the method for partial differential equations. Yanenko, [55] used it in the solution of the heat equation and Dyakonov [17] for a larger class of linear equations. Since then many scientific work have been devoted to examine the method from different aspects: Marchuk [46] and Strang [53], Csom´os et al.

[15], Farag´o et al [24] developed different version of splitting; Bjorhus [8] gave general results for linear problems with unbounded operators; [22] studied the convergence of the method. In practice the subproblems are solved numerically, thus it is crucial to analyze the interaction between splitting and numerical schemes in thesecombined methods, see e. g. [14] and [4]. Applications of operator splitting methods can be found—among many others —in [1], [2], [30], [45] for reaction–diffusion problems; and in [57] to the more complex problem of modeling of air pollution transport.

An important field of applications are reaction–diffusion problems. Reaction–diffusion equations arise in many area of natural sciences. They are important in chemistry, for example in the description of pattern–formation; in meteorology, the modeling of air pollution; in the modeling of reactors in chemical engineering; in chemical industry: the description of technological processes;

in mechanical engineering design: in combustion processes and flame propagation; in biochemistry:

in the study of metabolic processes; in water-quality modeling. Since in these cases the number of unknown functions is large and the number of spatial dimension can be two or three, therefore one needs efficient approximation methods to solve these systems of partial differential equation.

There is a natural way to apply operator splitting to reaction–diffusion equations. The original problem can be split into two subproblems corresponding to the two physical processes modeled by

(5)

the equation: reaction and diffusion. Thus one obtains a problem that is usually a system of ordinary differential equations, a kinetic differential equation; and a system of partial differential equations that describes the diffusion. These subproblems have different properties, different methods can be used to solve them efficiently.

One possible approach—which is commonly followed in the literature—to study reaction–

diffusion equations is to apply a spatial discretization on the problem, thereby replacing a system of partial differential equations by a system of ordinary differential equations. Although the obtained problem has substantially different characteristics, we can fruitfully use the results from its study for the original continuous problem.

As in the case of every approximation method the purpose of operator splitting is to generate an approximate solution of the given problem with a prescribed accuracy. The convergence of the method is a fundamental question and lies in the focus of many research papers. A typical scenario is to show the consistency of the method, then to find conditions under which the method is stable, finally based on the basic theorem ,,consistency+stability=convergence” the convergence of the method is proven.

In two main parts of this thesis results on the classical splittings and on iterative splittings are presented. In both cases after the discussion of linear problems the nonlinear problems are investigated and last the results for the combination of splitting with numerical methods are presented.

In every part the results are illustrated with examples, the numerical results strongly support the presented theoretical results of each part.

In Chapter 2 the basic definitions are given that are discussed and investigated in this work. I present the concept of operator semigroups and collect the propositions that are the basic tools for the investigations in the last large part of my work, Section 4.2. The basic concept of operator splitting is introduced at the end of the chapter.

The convergence of operator splitting is proven in Sections 3.2 and 4.2 for nonlinear problems. A scenario different from the one mentioned above (consistency+stability=convergence) is followed here: I give explicit estimates of the global error, thus convergence is proven in a direct manner. In Section 3.2 operator splitting is applied on a very general problem, namely the acting operator is only locally Lipschitz-continuous. This is an extension of the existing results in the literature, where at least multiple differentiability of the operator is assumed. Beside the strong restriction of multiple differentiability the evaluation of error bounds—given in the literature—is rather complicated in practice, in many cases it is impossible without knowing the exact solution of the given problem.

The explicit error estimates presented in Section 3.2 allows quantitative error evaluations, thus to choose a suitable time step to satisfy the prescribed accuracy. This section is based on the work [34].

In Sections 3.3.2 and 4.1 the local order of splitting is investigated, which is a significant step towards proving consistency. Section 3.3.2 contains the results of [33] on the local error of combined methods applied for semi-linear problems. Similar calculations can be found in [29], although the nonlinearity is very specific: a linear operator shifted with a known function. In [33] I consider a sum of a linear and a nonlinear operator of a much more general type, I calculate the local order of 12 (three types of splittings combined with four different numerical schemes) combined methods.

(6)

I Section 4.1 I present the possible extensions of iterative splitting applied to the case when the original problem is split into more then two subproblems. The local order is calculated for a large family of methods for bounded linear problems.

Section 4.2 discusses the waveform relaxation method. It can be introduced as a natural modifi- cation of iterative splitting, however it evolved independently from that. The method is applied to a class of semi-linear problems which contains reaction–diffusion and reaction–advection equations.

Beside the proof of convergence I give explicit error estimates, consider the windowing-technique and discuss the combination of the method with numerical schemes applied to solve the subproblems. These results are published in [35].

Beside the wide area of continuous problems the results of Section 3.2 and 3.3.2 can be applied to discretized reaction–diffusion problems. The main strength of the results of Section 4.2 is that they can be applied directly on the continuous problem, however the application to discretized problems is also presented in details.

(7)

2 Basic definitions

2.1 Initial value problem

Definition 1. LetXbe a Banach space,D⊂Xa nonempty open set,F:D⊂X →Xandu₀∈D(F).

Consider the autonomous initial value problem

u⁰(t) =F(u(t)), u(0) =u₀. (1)

In this work the problem is assumed to have solution on a finite time interval, t ∈[0,T], with someT ∈R⁺.IfX =Rⁿ then (1) formulates a system of ordinary differential equations. One can describe a system of partial differential equations withX =C^r(R^d,Rⁿ).An important example is the family of reaction–diffusion equationsU_t(t,x,y,z) =∆U(t,x,y,z) + f(U(t,x,y,z)), whereX = C²(R³,Rⁿ), u(t):(x,y,z)7→U(t,x,y,z) and f represents the chemical reactions, being usually a polynomial or a rational function.

Definition 2. Theclassical solutionof (1) is the continuously differentiable functionu, for which u(t)∈D(F),

h→0lim

u(t+h)−u(t)

h −F(u(t))

=0, for every t∈(0,T) and the limit fort=0 is to be understood as lim

h→0⁺

u(t+h)−u(t)

h .

The classical result ([12], [44], Picard [51] and Lindel¨of [43]) for ordinary differential equations says that ifF is a locally Lipschitz continuous function, then (1) has a unique solution.

Remark 1. IfF is a bounded linear operator onX, then the solution of (1) ise^tFu₀:=

∞ j=0

∑

t^j j!F^ju₀. In many practical applications the central issue is to find the solution of (1) or—due to the difficulty of this task—to provide a satisfactorily accurate approximate solution for it. Operator splitting is a family of methods that are used in the process to generate approximation solutions for (1). The main results of this thesis are related to the consistency and convergence of splitting methods applied to problem (1). The following definitions contain the basic concepts used in the analysis of approximation techniques.

Definition 3. A time-discretization method M is a family of operators M(τ):X →X,where τ ∈ (0,T]. An approximate solution u_N obtained by the time discretization method M of (1) at time t∈(0,T]is defined by

u_N(t) = (M(t/N))^N(u₀), for N∈N.

(8)

Definition 4. A time-discretization methodMisconsistentif for the solutionuof (1)

τ→0lim

u(t+τ)−M(τ)(u(t)) τ

=0 uniformly fort∈[0,T−τ].

Definition 5. Theorder of consistencyof a consistent time-discretization method M is the largest p∈Nfor which there exist positive constantscandτ^?∈(0,T]such that for the solutionuof (1)

sup

t∈[0,T−τ]

u(t+τ)−M(τ)(u(t)) τ

6cτ^p holds for everyτ ∈(0,τ^?].

In the literature of operator splitting thelocal splitting errorand thelocal orderare often examined in order to have an a priori guess concerning the consistency.

Definition 6. Thelocal errorof a time discretization method Misku(τ)−M(τ)(u₀)k,whereuis the solution of (1). Thelocal orderof the method is the largestq∈Nfor which there exist positive constantscandτ^?∈(0,T]such that for everyτ ∈(0,τ^?]

u(τ)−M(τ)(u₀) τ

6cτ^q holds for everyτ ∈(0,τ^?].

Consequence 1. Consider a consistent splitting method (as a time-discretization method) of order p and let1>τ∈(0,T],then for the order of the local splitting error q one has q> p.

Thus results on the order of the local order provides an upper bound on the order of consistency, furthermore the techniques of the analysis of the local order are often helpful in the study of consistency. I will study this property of splitting methods in Sections 3.3.2 and 4.1.

Definition 7. A time-discretization methodMisstableif there exist constantsKandτ^?such that k(M(τ))^Nk6K

holds for all suchN∈Nandτ ∈(0,τ^?]that satisfyNτ6T.

Definition 8. Letube the solution of (1). Theglobal errorof a time discretization method is E_N(t):=ku(t)−u_N(t)k=ku(t)−(M(t/N))^N(u₀)k, N∈N.

(9)

Definition 9. A time discretization method is calledconvergenton[0,T]if

N→∞lim E_N(t) =0 for everyt∈[0,T].

Definition 10. The order of convergence of a time discretization method is the largest p∈N for which

N→∞lim N^p−1E_N(t) =0 for everyt∈[0,T].

Remark 2. In Section 4.2 the approximate solutions are constructed in a different manner—without time discretization —, thus convergence is to be understood as uniform convergence there, that is,

m→∞lim max

t∈(0,T]

ku(t)−u_m(t)k=0.

2.2 Operator semigroups

An important set of problems of the form (1) can be investigated successfully with the concept of operator semigroup. Here only those basic definitions and relations are collected that are directly used in the results of this thesis. Operator semigroups are the primary topic of the works [18] and [50].

Definition 11. The set {S(t)}_t>₀ of bounded linear operators on a Banach space X is called a strongly continuous semigroupor aC₀semigroup if:

1. S(0) =I,

2. S(t+s) =S(t)S(s), for allt,s>0,

3. for everyx∈X the mapst7→S(t)xare continuous fromR⁺intoX.

A conventional notation for the semigroup is(S(t))_t_>₀.

Proposition 1. For every strongly continuous semigroup(S(t))_t_>₀there exist constants M>1and ω ∈Rsuch that for every x∈X

kS(t)xk6Me^ωtkxk for all t>0.

Definition 12. Let(S(t))_t_>₀be a strongly continuous semigroup on a Banach spaceX and letD(A) be the subspace ofX defined by

D(A):=

x∈X : lim

τ→0

S(τ)x−x τ exists

.

(10)

The operatorA:D(A)⊆X →X,

Ax:= lim

τ→0

S(τ)x−x τ

is called the generator of the semigroup(S(t))_t_>₀ or, in other words, A generates the semigroup (S(t))_t_>₀.

Proposition 2. If A generates a strongly continuous semigroup(S(t))_t_>₀then 1. A:D(A)⊆X→X is a linear operator.

2. If x∈D(A)then S(t)x∈D(A)and d

dtS(t)x=S(T)Ax=AS(t)x for every t>0.

Proposition 3. If A generates a strongly continuous semigroup(S(t))_t_>₀then for every x∈D(A) S(t)x−x=

Z _t

0

S(s)Axds.

Example 1. A bounded linear operatorA withD(A) =X generates the strongly continuous semi- groupe^tA:=

∞

∑

j=0

t^j

j!A^j.In fact it is a uniformly continuous semigroup.

Example 2. Consider the first-order differential operator on C₀(R^d)corresponding to the continuously differentiable vector fieldv:R^d→R^d

A f :=

d

∑

i=1

v_i∂if

for f ∈C¹_c(R^d):={f ∈C¹(R^d): f has compact support}andx∈R^d.The operatorAgenerates a strongly continuous semigroup.

Example 3. Consider the second-order differential operator A f :=

d

∑

i=1

∂_i²f

defined for every f in the Schwarz spaceS(R^d)—the space of infinitely many times differentiable, rapidly decreasing functions—, see [18, Section VI. 5.]. The operator Agenerates a strongly continuous semigroup.

Proposition 4. If F generates a strongly continuous semigroup(S(t))_t>₀, then the solution of (1)is u(t) =S(t)u₀,for every t∈(0,T].

(11)

2.3 Operator splitting

The basic idea of splitting is to decompose the operator on the right hand side of (1) into a sum of suboperators, and formulate a set of subproblems corresponding to the suboperators. In some sense the suboperators are simpler and the subproblems are easier to solve numerically or even symbolically in some cases. The most simple type of operator splitting is the sequential splitting.

The method of sequential splitting is as follows. LetN∈N,Nτ:={1, ...,N},τ:=T/Nand consider the time division{0,τ,2τ, ...,Nτ} ⊂[0,T].Suppose thatF is the sum ofAandB,approximate the solution of (1) by solving the following pair of subproblems fort∈[0,τ]

v⁰₁(t) =A(v₁(t)), v₁(0) =u₀, w⁰₁(t) =B(w₁(t)), w₁(0) =v₁(τ).

Let the approximation of u(τ) be defined as u^(N)_spl(τ):=w₁(τ). Next solve the subproblems again witht∈[τ,2τ]

v⁰₂(t) =A(v₂(t)), v₂(τ) =u^(N)_spl(τ), w⁰₂(t) =B(w₂(t)), w₂(τ) =v₂(2τ).

We define the approximation of u(2τ) as u^(N)_spl(2τ):=w₂(2τ), and so on. Four different splitting schemes will be introduced in the next chapter, the lower index ”spl” will be used as a general notation and it will change accordingly to the examined splitting scheme.

Definition 13. The method defined by solving the subproblems

v⁰_i(t) =A(v_i(t)), v_i((i−1)τ) =u^(N)_seq((i−1)τ), t∈[(i−1)τ,iτ]

w⁰_i(t) =B(w_i(t)), w_i((i−1)τ) =v_i(iτ), t∈[(i−1)τ,iτ] (2) u^(N)_seq(0):=u₀, u^(N)_seq(iτ):=w_i(iτ)

successively fori∈Nτ is calledsequential splitting, referred to as ”seq” in the followings.

Sequential splitting can be defined reversely as well

w⁰_i(t) =A(w_i(t)), w_i((i−1)τ) =u^(N)_seq((i−1)τ), t∈[(i−1)τ,iτ]

v⁰_i(t) =B(v_i(t)), v_i((i−1)τ) =w_i(iτ), t∈[(i−1)τ,iτ] (3) u^(N)_seq(0):=u₀, u^(N)_seq(iτ):=v_i(iτ), i∈Nτ.

The procedure provides a sequence of values{u^(N)_seq(0):=u₀,u^(N)_seq(τ), ...,u^(N)_seq(Nτ)}inX.Thus the sequential splitting is a time-discretization method which provides an approximation for the solution of (1) on the discrete set {0,τ, ...,Nτ}. The convergence of the method at an arbitrary time point t^?∈[0,T]can be investigated by defining the algorithm (2) on[0,t^?],with the time stepτ=t^?/N.

(12)

Example 4. IfA andB generate theC₀ semigroup(Q(t))_t_>₀ and(S(t))_t_>₀,then the approximate solution obtained by the sequential splitting (2) is

u^(N)_seq(iτ) = (S(τ)Q(τ))ⁱu₀, for everyi∈Nτ, and the splitting (3) gives

u^(N)_seq(iτ) = (Q(τ)S(τ))ⁱu₀, for everyi∈Nτ.

(13)

3 Classical splittings

Another splitting scheme, named after Marchuk and Strang who developed the new method independently, was first presented in the papers [46] and [53]. Consider (1) withF =A+B.

Definition 14. The algorithm defined as

v⁰_i(t) =A(v_i(t)), v_i((i−1)τ) =u^(N)_MS((i−1)τ), t ∈[(i−1)τ,(i−1/2)τ],

w⁰_i(t) =B(w_i(t)), w_i((i−1)τ) =v_i((i−1/2)τ), t ∈[(i−1)τ,iτ], (4) z⁰_i(t) =A(z_i(t)), z_i((i−1/2)τ) =w_i(iτ), t ∈[(i−1/2)τ,iτ],

u^(N)_MS(0) =u₀, u^(N)_MS(iτ) =z_i(iτ), i∈Nτ

is calledMarchuk–Strang type splitting.

The method ofsymmetrically weighted splittingwas first mentioned in [52]. The first theoretical investigation of the method was done in [15], and later it was succesfully applied—among others—

in air pollution modeling [9].

Definition 15. The method ofsymmetrically weighted splittingis defined as

v⁰_i(t) =A(v_i(t)), v_i((i−1)τ) =u^(N)_sw ((i−1)τ), t∈[(i−1)τ,iτ]

w⁰_i(t) =B(w_i(t)), w_i((i−1)τ) =v_i(iτ), t∈[(i−1)τ,iτ]

y⁰_i(t) =B(y_i(t)), y_i((i−1)τ) =u^(N)_sw ((i−1)τ), t∈[(i−1)τ,iτ]

z⁰_i(t) =A(z_i(t)), z_i((i−1)τ) =y_i(iτ), t∈[(i−1)τ,iτ] (5) u^(N)_sw (0):=u₀ u^(N)_sw (iτ) = ¹₂w_i(iτ) +¹₂z_i(iτ), i∈Nτ.

Example 5. IfAandBgenerate theC₀semigroup(Q(t))_t_>₀and(S(t))_t_>₀,then one time step with the Marchuk–Strang type splitting and the symmetrically weighted splitting can be written as

u^(N)_MS(iτ) =Q(τ/2)S(τ)Q(τ/2)u^(N)_MS((i−1)τ), and

u^(N)_sw (iτ) = 1

2Q(τ)S(τ) +1

2S(τ)Q(τ)

u^(N)_sw ((i−1)τ, for everyi∈Nτ

respectively.

The sequential, the Marchuk–Strang type and the symmetrically weighted splittings are called classical splittings. Observe that for these methods to be correctly defined the sets D(A)andD(B) need to obey some requirements, which have to be considered in the investigations. Because of the structural similarity theadditive splittingis also discussed in this chapter. The method was first introduced and tested in [27]. It was further examined for bounded linear problems in [24].

(14)

Definition 16. The method defined as

v⁰_i(t) =A(v_i(t)), v_i((i−1)τ) =u^(N)_as ((i−1)τ), t∈[(i−1)τ,iτ], y⁰_i(t) =B(y_i(t)), y_i((i−1)τ) =u^(N)_as ((i−1)τ), t∈[(i−1)τ,iτ], (6)

u^(N)_as (0):=u₀ u^(N)_as (iτ):=v_i(iτ) +y_i(iτ)−u^(N)_as ((i−1)τ), i∈Nτ

to approximate the solution of (1) is calledadditive splitting.

Example 6. In the case ofAandBbeing generators of theC₀ semigroups(Q(t))_t_>₀and(S(t))_t_>₀ one step with additive splitting can be written as

u^(N)_as (iτ) = (Q(τ) +S(τ)−I)u^(N)_as ((i−1)τ), for everyi∈Nτ, whereI:X →X is the identity operator.

3.1 Linear problems

The main results on operator splitting for linear problems is summarized briefly in this section. Con- sider (1) and suppose thatF is linear and it can be decomposed asF=A+B,whereAandBare linear, as well. A fundamental work is [8], where the splitting is examined for a more general problem, namely the right hand side of (1) is extended with an inhomogeneous part,u⁰(t) =F(u(t)) +f(t).

Here the results for the case f ≡0 are recalled.

Theorem 1(Corollary 4.4 in Bjørhus [8]). Let F (resp. A,B) be an infinitesimal generator of a C₀ semigroup(U(t))_t_>₀(resp. (Q(t))_t_>₀,(S(t))_t_>₀). Let D(A) =D(B) =D(F)and D(A²) =D(B²) = D(F²)be satisfied and let T >0.Then we have the approximation property

sup

06t6T−τ

kQ(τ)S(τ)u(t)−U(τ)u(t)k6τ²C(T), 06τ 6T whenever u₀∈D(A²),where C(T)is a constant independent ofτ.

This theorem means that the sequential splitting is consistent of order one, under the given conditions.

In [23] the argumentation of [8] is followed and it is shown that the Marchuk–Strang type and the symmetrically weighted splittings are consistent of order two. They requireD_k:=D(A^k)∩D(B^k)∩ D((A+B)^k) to be dense inX fork=1,2,3 and A^k|_D_k, B^k|_D_k, (A+B)^k|_D_k be closed operators for k=1,2,3.

Theorem 2(Thm. 3.4 and 4.2 in Farag´o and Havasi [23]). There exist constants C₁(T)and C₂(T) independently ofτsuch that

kQ(τ/2)S(τ)Q(τ/2)u(t)−U(τ)u(t)k6τ³C₁(T),

(15)

1

2Q(τ)S(τ) +1

2S(τ)Q(τ)

u(t)−U(τ)u(t)

6τ³C₂(T), for all u₀∈D,where D=

3

\

k=1

D_k.

Consequence 2. If A and B are bounded linear operators on X then the conditions of the above theorems are fulfilled, furthermore stability is also guaranteed since

k(S(τ)Q(τ))ⁿk=k(e^τBe^τA)ⁿk6e^nτkAke^nτkBk6e^T^kAke^T^kBk, k(Q(τ/2)S(τ)Q(τ/2))ⁿk6e^T^/2kAke^T^kBke^T^/2kAk

and

Q(τ)S(τ) +S(τ)Q(τ) 2

n

6kQ(τ)kⁿkS(τ)kⁿ6e^T^(kAk+kBk).

Consequently for bounded linear problems the classical splitting methods are convergent.

Theorem 3 (Theorem 2.1. in [24]). Let A and B be bounded linear operators on some Banach space X and consider the Cauchy problem(1). Then the local order of additive splitting(6)is1.

The work [4] is also to be mentioned here as a detailed study on the convergence of classical splittings for linear problems.

3.2 Nonlinear problems

Nonlinear problems are usually discussed in the context of Lie-operator formalism—e.q. [38], [22], [16], [32], [28]—, that allows brief formulation of the splitting error. Nevertheless, an actual representation of the Lie-operator in practice is complicated. It is more important that this technique demands the suboperators to be infinitely many times differentiable and the error estimations are expressed in terms of higher order derivatives. Let us mention just three examples, where a different approach is followed: in [5] the Schr¨odinger-equation is investigated with the sum of the differential operatori∆and a globally Lipschitz functionF, for whichF(0) =0 applies and its derivatives are bounded up to order four; in [53] nonlinear problems of a quite general form are shortly discussed assuming multiple differentiability of the operators acting on the right hand side; [41] also assumes multiple differentiability and uses series expansion for its analysis.

Beside the strong restriction of multiple differentiability the evaluation of error bounds—given in the literature—is rather complicated in practice, in many cases it is impossible without knowing the exact solution of the given problem.

The content of this section can be found in the work [34]. Operator splitting is analyzed in the general framework of Banach spaces and the acting operators are only assumed to be locally Lipschitz-continuous only. The convergence of the classical and additive splittings is proven in this very general case. The area of validity of the existing results for nonlinear problems is extended, on

(16)

the other hand the error bounds are expressed by the suboperators acting on the initial value, that is, in terms that are relatively easy to calculate in practice. Based on the results of [34] having two approximations of the unknown solution an a posteriori estimate can be calculated for the Lipschitz- constant. From that the global error can be estimated and a suitable time step can be chosen, thus one can ensure the approximation to fulfill the required accuracy. In the case of sequential splitting, based on the error estimates we propose an optimal decomposition which generates the lowest splitting error. This result is more of a heuristical nature and a deeper analysis would be necessary to clarify this question, although numerical tests on two relatively simple problems support our conjecture.

Operator splitting is frequently used in modeling large scale air pollution transport [57], among many other areas of applications, such as water quality modeling [20] or reaction kinetics [37]. The results of the present work form the general theoretical base for the application of operator splitting in a large class of problems.

Consider (1) with F = f+g, where D(f),D(g)⊂X and suppose that D:=D(f)∩D(g) is a nonempty open connected set. Let f :D(f)−→X and g:D(g)−→X be continuous operators.

Thus (1) can be written as

u⁰(t) = (f+g)(u(t)), u(t₀) =u₀, (7) withu₀∈D,t₀∈I,I⊂Rbeing an open interval.

Remark 3. If f andgare Lipschitz-continuous operators onD, —that is, there are constantsL_f,L_g such that for anyw,z∈Dthe relations

kf(w)−f(z)k6L_fkw−zk and kg(w)−g(z)k6L_gkw−zk

hold,— then (7) has a unique solution in D. In practice to demand the global Lipschitz property to hold on Dis too strict, therefore beside the existence of the solution of (7) we only assume the operators f andgto be locally Lipschitz continuous.

The following analysis is built on two assumptions.

1. Lett₀=0 and suppose that (7) has a unique solution[0,T]3t7→u(t),with someT >0.

2. Consider the following neighborhood of the unique solutionuof (7) S^r := ^[

t∈[0,T]

Sr(t), where Sr(t):={η :kη−u(t)k6r}, t∈[0,T],

whereris a given positive real number. Suppose that there existsLsuch that for anyt∈[0,T] and anyw,z∈Sr(t)the relations kf(w)− f(z)k6Lkw−zk andkg(w)−g(z)k6Lkw−zk hold.

In practice the constantL depends onras it will be illustrated in the numerical examples. Here L is chosen to the fixed value ofr. Lalso depends on f andg, although for our theoretical results no

(17)

real advantage is taken of this fact,Lcan be considered as max{L_f,L_g}. It is important to note that given this assumption the initial value problems

v⁰(s) = f(v(s)), v(0) =v₀ and w⁰(s) =g(w(s)), w(0) =w₀

have a unique solution inSr(t)for anyv₀,w₀∈Sr(t).The solutions of the sub-problems in (2) and (6) satisfy the relations

v_i(t) =u^(N)_seq((i−1)τ) + Z _t

0

f(v_i(s))ds, w_i(t) =v_i(τ) + Z _t

0

g(w_i(s))ds and

y_i(t) =u^(N)_as ((i−1)τ) + Z t

0

g(y_i(s))ds, v_i(t) =u^(N)_as ((i−1)τ) + Z t

0

f(v_i(s))ds, respectively for everyi∈Nτ.Furthermore the solution of (7) satisfies

u(t) =u₀+ Z _t

0

(f+g)(u(s))ds for every t∈[0,T]. (8) Note that the operator f+gis locally Lipschitz continuous in everySr(t)with the Lipschitz constant 2L.

In the following sections the convergence of splittings (2) and (6) is proven in a direct way, namely, for everyi∈Nτ an upper bound for the global error will be formulated as a function ofτ.

Then the limitN→∞can easily be obtained fromτ→0.Given the properties of f andgas above the results derived for (2) can be obtained in a similar way for (3) as well.

3.2.1 Preliminary estimations

The main challenge of the splitting procedure is to remain in the sufficiently small neighborhood of the real solution, that is, to ensureu^(N)_spl(iτ)∈Sr(iτ). The error estimations—as functions of the time step τ—are expressed under the assumption that every term lies inside S^r. However the obtained error bounds lead out of this set. The following lemma forms the base for us to show the existence of a thresholdτ^?,such that, the error bounds are valid for every time stepτ<τ^?.On the other hand it allows a constructive error estimation, that is, an appropriate time stepτ can be calculated.

Lemma 1. Consider the continuous functions a :R⁺₀ →R and b: R⁺₀ →R, where b is strictly monotone increasing and unbounded. Assume that a(t₀)6b(t₀)<K with some K∈Rand t₀∈R⁺₀, furthermore for all t>t₀such that a(t)6K we have a(t)6b(t).Then

a(t)6K for every t6b⁻¹(K).

Obviously a(t)6K holds for anyt >t₀that is close enough tot₀, the purpose of this lemma is to show that it holds on[t₀,b⁻¹(K)].

(18)

Proof. Indirect: suppose that ˆt<b⁻¹(K),where

tˆ:=max{t_K :t_K6b⁻¹(K)anda(t)6K,∀t∈[t₀,t_K]}.

Thena(ˆt)<Ksince

a(t)ˆ 6b(t)ˆ <b(b⁻¹(K)) =K.

Sinceais continuous there exist ˜t ∈(ˆt,b⁻¹(K))such thata(˜t) = a(ˆt) +K

2 <K.Thus we arrived at a contradiction in the definition of ˆt,consequently ˆt is not less thenb⁻¹(K).

Lemma 2. Suppose that h:D(h)−→X is Lipschitz-continuous on B_ρ(x₀):={η :kη−x₀k6ρ}

with the Lipschitz constant L_ρ.Consider the initial value problem x⁰(t) =h(x(t)), x(0) =x₀.As far as x(t)∈B_ρ(x₀),for t>0we have

kx(t)−x₀k6e^L^ρ^t−1

L_ρ kh(x₀)k.

Proof.

kx(t)−x₀k=kx₀+ Z _t

0

h(x(s))ds−x₀k6 Z _t

0

kh(x(s))−h(x₀)k+kh(x₀)kds6 6L_ρ

Z _t

0

kx(s)−x₀kds+tkh(x₀)k

is an integral-inequality for the functiont7→ kx(t)−x₀k,for which Gronwall’s lemma gives kx(t)−x₀k6

Z _t

0

skh(x₀)kL_ρe^L^ρ^(t−s)ds+tkh(x₀)k.

By calculating the integral and some simplifications we get the statement.

Consequence 3. Apply Lemma 1 with t₀:=0, a(t):=kx(t)−x₀k, b(t):= e^L^ρ^t−1

L_ρ kh(x₀)k and K:=ρ.Then for any06t6b⁻¹(ρ) =_L¹

ρln^ρL^ρ_kh(x^+kh(x⁰^)k

0)k the relationkx(t)−x₀k6e^L^ρ^t−1

L_ρ kh(x₀)k holds.

Consequence 4. Applying Lemma 2 for(7)with h:= f+g,x₀:=u(t),ρ:=r,B_ρ(x₀):=Sr(t)and L_ρ :=L,then one gets

ku(t+s)−u(t)k6 e^2Ls−1

2L k(f+g)(u(t))k,

(19)

for any s∈[t,t+t_r], with t_r:= 1

2Llnr2L+k(f+g)(u(t))k

k(f+g)(u(t))k . Considering(2)and(6) kv_i(t)−v_i(0)k6e^Lt−1

L kf(u^(N)_spl((i−1)τ))k, for any t ∈[0,t_v]and i∈Nτ,with t_v:= 1

LlnrL+kf(u^(N)_spl((i−1)τ))k kf(u^(N)_spl((i−1)τ))k

. Lemma 3. Consider the solution u of (7). Then for any t∈[0,T]

k(f+g)(u(t))k6e^2Ltk(f+g)(u₀)k and

kg(u(t))k6e^2Lt−1

2 k(f+g)(u₀)k+kg(u₀)k.

Proof. Chooset such thatSr(0)∩Sr(t)6= /0 andu(t)∈/Sr(0).Furthermore let ˜t<t be such that u(t)˜ ∈Sr(0)∩Sr(t).Then

k(f+g)(u(t))k6k(f+g)(u(t))−(f+g)(u(˜t))k+k(f+g)(u(˜t))k6 62Lku(t)−u(t)k˜ +k(f+g)(u(t))k.˜

Applying Consequence 4 one obtains k(f+g)(u(t))k62Le^2L(t−˜^t)−1

2L k(f+g)(u(˜t))k+k(f+g)(u(˜t))k6e^2L(t−˜^t)k(f+g)(u(˜t))k.

With similar stepsk(f+g)(u(˜t))k6e^2L˜^tk(f+g)(u₀)kcan be obtained, thus k(f+g)(u(t))k6e^2L(t−^t)^˜ k(f+g)(u(t˜))k6e^2Ltk(f+g)(u₀)k.

The second statement can be derived analogously

kg(u(t))k6kg(u(t))−g(u(t))k˜ +kg(u(t))k˜ 6Lku(t)−u(t)k˜ +kg(u(t))k˜ 6 6Le^2L(t−˜^t)−1

2L k(f+g)(u(˜t))k+kg(u(˜t))k6 6 e^2L(t−˜^t)−1

2 e^2L^t^˜k(f+g)(u₀)k+kg(u(t))˜ −g(u₀)k+kg(u₀)k6 6e^2L(t−^t)^˜ −1

2 e^2L˜^tk(f+g)(u₀)k+Lku(˜t)−u₀k+kg(u₀)k6

(20)

u₀

S

_r

H0L

uHt L uHtL

S

_r

HtL

r

Figure 1: Illustration of Lemma 3.

6 e^2L(t−^t)^˜ −1

2 e^2L˜^tk(f+g)(u₀)k+e^2L^t^˜−1

2 k(f+g)(u₀)k+kg(u₀)k=

= e^2Lt−1

2 k(f+g)(u₀)k+kg(u₀)k.

For any larger timetthe extension can be done by repeating the above argumentation. See Figure 1

for illustration.

At last recall that the remainder of the Taylor series of the exponential function can be expressed as

e^at−

n−1 i=0

∑

(at)ⁱ

i! = (at)ⁿ n! e^a^t^ˆ, where ˆt∈[0,t].Thus e^at−1

a =te^aˆ^t6te^at. 3.2.2 Convergence of sequential splitting

First a threshold for the time step of splitting is defined for which it can be shown that the following error estimates are valid. After that the convergence of splitting follows.

(21)

Definition 17. Considert^?∈(0,T]and the function G_seq(t):=e^Lt e^3Lt Lt²

2

e^2Lt^?−1

e^2Lt−1kg(u₀)k+Lt

4 t^?e^2Lt^?−te^2Lt^?−1 e^2Lt−1

!

k(f+g)(u₀)k

! +

+e^LtLt² 2

e^2Lt^?−1

2 k(f+g)(u₀)k+kg(u₀)k

!

+t e^2Lt^?−1

2 k(f+g)(u₀)k+kf(u₀)k+kg(u₀)k

!!

fort >0.As it can be shown to be strictly monotone increasing one can introduce τ^?:=G⁻¹_seq(r), obviouslyG_seq(t)6rfor everyt∈[0,τ^?].

It will turn out that an appropriate time step of sequential splitting is any numberτ 6τ^?. Theorem 4. Consider u, the solution of (48)and the sequential splitting (2), let0<τ 6τ^? such that t^?/τ=:N∈N.Then

E_seq^(N)(iτ)6e^2Lτ Lτ²

2

e^2Liτ−1

e^2Lτ−1kg(u₀)k+Lτ 4

iτe^2Liτ−τe^2Liτ−1 e^2Lτ−1

k(f+g)(u₀)k

, (9)

for every i∈ {0} ∪Nτ.

Proof. By definitionE_seq^(N)(0) =0.Suppose that the statement holds forE_seq^(N)((i−1)τ).Consider ku((i−1)τ+t)−w_i(t)k=

=

u((i−1)τ) + Z t

0

(f+g)(u((i−1)τ+s))ds−v_i(τ)− Z t

0

g(w_i(s))ds

=

u((i−1)τ) + Z _t

0

(f+g)(u((i−1)τ+s))ds−u^(N)_seq((i−1)τ)− Z _τ

0

f(v_i(s))ds− Z _t

0

g(w_i(s))ds 6 6E_seq^(N)((i−1)τ) +

Z t 0

kf(u((i−1)τ+s))− f(v_i(s))kds+ Z _τ

t

kf(v_i(s))kds+

+ Z _t

0

kg(u((i−1)τ+s))−g(w_i(s))kds6 6E_seq^(N)((i−1)τ) +L

Z _t

0

ku((i−1)τ+s)−v_i(s)kds+

+ Z _τ

t

kf(v_i(s))kds+ Z _t

0

Lku((i−1)τ+s)−w_i(s)kds. (10) The last estimation requires v_i(s) andw_i(s)to be in Sr(s), that is, ku(s)−v_i(s)k6r andku(s)− w_i(s)k6rfors∈[(i−1)τ,t].Next we show that these requirements are fulfilled.

(22)

The first integrand ku(s)−v_i(s)k in (10) is a continuous function of s and E_spl^(N)((i−1)τ)) = ku((i−1)τ)−v_i((i−1)τ)k<r, therefore for s∈[(i−1)τ,τ] that is close enough to (i−1)τ we haveku(s)−v_i(s)k6r.For simplicity let us change the variablestot,then

ku(t)−v_i(t)k=

u((i−1)τ) + Z _t

(i−1)τ(f+g)(u(s))ds−u^(N)_seq((i−1)τ)− Z _t

(i−1)τ

f(v_i(s))ds 6 6E_seq^(N)((i−1)τ) +

Z _t

(i−1)τkf(u(s))−f(v_i(s))kds+ Z _t

(i−1)τkg(u(s))kds.

Considering the Lipschitz-property and applying Lemma 3 for the last integrand one gets ku(t)−v_i(t)k6E_seq^(N)((i−1)τ) +

Z _t

(i−1)τ

Lku(s)−v_i(s)kds+

+ Z _t

(i−1)τ

e^2Ls−1

2 k(f+g)(u₀)k+kg(u₀)kds6 6E_seq^(N)((i−1)τ) +

Z t (i−1)τ

Lku(s)−v_i(s)kds+ (t−(i−1)τ)

e^2Lt−1

.

Applying Gronwall’s inequality one gets

ku(t)−v_i(t)k6 (11)

6eL(t−(i−1)τ)

E_seq^(N)((i−1)τ) + (t−(i−1)τ)

e^2Lt−1

6 6eL(t−(i−1)τ) e^2Lτ Lτ²

2

e^2L(i−1)τ−1

e^2Lτ−1 kg(u₀)k+

+Lτ

4 (i−1)τe^2L(i−1)τ−τe^2L(i−1)τ−1 e^2Lτ−1

!

k(f+g)(u₀)k

! +

+(t−(i−1)τ)

e^2Lt−1

.

Let us define the strictly monotone increasing function b_v(t):=eL(t−(i−1)τ) e^2Lτ Lτ²

2

e^2Lt^?−1

e^2Lτ−1kg(u₀)k+Lτ

4 t^?e^2Lt^?−τe^2Lt^?−1 e^2Lτ−1

!

k(f+g)(u₀)k

! +

+(t−(i−1)τ) e^2Lt^?−1

!!

,

(23)

then ku(t)−v_i(t)k 6b_v(t) holds for t close enough to (i−1)τ. Lemma 1 can be applied with t₀:= (i−1)τ,a(t):=ku(t)−v_i(t)k,b(t):=b_v(t)andK:=r.Fromb_v(iτ)<G_seq(τ)6rit follows thatku(t)−v_i(t)k6b_v(t)holds on[(i−1)τ,b⁻¹_v (r)], consequently it holds on[(i−1)τ,iτ]for every i∈Nτ.

Since in (10) ku(s)−w_i(s)k is a continuous function ofs and ku((i−1)τ)−w_i((i−1)τ)k= ku((i−1)τ)−v_i(iτ)k6b_v(iτ)<r,therefore fors∈[(i−1)τ,iτ]that is close enough to(i−1)τwe haveku(s)−w_i(s)k6r,thus the Lipschitz-condition holds and (10) is valid. Applying Gronwall’s inequality for (10) one gets

ku(t)−w_i(t)k6 6eL(t−(i−1)τ)

E_seq^(N)((i−1)τ) +L Z _t

(i−1)τku(s)−v_i(s)kds+ Z _iτ

t

kf(v_i(s))kds

. (12)

The last term can be estimated as Z iτ

t

kf(v_i(s))kds6 Z iτ

t

kf(v_i(s))−f(u(s))k+kf(u(s))kds6 6

Z _iτ

t

Lkv_i(s)−u(s)kds+ Z _iτ

t

kf(u(s))kds, thus (12) becomes

ku(t)−w_i(t)k6eL(t−(i−1)τ)

E_seq^(N)((i−1)τ) +L Z _iτ

(i−1)τku(s)−v_i(s)kds+ Z _iτ

t

kf(u(s))kds

. (13) Integrate (11) to get the second term in the parenthesis

L Z _iτ

(i−1)τku(s)−v_i(s)kds6 6L

Z _iτ

(i−1)τ

eL(s−(i−1)τ)

E_seq^(N)((i−1)τ) + (s−(i−1)τ)

e^2Ls−1

ds6 6(e^Lτ−1)E_seq^(N)((i−1)τ) +L

Z _iτ

(i−1)τ

e^Lτ(s−(i−1)τ)

e^2Liτ−1

ds=

= (e^Lτ−1)E_seq^(N)((i−1)τ) +e^LτLτ² 2

e^2Liτ−1

. (14)

Via Lemma 3 the third term in the parenthesis becomes Z iτ

t

kf(u(s))kds6 Z iτ

t

e^2Ls−1

2 k(f+g)(u₀)k+kf(u₀)kds6 6(iτ−t)

e^2Liτ−1

2 k(f+g)(u₀)k+kf(u₀)k

.

(24)

Thus (13) becomes

ku(t)−w_i(t)k6eL(t−(i−1)τ)

e^LτE_seq^(N)((i−1)τ) +e^LτLτ² 2

e^2Liτ−1

+

+(iτ−t)

e^2Liτ−1

,

substitutingE_seq^(N)((i−1)τ)and considering thatiτ−t6τ fort∈[(i−1)τ,iτ]one gets ku(t)−w_i(t)k6eL(t−(i−1)τ) e^3Lτ Lτ²

2

e^2L(i−1)τ−1

+Lτ

4 (i−1)τe^2L(i−1)τ−τe^2L(i−1)τ−1 e^2Lτ−1

!

k(f+g)(u₀)k

! + +e^LτLτ²

2

e^2Liτ−1

+τ

e^2Liτ−1

.

Let us define the strictly monotone increasing function b_w(t):=eL(t−(i−1)τ) e^3Lτ Lτ²

2

e^2Lt^?−1

4 t^?e^2Lt^?−τ

e^2Lt^?−1 e^2Lτ−1

!

k(f+g)(u₀)k

! +

+e^LτLτ² 2

e^2Lt^?−1

!

+τ e^2Lt^?−1

!!

,

thenku(t)−w_i(t)k6b_w(t),sinceiτ6t^?and ie^2Liτ−e^2Liτ−1

e^2Lτ−1 = (e^2Lτ−1)

i k=1

∑

ke^2L(k−1)τ

is increasing with i. Lemma 1 can be applied witht₀ := (i−1)τ, a(t):=ku(t)−w_i(t)k, b(t):=

b_w(t)andK :=r.Sinceb_w(iτ)<G_seq(τ)6r the estimationku(t)−w_i(t)k6b_w(t)holds on[(i− 1)τ,b⁻¹_i (r)], consequently it holds on[(i−1)τ,iτ]for everyi∈Nτ.

Finally forE_seq^(N)(iτ)start from relation (12), E_seq^(N)(iτ) =ku(iτ)−w_i(iτ)k6e^Lτ

E_seq^(N)((i−1)τ) +L Z _iτ

(i−1)τku(s)−v_i(s)kds

.

Using (14) one gets E_seq^(N)(iτ)6e^2Lτ

E_seq^(N)((i−1)τ) +Lτ² 2

e^2Liτ−1

2 k(f +g)(u₀)k+kg(u₀)k

.

(25)

Consider the assumption forE_seq^(N)((i−1)τ) E_seq^(N)(iτ)6e^2Lτ Lτ²

2 e^2Lτe^2L(i−1)τ−1

+Lτ²

4 e^2Lτ (i−1)e^2L(i−1)τ−e^2L(i−1)τ−1 e^2Lτ−1

!

k(f+g)(u₀)k+

+Lτ² 2

e^2Liτ−1

2 k(f+g)(u₀)k+Lτ²

2 kg(u₀)k

,

the first and the last term can be merged E_seq^(N)(iτ)6e^2Lτ

Lτ² 2

e^2Liτ−1

e^2Lτ−1kg(u₀)k+

+Lτ² 4

(i−1)e^2Liτ−e^2Liτ−e^2Lτ

e^2Lτ−1 +e^2Liτ−1

k(f+g)(u₀)k

=

=e^2Lτ Lτ²

2

e^2Liτ−1

e^2Lτ−1kg(u₀)k+Lτ² 4

ie^2Liτ−e^2Liτ−1 e^2Lτ−1

k(f+g)(u₀)k

,

which is the desired relation (9). Having E_seq^(N)(iτ)<G_seq(τ)6r it follows by induction that (9)

holds for everyi∈Nτ.

Theorem 5. Operator splitting(2)is convergent at every t^?∈[0,T]of order1.

Proof. Since the right side of (9) is an increasing function ofithe global error satisfies E_seq^(N)(t^?)6e^2Lτ Lτ²

2

e^2Lt^?−1

4 t^?e^2Lt^?−τ

e^2Lt^?−1 e^2Lτ−1

!

k(f+g)(u₀)k

! .

The right side is an increasing function oft^?thus E_seq^(N)(t^?)6e^2LτLτ

4

2τe^2LT−1

e^2Lτ−1kg(u₀)k+

Te^2LT−τe^2LT−1 e^2Lτ−1

k(f+g)(u₀)k

for everyt^?∈[0,T].The limit of the right side is zero as τ →0, which implies that it is zero as N→∞. One can easily check that the limit is not zero when it is multiplied byN,since

τ→0lim τ

e^2Lτ−1 = 1 2L.

(26)

Remark 4. The form of the global error (9) suggests that the decomposition where g(u₀) = 0 generates the smallest error, since the second term(f+g)(u₀)does not change. Thus an optimal splitting can be defined. Although, (9) is only an upper estimation for the error, thus this observation is more of a heuristical nature. This question will be examined later in numerical tests.

Consider an arbitrary decomposition of the right hand side of (48): f+g,then the decomposition f˜+g˜ with ˜f(v):= f(v) +g(u₀) and ˜g(w):=g(w)−g(u₀) has the above property. Refined error estimates can be derived with L_f and L_g, taking the dependence of L on the suboperators into account. A simple shift of the suboperators withg(u₀)however does not change the corresponding Lipschitz-constants.

3.2.3 Convergence of additive splitting

To show the convergence of additive splitting we follow the same scenario as in the previous section.

Definition 18. Considert^?∈(0,T]and the strictly monotone increasing function G_as(t):=e^Lt e^3LtLt²

2

e^2Lt^?−1

e^2Lt−1 +te^2Lt^?+1 2

!

(kf(u₀)k+kg(u₀)k) fort>0.Let us defineτ^?:=G⁻¹_as (r).

Theorem 6. Consider u, the solution of (7)and the splitting(6). Then forτ,such that t^?/τ=:N∈N and0<τ 6τ^?,

E_as^(N)(iτ)6e^3LτLτ² 2

e^2Liτ−1

e^2Lτ−1(kf(u₀)k+kg(u₀)k), (15) for every i∈ {0} ∪Nτ.

Proof. By definition E_as^(N)(0) =0.Suppose that the statement holds for E_as^(N)((i−1)τ).Then E_as^(N)((i−1)τ)<rand consider

E_as^(N)(iτ) =ku(iτ)−(v_i(iτ) +y_i(iτ)−u^(N)_as ((i−1)τ))k=

=

u((i−1)τ) + Z _iτ

(i−1)τ

(f+g)(u(s))ds−u^(N)_as ((i−1)τ)− Z _iτ

(i−1)τ

f(v_i(s))ds−

−u^(N)_as ((i−1)τ)− Z _iτ

(i−1)τg(y_i(s))ds+u^(N)_as ((i−1)τ) 6 6E_as^(N)((i−1)τ) +

Z _iτ

(i−1)τ

kf(u(s))−f(v_i(s))kds+ Z _iτ

(i−1)τ

kg(u(s))−g(y_i(s))kds6 6E_as^(N)((i−1)τ) +L

Z _iτ

(i−1)τku(s)−v_i(s)kds+L Z _iτ

(i−1)τku(s)−y_i(s)kds. (16)

2.1 Initial value problem

Contents

1 Introduction

2 Basic definitions

2.1 Initial value problem

∑

2.2 Operator semigroups

∑

∑

∑

2.3 Operator splitting

3 Classical splittings

3.1 Linear problems

3.2 Nonlinear problems

S

H0L

S

HtL

∑

∑