Convergence in the maximum norm of ADI-type methods for parabolic problems
aa r X i v : . [ m a t h . NA ] F e b Convergence in the maximum norm of ADI-type methods for parabolicproblems
S. Gonz´alez-Pinto and D. Hern´andez-Abreu a,1 a Departamento de An´alisis Matem´atico. Universidad de La Laguna. 38071. La Laguna, Spain.email: [email protected], [email protected]
Abstract
Results on unconditional convergence in the Maximum norm for ADI-type methods, such as the Douglasmethod, applied to the time integration of semilinear parabolic problems are quite difficult to get, mainlywhen the number of space dimensions m is greater than two. Such a result is obtained here under quitegeneral conditions on the PDE problem in case that time-independent Dirichlet boundary conditions areimposed. To get these bounds, a theorem that guarantees, in some sense, power-boundeness of the stabilityfunction independently of both the space and time resolutions is proved. Keywords:
Parabolic PDEs, time integration, stability, power boundedness, convergence, maximum norm,Approximate Matrix Factorization, W-methods, Alternating Direction Implicit schemes.
AMS subject classifications: 65M12, 65M20.
1. Introduction
The present article considers the numerical solution of ODE systems˙ U = D U + g ( t ) , U (0) = U , t ∈ [0 , t ∗ ] , D := m X j =1 D j , g ( t ) := m X j =1 g j ( t ) , (1)stemming from the spatial discretization by using finite differences (or finite volumes) of semilinear parabolicPDEs with constant diffusion coefficients β j > ∂ t u ( t, ~x ) = m X j =1 β j ∂ x j x j u ( t, ~x ) + c ( t, ~x ) , t ∈ [0 , t ∗ ] , ~x = ( x , . . . , x m ) ⊤ ∈ I := (0 , m ,u (0 , ~x ) = u ( ~x ) , ~x ∈ ∂ I ; u ( t, ~x ) = h ( t ) , t ∈ (0 , t ∗ ] , ~x ∈ ∂ I . (2)Here, c ( t, ~x ) is a source term and we assume that its discretization is entirely included in g ( t ), so that g j ( t ) , j = 2 , . . . , m , consist only of contributions from the boundary conditions in the j -direction. Inparticular, we shall be concerned with time-independent Dirichlet boundary conditions, in which case thevectors g ( t ) , . . . , g m ( t ) are constant.To prove convergence in the maximum norm for many numerical methods of splitting type applied to (1), itis customary to get uniform bounds for k R ( τ D , . . . , τ D m ) n k ∞ , n = 1 , , . . . , t ∗ /τ , where τ > R ( τ D , . . . , τ D m ) is a rational mapping acting on the matrices τ D j , 1 ≤ j ≤ m . Typically,we have D j = β j ( I N m ⊗ . . . ⊗ L j ⊗ . . . ⊗ I N ), with L j = tridiag(1 , − , / ∆ x j when second order central This work has been partially supported by the Spanish Project MTM2016-77735-C3-3-P of Ministerio de Econom´ıa, In-dustria y Competitividad.
Preprint submitted to Elsevier February 25, 2021 ifferences are considered in the spatial discretization of (2). Here, we denote the spacing ∆ x j = 1 / ( N j + 1),where N j is the number of equidistant grid-points on the j -direction, and ⊗ stands for the Kronecker productof matrices. It should be observed that the matrices D j pairwise commute. Such methods of splitting typewhen applied to (1) typically produce a recursion for the global errors E n := U n − U ( t n ), n ≥
0, of the form E n +1 = RE n + S n , where U n stands for the numerical solution at t n = nτ , S n denotes the local error and R is the stability matrix associated to the numerical integrator. For ADI-type integrators the stability matrixdepends on D , . . . , D m (see, e.g., [10, Sec. II.2.3]). A relevant example is R ( τ D , . . . , τ D m ) = I + Π( θ ) − τ D, D = m X j =1 D j , (3)where Π( θ ) = ( I − θτ D ) · · · ( I − θτ D m ), which has the associated stability function of m complex variables R ( z , z , . . . , z m ) = 1 + zQ , z := m X j =1 z j , Q := m Y j =1 (1 − θz j ) . (4)For the choice θ = 1 /
2, this is the stability matrix of the Peaceman-Rachford method (when m = 2), alsothe one of the Douglas scheme ([2], [7], [10, p. 373]) and the one of the one-stage AMF-W-method [4].Furthermore, the stability matrix of the so-called Hundsdorfer–Verwer scheme [10, Section IV.5.2], which isa 2-stage W-method of order 2 in general, and of order 3 for θ = (3 + √ /
6, is given by R ( τ D , . . . , τ D m ) = I + 2Π( θ ) − τ D − Π( θ ) − τ D + 12 (cid:0) Π( θ ) − τ D (cid:1) . In this case the stability function is given by R ( z , z , . . . , z m ) = 1 + zQ + z − z Q (with z and Q defined in(4)). The power boundedness in the maximum norm of some R -stability functions was already consideredin [5]. However, for m ≥ n or N [5], i.e., k R ( τ D , . . . , τ D m ) n k ∞ ≤ C m min { (ln( n +1)) m , (ln( N +1)) m } , τ > , nτ ≤ t ∗ , N = max { N j , j = 1 , . . . , m } . With this power bound, convergence results in the maximum norm of size O ( τ p | log( τ ) | m ), when the localerrors are of size O ( τ p +1 ) can be obtained. However, with power bounds of the stability matrix as the onein (5), it can be shown convergence of size O ( τ p ) in case of time independent BCs in (2). In Section 2, weprove a result related to the power boundedness for rational functions. This result is applied in Section 3to show unconditional convergence in the maximum norm for some ADI-methods. In Section 4, numericalexperiments are included to illustrate the orders of convergence regarding the PDE solution for some relevantADI-type methods.
2. Bounds in the maximum norm for rational functions
We look for bounds in the maximum norm of the form k D − R ( τ D , τ D , . . . , τ D m ) n k ∞ ≤ C < ∞ , t n = nτ, n = 0 , , , . . . , τ > ,D = m X j =1 D j , D j = β j ( I m ⊗ · · · ⊗ I j +1 ⊗ L j ⊗ I j − ⊗ · · · ⊗ I ) ,L j = x j tridiag(1 , − , N j , I j is the identity matrix of dimension N j , ∆ x j = N j +1 , (5)where R ( z , . . . , z m ) is a rational function (or a mapping when acting on the matrices D j ) of m complexvariables that is A m ( α )-stable, i.e. | R ( z , . . . , z m ) | ≤ , ∀ z j ∈ W ( α ) := { z ∈ C : | arg( − z ) | ≤ α } ∪ { } , for some 0 < α ≤ π/ . (6)Of course, if (6) holds true for some ¯ α ∈ (0 , π/
2] then it also holds for α = min { ¯ α, π/ } .2 heorem 1. If R ( z , z , . . . , z m ) is a rational function that satisfies (6), then there exists a constant K only depending on α and m such that (5) holds for C = K/ ˆ β, ˆ β =: min ≤ j ≤ m β j > . The proof of this theorem is given below and makes use of the following two lemmas.
Lemma 1.
For any matrix L j (5) it holds that k L − j k ∞ ≤ , (7) k ( zI j − τ L j ) − k ∞ ≤ sec( θ/ | z | , ∀ τ > z = | z | e i θ , | θ | < π, (8) k ( zI j − τ L j ) − k ∞ ≤ τ − | z | , ∀| z | < τ, (9) and k ( zI j − τ L j ) − k ∞ ≤ | z | − τ / ∆ x j , ∀| z | > τ / ∆ x j > . (10) Proof.
The formula in (7) is well known in the literature (see, e.g., [12, formula (4.10)] or [11, p. 43-45]). The formula (8) is an immediate consequence of Lemma 4.1 in [5] (see also [3, formula (5)]), with µ = τ / (∆ x j ) and T j = tridiag(1 , − , N j , since k ( zI j − τ L j ) − k ∞ = µ − k ( µ − zI j − T j ) − k ∞ ≤ µ − sec(arg( µ − z ) / µ − | z | = sec( θ/ | z | . To show (9), by considering | z | < τ , we have that k ( zI j − τ L j ) − k ∞ = k ( τ L j ) − ( I j − zτ − L − j ) − k ∞ ≤ τ − k L − j k ∞ − | z | τ − k L − j k ∞ ≤ τ − − − | z | τ − − = 18 τ − | z | . To show (10), for | z | > τ / ∆ x j it holds k ( zI j − τ L j ) − k ∞ = k z − ( I j − z − τ L j ) − k ∞ ≤ | z | − − | z | − τ k L j k ∞ ≤ | z | − τ / ∆ x j . (cid:3) Lemma 2.
Assume that for any positive integer m we have that z j = − r j e i θ j , r j > , − α ≤ θ j ≤ α, ( j = 1 , , . . . , m ) , ≤ α ≤ π/ . (11) Then | m X j =1 z j | ≥ vuut m X j =1 | z j | ≥ √ m m X j =1 | z j | ≥ √ m m Y j =1 | z j | /m . (12) Proof.
The last two inequalities in (12) follow from the fact that for positive numbers the Quadratic Meanis greater or equal than the Arithmetic Mean and this is greater or equal than the Geometric Mean. Toshow the first inequality, we observe that for complex numbers z , z satisfying (11) it holds that | z + z | = | z | + | z | + 2Re( z z ) ≥ | z | + | z | . Hence, s := z + z fulfils | s | ≥ p | z | + | z | and ithas an angle | ˆ θ | ≤ α ( ≤ π/
4) with the negative x − axis. In particular, s takes the form (11). Then,adding a new complex number z (11) and using the same argument we deduce that s := s + z fulfils | s | ≥ p | s | + | z | ≥ p | z | + | z | + | s | and it has an angle | ˆ θ | ≤ α ( ≤ π/
4) with the negative x − axis.The application of the induction principle concludes the proof. (cid:3) roof of Theorem 1. We define β := max { β , . . . , β m } and use below the following notation for theKronecker product of matrices [ ⊗ A l ] kl =1 := A k ⊗ A k − ⊗ . . . ⊗ A . Consider the positively oriented boundaryΓ = γ ∪ γ ∪ γ of the open domain Ω ⊂ { z : Re z ≤ } , which is symmetric with respect to the negativereal axis in the complex plane, γ := {− r e i θ ; − α ≤ θ ≤ α } , γ := {− r ∗ e i θ ; − α ≤ θ ≤ α } ,γ := {− ρ e − i α , r ≤ ρ ≤ r ∗ } ∪ {− ( r + r ∗ − ρ )e i α , r ≤ ρ ≤ r ∗ } ,r := 4 ˆ βτ, r ∗ := 8 βτ / (∆ x ) , ∆ x := min { ∆ x , . . . , ∆ x m } < , since each ∆ x j = 1 / ( N j + 1) < . (13)Observe that 0 < r < r ∗ . Let us define the rational function (and the associated mapping when acting onmatrices) φ ( z , z , . . . , z m ) = ( z + z + . . . + z m ) − R ( z , z , . . . , z m ) n . (14)Taking into account that φ ( z , z , . . . , z m ) is analytic if ( z , z , . . . , z m ) ∈ ¯Ω m , from the Cauchy’s integralformula applied on each variable we get the following formula by using iterated integrals φ ( z ∗ , z ∗ , . . . , z ∗ m ) = 1(2 π i) m I Γ · · · I Γ φ ( z , z , . . . , z m ) m Y j =1 dz j z j − z ∗ j , ∀ ( z ∗ , z ∗ , . . . , z ∗ m ) ∈ Ω m . (15)By considering the mapping acting on the matrices D j we deduce that φ ( τ D , . . . , τ D m ) = 1(2 π i) m I Γ · · · I Γ φ ( z , z , . . . , z m )[ ⊗ ( z j I j − τ β j L j ) − dz j ] mj =1 . (16)Observe that the eigenvalues of each matrix τ β j L j are λ ( j ) i = − β j τ (∆ x j ) sin (cid:0) π i ∆ x j (cid:1) , i = 1 , . . . , N j , where − r ∗ = − βτ (∆ x ) < λ ( j ) N j < λ ( j )1 = − β j τ (cid:18) sin( π/ · ∆ x j )∆ x j (cid:19) < − β j τ ≤ − r. Hence, the spectrum of τ β j L j falls in ( − r ∗ , − r ) ⊂ Ω.At this point we should notice the identity D − R ( τ D , τ D , . . . , τ D m ) n = τ φ ( τ D , . . . , τ D m ) . From here,taking the maximum norm and using that k M ⊗ M k ∞ = k M k ∞ k M k ∞ for two matrices M and M ,we get that k D − R ( τ D , τ D , . . . , τ D m ) n k ∞ ≤ A := τ (2 π ) m I Γ · · · I Γ | φ ( z , z , . . . , z m ) | m Y j =1 k ( z j I j − τ β j L j ) − k ∞ | dz j | . (17)Next we bound k ( z j I j − τ β j L j ) − dz j k ∞ when z j ∈ Γ. We distinguish three cases.1. z j ∈ γ , then z j = − re i θ j , − α ≤ θ j ≤ α , and | dz j | = rdθ j . From (9) it follows that k ( z j I j − τ β j L j ) − k ∞ ≤ τ β j − r ≤ τ ˆ β − r = 1 r and we deduce that z j ∈ γ = ⇒ k ( z j I j − τ β j L j ) − k ∞ | dz j | ≤ rr dθ j = dθ j . (18)2. z j ∈ γ , then z j = − ρ j e ± i α j , r ≤ ρ j ≤ r ∗ , and | dz j | = dρ j . From (8) it follows that k ( z j I j − τ β j L j ) − k ∞ ≤ sec (( π − α ) / | z j | = 1 ρ j sin( α/ z j ∈ γ = ⇒ k ( z j I j − τ β j L j ) − k ∞ | dz j | ≤ ρ j sin( α/ dρ j . (19)4. z j ∈ γ , then z j = − r ∗ e i θ j , − α ≤ θ j ≤ α , and | dz j | = r ∗ dθ j . From (10) it follows that k ( z j I j − τ β j L j ) − k ∞ ≤ | z j | − τ β j / ∆ x j ≤ r ∗ − r ∗ / r ∗ and we get that z j ∈ γ = ⇒ k ( z j I j − τ β j L j ) − k ∞ | dz j | ≤ dθ j . (20)From the A( α )-stability of R ( z , . . . , z m ) we deduce that | φ ( z , z , . . . , z m ) | ≤ | z + . . . + z m | , ∀ ( z , . . . , z m ) ∈ ( ¯Ω) m . (21)From (12) in Lemma 2, we have that | z + . . . + z m | ≥ √ m Q mj =1 | z j | /m , ∀ ( z , . . . , z m ) ∈ ( ¯Ω) m . Conse-quently, this together with (21) yields | φ ( z , z , . . . , z m ) | ≤ m − / m Y j =1 | z j | − /m , ∀ ( z , . . . , z m ) ∈ ( ¯Ω) m . (22)Now, by considering (17) and (22) we deduce that A ≤ B := τ √ m (2 π ) m X ≤ i ,i ,...,i m ≤ Z γ i · · · Z γ im m Y j =1 (cid:16) k ( z j I j − τ β j L j ) − k ∞ | z j | − /m | dz j | (cid:17) . (23)Taking account that all these iterated integrals can be transformed into products of integrals in one variable,we get B ≤ τ √ m (2 π ) m X m m m mmj ≥ m ! m ! m ! m ! ( A ) m ( A ) m ( A ) m , where, using (18), (19) and (20) respectively, A := Z α − α r − /m dθ = 2 αr − /m ,A := Z r ∗ r ρ − − /m sin( α/ dρ = m sin( α/ (cid:16) r − /m − ( r ∗ ) − /m (cid:17) < m sin( α/ r − /m , and A := Z α − α ( r ∗ ) − /m dθ = 4 α ( r ∗ ) − /m . (24)Then, we have for m + m + m = m , τ ( A ) m ( A ) m ( A ) m ≤ τ C m ,m ,m r − ( m + m ) /m ( r ∗ ) − m /m , with C m ,m ,m := (2 α ) m (cid:18) m sin( α/ (cid:19) m (4 α ) m , and r − ( m + m ) /m ( r ∗ ) − m /m = r − (cid:16) rr ∗ (cid:17) m /m ≤
14 ˆ βτ (cid:18) (cid:19) m /m . Hence, each term τ ( A ) m ( A ) m ( A ) m is bounded since τ ( A ) m ( A ) m ( A ) m ≤ ˆ β − C m ,m ,m · m /m . This concludes the proof. (cid:3) . Convergence in the uniform norm of some ADI-type methods The first goal of this section is to show unconditional convergence of order two in the maximum norm forsemilinear parabolic problems with constant diffusion coefficients (and a time dependent source term) andtime-independent Dirichlet boundary conditions (1)-(2), when the one-step AMF-W-method (henceforthdenoted as
AMF-W1 ) in [4, 6] is considered with the parameter choice θ = 1 / K (0)1 = τ D U n + τ g ( t n ) , ( I − θτ D j ) K ( j )1 = K ( j − + θτ ˙ g j ( t n ) , j = 1 , . . . , m,U n +1 = U n + K ( m )1 , (25)where ˙ v ( t ) stands for the derivative of a function v ( t ) regarding t . The following discussion can be appliedin similar terms to the Douglas method [10, p. 373]. We use the same notations as in [4]. The global errorat the time step t n = nτ is denoted as in [4, formula (2.3)] by E n = U n − U ( t n ) , where U n is the solution ofthe numerical method and U ( t ) = u ( t, ~x G ) is at the same time the exact solution of the (1) and the exactsolution of the PDE on the set of discrete points ~x G of the spatial mesh-grid G . Observe that we will notconsider in our analysis the truncation errors introduced in the spatial discretization of the PDE, since whenusing central differences we get a stable space discretization and the truncated spatial errors do not playany important role in the analysis of global errors (space truncation errors plus time integration errors) as itcan be seen e.g. in [10, Chapt. IV]. It should be remarked that the discretization of the source term c ( t, ~x )is entirely included in g ( t ) [4, Sect. 1]. Besides, the terms ϕ i ( t ) := D i U ( t ) + g i ( t ) , t ∈ [0 , t ∗ ] , ( i = 1 , . . . , m ) , satisfy m X i =1 ϕ i ( t ) = ˙ U ( t ) (26)and they are smooth (i.e. they have bounded first and second derivatives independently of the spatialresolution), since U ( t ) = u ( t, ~x G ) is a smooth function and we have (below δ i,j = 0 for i = j and δ i,i = 1) ϕ i ( t ) := D i U ( t ) + g i ( t ) = β i ∂ x i x i u ( t, ~x G ) + δ i · c ( t, ~x G ) + O (cid:0) (∆ x i ) ∂ x i x i x i x i u ( t, ~x G ) (cid:1) . Additionally, when time independent boundary conditions are assumed in the PDE problem (2), we have˙ g i ( t ) = 0 , i = 2 , . . . , m, and [4, Sect. 1 and 4], D i D i . . . D i r ˙ ϕ j ( t ) = β i β i · · · β i r ∂ r +2 ˙ u ( t, ~x G ) ∂x i ∂x i · · · ∂x i r ∂x j + O (cid:0) (∆ x i ) + · · · + (∆ x i r ) + (∆ x j ) (cid:1) , i < i < · · · < i r < j.D i D i . . . D i r ¨ ϕ j ( t ) = β i β i · · · β i r ∂ r +2 ¨ u ( t, ~x G ) ∂x i ∂x i · · · ∂x i r ∂x j + O (cid:0) (∆ x i ) + · · · + (∆ x i r ) + (∆ x j ) (cid:1) , i < i < · · · < i r < j. (27) Theorem 2.
Assume that the exact solution of the discretized problem (1) satisfies the following uniformbounds k U ( j ) ( t ) k ∞ ≤ C, k g ( i ) ( t ) k ∞ ≤ C, i = 0 , . . . , , j = 0 , . . . , , t ∈ [0 , t ∗ ] , that ˙ g j ( t ) = 0 , j = 2 , . . . , m, and that (27) holds. Then, the global errors, with DE = O ( τ ) , for theAMF-W-method (25) with θ = 1 / fulfill k E n k ≤ C ′ τ , n = 1 , . . . , n ∗ = t ∗ /τ, where the constant C ′ onlydepends on m and C . Proof.
According to [4, formula (2.11)] the global errors of (25) follow the recursion E n = ( R n D − ) DE + n − X j =0 R n − − j S j , n = 1 , , . . . , n ∗ = t ∗ /τ, (28)6here the matix R is given by (3) and the discretization local errors are given by [4, formula (2.10)] S n = Π( θ ) − (cid:0) τ ˙ U ( t n ) + θτ ˙ G ( t n ) (cid:1) − (cid:0) U ( t n + τ ) − U ( t n ) (cid:1) , (29)with G ( t ) = P mi =1 (cid:16)Q i − j =1 ( I − θτ D j ) (cid:17) g i ( t ), and the convention Q j =1 ( I − θτ D j ) = I (see [4, formula (2.4)]).We also make use of other expression for the global errors (see [4, formula (4.11)]), obtained by partialsummation in (28), E n = ( I − R n )( I − R ) − S + n − X j =0 ( I − R n − − j )( I − R ) − ( S j +1 − S j ) , n = 1 , , . . . , (30)and of a simplified expression for the local errors given in [4, formula (4.7)] S n = S (1) n + S (2) n ,S (1) n = τ − m X i =1 (Π i − Π) ˙ ϕ i ( t n ) , S (2) n = − τ Z (1 − s ) ... U ( t n + sτ ) ds, Π = ( I − θτ D ) · · · ( I − θτ D m ) , Π i := ( I − θτ D ) · · · ( I − θτ D i − ) , ( i > , Π := I. (31)Since, from (26), ϕ ( t ) = ˙ U ( t ) − P mi =2 ϕ i ( t ), we can split the term S (1) n of the local error in two parts as S (1) n = S (1 ,a ) n + S (1 ,b ) n , S (1 ,a ) n := τ − ( I − Π) ¨ U ( t n ) , S (1 ,b ) n := τ − m X i =2 (Π i − I ) ˙ ϕ i ( t n ) . (32)Now, to bound the global errors generated by the contributions of each term of the local error we make useof the following bounds, which are a consequence of Theorem 1. In this case, we can apply Theorem 1 dueto the result in [8] where is guaranted for R ( z , . . . , z m ) that α = α m = π m − , m ≥
2. Consequently itholds that k D − R n k ∞ ≤ C , n = 0 , , . . . , n ∗ = t ∗ /τ, τ > , N j ∈ N . (33)From the assumption ˙ g j ( t ) = 0 , j = 2 , . . . , m, on time-independent boundary conditions (see (27) and [4,Example 4.6]) k ( I − Π) ¨ U ( t ) k ∞ ≤ C τ, k ( I − Π)... U ( t ) k ∞ ≤ C τ, k ( I − Π i ) ˙ ϕ j ( t ) k ∞ ≤ C τ, k ( I − Π i ) ¨ ϕ j ( t ) k ∞ ≤ C τ, when j ≥ i. (34) (1) We start with the global errors E (1 ,a ) n generated by the local errors S (1 ,a ) n . To bound them, from (30)and (3) E (1 ,a ) n = ( I − R n )( I − R ) − S (1 ,a )0 + n − X j =0 ( I − R n − − j )( I − R ) − ( S (1 ,a ) j +1 − S (1 ,a ) j )= − τ I − R n ) D − ( I − Π) ¨ U ( t ) − τ n − X j =0 ( I − R n − − j ) D − ( I − Π) (cid:16) ¨ U ( t j + τ ) − ¨ U ( t j ) (cid:17) = O ( τ ) . (2) For the global errors E (1 ,b ) n generated by the local errors S (1 ,b ) n , we first take into account that S (1 ,b ) n = τ − m X i =2 (Π i − I ) ˙ ϕ i ( t n ) = τ D − ( R − I ) m X i =2 (Π i − I ) ˙ ϕ i ( t n ) . (35)7hen, it follows from (30) and (3) that E (1 ,b ) n = ( I − R n )( I − R ) − S (1 ,b )0 + n − X j =0 ( I − R n − − j )( I − R ) − ( S (1 ,b ) j +1 − S (1 ,b ) j )= − τ I − R n ) D − m X i =2 (Π i − I ) ˙ ϕ i ( t ) − τ m X i =2 n − X j =0 ( I − R n − − j ) D − (Π i − I ) ( ˙ ϕ i ( t j + τ ) − ˙ ϕ i ( t j )) = O ( τ ) . (3) For the global errors E (2) n generated by the local errors S (2) n , we use the formula (28). In this case wedefine V ( t ) := D (cid:16) − R (1 − s ) ... U ( t + sτ ) ds (cid:17) . (36)so that S (2) n = τ D − V ( t n ) . Besides, V ( t ) = − R (1 − s ) D ... U ( t + sτ ) ds = − R (1 − s ) (cid:0) U (4) ( t + sτ ) − ... g ( t + sτ ) (cid:1) ds. (37)From here we deduce (under the regularity assumption in the exact solution) that k V ( t ) k ∞ ≤ C . Now, thebound for the global errors follows from E (2) n = n − X j =0 R n − − j S (2) j = τ n − X j =0 R n − − j D − V ( t j ) = O ( τ ) . (38) (cid:3) Remark 1.
Second order of convergence in the maximum norm for the Douglas method with θ = isproved in [1, Theorem 3.1] under the assumption of power-boundedness for the stability matrix R in (3) andassuming that [1, (3.16b), p. 271] τ k − D − D l D l · · · D l k v ( t n ) = O (1) , ≤ l < . . . < l k < i ≤ m, ( v = ˙ ϕ i , ¨ ϕ i ) . (39)This assumption was also useful in [9, Theorem 3.2] in order to prove convergence for linear multistepmethods with stabilizing corrections applied to split ODEs. Although (39) is closely related to (27), theproof of convergence presented in [1] does require the assumption on power-boundedness for the stabilitymatrix R , which, as far as we are aware, has not been shown for m ≥ (cid:3) Convergence of order two in the maximum norm for the Douglas method [10, p. 373] and time independentboundary conditions can also be shown following similar steps as in the proof of Theorem 2. To this aim,let us consider the Douglas method applied to (1): v = U n + τ ( DU n + g ( t n )) v i = v i − + θτ (( D i v i + g i ( t n +1 )) − ( D i U n + g i ( t n ))) , i = 1 , . . . , m,U n +1 = v m . (40) Theorem 3.
Under the same assumptions of Theorem 2, the global errors, with DE = O ( τ ) , for theDouglas method (40) with θ = 1 / fulfill k E n k ≤ C ′ τ , n = 1 , . . . , n ∗ = t ∗ /τ, where the constant C ′ onlydepends on m and C . roof. The global errors E n = U n − U ( t n ) for the method (40) fulfill the recursion E n +1 = RE n + S n , n ≥
0, where the stability matrix R is given by (3) and the local errors S n are obtained as given in [10,formula (3.15)] by S n = − Q − m · · · Q − ( r + r ) − Q − m · · · Q − r − . . . − Q − m r m , (41)with Q i = I − θτ D i , 1 ≤ i ≤ m , whereas, taking into account that in (40) v i ≈ U ( t n + 1), 0 ≤ i ≤ m , for r i = r i ( t n ) it holds that r ( t n ) = U ( t n +1 ) − U ( t n ) − τ ˙ U ( t n ) = τ U ( t n ) + τ Z (1 − s ) ... U ( t n + sτ ) ds, (42)and, for 1 ≤ i ≤ m , r i ( t n ) = − θτ ( ϕ i ( t n +1 ) − ϕ i ( t n )) = − θτ Z ˙ ϕ i ( t n + sτ ) ds = − θτ ˙ ϕ i ( t n ) − θτ Z (1 − s ) ¨ ϕ i ( t n + sτ ) ds. (43)With Π and Π i defined in (31), we can rewrite (41) as S n = − Π − (( r + r ) + Q r + Q Q r + . . . + Q · · · Q m − r m )= − Π − (cid:0) m X i =0 r i (cid:1) + (Π − I ) r + (Π − I ) r + . . . + (Π m − I ) r m ! , r i = r i ( t n ) , ≤ i ≤ m. (44)With θ = , using (42), (43) and (26), it holds that m X i =0 r i ( t n ) = τ Z ( s − s )... U ( t n + sτ ) ds. (45)Now, we split the local error in two parts as S n = S ( a ) n + S ( b ) n ,S ( a ) n := − Π − (cid:0) P mi =0 r i (cid:1) , S ( b ) n := − Π − ((Π − I ) r + (Π − I ) r + . . . + (Π m − I ) r m ) . (46)Partial summation in the global error recursion leads us to the relation (30), with ( I − R ) − Π − = − τ − D − . Now, we bound the global errors generated by the contributions of each term of the local error. (a)
For the global errors E ( a ) n generated by the local errors S ( a ) n , using (30), (45) and (46) we have that E ( a ) n = ( I − R n ) D − τ − (cid:0) m X i =0 r i ( t ) (cid:1) + n − X j =0 ( I − R n − − j ) D − τ − (cid:0) m X i =0 r i ( t j +1 ) − r i ( t j ) (cid:1) = τ I − R n ) D − Z ( s − s )... U ( t + sτ ) ds + τ n − X j =0 ( I − R n − − j ) D − Z ( s − s ) (cid:0) ... U ( t j +1 + sτ ) − ... U ( t j + sτ ) (cid:1) ds = O ( τ ) . (b) For the global errors E ( b ) n generated by the local errors S ( b ) n , using (30), (43), (46) and (34) it holds that9 ( b ) n = ( I − R n ) D − τ − m X i =2 (Π i − I ) r i ( t ) + n − X j =0 ( I − R n − − j ) D − τ − m X i =2 (Π i − I )( r i ( t j +1 ) − r i ( t j ))= τ I − R n ) D − m X i =2 Z ( I − Π i ) ˙ ϕ i ( t + sτ ) ds + τ m X i =2 n − X j =0 ( I − R n − − j ) D − Z ( I − Π i ) (cid:0) ˙ ϕ i ( t j +1 + sτ ) − ˙ ϕ i ( t j + sτ ) (cid:1) ds = O ( τ ) . (cid:3) Remark 2.
A modified one-stage AMF-W method (henceforth denoted as modified AMF-W1 ) K (0)1 = τ D U n + τ g ( t n )( I − θτ D j ) K ( j )1 = K ( j − + θτ ˙ g j ( t n + τ / , j = 1 , . . . , m,U n +1 = U n + K ( m )1 , (47)was introduced in [4]. For this method, second order convergence in the ℓ ∞ − norm for time independentboundary conditions can be shown as in Theorem 2 above considering that its local error can be expressedas (see [4, formula (5.4)] ) S n = τ Π − Z / s ... U ( t n + sτ ) s. + τ − (cid:16) m X i =2 (cid:16) Π i − I (cid:17) ˙ ϕ i ( t n + τ ) (cid:17) − τ Z k ( s ) ... U ( t n + sτ ) s. , k ( s ) := min (cid:0) s , (1 − s ) (cid:1) . (cid:3) Theorem 4.
Under the same assumptions of Theorem 2, the global errors, with DE = O ( τ ) , for themodified AMF-W method (47) with θ = 1 / fulfill k E n k ≤ C ′ τ , n = 1 , . . . , n ∗ = t ∗ /τ, where the constant C ′ only depends on m and C . Proof.
The proof follows along the lines of the proofs of Theorems 2 and 3 (cid:3)
Remark 3.
For m ≥ ℓ ∞ -norm (up to a logarithmic factor) when applied to (1) with time dependent Dirichletboundary conditions. When m = 2, the methods (47) and (40) have the advantage that their global errorin the ℓ ∞ -norm is min {O ( τ | log h | ) , O ( τ ) } (see [4, Section 5]), whereas the method in (25) only hasconvergence of size O ( τ ) in the maximum norm . (cid:3)
4. Numerical Illustration
We first consider the linear diffusion partial differential equation (2) in three and four spatial dimensions,with diffusion coefficients β j = 1, 1 ≤ j ≤ m . Our aim is to illustrate numerically the second orderconvergence in the maximum norm for the one-stage AMF-W method (25) (and its modified version (47))and the Douglas method (40), both with parameter θ = , when time independent boundary conditions areimposed on the PDE. For time dependent boundary conditions, order one (up to a logarithmic factor) is10ttained by both methods when the spatial dimension is m ≥
3. For our numerical experiments we considerthat c ( t, ~x ) is selected in such way that u ( t, ~x ) = u e ( t, ~x ) := e t (cid:18) m m Y j =1 x j (1 − x j ) + κ m X j =1 (cid:16) x j + 1 j + 2 (cid:17) (cid:19) (48)is the exact solution of (2). We impose the initial condition u (0 , ~x ) = u e (0 , ~x ) and Dirichlet boundaryconditions. Here, we consider the cases m = 3 ,
4. If κ = 0 we have homogeneous boundary conditions, butwhen κ = 1 we get non-homogeneous time-dependent Dirichlet conditions.We apply the MOL approach on a uniform grid with meshwidth h = ∆ x i = 1 / ( N + 1), 1 ≤ i ≤ m , where N = 2 j − j = 2 , . . . , j max , with j max = 7 for m = 3 and j max = 5 if m = 4. Hence, a semi-discretizedsystem with corresponding dimension N m of the form (1) is obtained, where D is given in (5) and g ( t )includes the discretization of the term c ( t, ~x ) and the terms due to non-homogeneous boundary conditions.Observe that the exact solution (48) is a polynomial of degree 2 in each spatial variable so that the globalerrors come only from the time discretization. The methods (25), (47) and (40) are then applied to (1) withfixed step size τ = h = 2 − j , 2 ≤ j ≤ j max , and the corresponding global errors regarding the PDE solutionversus the stepsize are displayed below in Figure 1 in case of time independent boundary conditions ( κ = 0)and in Figure 2 in case of time dependent boundary conditions ( κ = 1). In the first case, all methods displaysecond order convergence in the ℓ ∞ -norm in both dimensions m = 3 and m = 4. In the second situation withtime dependent boundary conditions, all methods suffer an order reduction and the corresponding orders ofconvergence are at most one. −2 −1 −4 −3 −2 −1 stepsize ( τ =h) l ∞ − e rr o r Douglas ( θ =0.5)AMF−W1 ( θ =0.5)modified AMF−W1 ( θ =0.5) τ −3 −2 −1 stepsize ( τ =h) l ∞ − e rr o r Douglas ( θ =0.5)AMF−W1 ( θ =0.5)modified AMF−W1 ( θ =0.5) τ Figure 1: Error in the ℓ ∞ − norm vs stepsize on the linear model (2)-(48) with time independent boundary conditions ( κ = 0)and τ = h = ∆ x i , 1 ≤ i ≤ m . Spatial dimension m = 3 (left) and m = 4 (right). A dashed straight line with slope two isincluded to compare the PDE order of convergence. A second numerical experiment is included below in Figure 3 for the case of variable diffusion coefficients β j = β j ( ~x ), 1 ≤ j ≤ m . Although a theoretical analysis for such a case lies beyond the scope of this paper,similar orders of convergence are observed. To illustrate this assertion, we consider m = 3 spatial dimensionsand diffusion coefficients β = β ( x, y, z ) = (1 + xyz ) , β = β ( x, y, z ) = e x − y +3 z and β = β ( x, y, z ) = (1 + x ) e − y z . (49)Again c ( t, ~x ) is selected in such way that (48) is the exact solution of (2), with homogeneous boundaryconditions if κ = 0 and time-dependent Dirichlet conditions when κ = 1. The MOL approach is applied ona uniform grid with meshwidth h = ∆ x i = 1 / ( N + 1), 1 ≤ i ≤
3, where N = 2 j − j = 2 , . . . ,
7. The11 −2 −1 −4 −3 −2 −1 stepsize ( τ =h) l ∞ − e rr o r Douglas ( θ =0.5)AMF−W1 ( θ =0.5)modified AMF−W1 ( θ =0.5) τ −3 −2 −1 stepsize ( τ =h) l ∞ − e rr o r Douglas ( θ =0.5)AMF−W1 ( θ =0.5)modified AMF−W1 ( θ =0.5) τ Figure 2: Error in the ℓ ∞ − norm vs stepsize on the linear model (2)-(48) with time dependent boundary conditions ( κ = 1)and τ = h = ∆ x i , 1 ≤ i ≤ m . Spatial dimension m = 3 (left) and m = 4 (right). A dashed straight line with slope one isincluded to compare the PDE order of convergence. results displayed in Figure 3 (left) show that all methods provide second order of convergence in the ℓ ∞ norm when time-independent boundary conditions are considered ( κ = 0). For the case of time-dependentboundary conditions ( κ = 1), Figure 3 (right) show an order reduction to at most order one for the threemethods considered. −2 −1 −5 −4 −3 −2 −1 stepsize ( τ =h) l ∞ − e rr o r Douglas ( θ =0.5)AMF−W1 ( θ =0.5)modified AMF−W1 ( θ =0.5) τ −2 −1 −4 −3 −2 −1 stepsize ( τ =h) l ∞ − e rr o r Douglas ( θ =0.5)AMF−W1 ( θ =0.5)modified AMF−W1 ( θ =0.5) τ Figure 3: Error in the ℓ ∞ − norm vs stepsize on the linear model (2)-(48) with variable diffusion coefficients (49), τ = h = ∆ x i ,1 ≤ i ≤ m , and spatial dimension m = 3. Time independent boundary conditions κ = 0 (left) and time dependent boundaryconditions κ = 1 (right). Dashed straight lines with slopes two and one, respectively, are included to compare the PDE orderof convergence. References [1]
A. Arrar´as, K.J. in ’t Hout, W. Hundsdorfer and L. Portero . Modified Douglas splitting methods for reaction-diffusionequations. BIT, 57(2):261–285, 2017.[2]
J. Douglas Jr . Alternating direction methods for three space variables. Numer. Math. 4, 41–63, 1962.[3]
I. Farag´o and C. Palencia . Sharpening the estimate of the stability constant in the maximum-norm of the Crank-Nicolsonscheme for the one-dimensional heat equation.
Appl. Numer. Math. , 42(1-3):133–140, 2002. Ninth Seminar on NumericalSolution of Differential and Differential-Algebraic Equations (Halle, 2000). S. Gonz´alez-Pinto, E. Hairer and D. Hern´andez-Abreu . Convergence in ℓ and ℓ ∞ norm of one stage AMF-W-methodsfor parabolic problems. SIAM J. Numer. Anal. , 58 (2), 1117–1137, 2020.[5]
S. Gonz´alez-Pinto, E. Hairer and D. Hern´andez-Abreu . Power boundedness in the maximum norm ofstability matrices for ADI methods.
Accepted in BIT Numerical Math., 2021 ∼ hairer/preprints.html.[6] S. Gonz´alez-Pinto, E. Hairer, D. Hern´andez-Abreu and S. P´erez-Rodr´ıguez . AMF-type W-methods for parabolic problemswith mixed derivatives.
SIAM J. Sci. Comput. , 40(5):A2905–A2929, 2018.[7]
W. Hundsdorfer . A note on stability of the Douglas splitting.
Math. Comput.
67, 183–190, 1998.[8]
W. Hundsdorfer . Stability of approximate factorizations with θ -methods. BIT
39, 473–483, 1999.[9]
W. Hundsdorfer and K.J. in ’t Hout.
On multistep stabilizing correction splitting methods with applications to the Hestonmodel. SIAM J. Sci. Comput., 40(3), A1408–A1429, 2018.[10]
W. Hundsdorfer and J.G. Verwer . Numerical solution of time-dependent advection diffusion reaction equations.
Springerseries in comput. math. , Springer, 2003.[11]
S. Larsson and V. Thom´ee . Partial Differential Equations with Numerical Methods, Springer, 2009.[12]
R.M.M. Mattheij and M.D. Smooke , Estimates for the inverse of tridiagonal matrices arising in BVPs,
Linear Alg. Appl.
73, 33–57, 1986.73, 33–57, 1986.