MATHEMATICS OF COMPUTATION Volume 71, Number 240, Pages 1529–1543 S 0025-5718(01)01362-X Article electronically published on November 19, 2001
AVOIDING THE ORDER REDUCTION OF RUNGE-KUTTA METHODS FOR LINEAR INITIAL BOUNDARY VALUE PROBLEMS
M. P. CALVO AND C. PALENCIA
Abstract. A new strategy to avoid the order reduction of Runge-Kutta methods when integrating linear, autonomous, nonhomogeneous initial boundary value problems is presented. The solution is decomposed into two parts. One of them can be computed directly in terms of the data and the other satis?es an initial value problem without any order reduction. A numerical illustration is given. This idea applies to practical problems, where spatial discretization is also required, leading to the full order both in space and time.
1. Introduction Let us consider an abstract initial value problem (IVP) (1.1) u (t) = Au(t) + f (t), u(0) = u0 , 0 ≤ t ≤ T,
where A : D(A) ? X → X is the in?nitesimal generator of a C0 -semigroup in a Banach space X , f : [0, T ] → X and u0 ∈ X . It is well known that many evolutionary problems of practical interest, either hyperbolic or parabolic, can be written in this form [28]. Problem (1.1) is integrated in time by using a Runge-Kutta (RK) method with Butcher array (1.2) c A , bT
s s s s s×s . Given an with c = (ci )s i=1 ∈ R , b = (bi )i=1 ∈ R , A = (aij )i,j =1 ∈ R approximation un ∈ X to u(tn ), 0 ≤ tn < T , the numerical approximation un+1 to u(tn+1 ) at tn+1 = tn + kn ≤ T is de?ned by
(1.3)
un+1 = un + kn (bT ? A)Un + kn (bT ? I )F (tn ), Un = e ? un + kn (A ? A)Un + kn (A ? I )F (tn ),
s where F (tn ) = (f (tn + kn ci ))s i=1 and Un ∈ X is the solution of
(1.4)
Received by the editor January 14, 2000 and, in revised form, November 30, 2000. 2000 Mathematics Subject Classi?cation. Primary 65M12, 65M20. Key words and phrases. Abstract initial boundary value problems, Runge-Kutta, order reduction. This research has been supported by DGICYT under project PB95-705 and by Junta de Castilla y Le? on under project VA36/98.
c 2001 American Mathematical Society
1529
1530
M. P. CALVO AND C. PALENCIA
with e = [1, . . . , 1]T ∈ Rs . Here, and throughout the paper, I stands for the identity operator in any space. It is known that for su?ciently small step-sizes kn , (1.4) possesses a unique solution [5]. Let u ∈ C p+1 ([0, T ], X ) be a solution of (1.1), p being the classical order of the method. Under natural stability requirements, for ?xed step-sizes kn = k and homogeneous problems (i.e., f = 0), the error u(tn ) ? un behaves like O(k p ) [11]. However, for nonhomogeneous problems and unbounded operators A, it turns out that the classical order p is not attained [10, 24, 32]. This is the so-called order reduction phenomenon. In fact, for nonhomogeneous problems the order of convergence is rather governed by the stage order q of the method. Set wl = cl ? lAcl?1 , 1 ≤ l ≤ p,
s 0 where cl = (cl i )i=1 and c = e. Recall that q is de?ned as the largest integer for which wl = 0, 1 ≤ l ≤ q . Typically the order of convergence is min{p, q + 1 + θ} with 0 ≤ θ < 2 (possibly θ is not integer) [24]. Let us point out that, in the context of PDEs, this order reduction does not occur in the interior of the spatial domain [21]. IVPs of the form (1.1) are a particular instance of the more general initial boundary value problems (IBVPs) ? 0 ≤ t ≤ T, ? ? u (t) = Au(t) + f (t), (1.5) u(0) = u0 , ? ? ?u(t) = g (t), 0 ≤ t ≤ T,
where A : D(A) ? X → X is a linear extension of A, ? is a linear mapping from D(A) to another Banach space Y and g : [0, T ] → Y . The precise conditions to be ful?lled by these operators are described in Section 2 (see also [5]). Notice that the IVPs (1.1) arising in the applications always correspond to certain natural IBVPs (1.5) with g = 0. When (1.5) is integrated in time with the RK method (1.2), the order of convergence is min{p, q + θ} with 0 ≤ θ < 2 [5]. Thus, for IBVPs the order reduction is more severe than for IVPs. Di?erent remedies have been proposed in the literature to avoid the order reduction [1, 2, 3, 5, 12, 19, 20, 30, 32]. In the present paper a new strategy yielding the full order p in the time integration of (1.5) is introduced. This strategy cannot be applied directly to the IVP (1.1). However, this is not a serious drawback since in practice, as we mentioned above, (1.1) can be written in the format (1.5). The basic idea is as follows. The solution of (1.5) is decomposed as u(t) = v (t) + u? (t), where v (t) can be computed directly in terms of the data and u? (t) is the solution of a suitable IVP of the form (1.1). The key point of this construction is that u? (t) ∈ D(Ap?q ), so that there is no order reduction in the RK approximation to u? (t). This strategy is studied in Section 3. Since the implementation of our strategy demands the solution of additional stationary problems at each step (see Remark 3.6 below), a high cost may result. Thus, though the strategy is interesting from a theoretical point of view, the standard implementation of an RK method with appropriate stage order could be advantageous in some instances. Numerical experiments are presented in Section 4, illustrating the theoretical results. In the ?nal Section 5 it is shown that the remedy proposed in the present paper can also be accommodated in practical situations, where full discretizations of (1.5) are
AVOIDING THE ORDER REDUCTION FOR IBVPS
1531
required. It is worth mentioning here that, as di?erent from previous approaches [2, 3], standard spatial discretizations can be used. 2. Abstract initial boundary value problems In this section we describe brie?y the setting of [28]. Let X, Y be two Banach spaces. The norms in X and Y and the operator norms are denoted by · . Let A : D(A) ? X → X , ? : D(A) ? X → Y be linear operators. We assume that H1. The operator (A, ? ) : D(A) ? X → X × Y is closed. Moreover, ? : D(A) ? X → Y is onto. H2. Set D(A) = ker (? ) = { x ∈ D(A) / ?x = 0 }. Then the restriction A = A|D(A) : D(A) ? X → X is the in?nitesimal generator of a C0 -semigroup S (t), t ≥ 0, in X . It is well known that there exist M > 0 and ω ∈ R such that S (t) ≤ M eωt , t ≥ 0 (see, e.g., [31]). In [28] it is proved that H1–H2 imply the following condition: For each λ > ω there exists a linear and bounded operator K (λ) : Y → X such that, for each y ∈ Y , x = K (λ)y ∈ X is the unique solution of the problem (2.1) (A ? λ)x = 0, ?x = y.
Let us denote, for an integer r ≥ 0, Wr = D(Ar ) and Wr = D(Ar ). It turns out that for each λ > ω the operator K (λ)? : Wr → Wr is bounded, where Wr is endowed with the norm x
Wr
= x + Ax + · · · + Ar x ,
x ∈ Wr ,
which is equivalent to the graph-norm of Ar . Let T > 0 and let f : [0, T ] → X be continuous and g : [0, T ] → Y a C 1 -mapping. Assume that u0 ∈ W1 satis?es ?u0 = g (0). Then the IBVP (1.5) has got a unique solution u : [0, T ] → X [28]. Next we give an example to illustrate these concepts (see [4, 28] for additional examples). Example 2.1. Let us consider the initial boundary value problem for the onedimensional heat equation ? x ∈ [0, 1], 0 ≤ t ≤ T, ut (x, t) = uxx (x, t) + f (x, t), ? ? ? ? u(x, 0) = u (x), x ∈ [0, 1], 0 (2.2) ? 0 ≤ t ≤ T, u(0, t) = g0 (t), ? ? ? 0 ≤ t ≤ T, u(1, t) = g1 (t), where f : [0, 1] × [0, T ] → R, u0 : [0, 1] → R, g0 , g1 : [0, T ] → R are given functions. Problem (2.2) can be written in the abstract form (1.5) as follows: X = Lp ([0, 1]), 1 ≤ p < +∞, or X = C0 ([0, 1]) if p = +∞, Y = R × R, D(A) = W 2,p ([0, 1]) = { ? ∈ X / ? ∈ X }, A? = ? and ?? = [?(0), ?(1)]T , ? ∈ D(A). In this example 1,p ([0, 1]) = { ? ∈ W 2,p ([0, 1]) / ?(0) = ?(1) = 0 } and D(A) = W 2,p ([0, 1]) ∩ W0 T for ξ = [ξ0 , ξ1 ] ∈ Y , K (0)ξ is the mapping sending x ∈ [0, 1] to (1 ? x)ξ0 + xξ1 . Finally, for r ≥ 1, Wr = W 2r,p ([0, 1]) and Wr = { ? ∈ W 2r,p ([0, 1]) / ?(2j ) (0) = ?(2j ) (1) = 0, 0 ≤ j ≤ r }.
1532
M. P. CALVO AND C. PALENCIA
Assume that g ∈ C m ([0, T ], Y ) and f (j ) ∈ C m?j ?1 ([0, T ], Wm?j ?1 ), 0 ≤ j ≤ m ? 1, for some m ≥ 1. Notice that in practice this is a natural requirement since these conditions only demand regularity of f (j ) but they do not impose any restriction on the values of f (j ) at the boundary of the spatial domain. This is opposite to the unnatural assumptions f (j ) (t) ∈ Wm?j ?1 , 0 ≤ j ≤ m ? 1, 0 ≤ t ≤ T (see example above). Under certain compatibility conditions among the data f , g and u0 (see [2]) the solution u : [0, T ] → X of (1.5) belongs to C m ([0, T ], X ) ∩ C ([0, T ], W1 ). As long as H1 is satis?ed, we can di?erentiate in equation (1.5) and get that u (t) ∈ W1 and u (t) = Au (t) + f (t), ?u (t) = g (t), 0 ≤ t ≤ T, 0 ≤ t ≤ T.
Moreover, since u (t) = Au(t) + f (t) and f (t) ∈ W1 , we have u(t) ∈ W2 and then u (t) = A2 u(t) + Af (t) + f (t). In fact, by induction, it is easy to prove that, for 1 ≤ l ≤ m and 0 ≤ t ≤ T , u(t) ∈ Wl , u(l?1) (t) ∈ W1 and ? l? 1 ? ? u(l) (t) = Al u(t) + Al?j ?1 f (j ) (t), (2.3) j =0 ? ? ?u(l) (t) = g (l) (t). This also yields
l? 1
Al u(t) = u(l) (t) ?
j =0
Al?j ?1 f (j ) (t),
so that
l? 1
(2.4)
? Al u(t) = g (l) (t) ?
j =0
? Al?j ?1 f (j ) (t).
In summary, we have shown the way to obtain the boundary values of Al u in terms of the data f and g and their derivatives. In the next section we will use the related expression for λ > ω (2.5)
l j ?1 i=0
? (λ ? A)l u(t) = λl g (t) +
j =1
l j
λl?j (?1)j g (j ) (t) ?
? Aj ?i?1 f (i) (t) ,
which is a direct consequence of (2.4). 3. The main result Let u : [0, T ] → X be a solution of (1.5). We try to write u = v + u? , where u : [0, T ] → X is the solution of a suitable IVP
?
(3.1)
u? (t) = Au? (t) + f ? (t), u? (0) = u? 0,
AVOIDING THE ORDER REDUCTION FOR IBVPS
1533
? and where v : [0, T ] → X , u? 0 and f : [0, T ] → X can be computed from the data u0 , f and g in such a way that the RK method applied to (3.1) does not su?er from order reduction. This idea is frequently used with v (t) = Eg (t), where E : Y → X is a linear operator such that ?Ey = y for y ∈ Y . This choice leads to the initial value problem
(3.2)
u? (t) = Au? (t) + f (t) ? Eg (t) + AEg (t), u? (0) = u0 ? Eg (0).
Since u? ∈ D(A), then the order of the RK discretization of (3.2) is k q+1 (see [5, 24]). Thus, for q < p ? 1 and this choice of v , the order reduction still appears for u? . Our goal is to construct a more appropriate correcting term v . To this end we de?ne the operators Cr (λ) : Wr → X as
r ?1
(3.3)
Cr (λ) =
s=0
(λ ? A)?s K (λ)? (λ ? A)s ,
λ > ω,
r ≥ 1.
Before de?ning v we give the following technical lemma which will be used in the proof of Theorem 1. Lemma 3.1. Let r ≥ 1, λ > ω . Then (i) Cr (λ) is a bounded operator from Wr to Wr . (ii) For 0 ≤ j ≤ r ? 1 and x ∈ Wr it holds that ? (λ ? A)j Cr (λ)x = ? (λ ? A)j x. Proof. Since K (λ)? is a bounded operator from Wr?s to Wr?s , 0 ≤ s ≤ r ? 1, it is clear that (λ ? A)?s K (λ)? (λ ? A)s is a bounded operator from Wr to Wr . This proves (i). Let x ∈ W1 . Since (λ ? A)K (λ)?x = 0 and (λ ? A)(λ ? A)?1 x = x, it turns out that, for x ∈ Wr , r ≥ 2,
r ?1
(λ ? A)Cr (λ)x (3.4)
= =
(λ ? A)K (λ)?x +
s=1
(λ ? A)(λ ? A)?s K (λ)? (λ ? A)s x
Cr?1 (λ)(λ ? A)x.
We now proceed by induction in the proof of (ii). For x ∈ W1 , C1 (λ)x = K (λ)?x and (ii) for r = 1 is a direct consequence of the de?nition of K (λ). Let us now assume that (ii) holds for r ? 1, r ≥ 2. Let x ∈ Wr . By (3.4), for 1 ≤ j ≤ r ? 1 ? (λ ? A)j Cr (λ)x = ? (λ ? A)j ?1 Cr?1 (λ)(λ ? A)x. Then using (ii) with r ? 1, we conclude that ? (λ ? A)j Cr (λ)x = ? (λ ? A)j ?1 (λ ? A)x = ? (λ ? A)j x. Moreover, (ii) is clear for j = 0 since ? (λ ? A)?1 = 0. We are now in a position to de?ne the correcting term. Assume that the solution u of (1.5) belongs to C p+1 ([0, T ], Wr ) for some 1 ≤ r ≤ p ? q . We want to stress
1534
M. P. CALVO AND C. PALENCIA
again that, in practice, this condition is natural, since it only demands the solution be smooth both in time and space. We de?ne (3.5) v (t) = Cr (λ)u(t), 0 ≤ t ≤ T,
where Cr (λ) is given by (3.3) with λ > ω ?xed. Observe that r = 1 corresponds to the usual correction mentioned above with E = K (0). Notice that, due to (2.5), v (t) can certainly be computed in terms of f and g and their derivatives. For this choice of v it turns out that f ? (t) u? 0 = = f (t) ? Cr (λ)u (t) + ACr (λ)u(t), [I ? Cr (λ)]u0 .
Assume also that f ∈ C p ([0, T ], X ), g ∈ C p+1 ([0, T ], Y ). Taking derivatives in equation (1.5) and using H1 as in Section 2 prove that f (j ) (t) ∈ Wr?j ?1 , 0 ≤ j ≤ r ? 1, 0 ≤ t ≤ T and then f ? (t) can again be computed in terms of the data and their derivatives. More precisely, f ? (t) = [I ? Cr (λ)]f (t) + (λ ? A)?(r?1) K (λ)? (λ ? A)r u(t), where ? (λ ? A)r u(t) is given by (2.5). If ω < 0, we can choose λ = 0 and the simpler expression ? ? f ? (t) = [I ? Cr (0)]f (t) ? A?(r?1) K (0) ?g (r) (t) ?
r ?1 j =0
? Ar?j ?1 f (j ) (t)?
results. Recall that the stability function r(z ) of the RK method is de?ned by r(z ) = ? > 0 such that, for 1 + z bT (I ? z A)?1 e. Hereafter, we suppose that there exists k s s s ? 0 < k ≤ k , the operator I ? A ? kA : D(A) ? X → X is invertible (see, e.g., ? the equations de?ning the RK approximations are [5]). Therefore, for 0 < k ≤ k solvable and the operator r(kA) is well de?ned in X . ? . Let us denote by u? the In the following theorem we assume that max km ≤ k n RK approximation to the solution of (3.1) at time level tn , 0 ≤ n ≤ N . Thus, we adopt un = u? n + v (tn ) as the numerical approximation to u(tn ). Theorem 3.2. Let u ∈ C p+1 ([0, T ], Wr ), r = p ? q , be the solution of the IBVP (1.5). Assume also that f ∈ C p ([0, T ], X ) and g ∈ C p+1 ([0, T ], Y ). Let v : [0, T ] → X be the correcting term de?ned in (3.5). Then
n
u(tn ) ? un = u? (tn ) ? u? n ≤ Cρn
m=1
p+1 km ?1 Im (u),
where C > 0 is a constant and
n
ρn
=
m=1
r(km?1 A) ,
q+1≤l≤p+1
Im (u) =
sup
u ( l)
L∞ ([tm?1 ,tm ],Wr )
.
The constant C depends on T , A and the RK method but it is independent of u, f and g .
AVOIDING THE ORDER REDUCTION FOR IBVPS
1535
Proof. By Lemma 1, I ? Cr (λ) : Wr → Wr is bounded. The norm of this operator is denoted by L. Since u? (t) = [I ? Cr (λ)]u(t), 0 ≤ t ≤ T , then u? ∈ C p+1 ([0, T ], Wr ) and u?(l) (t)
Wr
≤ L u(l) (t)
Wr
,
0 ≤ l ≤ p + 1. 0 ≤ j ≤ r ? 1,
Moreover, again by Lemma 1, we have ? (λ ? A)j u(l) (t) = ? (λ ? A)j Cr (λ)u(l) (t), whence ? (λ ? A)j u?(l) (t) = 0, Aj u?(l) (t) ≤ L u(l) (t) 0 ≤ j ≤ r ? 1, 0 ≤ j ≤ r, 0 ≤ l ≤ p + 1. so that u?(l) (t) ∈ Wr for 0 ≤ l ≤ p + 1, 0 ≤ t ≤ T . Thus, u? ∈ C p+1 ([0, T ], Wr ) and (3.6)
Wr
,
Now, the proof is straightforward by well-known results concerning RK methods for nonhomogeneous problems [5, 10] and by using (3.6). Remark 3.3. It makes sense to carry out corrections with 1 ≤ r < p ? q . Then the arguments used in the proof of Theorem 1 show that
n
(3.7)
u(tn ) ? un ≤ Cρn
m=1
q+r +1 km ?1 Im (u).
Furthermore the results in [5, 24] imply that if r(∞) = 1 and under natural stability requirements, then for ?xed step-sizes km = k
n
(3.8)
u(tn ) ? un ≤ Cρn k q+r+1
m=1
Im (u).
Thus, in this situation the full order p is achieved with r = p ? q ? 1. Remark 3.4. Fractional orders of convergence are also possible. Let us assume that there exists ? ∈ (0, 1) such that W1 ? W? , where W? is D((ω ? ? A)? ), for some ω ? > ω . Let us point out that this is the case indeed in many examples (see [4]). For instance in the example of equation (2.2), it turns out ? = 1/(2p), 1 ≤ p ≤ +∞. It is easy to prove that I ? Cr (λ) maps Wr+1 to Wr+? continuously. Then, according to the results in [5], an extra order ? is present in (3.7) and (3.8). Remark 3.5. If A is the in?nitesimal generator of a C0 -semigroup of contractions in a Hilbert space and the RK method is A-stable, then [16, 35] ρn ≤ etn max{0,ω} . For general C0 -semigroups in Banach spaces and ?xed step-sizes the bound √ ρn ≤ etn max{0,ω} O( n) is optimal [11, 14, 15] and half an order of convergence could be lost. In special situations this loss might not occur (see [11]). On the other hand, if A is the in?nitesimal generator of an analytic semigroup (see, e.g., [6, 14, 25, 26]), then ρn ≤ etn max{0,ω} O(1). This bound is valid even for variable step-sizes if |r(∞)| < 1 [8, 26].
1536
M. P. CALVO AND C. PALENCIA
Remark 3.6. From the de?nition of Cr (λ) it is clear that the computation of the ? initial condition u? 0 and each evaluation of f and the correcting term v require the solution of 2r ? 1 stationary problems as (2.1). In these stationary problems, derivatives of f and g are involved which, in the evaluation of v and f ? , may be replaced by appropriate approximations of order at least p. Notice also that the computation of v (t) is only necessary at the time levels where output is required. 4. Numerical illustration In this section we present two simple experiments to illustrate the results in Section 3. Let us consider the IBVP for the one-dimensional heat equation (2.2) in the example, with p = +∞, where f , u0 , g0 and g1 are chosen for u(x, t) = sin x/(1 + t) to be its exact solution. As mentioned at the beginning of Section 3, introducing the usual correcting term v (t) = Eg (t) = g0 (t) + x[g1 (t) ? g0 (t)], and de?ning u = u ? v , the original problem is transformed into the IBVP with homogeneous boundary conditions ? ? ut (x, t) = u? ? xx (x, t) + f (x, t) ? g0 (t) ? x[g1 (t) ? g0 (t)], ? ? ? u? (x, 0) = u (x) ? g (0) ? x[g (0) ? g (0)], 0 0 1 0 (4.1) ? ? (0 , t ) = 0 , u ? ? ? ? u (1, t) = 0. (x ∈ [0, 1], 0 ≤ t ≤ T ). For the spatial discretization of (4.1), given an integer J ≥ 2, we de?ne the uniform grid xj = jh, 1 ≤ j ≤ J , in the interval [0, 1], where h = 1/(J + 1). The second derivative is approximated using the standard three-point ?nite di?erences. For the time integration we use the singly diagonally implicit RK method proposed in [13] with Butcher array γ 1/2 1?γ where 1 π 1 + , γ = √ cos 18 2 3 δ= 1 . 6(2γ ? 1)2 γ 1/2 ? γ 2γ δ γ 1 ? 4γ 1 ? 2δ
?
γ δ
,
This method has order p = 4, stage order q = 1, is A-stable and |r(∞)| < 1. According to the results included in [5, 24] and, since u? ∈ D(A) and ? = 0 (see Remark 2), the expected order of the error of the time integrator is three. In our numerical experiments, to clearly observe the error due to the time integrator, we integrate the di?erential system (4.2) where ?? (t) = (Ph A ? Ah Ph )u? (x, t) + Ph f ? (t). f h
? ?? 0 ≤ t ≤ T, u? h (t) = Ah uh (t) + fh (t), ? ? ? uh (0) = [u (x1 , 0), . . . , u (xJ , 0)]T ,
AVOIDING THE ORDER REDUCTION FOR IBVPS
1537
Table 1. Errors for r = 1 h=k Error Order 5.000e-02 2.500e-02 1.250e-02 6.250e-03 3.125e-03 5.9752e-07 7.0579e-08 8.2533e-09 9.8621e-10 1.2009e-10 3.08 3.10 3.10 3.04 Table 2. Errors for r = 2 h=k Error Order 5.000e-02 2.500e-02 1.250e-02 6.250e-03 3.125e-03 5.5893e-08 5.5342e-09 4.6049e-10 3.4250e-11 2.3790e-12 3.34 3.59 3.75 3.85
Here for ? ∈ D(A), Ph ?(x, t) = [?(x1 , t), . . . , ?(xJ , t)]T and Ah is the three-point approximation to A. The exact solution of (4.2) is the restriction of the exact solution of (4.1) to the grid points, so that there is no spatial error. In Table 1 the global errors when h = k have been reported. We observe that when halving k the global errors are approximately divided by 8, as corresponding to the third order. This means that the order reduction is still present because the classical order of the RK method is four. In our second experiment we transform the original IBVP (2.2) into the IVP (3.1) by introducing the correcting term with r = 2 v = C2 (0)u = Eg + A?1 E [g ? ?f ]. Now u? = u ? v ∈ D(A2 ) and, according to Remark 1, the global error of the time integrator should exhibit the full order four. In Table 2 the global errors for this second experiment are presented. The quantities in the last row approach the full order four as k goes to zero, as expected. 5. The fully discrete problem When solving IBVPs arising in PDEs using the method of lines, ?rst it is necessary to discretize with respect to the spatial variable and then to time integrate the resulting system of ordinary di?erential equations. Then, in practice, we do not time integrate the IVP (3.1) but a discrete version of it. In the numerical ? experiments shown in Section 4 we have computed u? 0 and f exactly and then we have discretized with respect to the spatial variable. However, in general, it is not possible to compute (3.1) exactly because it would require the solution of certain stationary problems which, in most cases, must be solved numerically. In this section we prove that the result in Section 3 also applies to the fully discrete problem. As we mentioned in the introduction, one advantage of our approach is that standard spatial discretizations can be used. Let us denote by 0 < h < h0 the parameter of the spatial discretization and let (Xh , · ) be a family of normed spaces approximating X [17, 29]. The link between X and Xh is given by certain linear mappings Ph : X → Xh . For x ∈ X , Ph x is the approximation to x in Xh . The operator A : D(A) → X is approximated by linear operators Ah : Xh → Xh . As usual [17, 29], we assume that the semigroups etAh are uniformly bounded, i.e., etAh ≤ M eωt , t ≥ 0, for some M and ω independent of h. In this situation, the consistency of the spatial discretization is measured by (5.1) ?s (h) := (ω ? ? Ah )?1 (Ah Ph ? Ph A)
? s Ws →Xh ,
for some s > 0 and ω > ω . Recall that Ws = D(A ).
1538
M. P. CALVO AND C. PALENCIA
For instance, in the case of ?nite elements of degree 2s ? 1 and a second-order elliptic operator A, it turns out that ?s (h) = O(h2s ) in the L2 -norm [36]. If a ?nite di?erence scheme of order 2s ? 2 is used instead, again for an elliptic operator A of second order, we have (Ah Ph ? Ph A)
Ws →Xh
= O(h2s?2 )
and, since (ω ? ? Ah )?1 ≤ M/(ω ? ? ω ), then ?s (h) = O(h2s?2 ). We also assume that ?s (h) balances the behaviour of Ah in the sense that (5.2)
?1 (Ah Ph ? Ph A) Am h Wm →Xh
≤ λ,
1 ≤ m ≤ s,
for some λ > 0 independent of h. This condition is clearly satis?ed in the context of either ?nite elements or ?nite di?erences mentioned above, where Ah = O(h?2 ). In the present context of IBVPs also the operators Cr (λ) must be approximated. In practice, this can be done by using either boundary elements [9, 18] or ?nite di?erences [22, 23] for discretizing K (λ)? . Thus, we assume that Ph f ? (t), Ph v (t) ? ? and Ph u? 0 are approximated by fh (t), vh (t) and uh,0 in Xh . For solutions u taking values in Ws we denote by ψs (h) the error of these approximations, i.e.,
? (t) Ph f ? (t) ? fh
= O(ψs (h)), = O(ψs (h)), = O(ψs (h)),
Ph u? 0
?
u? h,0
Ph v (t) ? vh (t)
as h → 0+, uniformly in 0 ≤ t ≤ T . After the spatial discretization of (3.1), we have to time integrate the di?erential system (5.3)
? ? u? h (t) = Ah uh (t) + fh (t), ? ? uh (0) = uh,0 .
Let us denote by u? h,n the RK approximation to the solution of (5.3) at time level tn = nk . The fully discrete approximation to u(tn ) is thus de?ned by uh,n = u? h,n + vh (tn ). ? > 0 and ρ = ρ(T ) > 0 The fully discrete scheme is said to be stable if there exist k ? such that, for 0 < k ≤ k, the operator I ? A ? kAh : Xh → Xh is invertible, 0 < h < h0 , and r(kAh )n ≤ ρ, 0 ≤ nk ≤ T. Under these assumptions we can establish the result for the full discretization. Theorem 5.1. Let u ∈ C p+1 ([0, T ], Wr ), r = p ? q , be the solution of the IBVP (1.5). Assume also that f ∈ C p ([0, T ], X ), g ∈ C p+1 ([0, T ], Y ). Assume that (5.2) is ?, satis?ed for s = r and that the fully discrete scheme is stable. Then, for 0 < k ≤ k Ph u(tn ) ? uh,n ≤ Ctn [k p + ?r (h) + ψr (h)], where C = C (u, λ, ρ). Proof. For simplicity in the proof we suppose ω < 0. ? ? ? Set e? h,n = Ph u (tn ) ? uh,n . From the de?nition of uh,n and the relation u (t) = u(t) ? v (t) it is clear that Ph u(tn ) ? uh,n ≤ e? h,n + Ph v (tn ) ? vh (tn ) .
AVOIDING THE ORDER REDUCTION FOR IBVPS
1539
The second term is bounded by O(ψr (h)) and only e? h,n must be analyzed. Let u ?? be the RK approximation to the solution of h,n+1 (5.4)
? (t) xh (t) = Ah xh (t) + fh
at tn+1 , starting from xh (tn ) = Ph u? (tn ). Then,
? ? ?? ?? e? h,n+1 = u h,n+1 ? uh,n+1 + Ph u (tn+1 ) ? u h,n+1 .
The ?rst two terms on the right-hand side are numerical solutions of (5.4) after one step with initial conditions Ph u? (tn ) and u? h,n respectively, so that (5.5)
? ? ?? e? h,n+1 = r(kAh )eh,n + Ph u (tn+1 ) ? u h,n+1 .
Next, in order to use the variation-of-constants formula, it is su?cient to study ? ? s ? ?? the truncation term Ph u? (tn+1 ) ? u h,n+1 . We set F (t) = (f (t + ci k ))i=1 , U (t) = ? s ? ? s (u (t + ci k ))i=1 and U (t) = (u (t + ci k ))i=1 for the di?erential equation (3.1) ? ? (t) = (fh (t + ci k ))s and Fh i=1 for (5.4). Direct substitution of Ph u? into the equations de?ning the stages of the RK approximation to the solution of (5.4) de?nes the defects τh,n by (5.6) ? (tn ) + τh,n , Ph U ? (tn ) = (e ? I )Ph u? (tn ) + (A ? kAh )Ph U ? (tn ) + k (A ? I )Fh which can also be written as Ph U ? (tn ) = + Ph (e ? I )u? (tn ) + k A ? U ? (tn )
? ? Ph F ? )(tn ) + τh,n . k A ? (Ah Ph ? Ph A)U ? (tn ) + k (A ? I )(Fh
Taylor expansions of U ? and U ? at tn show
p
U ? (tn ) ? (e ? I )u? (tn ) ? k A ? U ? (tn ) =
l=q+1
kl (wl ? I ) u?(l) (tn ) + O(k p+1 ) l!
(recall that wl = cl ? lAcl?1 , for q + 1 ≤ l ≤ p). Then
p
(5.7)
τh,n
=
l=q+1
kl (wl ? I ) Ph u?(l) (tn ) ? k A ? (Ah Ph ? Ph A)U ? (tn ) l!
?
? ? Ph F ? )(tn ) + O(k p+1 ). k (A ? I )(Fh
? ? for the RK approximation to the On the other hand, the internal stages U h,n solution of (5.4) starting from Ph u? (tn ) satisfy (5.8) ? ? = (e ? I )Ph u? (tn ) + (A ? kAh )U ? ? + k (A ? I )F ? (tn ). U h,n h,n h
? ? ?h,n ?h,n = (A ? kAh ) Ph U ? (tn ) ? U + τh,n , Ph U ? (tn ) ? U
Substracting (5.8) from (5.6) we get
which using (5.7) leads to
p ? ?h,n Ph U ? (tn ) ? U
= (I ? A ? kAh )?1
l=q+1
kl (wl ? I ) Ph u?(l) (tn ) l!
(5.9)
? (I ? A ? kAh )?1 k A ? (Ah Ph ? Ph A)U ? (tn )
? ? Ph F ? )(tn ) + O(k p+1 ). ? (I ? A ? kAh )?1 k (A ? I )(Fh
1540
M. P. CALVO AND C. PALENCIA
Let us now consider the equation that de?nes the RK approximation to the solution of (5.4) at time level tn+1 . Direct substitution of Ph u? into this equation de?nes the defects δh,n+1 by (5.10) ? (tn ) + δh,n+1 , Ph u? (tn+1 ) = Ph u? (tn ) + k (bT ? Ah )Ph U ? (tn ) + k (bT ? I )Fh which can be written as Ph u? (tn+1 ) = + Ph u? (tn ) + k (bT ? I )U ? (tn ) + k (bT ? I )(Ah Ph ? Ph A)U ? (tn )
? ? Ph F ? )(tn ) + δh,n+1 . k (bT ? I )(Fh
As u? (tn+1 ) ? u? (tn ) ? k (bT ? I )U ? (tn ) = O(k p+1 ), then (5.11) δh,n+1 = ? ?k (bT ? I )(Ah Ph ? Ph A)U ? (tn )
? ? Ph F ? )(tn ) + O(k p+1 ). k (bT ? I )(Fh
? For the numerical solution u ?? h,n+1 of (5.4) starting from Ph u (tn ) we have
(5.12)
? T T ? ?? u ?? h,n+1 = Ph u (tn ) + k (b ? Ah )Uh,n + k (b ? I )Fh (tn ).
Substracting (5.12) from (5.10) we get
T ? ?? ?? Ph u? (tn+1 ) ? u h,n+1 = (b ? kAh ) Ph U (tn ) ? Uh,n + δh,n+1 ,
which, after using (5.9) and (5.11) and some manipulations, leads to (5.13) where
p n Kh n n n n p+1 ?? ), Ph u? (tn+1 ) ? u h,n+1 = Kh + Lh + Mh + Nh + O(k
= = =
n Nh
(bT ? kAh )(I ? k A ? Ah )?1
l=q+1
kl (wl ? I ) Ph u?(l) (tn ), l!
Ln h
n Mh
?k b (I ? k A ? Ah )
T
?1
? (Fh
? Ph F ? )(tn ),
?k bT (I ? k A ? Ah )?1 [I ? (Ah Ph ? Ph A)] (U ? (tn ) ? (e ? I )u? (tn )) , = = ?k bT (I ? k A ? Ah )?1 [I ? (Ah Ph ? Ph A)](e ? I )u? (tn )
1 ? (I ? r(kAh )) A? h (Ah Ph ? Ph A)(e ? I )u (tn ).
and
Now the variation-of-constants formula applied to (5.5) reads (5.14)
n n ? ? e? h,n = r(kAh ) [Ph u0 ? uh,0 ] + j =1 j j j p r(kAh )n?j [Kh + Lj h + Mh + Nh ] + O(k ).
The ?rst term on the right-hand side of (5.14) can be bounded by O(ρψr (h)). j turns out to be O(k p+1 ). To see this, notice that the equality The term Kh
l ?( l ) = Ph Al u?(l) + Al h Ph u m=1 ?( l ) for 1 ≤ l ≤ r. together with (5.1) and (5.2), implies the boundedness of Al h Ph u j p+1 ). Then, the contribution The proof of Theorem 1 in [5] shows that Kh = O(k of these terms to the error is O(ρtn k p ). ?1 l ? m ?( l ) Am u , h Ah (Ah Ph ? Ph A)A
AVOIDING THE ORDER REDUCTION FOR IBVPS
1541
On the other hand, Lj h behaves as kO(ψr (h)) and, since
1 ? ? [I ? A? h (Ah Ph ? Ph A)] (U (tn ) ? (e ? I )u (tn )) = kO(?r (h)) j behaves as and the operator bT (I ? A ? kAh )?1 (kAh ) is bounded, the term Mh k?r (h). Therefore, these terms contribute to the error with
ρtn O(ψr (h)) + O(?r (h)) . It remains to analyze
n
(5.15)
j =1
j r(kAh )n?j Nh
which, after summation by parts, can be written as
n?1 j =1 1 ? n ?1 ? +A? h (Ah Ph ? Ph A)u (tn ) ? r(kAh ) Ah (Ah Ph ? Ph A)u (t1 ). 1 ? ? r(kAh )n?j A? h (Ah Ph ? Ph A)(e ? I ) (u (tj ) ? u (tj +1 ))
As long as (5.1) is satis?ed and u? (tj ) ? u?(tj +1 ) = O(k ), then the norm of (5.15) is ρtn O(?r (h)). Remark 5.2. When using ?nite elements in the context of the maximum norm, it ? ), where s ? = s + 1, turns out that the norms in (5.2) behave like λ(h) = O(|lnh|s for piecewise linear elements, or s ? = s, for higher order elements [33]. This is not a serious problem since λ = λ(h) enters as a factor in the error constant C in Theorem 5.1. Remark 5.3. For simplicity we have assumed that the fully discrete scheme is stable. In practice, ρ might depend on h and k . This dependence has been widely studied in the literature (see references in Remark 3.5 and, for second-order elliptic operators and the maximum norm, [7] for ?nite di?erences and [27, 34, 36] for ?nite elements). Typically ρ exhibits at most weak singularities when h and k go to 0. Thus, as ρ enters again as a factor in C , the full order of convergence could be slightly spoiled, due to stability reasons. Remark 5.4. If r(∞) = 1, then the optimal order p may be achieved with r = p ? q ? 1 (see Remark 3.3). Remark 5.5. Fixed step-sizes are required for the treatment of (5.15), when summation by parts applies. However, in cases for which (Ah Ph ? Ph A)u? (t) = O(?r (h))
n n + Nh (for instance when ?nite di?erences are used), the joint contribution Mh can be analyzed without using summation by parts. In these situations variable step-sizes are also allowed.
Remark 5.6. If W1 ? W? , 0 < ? < 1, and r < p ? q , then a fractional order in time occurs as in Remark 3.4. In fact, by interpolating the two inequalities Ah Ph ξ ≤ λ ξ + Ph Ph ξ ≤ Ph ξ , Aξ ,
1542
M. P. CALVO AND C. PALENCIA
ξ ∈ W1 , we deduce that when Ph = O(1) (a standard situation), (ω ? Ah )? Ph ξ ≤ C? ξ
?,
ξ ∈ W? , for all ? < ?. Thus, an extra order ? is present for all ? < ?. This explains why the extra order ? is observed in full discretizations of (3.2) (see, e.g., the numerical experiments in [24]). References
1. S. Abarbanel, D. Gottlieb & M. H. Carpenter, On the removal of boundary errors caused by Runge-Kutta integration of nonlinear partial di?erential equations, SIAM J. Sci. Comput. 17 (1996) 777-782. MR 96m:65071 2. I. Alonso-Mallo, Rational methods with optimal order of convergence for partial di?erential equations, Appl. Numer. Math. 35 (2000) 265–292. CMP 2001:05 3. I. Alonso-Mallo, Explicit single step methods with optimal order of convergence for partial di?erential equations, Appl. Numer. Math. 31 (1999) 117–131. MR 2000d:65186 4. I. Alonso-Mallo & C. Palencia, On the convolution operators arising in the study of abstract initial boundary value problems, Proc. Royal Soc. Edinburgh 126A (1996) 515–539. MR 97d:34068 5. I. Alonso-Mallo & C. Palencia, Optimal orders of convergence in Runge-Kutta methods for linear, non-homogeneous PDEs with singular source terms, preprint. 6. N. Yu. Bakaev, On the stability of some general discretization methods, Dokl. Akad. Nauk. SSSR 309 (1989) 11–15; English transl. in Soviet Math. Dokl. 40 (1990). MR 91e:65104 7. N. Yu. Bakaev, Estimates for the resolvent of a multi-dimensional elliptic di?erence operator, Function.-Di?. Eqs., Perm’ PPI (1991) 118–126. MR 94e:39001 8. N. Yu. Bakaev, On variable stepsize Runge-Kutta approximations of a Cauchy problem for the evolution equation, BIT 38 (1998) 462–485. MR 99i:65069 9. C. A. Brebbia, J. C. F. Telles & L. C. Wrobel, Boundary element techniques. Theory and applications in engineering, Springer, Berlin, 1984. MR 89i:65002 10. P. Brenner, M. Crouzeix & V. Thom? ee, Single step methods for inhomogeneous linear di?erential equations in Banach space, R.A.I.R.O. Anal. Numer. 16 (1982) 5–26. MR 83d:65268 11. P. Brenner & V. Thom? ee, On rational approximations of semigroups, SIAM J. Numer. Anal. 16 (1979) 683–694. MR 80j:47052 12. M. H. Carpenter, D. Gottlieb, S. Abarbanel & W. S. Don, The theoretical accuracy of RungeKutta time discretizations for the initial boundary value problem: a study of the boundary error, SIAM J. Sci. Comput. 16 (1995) 1241–1252. MR 96h:65088 13. M. Crouzeix & P. A. Raviart, M? ethodes de Runge-Kutta, Unpublished Lecture Notes, Universit? e de Rennes (1980). 14. M. Crouzeix, S. Larsson, S. Piskarev & V. Thom? ee, The stability of rational approximations of holomorphic semigroups, BIT 33 (1993) 74–84. MR 96f:65069 15. J. L. M. van Dorsselaer, J. F. B. M. Kraaijenvanger & M. N. Spijker, Linear stability analysis in the numerical solution of initial value problems, in Acta Numer. 1993, 199–237. MR 94e:65051 16. E. Hairer & G. Wanner, Solving ordinary di?erential equations II, Sti? and di?erentialalgebraic problems, 2nd Ed., Springer, Berlin, 1996. MR 97m:65007 17. K. Ito & F. Kappel, The Trotter-Kato theorem and approximations of PDEs, Math. Comp. 67 (1998) 21–44. MR 98e:47060 18. C. Johnson, Numerical solution of partial di?erential equations by the ?nite element method, Cambridge University Press, Cambridge, 1987. MR 89b:65003a 19. S. L. Keeling, Galerkin/Runge-Kutta discretizations for parabolic equations with timedependent coe?cients, Math. Comp. 52 (1989) 561–586. MR 90a:65239 20. S. L. Keeling, Galerkin/Runge-Kutta discretizations for semilinear parabolic equations, SIAM J. Numer. Anal. 27 (1990) 394–418. MR 91d:65140 21. Ch. Lubich & A. Ostermann, Interior estimates for time discretizations of parabolic equations, Appl. Numer. Math. 18 (1995) 241–251. MR 96f:65124 22. A. R. Mitchell & D. F. Gri?ths, The ?nite di?erence method in partial di?erential equations, Wiley-Interscience, New York, 1980. MR 82a:65002
AVOIDING THE ORDER REDUCTION FOR IBVPS
1543
23. K. W. Morton & D. F. Mayers, Numerical solution of partial di?erential equations, Cambridge University Press, Cambridge, 1994. MR 96b:65001 24. A. Ostermann & M. Roche, Runge-Kutta methods for partial di?erential equations and fractional orders of convergence, Math. Comp. 59 (1992) 403–420. MR 93a:65125 25. C. Palencia, A stability result for sectorial operators in Banach spaces, SIAM J. Numer. Anal. 30 (1993) 1373–1384. MR 94j:65109 26. C. Palencia, On the stability of variable stepsize rational approximations of holomorphic semigroups, Math. Comput. 62 (1994) 93–103. MR 94c:47066 27. C. Palencia, Maximum-norm analysis of completely discrete ?nite element methods for parabolic problems, SIAM J. Numer. Anal. 33 (1996) 1654–1668. MR 97e:65099 28. C. Palencia & I. Alonso-Mallo, Abstract initial boundary value problems, Proc. Royal Soc. Edinburgh 124A (1994) 879–908. MR 95k:35090 29. C. Palencia & J. M. Sanz-Serna, An extension of the Lax-Richtmyer theory, Numer. Math. 44 (1984) 279–283. MR 86c:65096 30. D. Pathria, The correct formulation of intermediate boundary conditions for Runge-Kutta time integration of initial boundary value problems, SIAM J. Sci. Comput. 18 (1997) 1255– 1266. MR 98d:65100 31. A. Pazy, Semigroups of linear operators and applications to partial di?erential equations, Springer, Berlin, 1983. MR 85g:47061 32. J. M. Sanz-Serna, J. G. Verwer & W. H. Hundsdorfer, Convergence and order reduction of Runge-Kutta schemes applied to evolutionary problems in partial di?erential equations, Numer. Math. 50 (1986) 405–418. MR 88f:65146 33. A. H. Schatz & L. B. Wahlbin, On the quasi-optimality in L∞ of the H1 -projection into ?nite element spaces, Math. Comp. 38 (1982) 1–22. MR 82m:65106 34. A. H. Schatz, V. Thom? ee & L. B. Wahlbin, Stability, analyticity, and almost best approximation in maximum norm for parabolic ?nite element equations, Comm. Pure Appl. Math. 51 (1998) 1349–1385. MR 99h:65171 35. B. Sz.-Nagy & C. Foias, Harmonic Analysis of Operators on Hilbert Spaces, North-Holland, Amsterdam, 1970. MR 43:947 36. V. Thom? ee, Galerkin ?nite element methods for parabolic problems, Springer, Berlin, 1997. MR 98m:65007 ? tica Aplicada y Computacio ? n, Universidad de Valladolid, Departamento de Matema Valladolid, Spain E-mail address : maripaz@mac.cie.uva.es ? tica Aplicada y Computacio ? n, Universidad de Valladolid, Departamento de Matema Valladolid, Spain E-mail address : palencia@mac.cie.uva.es