OVERCOMING ORDER REDUCTION IN DIFFUSION-REACTION SPLITTING. PART 2: OBLIQUE BOUNDARY CONDITIONS LUKAS EINKEMMER∗ AND ALEXANDER OSTERMANN∗ Abstract. Splittingmethodsconstituteawell-establishedclassofnumericalschemesforthetime integrationofpartialdifferentialequations. Theirmainadvantagesovermoretraditionalschemesare 6 computational efficiencyandsuperiorgeometricproperties. Inthepresenceofnon-trivialboundary 1 conditions, however, splittingmethods usuallysufferfromorderreductionandsomeadditional loss 0 of accuracy. For diffusion-reaction equations with inhomogeneous oblique boundary conditions, a 2 modification of the classic second-order Strang splitting is proposed that successfully resolves the problemoforderreduction. Thesamecorrectionalsoimprovestheaccuracyoftheclassicfirst-order n Lie splitting. The proposed modification only depends on the available boundary data and, in the a case ofnonDirichletboundary conditions, onthe computed numericalsolution. Consequently, this J modification can be implemented in an efficient way, which makes the modified splitting schemes 1 superior to their classic versions. The framework employed in our error analysis also allows us to 1 explain the fractional orders of convergence that are often encountered for classic Strang splitting. Numericalexperiments thatillustratethetheoryareprovided. ] A Key words. Strang splitting, diffusion-reaction equation, oblique boundary conditions, Neu- N mannboundaryconditions,Robinboundaryconditions,mixedboundaryconditions,orderreduction . h AMS subject classifications. 65M20,65M12,65L04 t a m 1. Introduction. For the numerical solution of evolution equations, splitting [ methods haveattractedmuchattentioninrecentyears. The mainadvantageofsplit- ting methods over more traditionaltime integrationschemes is the fact that the par- 1 tial vector fields in the splitting can usually be integratedmore efficiently, sometimes v much more efficiently, comparedto applying a numericalmethod to the full problem. 8 8 Moreover, splitting methods have, in general, better geometric properties [10]. As a 2 typical example, we mention splitting of the Schr¨odingerequation [6] into the kinetic 2 and the potential part. The former is integrated in frequency space using the fast 0 Fourier transform, whereas the latter can often be integrated exactly. This requires . 1 considerablylesscomputationaleffortcomparedtodirectlyapplying,forexample,an 0 implicit Runge–Kutta method to the Schr¨odinger equation. 6 In this paper, we will consider a splitting approach for diffusion-reaction equa- 1 : tions. Inourproblem,thediffusionismodelledbyalinearellipticdifferentialoperator v andthereactionbyanonlinearsmoothfunction. Suchasplittingintoalinearelliptic i X problem, which can be solved efficiently by fast Poisson techniques [17], and a non- r linear ordinary differential equation, which is simply solved pointwise in space, has a attracted much attention in the literature (see, e.g., the articles [1, 3, 5, 8, 11] and the monograph [14]). Although splitting methods are widely used in multiphysics problems, a rigorous error analysis is still missing for many important applications. For instance, the con- vergenceanalysisofsplittingmethodsisoftencarriedouteitherforperiodicboundary conditions or for full space problems. On the other hand, it is known from numerical experiments that non-trivial boundary conditions lead to order reduction (see, e.g., the experiments in [4, 12]). In order to avoidsuch an order reduction, the splitting must be modified in such a way that the internal steps of the splitting become compatible with the prescribed ∗Department of Mathematics, University of Innsbruck, Technikerstraße 13, Innsbruck, Austria ([email protected], [email protected]). 1 2 L.EINKEMMER,A.OSTERMANN boundary conditions. For reaction-diffusion equations with Dirichlet boundary con- ditions, we presented such a correction in our recent paper [4]. Here, we propose a modification of that idea which, on the one hand, is easier to implement, and on the other hand canbe extendedto more generaloblique boundary conditions. The latter comprise Dirchlet, Neumann, Robin, and mixed boundary conditions. To illustrate the underlying ideas of our approach,we consider for a moment the semilinear heat equation (1.1) ∂ u=∆u+f(u), αu+β∂ u| =b. t n ∂Ω Depending on the choice of α and β, it is equipped with Dirichlet, Neumann, Robin, or a mixture of these boundary conditions. The problem of order reduction in splitting methods has the following origin. When the right-hand side of the differential equation (1.1) is split into the two parts ∆u and f(u), the nonlinearity f(u) does not satisfy, in general, the homogeneous boundary conditions αu + β∂ u| = 0 with f(u) in place of u. However, these n ∂Ω boundary conditions are essential as they characterize the domain of the Laplacian ∆, which plays a paramount role in the error analysis. The noncompliance of the nonlinearity with the homogeneous boundary conditions is the main source of order reduction. We therefore look for a function q that satisfies the boundary conditions of the nonlinearity f(u), that is ′ (1.2) αq+β∂ q| =αf(u)+βf (u)∂ u| , n ∂Ω n ∂Ω and define the splitting ∂ v =∆v+q, αu+β∂ u| =b t n ∂Ω (1.3) ∂ w =f(w)−q t with the help of this correction. Then ∂ w satisfies the desired boundary conditions. t Note that q still depends on the value of u and ∂ u on the boundary. We can n eliminate ∂ u from the relation by using β∂ u| =b−αu| . This then yields n n ∂Ω ∂Ω ′ αq+β∂ q| =αf(u)+f (u)(b−αu)| . n ∂Ω ∂Ω It is instructive to consider a few special cases. For Dirichlet boundary conditions (α=1, β =0) we obtain the familiar relation q| =f(b) ∂Ω which, in the notation of our previous paper [4] constitutes two steps. First, we compute anextensionz ofthe boundary data,i.e. z| =b,andthen q =f(z)which ∂Ω (onthe boundary)yields the same result. Note thatq isindependent ofu andcanbe precomputed (for time independent boundary data b) in this case. For Neumann boundary conditions (α=0, β =1) we obtain ′ ∂ q| =f (u)b. n ∂Ω Note that in this case q depends on u| . However, we can substitute u by the ∂Ω numerical solution at the beginning of the time step and still obtain a second-order approximation, see section 4.2. OVERCOMINGORDERREDUCTIONINSPLITTINGMETHODSII 3 In the general case, the continuation q will thus be constraint by ′ (1.4) αq(t )+β∂ q(t )| =αf(u )+f (u )(b−αu )| , n n n ∂Ω n n n ∂Ω where u denotes the numerical solution at time t . n n Theoutlineofthepaperisasfollows. Insection2wepresentourmodelproblem, a semilinear parabolicequationsubject to oblique boundary conditions. The numeri- calschemeisexplainedinsection3. Withthehelpofacorrectionqthathastosatisfy acertainboundarycondition,wedefineaso-calledboundary-correctedsplitting. Sec- tion4addressestheerroranalysisofourmodifiedsplittingschemes. Theconvergence properties of the Lie splitting are analyzedin section 4.1, that of the Strang splitting in section 4.2. The main result in this context is Theorem 4.5 which proves second order convergence of the modified Strang splitting under natural assumptions on the data. The order reduction, which is observed for the classic Strang splitting, is fur- ther discussedandrigorouslyexplainedinsection4.3. Some representativenumerical results in one and two space dimensions are given in sections 5 and 6, respectively. Issues of efficient implementation are finally discussed in section 7. Weconcludethissectionbysomeremarksontheemployednotation. TheletterC alwaysdenotesagenericconstant. Itmaytakeondifferentnumericvaluesatdifferent occurrences. Most functions in our paper depend on time t and space x. Whenever convenient,however,wesuppressthesecondvariableinthenotationandsimplywrite g(t) instead of g(t,·). Further, we denote g (s) = g(t +s) for t = nτ, where τ is n n n the time step size of the numerical method. The very exception to this rule is the exactsolutionof the problem,denoted by u. Here, the symbol u alwaysdenotes the n numerical approximation to u(t) at time t=t . n 2. Model problem. The purpose of this section is to present a model problem for which a rigorous error analysis of Lie and Strang splitting will be provided. On the one hand, the problem shouldbe simple enoughso that all steps in the proof can be presentedin detail. On the other hand, it shouldbe sufficiently generalto contain interestingproblemclasses. Wechooseaclassofsemilinearparabolicproblemssubject to oblique boundary conditions. LetΩ⊂Rd be anopen,boundedsubsetwith smoothboundary∂Ω. We consider a second-order elliptic differential operator d d (2.1) D = d (x)∂ + d (x)∂ +d (x)I ij ij i i 0 X X i,j=1 i=1 with smooth coefficients. The matrix (d (x)) is assumed to be symmetric and uni- ij formly positive definite on Ω. Further, we consider the first-order boundary operator d (2.2) B = β (x)∂ +α(x)I i i X i=1 withsufficientlysmoothcoefficients. TheoperatorB isassumedtosatisfytheuniform non tangentiality condition d (2.3) inf β (x)n (x) >0. x∈∂Ω(cid:12) i=1 i i (cid:12) (cid:12)P (cid:12) (cid:12) (cid:12) Here, n(x) denotes the outer normal of Ω at x ∈ ∂Ω. We further assume that f is a sufficiently smooth real function. For more details concerning the functional analytic framework, we refer to [16, Sect. 3.1]. 4 L.EINKEMMER,A.OSTERMANN Henceforth, we consider the following semilinear parabolic problem with oblique boundary conditions (2.4a) ∂ u=Du+f(u) t (2.4b) Bu| =b ∂Ω (2.4c) u(0)=u 0 as the model problem in our analysis. The nonlinearity f, the inhomogeneity b and the initial data u are assumedto be sufficiently smooth. Depending onthe choiceof 0 coefficientsinB,differentboundaryconditionsaremodelled. Forexample,thechoice β =...=β =0 1 s gives a Dirichlet problem, whereas the choice d α=0, β (x)= d (x)n (x), 1≤i≤d i j=1 ij j P correspondstoaNeumannproblem. Notethat(2.3)issatisfiedinbothofthesecases. 3. Description of the numerical method. We next describe how to carry out one splitting step in the time interval [t ,t ] of length τ. As motivated in the n n+1 introduction, wechooseas correctiona smoothfunction q thatsatisfies the boundary conditions of f(u). Actually, it is sufficient to require (as we will see later) (3.1) Bq (0)| =Bf(u(t ))| +O(τ). n ∂Ω n ∂Ω Since our numerical methods converge at least with order one, we can simply take ′ (3.2) Bq | =αf(u )+f (u ) b −αu | , n ∂Ω n n n n ∂Ω (cid:0) (cid:1) as this choice then satisfies (3.1). With this correction q at hand, we consider the n boundary-correctedsplitting (3.3a) ∂ v =Dv +q , Bv | =b t n n n n ∂Ω n (3.3b) ∂ w =f(w )−q , t n n n and solve it on the time interval [t ,t ] by the standard Lie or Strang approach. n n+1 In the caseof Lie splitting, for a giveninitial value u , we first1 solve (3.3b) with n initial value w (0)=u to obtainw (τ). Then, we integrate (3.3a) with initial value n n n v (0)=w (τ) and finally define the numerical solution u at time t as n n n+1 n+1 (3.4) u =L u =v (τ). n+1 τ n n Henceforth, the numerical scheme (3.4) with operator L will be referred to as the τ modified Lie splitting step of size τ. In the case of Strang splitting, for a given initial value u , we first solve (3.3a) n with initial value v (0) = u to obtain v (τ). Next, we integrate (3.3b) with initial n n n 2 valuew (0)=v (τ)toobtainw (τ). Finally,weintegrateoncemore(3.3a),butthis n n 2 n time withinitialvalue v (0)=w (τ), anddefine the numericalsolutionu attime n n n+1 t as n+1 (3.5) u =S u =v (τ). n+1 τ n n 2 1Wenotethattheorderofthetwopartialflowscanalsobeinterchanged. Theresultingscheme canbeanalyzedinasimilarway,seealsoLemma4.6below. OVERCOMINGORDERREDUCTIONINSPLITTINGMETHODSII 5 Thenumericalscheme(3.5)withoperatorS willbereferredtoasthemodified Strang τ splitting step of size τ henceforth. Note that the classic Lie and Strang splitting schemes can formally also be an- alyzed in the above framework, simply by setting q = 0. We will make use of that n later on. For the choice q =0, however,condition (3.1) is not satisfied, in general. n 4. Error analysis. The purposeofthis sectionis to give athorougherroranal- ysis of Lie and Strang splitting applied to (2.4). Our analysis will be based on the framework of analytic semigroups, see [18, 13, 16]. For the purpose, we first have to transformthepartialdifferentialequation(2.4)tohomogeneousboundaryconditions. Note, however, that this transformation is only employed in the proof and will not affect our numerical methods. Let z denote a smooth function that fulfills the boundary conditions (4.1) Bz| =b, ∂Ω and let u=u−z. Then u satisfies the abstract initial value problem (4.2) b ∂ u=Au+b f(u+z)+Dz−∂ z, u(0)=u(0)−z(0), t t where A is the inbfinitesbimal gbenerator of an analytibc semigroup on the Banach space Lp(Ω). Its domain D(A)={ψ∈W2,p(Ω)|Bψ| =0}⊂Lp(Ω) ∂Ω incorporates the homogeneous boundary conditions of the problem. For functions ψ ∈D(A), the action of A is defined by Aψ =Dψ. For details, we again refer to [16, Sect. 3.1]. For later use we recall that analytic semigroups enjoy the parabolic smoothing property, i.e. (4.3) ketA(−A)γk≤Ct−γ, γ ≥0 for all t∈(0,T]. Without any loss of generality,we assumed here that A is invertible with a bounded inverse. This can always be achieved by an appropriate scaling of u. Let ϕ denote the entire functions, recursively defined by k ϕ (z)− 1 ϕ (z)= k k!, k ≥0, ϕ (z)=ez. k+1 0 z From the integral representation 1 sk−1 (4.4) ϕ (z)= e(1−s)z ds, k ≥1, k Z (k−1)! 0 we infer that the operators ϕ (tA) are uniformly bounded for 0 ≤ t ≤ T. Moreover, k the following relation holds (4.5) eτA−ϕ (τA)=τA ϕ (τA)−ϕ (τA) , 1 1 2 (cid:0) (cid:1) which will be employed later. We still have to formulate an appropriate framework to include the nonlinearity f. Forsimplicity, wewanttouse the same normforthe analysisof (3.3a) and(3.3b). 6 L.EINKEMMER,A.OSTERMANN This, however, implies that f must be considered as a smooth function on Lp(Ω), defined in a neighborhood of the exact solution with values in Lp(Ω). This seems to be quite a restrictive assumption. However, we can use the well-known fact that the norm is differentiable in Lp(Ω) \ {0} for 1 < p < ∞ and twice differentiable for 2 ≤ p < ∞, see [19, Sect. 2.2]. Therefore, f(u) can be replaced in the analysis by χ(kuk)f(u), where χ is a suitably chosen real smooth cut-off function. Without further mention, we will follow this approach. Still another possibility would be to consider f : V → Lp(Ω), where V denotes the fractional power space V = D((−A)γ) for some 0 ≤ γ < 1. This framework is used often in the literature, for example in [13]. For simplicity, however, we will not consider it further here. Theexactsolutionoftheabstractproblem(4.2)canbeexpressedbythevariation- of-constants formula t u(t)=eτAu(0)+ e(t−s)A f(u(s)+z(s))+Dz(s)−∂ z(s) ds. Z0 (cid:16) t (cid:17) b b b This allows us to express the exact solution of (2.4) at time t = t +τ in the n+1 n following way: u(t )=z (τ)+eτA u(t )−z (0) n+1 n n n (4.6) τ (cid:0) (cid:1) + e(τ−s)A f(u(t +s))+Dz (s)−∂ z (s) ds, Z0 (cid:16) n n t n (cid:17) wherez (s)isasmoothfunctionsatisfyingtheboundaryconditionBz (0)| =b(t ). n n ∂Ω n For example, we can take z (s)=z(t +s). n n The convergence analysis is quite similar to that one carried out in our previous paper [4] in the Dirichlet case. The main difference comes from the fact that the correction q is now solution dependent, in general. 4.1. Lie and modified Lie splitting. We commence this section by studying thelocalerrorofLiesplittingappliedto(2.4). Forthispurpose,westartthenumerical solution at time t with the initial value u = u(t ) on the exact solution. We first n n n carry out a substep of size τ with the vector field (3.3b), which gives e (4.7) w (τ)=u +τ f(u )−q (0) +O(τ2), n n n n (cid:0) (cid:1) and then conclude the Lie stepeby integraeting (3.3a) with initial value v (0)=w (τ) n n up to time τ. This provides the sought-after numerical solution L u =z (τ)+eτA w (τ)−z (0) τ n n n n (4.8) τ (cid:0) (cid:1) e + e(τ−s)A q (s)+Dz (s)−∂ z (s) ds. Z0 (cid:16) n n t n (cid:17) The local error δ = L u −u(t ) is obtained by inserting (4.7) into (4.8) and n+1 τ n n+1 subtracting (4.6). This gives the following representation of the local error e δ =τeτA f(u )−q (0) n+1 n n (4.9) (cid:0) τ (cid:1) + ee(τ−s)A q (s)−f(u(t +s)) ds+O(τ2). Z0 (cid:16) n n (cid:17) Expanding q (s)−f(u(t +s))=q (0)−f(u )+O(τ) n n n n e OVERCOMINGORDERREDUCTIONINSPLITTINGMETHODSII 7 we get δ =τ eτA−ϕ (τA) q (0)−f(u ) +O(τ2) n+1 1 n n (4.10) =τ(cid:0)2A ϕ (τA)−ϕ(cid:1)(cid:0)(τA) q (0)−(cid:1)f(u ) +O(τ2), 1 2 n e n (cid:0) (cid:1)(cid:0) (cid:1) where we used (4.4) and (4.5). Henceforth, we will emploey the following assumption on the data of (2.4). Assumption 4.1. Let the domain Ω, the differential operators D and B, and the inhomogeneity d satisfy the assumptions of section 2, let f be differentiable, and assume that u is smooth and satisfies the boundary conditions. 0 UndertheseassumptionsAgeneratesananalyticsemigroup[16,Sect.3.1]. More- over, Proposition 7.1.10 in [16] shows that the solution u of (2.4) is continuously differentiable. Theorem 4.2 (Convergenceofthe classicLiesplitting). Under Assumption 4.1, the classic Lie splitting is convergent of order τ|logτ|, i.e., the global error satisfies the bound ku −u(t )k≤Cτ(1+|logτ|), 0≤nτ ≤T, n n where the constant C depends on T but is independent of τ and n. Proof. The global error e =u −u(t ) satisfies the recursion n n n (4.11) e =L u −L u +δ n+1 τ n τ n n+1 which, by (4.7) and (4.8), can be brought intoethe following form: (4.12) e =eτA e +τ f(u )−f(u ) +δ +O(τ2). n+1 n n n n+1 (cid:16) (cid:0) (cid:1)(cid:17) e Solving this recursion, using the local Lipschitz continuity of f, the very form of the local errors (4.10) with q =0, and the parabolic smoothing property (4.3), we get n n−1 n−1 1 ke k≤Cke k+Cτ ke k+Cτ2 +Cτ. n 0 k kτ X X k=0 k=1 Note that the third term on the right-hand side can by estimated by Cτ(1+|logτ|). This, together with Gronwall’s inequality and ke k=0 gives the desired bound. 0 Under the slightly stronger assumption (4.13) q (0)−f(u(t ))∈D((−A)γ) for some 0<γ ≤1, n n the logτ term in the above theorem can be omitted. This follows at once from the representation of the local errors δ =−τ2(−A)1−γ ϕ (τA)−ϕ (τA) (−A)γ q (0)−f(u(t )) +O(τ2). n+1 1 2 n n (cid:0) (cid:1) (cid:0) (cid:1) For the oblique boundary conditions considered in this paper, we can choose any γ < 1 , see [9]. 2p In particular, for the modified Lie splitting (which is equivalent to setting γ = 1 in (4.13)) we get the following result. Theorem 4.3 (Convergence of the modified Lie splitting). Under Assump- tion 4.1, the modified Lie splitting with q satisfying (3.2) is convergent of order n one. 8 L.EINKEMMER,A.OSTERMANN 4.2. Strang and modified Strang splitting. In order to study the conver- gence properties of Strang splitting and its modifications, we again first analyze its local error when applied to (2.4). For this purpose, we consider one step of the nu- merical solution, starting at time t with the initial value u = u(t ) on the exact n n n solution. The firsthalf step of the modified Strang splitting (see (3.3a)) is then given e by vn(τ2)=zn(τ2)+eτ2A un−zn(0) (4.14) +Z0τ2 e(τ2−(cid:0)se)A(cid:16)qn(s)+(cid:1)Dzn(s)−∂tzn(s)(cid:17)ds, where z denotes a sufficiently smooth function that satisfies the boundary condition n Bz (0)| = b (0). We use here the same function z as in (4.6). The subsequent n ∂Ω n n full step with the vector field (3.3b) then has the representation (4.15) w (s)=v (τ)+t f(w (s))−q (s) +O(s3), 0≤s≤τ. n n 2 n 2 n 2 (cid:0) (cid:1) Note that we have employed a symmetric expansion based on the midpoint of the interval[0,s]. Therefore,theO(s2)termdoesnotshowupintheexpansion. Carrying outagainonehalfstepof (3.3a)withstartingvaluew (τ)providesonestepofStrang n splitting Shun =zn(τ)+eτ2A wn(τ)−zn(τ2) (4.16) e +Z0τ2 e(τ2−(cid:0)s)A(cid:16)qn(τ2 +s)(cid:1)+Dzn(τ2 +s)−∂tzn(τ2 +s)(cid:17)ds. Inserting (4.15) and(4.14) into (4.16) finally gives the following representationofthe numerical solution S u =z (τ)+eτA u −z (0) +τeτA f(w (τ))−q (τ) h n n n n n 2 n 2 (4.17) τ (cid:0) (cid:1) (cid:0) (cid:1) e + e(τ−se)A q (s)+Dz (s)−∂ z (s)) ds+O(τ3). Z0 (cid:16) n n t n (cid:17) The local error is given by δ = S u − u(t ). Taking the difference of the n+1 h n n+1 numerical (4.17) and the exact solution (4.6), we get e δn+1 =τeτ2A f(wn(τ2))−f(u(tn+ τ2)) +τeτ2A f(u(tn+ τ2))−qn(τ2) (cid:0)τ (cid:1) (cid:0) (cid:1) (4.18) + e(τ−s)A q (s)−f(u(t +s)) ds+O(τ3) Z0 (cid:16) n n (cid:17) =(cid:13)1 +(cid:13)2 +(cid:13)3 +O(τ3). This error representation is composed of three terms and a remainder. In order to bound the first term in (4.18), we use the Lipschitz continuity of f and the bound (4.19) kw (τ)−u(t + τ)k≤Cτ2, n 2 n 2 which holds uniformly in n and τ on compact time intervals (see Lemma 4.6 below). This shows that (cid:13)1 is of size O(τ3). Theremainingterms(cid:13)2 and(cid:13)3 aretreatedtogether. EmployingthePeanokernel representation of the error of the midpoint rule τ τ2 s2 τ (τ −s)2 (4.20) ψ(s)ds=τψ(τ)+ ψ′′(s)ds+ ψ′′(s)ds Z0 2 Z0 2 Zτ2 2 OVERCOMINGORDERREDUCTIONINSPLITTINGMETHODSII 9 with ψ(s)=e(τ−s)Aψ(s), ψ(s)=q (s)−f(u(t +s)) n n b b shows that we just have to estimate the two integral remainder terms on the right- hand side of (4.20). The term ψ′′(s) consists of three different terms ψ′′(s)=Ae(τ−s)AAψ(s)−2Ae(τ−s)Aψ′(s)+e(τ−s)Aψ′′(s). b b b The leftmost A in the first two terms can be compensated by parabolic smoothing in exactly the same way as for Lie splitting (see also the proof of Theorem 4.5 below). Therefore, the only term that requires attention is (4.21) χ(s)=e(τ−s)AAψ(s)=e(τ−s)AA q (0)−f(u ) +ρ . n n n (cid:0) (cid:1) Due to the parabolic smoothingbproperty, the rembainder ρ esatisfies the bound n kρ k≤Cτ(τ −s)−1 n uniformly in n. This together with (3.1) shows that χ is uniformly bounded, and proves the following representationof the local error (4.22) δ =Aδ +O(τ3), δ =O(τ3). n+1 n+1 n+1 b b We will employ the following assumption on the data of (2.4). Assumption 4.4. Let the domain Ω, the differential operators D and B, and the inhomogeneity d satisfy the assumptions of section 2, let f be twice differentiable, and assume that u and Du are smooth and satisfy the boundary conditions. 0 0 WearenowinthepositiontostatetheconvergenceresultforthemodifiedStrang splitting. Theorem 4.5 (Convergence of the modified Strang splitting). Let q satisfy n (3.2). Then, under Assumption 4.4 the modified Strang splitting scheme is convergent of order τ2|logτ|, i.e., the global error satisfies the bound ku −u(t )k≤Cτ2(1+|logτ|), 0≤nτ ≤T, n n where the constant C depends on T but is independent of τ and n. Under the assumption q (0)−f(u(t ))∈D((−A)1+γ) for some 0<γ ≤1, n n which is slightly stronger than (3.2), the logτ term in the above theorem can be omitted. This follows in the same way as the corresponding result for Lie splitting, discussed above. For the oblique boundary conditions considered in this paper, we can choose again any γ < 1 , see [9]. 2p Proof of Theorem 4.5. The global error e =u −u(t ) satisfies the recursion n n n (4.23) e =S u −S u +δ n+1 τ n τ n n+1 which, by (4.17) and (4.15), can be brought toethe following form (4.24) e =eτA e +τE(u ,u ) +δ +O(τ2), n+1 n n n n+1 (cid:0) (cid:1) e 10 L.EINKEMMER,A.OSTERMANN Table1 Expected orders of convergence for classic Strang splitting invarious norms. boundary type L1 L2 L∞ β =...=β =0 1.50 1.25 1.00 1 d ∃j with β 6=0 2.00 1.75 1.50 j where E satisfies the bound kE(u ,u )k ≤ Cke k uniformly in n. Using the repre- n n n sentation (4.22) of the local errors δ , we can proceed from here on literally as in n+1 e the proof of Theorem 4.2. The following Lemma was used for the proof of Theorem 4.5. It also shows that modified Lie splitting, carriedout with the flows interchanged,has the same order of convergence. Lemma 4.6. Let q satisfy (3.2). Then, under the assumptions of Theorem 4.2, n the following bound holds (4.25) w (τ)−u(t + τ) ≤Cτ2 n 2 n 2 (cid:13) (cid:13) with a constant C that can(cid:13)be chosen independ(cid:13)ently of n and τ on compact time intervals 0≤t =nτ ≤T. n Proof. From (4.15) with t= τ we infer that 2 w (τ)=v (τ)+ τ f(w (0))−q (0) +O(τ2) n 2 n 2 2 n n (cid:0) (cid:1) with v (τ) givenby (4.14). Using againthe notationψ(s)=q (s)−f(u(t +s)), we n 2 n n have b τ wn(τ2)−u(tn+ τ2)= τ2 f(wn(0))−qn(0) +Z 2 e(τ2−s)Aψ(s)ds (cid:0) (cid:1) 0 = τ2 f(wn(0))−f(un) +Z τ2Z se(τ2−ξb)A(Aψ(ξ)+ψ′(ξ))dξds. (cid:0) (cid:1) 0 0 e b b Since q satisfies the boundary conditions of f(u ), the above integrand fulfills the n n bound e (4.26) ke(τ2−ξ)A(Aψ(ξ)+ψ′(ξ))k≤C andthedoubleintegralisappropriatelybboundebd. Further,sincef islocallyLipschitz continuous, it remains to estimate the difference w (0)−u =v (τ)−u n n n 2 n e =zn(τ2)+eeτ2A un−zn(0) −un = τAϕ (τA) z(cid:0) (0)−u (cid:1)+O(τ). 2 1 2 ne n e (cid:0) (cid:1) But both, zn(0) and un satisfy the boundary conditioens. Therefore, their difference z (0)−u lies in the domain of A. This gives at once the required bound (4.25). n n e 4.3. Explanation of the encountered fractional orders of convergence. e We arenowinthepositiontoexplaintheorderreductionwhichisencounteredinour experimentswiththeclassicStrangscheme. Undertherelaxedassumption(4.13),we get instead of (4.26) the weaker bound (4.27) ke(τ2−ξ)A(Aψ(ξ)+ψ′(ξ))k≤C(τ −ξ)γ−1. 2 b