Order reduction and how to avoid it when Lawson methods integrate reaction-diffusion boundary value problems
aa r X i v : . [ m a t h . NA ] S e p Order reduction and how to avoid it when Lawsonmethods integrate reaction-diffusion boundary valueproblems
B. Cano ∗ IMUVA, Departamento de Matem´atica Aplicada,Facultad de Ciencias, Universidad de Valladolid,Paseo de Bel´en 7, 47011 Valladolid,Spain andN. Reguera † IMUVA, Departamento de Matem´aticas y Computaci´on,Escuela Polit´ecnica Superior, Universidad de Burgos,Avda. Cantabria, 09006 Burgos,Spain
Abstract
It is well known that Lawson methods suffer from a severe order reductionwhen integrating initial boundary value problems where the solutions are notperiodic in space or do not satisfy enough conditions of annihilation on theboundary. However, in a previous paper, a modification of Lawson quadraturerules has been suggested so that no order reduction turns up when integratinglinear problems subject to even time-dependent boundary conditions. In thispaper, we describe and thoroughly analyse a technique to avoid also order reduc-tion when integrating nonlinear problems. This is very useful because, given anyRunge-Kutta method of any classical order, a Lawson method can be constructedassociated to it for which the order is conserved.
There is an effort in the recent literature to understand the order reduction which turnsup when integrating initial boundary value problems subject to non-periodic boundary ∗ Corresponding author. Email: [email protected] † Email: [email protected] ≥
2, global order 2 is achieved if the summation-by-parts argument is applied. (For that, apart from enough regularity, we need theassumptions (22),(34),(41), which seem to be true but which proof is out of the scopeof this paper.) Moreover, by resorting to numerical differentiation, local order 3 canbe achieved and, whenever a CFL condition is satisfied (61), global error ≥ γ = 1 in (61)and therefore it just means that the timestepsize is not too big with respect to thespace grid.) Finally, although the formulas get more complicated, by using numericaldifferentiation, local order 4 can also be achieved if the underlying RK method is of2rder ≥
3. We have not described the formulas to get local order ≥ Let X and Y be Banach spaces and let A : D ( A ) ⊂ X → X and ∂ : X → Y be linear operators. Our goal is to avoid order reduction when integrating in timethrough Lawson methods the nonlinear abstract non homogeneous initial boundaryvalue problem u ′ ( t ) = Au ( t ) + f ( t, u ( t )) , ≤ t ≤ T,u (0) = u ∈ X,∂u ( t ) = g ( t ) ∈ Y, ≤ t ≤ T. (1)We will assume that the functions f : [0 , T ] × X → X (in general nonlinear) and g : [0 , T ] → Y are regular enough.The abstract setting (1) permits to cover a wide range of nonlinear evolutionaryproblems governed by partial differential equations. We use the following hypotheses,similar to the ones in [4] when avoiding the same kind of problems with exponentialsplitting methods.(A1) The boundary operator ∂ : D ( A ) ⊂ X → Y is onto and g ∈ C ([0 , T ] , Y ).(A2) Ker( ∂ ) is dense in X and A : D ( A ) = ker( ∂ ) ⊂ X → X , the restriction of A toKer( ∂ ), is the infinitesimal generator of a C - semigroup { e tA } t ≥ in X , whichtype ω is assumed to be negative.(A3) If z ∈ C satisfies ℜ ( z ) > v ∈ Y , then the steady state problem Ax = zx, (2) ∂x = v, (3)3ossesses a unique solution denoted by x = K ( z ) v . Moreover, the linear operator K ( z ) : Y → D ( A ) satisfies k K ( z ) v k ≤ C k v k , (4)where the constant C holds for any z such that Re ( z ) ≥ ω > ω .(A4) The nonlinear source f belongs to C ([0 , T ] × X, X ).(A5) The solution u of (1) satisfies u ∈ C ([0 , T ] , X ), u ( t ) ∈ D ( A ) for all t ∈ [0 , T ] and Au ∈ C ([0 , T ] , X ).In the remaining of the paper, we always suppose that (A1)-(A5) are satisfied.However, we notice that we also assume more regularity in certain results. As for thewell-posedness of problem (1), the same remark as in [4] can be made. We repeat ithere for the sake of clarity. Remark 1.
From hypotheses (A1)-(A4), the problem (1) with homogeneous boundaryconditions has a unique classical solution for small enough time intervals (see Theorem6.1.5 in [15]).Regarding the nonhomogeneous case, as in (A1) we also assume that g ∈ C ([0 , T ] , Y ) ,we can look for a solution of (1) given by: u ( t ) = v ( t ) + K ( z ) g ( t ) , t ≥ , for some fixed ℜ ( z ) > ω . Then, v is solution of an IBVP with vanishing boundaryvalues similar to the one in [15] and the well-posedness for the case of nonhomoge-neous boundary values is a direct consequence if we take the abstract theory for initialboundary value problems in [5, 14] into account.However, condition (A4) may be very strong. When X is a function space with anorm L p , ≤ p < + ∞ , and f is a Neminskii operator, u → f ( u ) = φ ( u ) , with φ : C → C , (A4) implies that φ is globally Lipschitz in C . This objection disappearsby considering the supremum norm, which is used in our numerical examples, wherethe nonlinear source is given by u → f ( t, u ) = φ ( u ) + h ( t ) , (5) with h : [0 , T ] → X , that is, f is the sum of a Neminskii operator and a linear term.In this way, problem (1) is well posed. Because of hypothesis (A2), { ϕ j ( tA ) } j =1 are bounded operators for t >
0, where { ϕ j } are the standard functions which are used in exponential methods [12] and whichare defined by ϕ j ( tA ) = 1 t j Z t e ( t − τ ) A τ j − ( j − dτ, j ≥ . (6)4t is well-known that they can be calculated in a recursive way through the formulas ϕ j +1 ( z ) = ϕ j ( z ) − /j ! z , z = 0 , ϕ j +1 (0) = 1( j + 1)! , ϕ ( z ) = e z . (7)For the time integration, we will center on exponential Lawson methods whichare determined by an explicit Runge-Kutta tableau and which, when applied to thefinite-dimensional nonlinear problem like U ′ ( t ) = M U ( t ) + F ( t, U ( t )) , (8)where M is a matrix, read like this at each step K n,i = e c i kM U n + k i − X j =1 a ij e ( c i − c j ) kM F ( t n + c j k, K n,j ) , i = 1 , . . . , s, (9) U n +1 = e kM U n + k s X i =1 b i e (1 − c i ) kM F ( t n + c i k, K n,i ) , (10)where k > t n = t + nk .Following the example in Section 2 of [4], we take X = C (Ω) for a certain boundeddomain Ω ∈ R d . There, we consider the maximum norm and a certain grid Ω h (ofΩ) over which the approximated numerical solution will be defined. In this way, thisnumerical approximation belongs to C N , where N is the number of nodes in the grid,and the maximum norm k u h k h = k [ u , . . . , u N ] T k h = max ≤ i ≤ N | u i | is considered.Notice that, usually, when considering Dirichlet boundary conditions, nodes on theboundary are not considered while, when using Neumann or Robin boundary condi-tions, the nodes on the boundary are taken into account.In that sense, we consider the projection operator P h : X → C N , (11)which takes a function to its values over the grid Ω h . On the other hand, the operator A , when applied over functions which satisfy a certain condition on the boundary ∂u = g , is discretized by A h, U h + C h g, where A h, is the matrix which discretizes A and C h : Y → C N is another operator,which is the one which contains the information on the boundary.We also assume that the source function f has also sense as function from [0 , T ] × C N on C N and, for each t ∈ [0 , T ] and u ∈ X , P h f ( t, u ) = f ( t, P h u ) . (12)This fact is obvious when f is given by (5).In a similar way to [4], we consider the following hypotheses:5H1) The matrix A h, satisfies(a) k e tA h, k h ≤ A h, is invertible and k A − h, k h ≤ C for some constant C which does notdepend on h .(H2) We define the elliptic projection R h : D ( A ) → C N as the solution of A h, R h u + C h ∂u = P h Au. (13)We assume that there exists a subspace Z ⊂ D ( A ), such that, for u ∈ Z ,(a) A − u ∈ Z and e tA u, f ( t, u ) ∈ Z , for t ∈ [0 , T ].(b) for some ε h and η h which are both small with h , k A h, ( P h u − R h u ) k ≤ ε h k u k Z , k P h u − R h u k ≤ η h k u k Z . (14)(Although obviously, because of (H1b), η h could be taken as Cε h , for somediscretizations η h can decrease more quickly with h than ε h and that leadsto better error bounds in the following sections.)(c) k A − h, C h k h ≤ C ′′ for some constant C ′′ which does not depend on h . Thisresembles the continuous maximum principle which is satisfied because of(4) when z = 0.(H3) The nonlinear source f belongs to C ([0 , T ] × C N , C N ) and the derivative withrespect to the variable in C N is uniformly bounded in a neighbourhood of thesolution where the numerical approximation stays.As in [4], hypothesis (H1a) can be deduced in our numerical experiments by usingthe logarithmic norm of matrix A h, . It was already proved in [2] that Lawson methods, even when applied to linear problemswith vanishing boundary conditions, show in general just order 1 in time. For thatreason, with nonlinear problems and probably non-vanishing boundary conditions, ingeneral we cannot expect more than that order. In any case, in this section we generalisethe main result in [2] not only in the sense of considering more general problems butalso in the sense of taking the error coming from the space discretization into account.Notice that, by using the space discretization which is described in the previoussection, the following semidiscrete problem arises after discretising (1), U ′ h ( t ) = A h, U h ( t ) + C h g ( t ) + f ( t, U h ( t )) ,U h ( t ) = P h u ( t ) . (15)6hen, applying Lawson method (9)-(10) to this, the following formulas define one stepfrom U nh to U n +1 h : K nh,i = e c i kA h, U nh + k i − X j =1 a ij e ( c i − c j ) kA h, [ C h g ( t n + c j k ) + f ( t n + c j k, K nh,j )] , i = 1 , . . . , s,U n +1 h = e kA h, U nh + k s X i =1 b i e (1 − c i ) kA h, [ C h g ( t n + c i k ) + f ( t n + c i k, K nh,i )] . (16)We define the local error as ρ h,n +1 = ¯ U n +1 h − P h u ( t n +1 ), where u ( t ) is the solution of(1) and ¯ U n +1 h is deduced as U n +1 h but starting from P h u ( t n ) instead of U nh . Then,¯ U n +1 h = e kA h, P h u ( t n ) + k s X i =1 b i e (1 − c i ) kA h, [ C h g ( t n + c i k ) + f ( t n + c i k, ¯ K nh,i )] , (17)where, for i = 1 , . . . , s ,¯ K nh,i = e c i kA h, P h u ( t n ) + k i − X j =1 a ij e ( c i − c j ) kA h, [ C h g ( t n + c j k ) + f ( t n + c j k, ¯ K nh,j )] . (18) g ( t ) ≡ In this subsection, we will analyse the error with the classical approach under theassumption that the boundary conditions vanish, which is the only case which has beenstudied for other standard exponential methods when integrating nonlinear problemsand where just the error coming from the time integration has been analysed [11]. Asdistinct, here we also consider the error coming from the space discretization.
Theorem 2.
Under hypotheses (A1)-(A5) and (H1)-(H3), whenever u ∈ C ([0 , T ] , Z ) and ∂u ≡ , ρ h,n = ¯ U n +1 h − P h u ( t n +1 ) = O ( k ) , where the constant in Landau notation is independent of k and h . Moreover, if P si =1 b i =1 and u ∈ C ([0 , T ] , X ) , it happens that A − h, ρ h,n = O ( η h k + k ) . Proof.
Firstly we notice that, for i = 1 , . . . , s , ¯ K nh,i in (18) are uniformly bounded on h when g ( t ) ≡
0, which can be proved by induction on i by using (H1a), (A4) and (A5).Then, by using (7),¯ U n +1 h = P h u ( t n ) + kA h, ϕ ( kA h, ) P h u ( t n ) + k s X i =1 b i e (1 − c i ) kA h, f ( t n + c i k, ¯ K nh,i )= P h u ( t n ) + O ( k ) , A h, ϕ ( kA h, ) P h u ( t n ) is uniformly bounded on h because A h, ϕ ( kA h, ) P h u ( t n ) = ϕ ( kA h, ) A h, R h u ( t n ) + ϕ ( kA h, ) A h, ( P h − R h ) u ( t n )= ϕ ( kA h, ) P h Au ( t n ) + O ( ε h ) . (19)(Here, the second equality comes from the fact that A h, R h u ( t ) = P h Au ( t ) due to (13)with ∂u = 0, and also to (14) considering that u ∈ C ([0 , T ] , Z ).)As for the second result, notice that, by using (7) again and an argument similarto (19) , ¯ K nh,i can also be written as¯ K nh,i = P h u ( t n ) + c i kA h, ϕ ( c i kA h, ) P h u ( t n ) + O ( k ) = P h u ( t n ) + O ( k ) . (20)Then, using now (7) to expand e kA h, till ϕ ( kA h, ) and e (1 − c i ) kA h, till ϕ ((1 − c i ) kA h, )and again an argument similar to (19) for A h, ϕ ( kA h, ) P h u ( t n ), A − h, ρ h,n +1 = A − h, (cid:20) P h u ( t n ) + kA h, P h u ( t n ) + k A h, ϕ ( kA h, ) P h u ( t n )+ k s X i =1 b i [ f ( t n + c i k, ¯ K nh,i ) + (1 − c i ) kA h, ϕ ((1 − c i ) kA h, ) f ( t n + c i k, ¯ K nh,i )] − P h u ( t n ) − kP h ˙ u ( t n ) + O ( k ) (cid:21) = k [ P h u ( t n ) − A − h, P h ˙ u ( t n ) + A − h, f ( t n , P h u ( t n ))] + O ( k ) , (21)where (H1b) and (20) have been used as well as the fact that P si =1 b i = 1. Using nowthat A − h, P h ˙ u ( t ) = A − h, P h [ Au ( t ) + f ( t, u ( t ))] = R h u ( t ) + A − h, P h f ( t, u ( t )) , the bracket in (21) is O ( η h ) according to (14) and using (12), from what the resultfollows.Using the classical argument for the global error, the first result of the previoustheorem would not lead to convergence. However, by using the second result and afew more assumptions, a summation-by-parts argument very similar to that given in[4] for Strang method and parabolic problems when avoiding order reduction does leadto the following result. Theorem 3.
Under the hypotheses of Theorem 2, and assuming also that u ∈ C ([0 , T ] , X ) ∩ C ([0 , T ] , Z ) , with ˙ u ( t ) ∈ D ( A ) for t ∈ [0 , T ] , ˙ u ∈ C ([0 , T ] , Z ) , A ˙ u ∈ C ([0 , T ] , X ) and k kA h, n − X r =1 e rkA h, k h ≤ C, ≤ nk ≤ T, (22) it happens that e h,n = U nh − P h u ( t n ) = O ( η h + k ) . .2 Local and global error when g ( t ) = 0 In this subsection we study the more general case in which g ( t ) = 0. For the sake ofbrevity, we will consider just the case in which all the nodes c i are strictly increasingwith the index i , but we will distinguish between the case in which c i = 1 for i = 1 , . . . , s and the case in which some of those c i is equal to 1. Both types of results show a verypoor behaviour, but our aim is just to explain the differences between them. Theorem 4.
Under hypotheses (A1)-(A5) and (H1)-(H3), assuming also that ∂u = 0 , c < c < · · · < c s , c i = 1 for i = 1 , . . . , s and that there exists C ′ , h such that k τ A h, e τA h, k h ≤ C ′ , τ ≥ , h ≤ h , (23) it happens that ρ h,n = O (1) , A − h, ρ h,n = O ( k ) , where the constant in Landau notation is independent of k and h .Proof. Firstly we notice that now the stages ¯ K nh,i in (18) are also uniformly boundedbecause, as the nodes are different, ke ( c i − c j ) kA h, C h g ( t n + c j k ) = 1( c i − c j ) [( c i − c j ) kA h, e ( c i − c j ) kA h, ] A − h, C h g ( t n + c j k ) . Then, the bracket is bounded because of (23) and the last factor because of (H2c).With the same argument, as c i = 1, ¯ U n +1 h in (17) is also uniformly bounded, whichimplies that ρ h,n is just O (1). On the other hand, by using (7) in the first term of (17), A − h, ρ h,n +1 = A − h, (cid:20) P h u ( t n ) + kA h, ϕ ( kA h, ) P h u ( t n )+ k s X i =1 b i e (1 − c i ) kA h, [ f ( t n + c i k, ¯ K nh,i ) + C h g ( t n + c i k )] − P h u ( t n ) + O ( k ) (cid:21) , which again is O ( k ) because of (H2c).Now, with the same summation-by-parts argument as in [4], the following resultfollows for the global error. Theorem 5.
Under the same hypotheses of Theorem 4, if u ∈ C ([0 , T ] , X ) and (22)holds, then e h,n = U nh − P h u ( t n ) = O (1) . Let us now consider the case in which some c i = 1. Theorem 6.
Under hypotheses (A1)-(A5) and (H1)-(H3), assuming also that ∂u = 0 , c < c < · · · < c s , c i = 1 for some i ∈ { , . . . , s } and (23), it happens that ρ h,n = O (1 + k k C h k ) , A − h, ρ h,n = O ( k ) , where the constant in Landau notation is independent of k and h . roof. The proof is the same as that of Theorem 4, with the difference that now one ofthe terms in (17) can just be bounded by k k C h k because c i = 1 for some i . However,the result for A − h, ρ h,n is the same as in that theorem because of (H2c).As we will see in the numerical experiments, this explains that the local errorbehaves very badly, because it grows when h diminishes. However, in spite of that, forfixed but small h , it behaves with order 1 in k because the term in k k C h k dominates.The same happens with the global error, as the following theorem states by using againa summation-by-parts argument. Theorem 7.
Under the same hypotheses of Theorem 6, if u ∈ C ([0 , T ] , X ) and (22)holds, then e h,n = U nh − P h u ( t n ) = O (1 + k k C h k ) . When discretizing firstly in time, some suggestion must be given for the exponential-type operators which turn up in Lawson formulas (9)-(10). If u ( t n ) vanished at theboundary, it could seem natural to substitute those exponential operators by the C -semigroup e τA for suitable scalar values τ . However, that may also lead to orderreduction since the solution of˙ u ( τ ) = A u ( τ ) , u (0) = u , cannot be accurately enough approximated by an expansion of the form u + τ A u + τ A u + · · · unless u ∈ D ( A l ) for a high enough l . Besides, as u ( t n ) does not vanish in gen-eral, suitable initial boundary value problems must be considered with the appropiateboundary values. We will consider three different cases depending on whether we wantto achieve local order 2, 3 or 4. For local order 2 and Dirichlet boundary conditions,those boundaries can be calculated directly in terms of the data f and g of the problem(1). For the same order, but with Neumann/Robin boundary conditions, the approxi-mation at the boundary given by the space discretization of the problem must be used.For higher orders, we will have to resort to numerical differentiation to calculate thoseboundaries.From now on, we will assume that the coefficients of Butcher tableau satisfy thefollowing standard equalities: s X i =1 b i = 1 , i − X j =1 a ij = c i , i = 1 , . . . , s. (24)10 .1 Searching for local order For the stages, starting from the continuous approximation u n at t = t n , we considerrecursively K n,i = v n ( c i k ) + k i − X j =1 a ij w n,j (( c i − c j ) k ) , i = 1 , . . . , s, (25)where ˙ v n ( s ) = Av n ( s ) ,v n (0) = u n ,∂v n ( s ) = ∂u ( t n ) , ˙ w n,j ( s ) = Aw n,j ( s ) ,w n,j (0) = f ( t n + c j k, K n,j ) ,∂w n,j ( s ) = 0 , (26)Then, we suggest as the numerical approximation at the next step u n +1 = ˜ v n ( k ) + k s X i =1 b i ˜ w n,i ((1 − c i ) k ) , (27)where ˙˜ v n ( s ) = A ˜ v n ( s ) , ˜ v n (0) = u n ,∂ ˜ v n ( s ) = ∂ [ u ( t n ) + sAu ( t n )] , ˙˜ w n,j ( s ) = A ˜ w n,j ( s ) , ˜ w n,j (0) = f ( t n + c j k, K n,j ) ,∂ ˜ w n,j ( s ) = ∂f ( t n , u ( t n )) . (28)After discretizing (26) and (28) in space, the following systems arise when startingfrom the discrete numerical approximation U n,h at the previous step, and denoting by K n,h,j to the the discretized stages, (cid:26) ˙ V n,h ( s ) = A h, V n,h ( s ) + C h ∂u ( t n ) ,V n,h (0) = U n,h , (cid:26) ˙ W n,h,j ( s ) = A h, W n,h,j ( s ) ,W n,h,j (0) = f ( t n + c j k, K n,h,j ) . ( ˙˜ V n,h ( s ) = A h, ˜ V n,h ( s ) + C h [ ∂u ( t n ) + sAu ( t n )] , ˜ V n,h (0) = U n,h , ( ˙˜ W n,h,j ( s ) = A h, ˜ W n,h,j ( s ) + C h ∂f ( t n , u ( t n )) , ˜ W n,h,j (0) = f ( t n + c j k, K n,h,j ) . (29)By using the variation-of-constants formula and the definition of ϕ j in (6), we have V n,h ( s ) = e sA h, U n,h + Z s e ( s − τ ) A h, C h ∂u ( t n ) dτ = e sA h, U n,h + sϕ ( sA h, ) C h ∂u ( t n ) , ˜ V n,h ( s ) = e sA h, U n,h + Z s e ( s − τ ) A h, C h [ ∂u ( t n ) + τ Au ( t n )] dτ = e sA h, U n,h + sϕ ( sA h, ) C h ∂u ( t n ) + s ϕ ( sA h, ) C h ∂Au ( t n ) , ˜ W n,h,j ( s ) = e sA h, f ( t n + c j k, K n,h,j ) + sϕ ( sA h, ) C h ∂f ( t n , u ( t n )) . K n,h,i = e c i kA h, U n,h + c i kϕ ( c i kA h, ) C h ∂u ( t n ) + k i − X j =1 a ij e ( c i − c j ) kA h, f ( t n + c j k, K n,h,j ) ,U n +1 ,h = e kA h, U n,h + kϕ ( kA h, ) C h ∂u ( t n ) + k ϕ ( kA h, ) C h ∂Au ( t n )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, f ( t n + c i k, K n,h,i )+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h ∂f ( t n , u ( t n )) (cid:21) . (30) Remark 8.
Notice that the three terms on the boundary ∂u ( t n ) , ∂Au ( t n ) and ∂f ( t n , u ( t n )) ,are necessary to consider this approximation. However, as ∂u ( t n ) = g ( t n ) , ∂Au ( t n ) = ˙ g ( t n ) − ∂f ( t n , u ( t n )) , all reduces to calculate ∂f ( t n , u ( t n )) . In the same way as it was stated in [4, 9], withDirichlet boundary conditions, when the nonlinear term is given by (5), that term canbe calculated exactly as ∂f ( t n , u ( t n )) = φ ( g ( t n )) + ∂h ( t n ) . With Neumann or Robin boundary conditions ∂u ( t ) = αu ( t ) | ∂ Ω + β∂ n u ( t ) | ∂ Ω = g ( t ) , β = 0 , (31) as ∂f ( t n , u ( t n )) = α [ φ ( u ( t n ) | ∂ Ω ) + h ( t n ) | ∂ Ω ] + β [ φ ′ ( u ( t n ) | ∂ Ω ) ∂ n u ( t n ) | ∂ Ω + ∂ n h ( t n ) | ∂ Ω ] ,u ( t n ) | ∂ Ω can be substituted by the numerical approximation which the space discretiza-tion of the problem necessarily gives in this case and ∂ n u ( t n ) | ∂ Ω by the result of solvingfrom (31). In any case, the error which comes from the approximation of the boundaryterms in (30) is given in Table 1. Remark 9.
We also notice that, when ∂u ( t ) = ∂Au ( t ) = 0 , because of (1), it neces-sarily happens that ∂f ( t, u ( t )) = 0 . Besides, in this case, (30) is equivalent to (16). Inthis way, through Theorems 10 and 11, we will be implicitly proving local order forthe classical approach in such a case. In order to define the local error of the time semidiscretization, we consider¯ K n,i = ¯ v n ( c i k ) + k i − X j =1 a ij ¯ w n,j (( c i − c j ) k ) , i = 1 , . . . , s, ¯ u n +1 = ¯˜ v n ( k ) + k s X i =1 b i ¯˜ w n,i ((1 − c i ) k ) , (32)12here ¯ v n and ¯˜ v n satisfy the same equation and boundary conditions as v n and ˜ v n ,but starting from u ( t n ) instead of u n . The same happens with ¯ w n,i , ¯˜ w n,i and w n,i ,˜ w n,i with the difference that the initial condition is now f ( t n + c i k, ¯ K n,i ) instead of f ( t n + c i k, K n,i ). Then, for the local error ρ n = ¯ u n +1 − u ( t n +1 ), we have the followingresult: Theorem 10.
Under hypotheses (A1)-(A5) and (H1)-(H3), assuming also that, forevery t ∈ [0 , T ] , f ( t, u ( t )) ∈ D ( A ) , and Af ( t, u ( t )) ∈ C ([0 , T ] , X ) , u ∈ C ([0 , T ] , D ( A )) ∩ C ([0 , T ] , X ) , (33) it follows that ρ n = O ( k ) . Moreover, if f ∈ C ([0 , T ] × X, X ) , u ∈ C ([0 , T ] , X ) , thereexists a constant C such that the following bound holds k A − f u ( t, u ( t )) A w k ≤ C, for every w ∈ D ( A ) and t ∈ [0 , T ] , (34) and the Runge-Kutta tableau corresponds to a method of classical order ≥ , it followsthat A − ρ n = O ( k ) .Proof. Using Lemma 3.1 in [3],¯ v n ( s ) = u ( t n ) + sϕ ( sA ) Au ( t n ) , (35)¯˜ v ( s ) = u ( t n ) + sAu ( t n ) + s ϕ ( sA ) A u ( t n ) , ¯ w n,i ( s ) = e sA f ( t n + c i k, ¯ K n,i ) , ¯˜ w n,i ( s ) = e sA [ f ( t n + c i k, ¯ K n,i ) − f ( t n , u ( t n ))] + f ( t n , u ( t n )) + sϕ ( sA ) Af ( t n , u ( t n )) . Then,¯ K n,i = u ( t n ) + c i kϕ ( c i kA ) Au ( t n ) + k i − X j =1 a ij e ( c i − c j ) kA f ( t n + c j k, ¯ K n,j ) (36)= u ( t n ) + O ( k ) , i = 1 , . . . , s, (37)¯ u n +1 = u ( t n ) + kAu ( t n ) + k ϕ ( kA ) A u ( t n )+ k s X i =1 b i (cid:20) e (1 − c i ) kA [ f ( t n + c i k, ¯ K n,i ) − f ( t n , u ( t n ))] + f ( t n , u ( t n ))+(1 − c i ) kϕ ((1 − c i ) kA ) Af ( t n , u ( t n )) (cid:21) (38)= u ( t n ) + k [ Au ( t n ) + f ( t n , u ( t n ))] + O ( k ) = u ( t n ) + k ˙ u ( t n ) + O ( k ) , where, for the last line, (33), (A4) together with (37) and the first condition of (24)have been used. From this, the first result on the local error follows.As for the second result, looking at the term in k in ¯ u n +1 and using that u ∈ ([0 , T ] , X ), we can notice that A − ρ n +1 = k ( kA ) − ( ϕ ( kA ) − I ) A u ( t n ) + k A − A u ( t n )+ k s X i =1 b i (1 − c i )((1 − c i ) kA ) − [ e (1 − c i ) kA − I ][ f ( t n + c i k, ¯ K n,i ) − f ( t n , u ( t n ))]+ k s X i =1 b i A − [ f ( t n + c i k, ¯ K n,i ) − f ( t n , u ( t n ))] (39)+ k s X i =1 b i (1 − c i ) ((1 − c i ) kA ) − ( ϕ ((1 − c i ) kA ) − I ) Af ( t n , u ( t n ))+ k s X i =1 b i (1 − c i ) A − Af ( t n , u ( t n )) − k u ( t n ) + O ( k ) . Using (7), (37), (A4) and (33), the first, third and fifth terms are O ( k ). As for thefourth one, considering (36), it can be written as k s X i =1 b i A − f u ( t n , u ( t n ))[ c i ϕ ( c i kA ) Au ( t n ) + i − X j =1 a ij e ( c i − c j ) kA f ( t n + c j k, ¯ K n,j )]+ k s X i =1 b i c i A − f t ( t n , u ( t n )) + O ( k )= k s X i =1 b i A − f u ( t n , u ( t n )) A c i ( kA ) − [ ϕ ( c i kA ) − I ] Au ( t n )+ k ( s X i =1 b i c i ) A − f u ( t n , u ( t n )) Au ( t n )+ k s X i =1 b i i − X j =1 a ij A − f u ( t n , u ( t n )) A ( kA ) − [ e ( c i − c j ) kA − I ] f ( t n , u ( t n ))+ k s X i =1 b i i − X j =1 a ij A − f u ( t n , u ( t n )) f ( t n , u ( t n ))+ k s X i =1 b i c i A − f t ( t n , u ( t n )) + O ( k )= k A − [ f u ( t n , u ( t n )) ˙ u ( t n ) + f t ( t n , u ( t n ))] + O ( k ) , where, for the last equality, we have used (7) again, (34), the second condition in (24)and the fact that P si =1 b i c i = 1 / A − ρ n +1 = k A − [ A u + f u ˙ u + f t + Af − ¨ u ] + O ( k ) = O ( k ) , We define the local error of the full discretization as ρ n +1 ,h = ¯ U n +1 ,h − P h u ( t n +1 ), where¯ U n +1 ,h is defined as U n +1 ,h but starting from P h u ( t n ). More precisely, in a similar wayto the derivation of (30), ¯ U n +1 ,h is defined through some stages ¯ K n,i,h in the followingway: ¯ K n,h,i = ¯ V n,h ( c i k ) + k i − X j =1 a ij ¯ W n,h,j (( c i − c j ) k ) , i = 1 , . . . , s, ¯ U n +1 ,h = ¯˜ V n,h ( k ) + k s X i =1 b i ¯˜ W n,h,i ((1 − c i ) k ) , where ¯ V n,h , ¯ W n,j,h , ¯˜ V n,h and ¯˜ W n,j,h are defined as in (29), but changing U n,h by P h u ( t n )and K n,h,i by ¯ K n,h,i . We then have the following result for the full discretized localerror: Theorem 11.
Under the same hypotheses of the first part of Theorem 10, and assumingalso that, for t ∈ [0 , T ] , A l u ( t ) ∈ Z, l = 0 , , , Af ( t, u ( t )) ∈ Z, (40) it happens that ρ n,h = O ( k + kε h ) where the constant in Landau notation is independentof k and h . Moreover, under the additional hypotheses of the second part of Theorem10, together with the following condition which is related to (34), k A − h, f u ( t, P h u ( t )) A h, k h ≤ C, t ∈ [0 , T ] , h ≤ h , (41) it follows that A − h, ρ n,h = O ( k + kη h + k ǫ h ) .Proof. For the first part of the theorem, as ρ n +1 ,h = ¯ U n +1 ,h − P h u ( t n +1 ) = ( ¯ U n +1 ,h − P h ¯ u n +1 ) + P h (¯ u n +1 − u ( t n +1 ))= ( ¯ U n +1 ,h − P h ¯ u n +1 ) + P h ρ n +1 , (42)because of Theorem 10, it suffices to prove that ¯ U n +1 ,h − P h ¯ u n +1 = O ( kε h ). Consideringthe definition of ¯ u n +1 (32),¯ U n +1 ,h − P h ¯ u n +1 = ¯˜ V n,h ( k ) − P h ¯˜ v n ( k ) + k s X i =1 b i [ ¯˜ W n,i,h − P h ¯˜ w n,i ((1 − c i ) k )] , (43)where ¯˜ v n and ¯˜ w n,i satisfy (28) with u n substituted by u ( t n ) and K n,j substituted by¯ K n,j . Then,˙¯˜ V n,h ( s ) − P h ˙¯˜ v n ( s ) = A h, ¯˜ V n,h ( s ) − A h, R h ¯˜ v n ( s )= A h, ( ¯˜ V n,h ( s ) − P h ¯˜ v n ( s )) + A h, ( P h − R h )¯˜ v n ( s ) , ¯˜ V n,h (0) − P h ¯˜ v n (0) = 0 . v n ( s ) ∈ Z , which implies,using (14) and the variation-of-constants formula, that¯˜ V n,h ( k ) − P h ¯˜ v n ( k ) = Z k e ( k − s ) A h, A h, ( P h − R h )¯˜ v n ( s ) = O ( kε h ) . (44)In a similar way, using also (12),˙¯˜ W n,h,i ( s ) − P h ˙¯˜ w n,i ( s ) = A h, ¯˜ W n,i,h ( s ) − A h, R h ¯˜ w n,i ( s )= A h, ( ¯˜ W n,h,i ( s ) − P h ¯˜ w n,i ( s )) + A h, ( P h − R h ) ¯˜ w n,i ( s ) , ¯˜ W n,h,i (0) − P h ¯˜ w n,i (0) = f ( t n + c i k, ¯ K n,h,i ) − f ( t n + c i k, P h ¯ K n,i ) , where ¯˜ w n,i ( s ) again belongs to Z because of the definition of ¯ K n,i (32), (35), (H2a) andthe conditions of regularity (40). Proceeding as before,¯˜ W n,h,i ((1 − c i ) k ) − P h ¯˜ w n,i ((1 − c i ) k )= e (1 − c i ) kA h, [ f ( t n + c i k, ¯ K n,h,i ) − f ( t n + c i k, P h ¯ K n,i )] + O ( kε h ) , (45)where¯ K n,h,i − P h ¯ K n,i = ¯ V n,h ( c i k ) − P h ¯ v n ( c i k ) + k i − X j =1 a ij [ ¯ W n,h,j (( c i − c j ) k ) − P h ¯ w n,j (( c i − c j ) k )] . (46)Now, with similar arguments as those for deducing (44) and (45),¯ V n,h ( c i k ) − P h ¯ v n ( c i k ) = O ( kε h )¯ W n,h,j (( c i − c j ) k ) − P h ¯ w n,j (( c i − c j ) k ) = e ( c i − c j ) kA h, [ f ( t n + c j k, ¯ K n,h,j ) − f ( t n + c j k, P h ¯ K n,j )] . Therefore, writing this in (46), it is inductively proved that¯ K n,h,i − P h ¯ K n,i = O ( kε h ) , (47)with which the first result of the theorem is proved considering that in (45) and thenin (43) together with (44).As for the second part of the theorem, notice that from (42), A − h, ρ n,h = A − h, ( ¯ U n +1 ,h − P h ¯ u n +1 ) + A − h, P h ρ n +1 . (48)For the first term in (48), we then notice that A − h, ( ¯˜ V n,h ( k ) − P h ¯˜ v n ( k )) = O ( kη h ) , (49)considering (44) multiplied by A − h, and the second formula of (14). On the other hand,considering (45), A − h, [ ¯˜ W n,h,i ((1 − c i ) k ) − P h ¯˜ w n,i ((1 − c i ) k )]= A − h, e (1 − c i ) kA h, [ f ( t n + c i k, ¯ K n,h,i ) − f ( t n + c i k, P h ¯ K n,i )] + O ( kη h ) , (50)16ut A − h, e (1 − c i ) kA h, [ f ( t n + c i k, ¯ K n,h,i ) − f ( t n + c i k, P h ¯ K n,i )]= e (1 − c i ) kA h, A − h, f u ( t n , P h u ( t n ))( ¯ K n,h,i − P h ¯ K n,i ) + O ( k ε h ) , taking (37) and (47) into account. Moreover, by similar arguments as above, A − h, ( ¯ K n,h,i − P h ¯ K n,i ) = O ( kη h ) . Therefore, using (41), A − h, e (1 − c i ) kA h, [ f ( t n + c i k, ¯ K n,h,i ) − f ( t n + c i k, P h ¯ K n,i )] = O ( kη h + k ε h ) , and so, inserting this in (50) and using also (49), A − h, [ ¯ U n +1 ,h − P h ¯ u n +1 ] = O ( kη h + k ε h ) . (51)For the second term in (48, we notice that, as ∂A − ρ n +1 = 0, applying (13), it followsthat A h, R h ( A − ρ n +1 ) = P h ρ n +1 , from what A − h, P h ρ n +1 = R h ( A − ρ n +1 ) = P h ( A − ρ n +1 ) + ( R h − P h )( A − ρ n +1 ) . (52)Then, P h A − ρ n +1 here is O ( k ) because of the second part of Theorem 10. Secondly,( R h − P h ) A − ρ n +1 is O ( kη h ) because A − ρ n +1 belongs to Z due to (38), (36), the condi-tions (40) and (H2a). Moreover, k A − ρ n +1 k Z = O ( k ) because u ( t n +1 ) = u ( t n ) + k ˙ u ( t ∗ n )for t ∗ n ∈ [ t n , t n +1 ], with ˙ u ( t ∗ n ) = Au ( t ∗ n ) + f ( t ∗ n , u ( t ∗ n )) ∈ Z . Therefore, A − h, P h ρ n +1 = O ( k + kη h ) and the result follows considering also (51) in (48). From the first part of Theorem 11, again the classical argument would lead to globalerror e n,h = U n,h − P h u ( t n ) = O ( k + ε h ). However, for parabolic problems, for which(22) is expected to hold, using the second part of the same theorem, a summation-by-parts argument leads to second order in time, as the following theorem states. Thetheorem is valid for both Dirichlet and Neumann/Robin boundary conditions, in spiteof the fact that, for the latter, ∂f ( t, u ( t )) must be approximated through the numericalsolution itself, as explained in Remark 8. Theorem 12.
Under hypotheses of Theorem 11, but assuming also (22) and that f ∈ C ([0 , t ] × X, X ) , Af ( t, u ( t )) ∈ C ([0 , T ] , X ) , u ∈ C ([0 , T ] , D ( A )) ∩ C ([0 , T ] , X ) ,A l ˙ u ( t ) ∈ Z, l = 0 , , , ddt A l f ( t, u ( t )) ∈ Z, l = 0 , , it follows that e n,h = U n,h − P h u ( t n ) = O ( k + kε h + η h ) . roof. Firstly notice that e n +1 ,h = [ U n +1 ,h − ¯ U n +1 ,h ] + [ ¯ U n +1 ,h − P h u ( t n +1 )] = [ U n +1 ,h − ¯ U n +1 ,h ] + ρ n +1 ,h . (53)Then, using (30), when considering Dirichlet boundary conditions, in which case ∂Au ( t n )and ∂f ( t n , u ( t n )) are calculated exactly in terms of data, as ¯ U n +1 ,h is the same as U n +1 ,h but starting from P h u ( t n ) instead of U n,h , U n +1 ,h − ¯ U n +1 ,h = e kA h, [ U n,h − P h u ( t n )]+ k s X i =1 b i e (1 − c i ) kA h, [ f ( t n + c i k, K n,h,i ) − f ( t n + c i k, ¯ K n,h,i )] , (54)where, recursively, for i = 1 , . . . , s , K n +1 ,h,i − ¯ K n +1 ,h,i = e c i kA h, [ U n,h − P h u ( t n )]+ k i − X j =1 a ij e ( c i − c j ) kA h, [ f ( t n + c j k, K n,h,j ) − f ( t n + c j k, ¯ K n,h,j )] . (55)In such a way, it is inductively proved that K n +1 ,h,i − ¯ K n +1 ,h,i = O ( e n,h ) and finally,using (H3), (54) and (53), e n +1 ,h = e kA h, e n,h + O ( ke n,h ) + ρ n +1 ,h , (56)from what the result follows from Theorem 11 by a summation-by-parts argument anda discrete Gronwall lemma in the same way than the proof of Theorem 22 in [4] forStrang method.On the other hand, when considering Robin/Neumann boundary conditions, as,according to Remark 8, ∂f ( t n , u ( t n )) is just calculated approximately with an errorwhich is O ( e n,h ), using (30) again, U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k ϕ ( kA h, ) C h O ( e n,h )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, [ f ( t n + c i k, K n,h,i ) − f ( t n + c i k, ¯ K n,h,i )]+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h O ( e n,h ) (cid:21) , where K n +1 ,h,i − ¯ K n +1 ,h,i is the same as in (55) because ∂u ( t n ) is given exactly in termsof data with this type of boundary conditions. Then, using (7) and (H3), U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k [ ϕ ( kA h, ) − I ] A − h, C h O ( e n,h )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, O ( e n,h ) + [ e (1 − c i ) kA h, − I ] A − h, C h O ( e n,h ) (cid:21) . Using now (H2c), it follows that U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + O ( ke n,h ), from what(56) applies again and the result follows in the same way as above.18 .2 Searching for local order 3 For the stages, we again consider (25), but where now v n and w n.,j are those in (28)instead of those in (26). On the other hand, u n +1 is calculated through (27) where now˜ v n , ˜ w n,j ( j = 1 , . . . , s ) satisfy ˙˜ v n ( s ) = A ˜ v n ( s ) , ˜ v n (0) = u n ,∂ ˜ v n ( s ) = ∂ [ u ( t n ) + sAu ( t n ) + s A u ( t n )] , ˙˜ w n,j ( s ) = A ˜ w n,j ( s ) , ˜ w n,j (0) = f ( t n + c j k, K n,j ) ,∂ ˜ w n,j ( s ) = ∂ [ f ( t n + c j k, u ( t n ) + c j k ˙ u ( t n )) + sAf ( t n , u ( t n ))] . (57)After discretizing in space and using the variation-of-constants formula as in Subsection4.1, the full discretized numerical solution after one step is given by K n,h,i = e c i kA h, U n,h + c i kϕ ( c i kA h, ) C h ∂u ( t n ) + c i k ϕ ( c i kA h, ) C h ∂Au ( t n )+ k i − X j =1 a ij (cid:20) e ( c i − c j ) kA h, f ( t n + c j k, K n,h,j )+( c i − c j ) kϕ (( c i − c j ) kA h, ) C h ∂f ( t n , u ( t n )) (cid:21) , i = 1 , . . . , s,U n +1 ,h = e kA h, U n,h + X l =1 k l ϕ l ( kA h, ) C h ∂A l − u ( t n )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, f ( t n + c i k, K n,h,i )+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n ))+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h ∂Af ( t n , u ( t n )) (cid:21) . (58) Remark 13.
Notice that, apart from the terms on the boundary which were alreadynecessary to achieve local order and which can be calculated according to Remark 8,now we also need ∂A u ( t n ) , ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) and ∂Af ( t n , u ( t n )) .With Dirichlet boundary conditions and f like in (5), it follows that ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) = φ ( g ( t n ) + c i k ˙ g ( t n )) + ∂h ( t n + c i k ) ,∂A u ( t n ) = ¨ g ( t n ) − ∂ ˙ h ( t n ) − φ ′ ( g ( t n )) ˙ g ( t n ) − ∂Af ( t n , u ( t n )) , and the only term which cannot be calculated exactly in terms of data is ∂Af ( t n , u ( t n )) .However, that can be approximated recurring to numerical differentiation. For example,in one dimension and assuming that A is the second spatial derivative, Af ( t n , u ( t n )) = φ ′′ ( u ( t n )) u x ( t n ) + φ ′ ( u ( t n )) u xx ( t n ) + h xx ( t n ) , rom what ∂Af ( t n , u ( t n )) ≈ φ ′′ ( g ( t n ))ˆ u x ( t n ) | ∂ Ω + φ ′ ( g ( t n ))( ˙ g ( t n ) − φ ( g ( t n )) − ∂h ( t n )) + ∂h xx ( t n ) , where ˆ u x ( t n ) | ∂ Ω is the result of applying numerical differentiation to approximate u x ( t n ) on the boundary. For that, both the exact values at the boundary and the approximatedvalues at the interior of the domain given by the numerical approximation must beused. As a result, ˆ u x ( t n ) − u x ( t n ) = O ( ν h + e n,h h ) , where ν h decreases with h and comesfrom the error of the numerical differentiation if the exact values of the solutions wereused. The second term e n,h h comes from the fact that the values at the interior are justthe approximations which are given by the numerical solution and to the necessity ofdividing by h when approximating a first derivative in space. For a general operator A ,we will assume that the error when approximating both ∂A u ( t n ) and ∂Af ( t n , u ( t n )) is as specified in Table 1 for some real value γ , where ν h comes from the numericalapproximation of the corresponding derivatives in space if the exact values had beentaken.As for Robin/Neumann boundary conditions (31), we notice that ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) = α [ φ ( u ( t n ) | ∂ Ω + c i k ˙ u ( t n ) | ∂ Ω ) + h ( t n + c i k ) | ∂ Ω ]+ β [ φ ′ ( u ( t n ) | ∂ Ω + c i k ˙ u ( t n ) | ∂ Ω )( ∂ n u ( t n ) | ∂ Ω + c i k∂ n ˙ u ( t n ) | ∂ Ω ) + ∂ n h ( t n + c i k ) | ∂ Ω ] . Here u ( t n ) | ∂ Ω and ∂ n u ( t n ) | ∂ Ω are approximated through the numerical solution, as inRemark 8, ˙ u ( t n ) | ∂ Ω is then approximated through numerical differentiation in time fromthe approximated values at the boundary and ∂ n ˙ u ( t n ) | ∂ Ω is then solved from the differ-entiation in time of (31). In such a way, if the error coming from the numericaldifferentiation in time from the exact values is O ( µ k, ) , with µ k, decreasing when k decreases, it happens that the error when approximating ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) is O ( kµ k, + e n,h ) , with the same argument as before for space numerical differentia-tion, and taking now into account the factor k which is multiplying the correspondingderivative. As for ∂A u ( t n ) , using (1), ∂A u = ¨ g − ∂ [ ˙ h + φ ′ ( u ) ˙ u ] − ∂ ( Af ) . Moreover, ∂ [ ˙ h + φ ′ ( u ) ˙ u ] = α [ ˙ h | ∂ Ω + φ ′ ( u | ∂ Ω ) ˙ u | ∂ Ω ] + β [ ∂ n ˙ h | ∂ Ω + φ ′′ ( u | ∂ Ω ) ∂ n u | ∂ Ω ˙ u | ∂ Ω + φ ′ ( u | ∂ Ω ) ∂ n ˙ u | ∂ Ω ] , and it is then necessary to approximate u | ∂ Ω , ∂ n u | ∂ Ω , ˙ u | ∂ Ω and ∂ n ˙ u | ∂ Ω , as before. Onthe other hand, when A is the second derivative in one dimension, ∂ ( Af ) = α (cid:20) φ ′′ ( u | ∂ Ω ) u x | ∂ Ω + φ ′ ( u | ∂ Ω )[ ˙ u | ∂ Ω − φ ( u | ∂ Ω ) − h | ∂ Ω ] + h xx | ∂ Ω (cid:21) + β (cid:20) φ ′′′ ( u | ∂ Ω ) u x | ∂ Ω + 3 φ ′′ ( u | ∂ Ω ) u x | ∂ Ω [ ˙ u | ∂ Ω − φ ( u | ∂ Ω ) − h | ∂ Ω ]+ φ ′ ( u | ∂ Ω )[ ˙ u x | ∂ Ω − φ ′ ( u | ∂ Ω ) u x | ∂ Ω − h x | ∂ Ω ] + h xxx | ∂ Ω (cid:21) . Therefore, in this particular case, approximating u | ∂ Ω , ˙ u | ∂ Ω , u x | ∂ Ω and ˙ u x | ∂ Ω as above,the error which comes from calculating ∂A u ( t n ) and ∂Af ( t n , u ( t n )) is O ( µ k, + e n,h k ) .However, for more general operators, space numerical differentiation may be also needed,and therefore we will assume that, in general, the error coming from the calculation ofthose boundaries is as specified in Table 1. emark 14. We notice that, if ∂u ( t ) = ∂Au ( t ) = ∂A u ( t ) = 0 , from (1) it follows that ∂f ( t, u ( t )) = ∂Af ( t, u ( t )) = 0 . In particular, this implies that ∂f ( t n + c i k, u ( t n + c i k )) =0 , which differs from ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) in O ( k ) . Then, in this case, eachstep in (58) differs from the classical approach (16) in O ( k ) since the difference is k s X i =1 b i (1 − c i ) kϕ ((1 − c i ) kA h, ) C h O ( k ) = k s X i =1 b i [ e (1 − c i ) kA h, − I ] A − h, C h O ( k ) , and, according to (H2c), A − h, C h is uniformly bounded. This justifies, through Theo-rems 15 and 16, that the local error with the classical approach, under these particularboundary conditions, behaves with order under the assumptions of those theorems. In a similar way to Theorem 10, we have the following result:
Theorem 15.
Under hypotheses (A1)-(A5) and (H1)-(H3), assuming also that, forevery t ∈ [0 , T ] , f ( t, u ( t )) ∈ D ( A ) , that, for small enough τ ( τ ≤ τ ), f ( t + τ, u ( t ) + τ ˙ u ( t )) , f t ( t + τ, u ( t ) + τ ˙ u ( t )) , f u ( t + τ, u ( t ) + τ ˙ u ( t )) ˙ u ( t ) ∈ D ( A ) and u ∈ C ([0 , T ] , D ( A )) ∩ C ([0 , T ] , X ) ,f ∈ C ([0 , T ] × X, X ) ,A l f ( · , u ( · )) ∈ C ([0 , T ] , X ) , l = 1 , ,Af t ( t + τ, u ( t ) + τ ˙ u ( t )) , A [ f u ( t + τ, u ( t ) + τ ˙ u ( t )) ˙ u ( t )] ∈ C ([0 , T ] × [0 , τ ] , X ) , (59) and that the Runge-Kutta tableau corresponds to a method of classical order ≥ , itfollows that ρ n = O ( k ) . Moreover, if f ∈ C ([0 , T ] × X, X ) , u ∈ C ([0 , T ] , X ) , (34)holds and the Runge-Kutta tableau corresponds to a method of classical order ≥ , itfollows that A − ρ n = O ( k ) .Proof. Using Lemma 3.1 in [3] again,¯ v n ( s ) = u ( t n ) + sAu ( t n ) + s ϕ ( sA ) A u ( t n ) , ¯˜ v n ( s ) = u ( t n ) + sAu ( t n ) + s A u ( t n ) + s ϕ ( sA ) A u ( t n ) , ¯ w n,i ( s ) = e sA [ f ( t n + c i k, ¯ K n,i ) − f ( t n , u ( t n ))] + f ( t n , u ( t n )) + sϕ ( sA ) Af ( t n , u ( t n )) , and with a similar argument to that used in that lemma,¯˜ w n,i ( s ) = e sA [ f ( t n + c i k, ¯ K n,i ) − f ( t n , u ( t n ) + c i k ˙ u ( t n ))] + f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n ))+ sAf ( t n , u ( t n )) + sϕ ( sA )[ Af ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) − Af ( t n , u ( t n ))]+ s ϕ ( sA ) A f ( t n , u ( t n )) . K n,i = u ( t n ) + c i kAu ( t n ) + c i k ϕ ( kA ) A u ( t n )+ k i − X j =1 a ij (cid:20) e ( c i − c j ) kA [ f ( t n + c j k, ¯ K n,j ) − f ( t n , u ( t n ))] + f ( t n , u ( t n ))+( c i − c j ) kϕ (( c i − c j ) kA ) Af ( t n , u ( t n )) (cid:21) = u ( t n ) + c i k ˙ u ( t n ) + O ( k ) , i = 1 , . . . , s, ¯ u n +1 = u ( t n ) + kAu ( t n ) + k A u ( t n ) + k ϕ ( kA ) A u ( t n )+ k s X i =1 b i (cid:20) f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) + (1 − c i ) kAf ( t n , u ( t n ))+ e (1 − c i ) kA [ f ( t n + c i k, ¯ K n,i ) − f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n ))]+(1 − c i ) kϕ ((1 − c i ) kA )[ Af ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) − Af ( t n , u ( t n ))]+(1 − c i ) k ϕ ((1 − c i ) kA ) A f ( t n , u ( t n )) (cid:21) = u ( t n ) + k [ Au ( t n ) + f ( t n , u ( t n ))]+ k A u ( t n ) + f t ( t n , u ( t n )) + f u ( t n , u ( t n )) ˙ u ( t n ) + Af ( t n , u ( t n ))] + O ( k ) , = u ( t n ) + k ˙ u ( t n ) + k u ( t n ) + O ( k ) , where (59) has been used, and so the first part of the theorem follows.For the second part of the theorem, in a similar way to the proof of Theorem 10,but looking now at the term in k and using the new stronger hypotheses, it can beproved that A − ρ n +1 = k A − (cid:20) A u + f tt + 2 f tu ˙ u + f uu ˙ u + f u A u + f u [ f t + f u ˙ u + Af ]+ Af t + A ( f u ˙ u ) + A f − ... u (cid:21) + O ( k ) , (60)and the term in bracket vanishes by differentiating (1) twice. Following a similar proof to that of Theorem 11, the following result turns up:
Theorem 16.
Under the same hypotheses of the first part of Theorem 15, and assumingalso that, for t ∈ [0 , T ] and τ ∈ [0 , τ ] , A l u ( t ) ∈ Z, l = 0 , , , , A l f ( t, u ( t )) ∈ Z, l = 1 , , Af ( t + τ, u ( t ) + τ ˙ u ( t )) ∈ Z, it happens that ρ n,h = O ( k + kε h ) where the constant in Landau notation is independentof k and h . Moreover, under the additional hypotheses of the second part of Theorem15, together with condition (41), it follows that A − h, ρ n,h = O ( k + kη h + k ǫ h ) . .2.3 Global error of the full discretization This subsection is different from 4.1.3 in the fact that, for both Dirichlet and Robin/Neumannboundary conditions, numerical differentiation must be used to approximate the cor-responding boundary values in (58). Because of that, in order to assure convergencewith this technique, we will have to ask that k is sufficiently small with respect to h γ ,where γ is the parameter which turns up when applying numerical differentiation, asstated in Remark 13. Again, the classical argument would lead to a worse bound forthe global error than in parabolic problems, when a summation-by-parts argumentscan be used. Theorem 17.
Under the hypotheses of the first part of Theorem 16, if there exists aconstant C such that kh γ ≤ C, (61) when considering Dirichlet boundary conditions, e n,h = O ( k + ε h + kν h ) and, withRobin/Neumann boundary conditions, e n,h = O ( k + ε h + kν h + kµ k, ) , where ν h and µ k, are the errors coming respectively from numerical differentiation in space and timeaccording to Remark 13. On the other hand, under the hypotheses of the second partof Theorem 16, but assuming also (22) and that f ∈ C ([0 , T ] × X, X ) , u ∈ C ([0 , T ] , D ( A )) ∩ C ([0 , T ] , X ) ,A l f ( t, u ( t )) ∈ C ([0 , T ] , X ) , l = 1 , , Af ( t + τ, u ( t ) + τ ˙ u ( t )) ∈ C ([0 , T ] × [0 , τ ] , X ) ,A l ˙ u ( t ) ∈ Z, l = 0 , , , , ddt A l f ( t, u ( t )) ∈ Z, l = 0 , , , t ∈ [0 , T ] ,ddt A l f ( t + τ, u ( t ) + τ ˙ u ( t )) ∈ Z, l = 0 , , τ ∈ (0 , τ ] , (62) it follows that, when considering Dirichlet boundary conditions, e n,h = O ( k + kε h + η h + kν h ) and, with Robin/Neumann boundary conditions, e n,h = O ( k + kε h + η h + kµ k, + kν h ) .Proof. For the proof, as in Theorem 12, we must consider the decomposition (53) where¯ U n +1 ,h is calculated as U n +1 ,h but starting from P h u ( t n ) and calculating the boundariesin (58) in an exact way. In contrast, according to Table 1, when considering U n +1 ,h ,the boundaries in (58) can just be calculated approximately.More precisely, with Dirichlet boundary conditions, the terms on the boundaryfor the stages in (58) can be calculated exactly. However, when calculating U n +1 ,h , ∂A u ( t n ) and ∂Af ( t n , u ( t n )) can just be calculated except for O ( ν h + e n,h h γ ). Because ofthis, U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k ϕ ( kA h, ) C h O ( ν h + e n,h h γ )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, [ f ( t n + c i k, K n,h,i ) − f ( t n + c i k, ¯ K n,h,i )]+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( ν h + e n,h h γ ) (cid:21) , K n,h,i − ¯ K n,h,i = O ( e n,h ) as in the proof of Theorem 12. Therefore, using (7) and(H2c), U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k ( ϕ ( kA h, ) − I ) O ( ν h + e n,h h γ ) + O ( ke n,h )+ k s X i =1 b i (1 − c i )( ϕ ( kA h, ) − I ) O ( ν h + e n,h h γ ) , from what, using condition (61), e n +1 ,h = e kA h, e n,h + O ( ke n,h ) + O ( k ν h ) + ρ n +1 ,h . The classical argument of convergence and the first part of Theorem 16 leads thento the first result of this theorem for Dirichlet boundary conditions. For the secondpart, the second part of Theorem 16 must be used, apart from (22) and the additionalregularity (62).On the other hand, with Robin/Neumann boundary conditions, there is some errorwhen approximating the boundaries for both the stages and the numerical solution.More precisely, using Table 1 and (58), K n,h,i − ¯ K n,h,i = e c i kA h, e n,h + c i k ϕ ( c i kA h, ) C h O ( e n,h )+ k i − X j =1 a ij [ O ( e n,h ) + ( c i − c j ) kϕ (( c i − c j ) kA h, ) C h O ( e n,h )]= e c i kA h, e n,h + c i k [ ϕ ( c i kA h, ) − I ] O ( e n,h )+ k i − X j =1 a ij (cid:20) O ( e n,h ) + [ e ( c i − c j ) kA h, − I ] O ( e n,h ) (cid:21) = e c i kA h, e n,h + O ( ke n,h ) = O ( e n,h ) , and then U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k ϕ ( kA h, ) C h O ( µ k, + e n,h k + ν h + e n,h h γ )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, [ f ( t n + c i k, K n,h,i ) − f ( t n + c i k, ¯ K n,h,i )]+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h O ( kµ k, + e n,h )+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( µ k, + e n,h k + ν h + e n,h h γ ) (cid:21) = e kA h, e n,h + k [ ϕ ( kA h, ) − I ] O ( µ k, + e n,h k + ν h + e n,h h γ )+ k s X i =1 b i (cid:20) O ( e n,h ) + [ e (1 − c i ) kA h, − I ] O ( kµ k, + e n,h )+(1 − c i ) k [ ϕ ( kA h, ) − I ] O ( µ k, + e n,h k + ν h + e n,h h γ ) (cid:21) . e n +1 ,h = e kA h, e n,h + O ( ke n,h + k µ k, + k ν h ) + ρ n +1 ,h , so that, using the first part of Theorem 16 and the classical argument of convergence, e n,h = O ( k + ε h + kν h + kµ k, ). Again, under the second set of hypotheses in Theorem16 and using (22) and the regularity (62), the finer result e n,h = O ( k + kε h + η h + kµ k, + kν h ) can be achieved. The idea is again to calculate the stages as in (25) but with v n , w n,j as in (57), and u n +1 through (27) but with ˜ v n , ˜ w n,j satisfying ˙˜ v n ( s ) = A ˜ v n ( s ) , ˜ v n (0) = u n ,∂ ˜ v n ( s ) = ∂ [ u ( t n ) + sAu ( t n ) + s A u ( t n ) + s A u ( t n )] , ˙˜ w n,j ( s ) = A ˜ w n,j ( s ) , ˜ w n,j (0) = f ( t n + c j k, K n,j ) ,∂ ˜ w n,j ( s ) = ∂ (cid:20) f (cid:18) t n + c j k, u ( t n ) + c j kAu ( t n ) + c j k A u ( t n )+ k P j − r =1 a j,r [ f ( t n + c r k, u ( t n ) + c r k ˙ u ( t n )) + ( c j − c r ) kAf ( t n , u ( t n ))] (cid:19) + sAf ( t n + c j k, u ( t n ) + c j k ˙ u ( t n )) + s A f ( t n , u ( t n )) (cid:21) . ∂u ( t n ) - - ∂Au ( t n ) /∂f ( t n , u ( t n )) - O ( e n,h ) ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) - O ( kµ k, + e n,h ) ∂A u ( t n ) /∂Af ( t n , u ( t n )) O ( ν h + e n,h h γ ) O ( µ k, + e n,h k + ν h + e n,h h γ ) ∂Af ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) O ( ν h + e n,h h γ + kµ k, h γ ) O ( kµ k, + kµ k, + e n,h k + ν h + e n,h h γ ) ∂A u ( t n ) /∂A f ( t n , u ( t n )) O ( ν h + e n,h kh γ + µ k, h γ ) O ( µ k, + µ k, + e n,h k + ν h + e n,h kh γ )Table 1: Errors which are committed at each step when approximating the correspond-ing boundary terms with the suggested technique to avoid order reduction, as justifiedin Remarks 8, 13 and 18.Again, after discretizing these problems in space and using the variation-of-constantsformula, the following full discretation formulas arise: K n,h,i = e c i kA h, U n,h + X l =1 c li k l ϕ l ( c i kA h, ) C h ∂A l − u ( t n )+ k i − X j =1 a ij (cid:20) e ( c i − c j ) kA h, f ( t n + c j k, K n,h,j )+( c i − c j ) kϕ (( c i − c j ) kA h, ) C h ∂f ( t n + c j k, u ( t n ) + c j k ˙ u ( t n ))+( c i − c j ) k ϕ (( c i − c j ) kA h, ) C h ∂Af ( t n , u ( t n )) (cid:21) ,U n +1 ,h = e kA h, U n,h + X l =1 k l ϕ l ( kA h, ) C h ∂A l − u ( t n )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, f ( t n + c i k, K n,h,i )+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h ∂f (cid:18) t n + c i k, u ( t n ) + c i kAu ( t n ) + c i k A u ( t n )+ k i − X j =1 a ij [ f ( t n + c j k, u ( t n ) + c j k ˙ u ( t n )) + ( c i − c j ) kAf ( t n , u ( t n ))] (cid:19) +(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h ∂Af ( t n + c i k, u ( t n ) + c i k ˙ u ( t n ))+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h ∂A f ( t n , u ( t n )) (cid:21) . (63) Remark 18.
Apart from the terms on the boundaries which were already necessary toachieve local order , now ∂A u ( t n ) , ∂A f ( t n , u ( t n )) and ∂Af ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) (64) must also be calculated. Using (1) and simplifying notation, ∂A u = ∂ [ ... u − ( f tt + 2 f tu ˙ u + f uu ˙ u + f u ¨ u )] − ∂ [ A ( f t + f u ˙ u )] − ∂A f, (65)26 here here everything is assumed to be evaluated either on t n or ( t n , u ( t n )) . Now, inorder to calculate ∂ [ A ( f t + f u ˙ u )] and ∂A f , we can see that, with A the second derivativein one dimension and f like in (5), A ( f t + f u ˙ u ) = φ ′′′ ( u ) u x ˙ u + φ ′′ ( u ) u xx ˙ u + 2 φ ′′ ( u ) u x ˙ u x + φ ′ ( u ) ˙ u xx + h txx , and it happens that u xx and ˙ u xx can be calculated through u xx = ˙ u − φ ( u ) − h, ˙ u xx = ¨ u − φ ′ ( u ) ˙ u − h t . (66) As for A f , A f = φ (4) ( u ) u x + 6 φ (3) ( u ) u x u xx + 3 φ ′′ ( u ) u xx + 4 φ ′′ ( u ) u x u xxx + φ ′ ( u ) u xxxx + h xxxx , where u xx can be calculated as in (66) and u xxx = ˙ u x − φ ′ ( u ) u x − h x , u xxxx = ¨ u − φ ′ ( u ) ˙ u − h t − φ ′′ ( u ) u x − φ ′ ( u )( ˙ u − φ ( u ) − h ) − h xx . Finally, Af ( t n + c i k, u ( t n )+ c i k ˙ u ( t n )) = φ ′′ ( u + c i k ˙ u )( u x + c i k ˙ u x ) + φ ′ ( u + c i k ˙ u )( u xx + c i k ˙ u xx )+ h xx ( t n + c i k ) , where, in the right-hand-side u and its derivatives are all evaluated at t = t n and, for u xx , (66) can again be used.Therefore, with Dirichlet boundary conditions, the boundary of every term is exactlycalculable in terms of data except for u x and ˙ u x . Then, u x can be approximated as inRemark 13 and so ˙ u x considering also space numerical differentiation over the exactvalues of ˙ u on the boundary and the approximated values of ˙ u in the interior of thedomain, which can again be approximated by numerical differentiation in time fromthe values of the numerical solution. It can thus be deduced that, in this case, theapproximation of ˙ u x at the boundary differs from the exact in O ( ν h + e n,h hk + µ k, h ) , where µ k, is the error which comes from the approximation of the first derivative if the valuesfrom which it is calculated were all exact. Because of this, the error of approximationof the first two terms in (64) behaves as O ( ν h + e n,h hk + µ k, h ) and, for the last one, as O ( ν h + e n,h hk + k µ k, h ) due to the factor k multiplying ˙ u x . For a general operator A , wewill assume all terms in (64) are calculated except for an error as that in Table 1 forsome real value γ .With Robin/Neumann boundary conditions, in one dimension and with A the secondderivative in space, u | ∂ Ω , u x | ∂ Ω , ˙ u | ∂ Ω , ˙ u x | ∂ Ω , ¨ u | ∂ Ω , ¨ u x | ∂ Ω will also be needed. (Notice that ... u just appears linearly in (65) and therefore is given directly in terms of data through ∂ ... u = ... g .) Then, u | ∂ Ω and u x | ∂ Ω are calculated except for O ( e n,h ) as in Remark 8; ˙ u | ∂ Ω and ˙ u x | ∂ Ω except for O ( µ k, + e n,h k ) as in Remark 13 and ¨ u | ∂ Ω and ¨ u x | ∂ Ω in a similar waythrough numerical differentiation except for O ( µ k, + e n,h k ) , where µ k, comes from theerror in the numerical approximation of the second derivative. For a general operator A , we thus assume that the error in the calculation of the first two boundaries of (64)is as written in the right-bottom part of Table 1 where the last two terms come fromthe possible error in the numerical approximation of some spatial derivatives of u and ˙ u . For the boundary of the last term in (64), the error in the calculation, due to thefactor k multiplying ˙ u , is also as specified in Table 1. emark 19. We notice that, if ∂u ( t ) = ∂Au ( t ) = ∂A u ( t ) = ∂A u ( t ) = 0 , from(1) it follows that ∂f ( t, u ( t )) = ∂Af ( t, u ( t )) = ∂A f ( t, u ( t )) = 0 . As in Remark 14, ∂f ( t n + c i k, u ( t n + c i k )) = 0 differs from ∂f ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) in O ( k ) . Then,the stages in (63) differ from those in the classical approach (16) in k i − X j =1 a ij ( c i − c j ) kϕ (( c i − c j ) kA h, ) C h O ( k ) = k i − X j =1 a ij [ e ( c i − c j ) kA h, − I ] A − h, C h O ( k ) = O ( k ) . On the other hand, ∂f ( t n + c i k, u ( t n + c i k )) = 0 also differs from the factor withthe big parenthesis in (63) in O ( k ) and ∂Af ( t n + c i k, u ( t n + c i k )) = 0 differs from ∂Af ( t n + c i k, u ( t n ) + c i k ˙ u ( t n )) in O ( k ) . Because of all this, the difference in thenumerical solution between (63) and the classical approach is k s X i =1 b i (cid:20) e (1 − c i ) kA h, O ( k ) + (1 − c i ) kϕ ((1 − c i ) kA h, ) C h O ( k )+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( k ) (cid:21) = k s X i =1 b i (cid:20) e (1 − c i ) kA h, O ( k ) + [ e (1 − c i ) kA h, − I ] A − h, C h O ( k )+(1 − c i ) k [ ϕ ((1 − c i ) kA h, ) − I ] A − h, C h O ( k ) (cid:21) = O ( k ) . This justifies, through Theorems 20 and 21, that the local error with the classical ap-proach, under these particular boundary conditions, behaves with order under theassumptions of those theorems. With a similar proof to that of Theorems 10 and 15, the following result follows:
Theorem 20.
Under hypotheses (A1)-(A5) and (H1)-(H3), assuming also that, forevery t ∈ [0 , T ] , f ( t, u ( t )) ∈ D ( A ) ; for small enough τ ( τ ≤ τ ), f ( t + τ, u ( t ) + τ ˙ u ( t )) , f t ( t + τ, u ( t ) + τ ˙ u ( t )) , f u ( t + τ, u ( t ) + τ ˙ u ( t )) ˙ u ( t ) ∈ D ( A ) ; for small enough k and σ ( k ≤ k , σ ≤ σ ) and i = 1 , . . . , s , f (cid:18) t + c i k, u ( t ) + c i kAu ( t ) + c i σ A u ( t )+ σ i − X j =1 a ij [ f ( t + c j k, u ( t ) + c j k ˙ u ( t )) + ( c i − c j ) kAf ( t, u ( t ))] (cid:19) ∈ D ( A ) ,f u (cid:18) t + c i k, u ( t ) + c i kAu ( t ) + c i σ A u ( t )+ σ i − X j =1 a ij [ f ( t + c j k, u ( t ) + c j k ˙ u ( t )) + ( c i − c j ) kAf ( t, u ( t ))] (cid:19) ∈ D ( A );28 ∈ C ([0 , T ] , D ( A )) ∩ C ([0 , T ] , X ) ,f ∈ C ([0 , T ] × X, X ) ,A l f ( t, u ( t )) ∈ C ([0 , T ] , X ) , l = 1 , , ,A l f t ( t + τ, u ( t ) + τ ˙ u ( t )) , A l [ f u ( t + τ, u ( t ) + τ ˙ u ( t )) ˙ u ( t )] ∈ C ([0 , T ] × [0 , τ ] , X ) , l = 1 , ,A (cid:20) f u (cid:18) t + c i k, u ( t ) + c i k ˙ u ( t ) + c i σ A u ( t )+ σ i − X j =1 a ij [ f ( t + c j k, u ( t ) + c j k ˙ u ( t )) + ( c i − c j ) kAf ( t, u ( t ))] (cid:19) · [ c i σA u ( t ) + i − X j =1 a ij [ f ( t + c j k, u ( t ) + c j k ˙ u ( t )) + ( c i − c j ) kAf ( t, u ( t ))] (cid:21) ∈ C ([0 , T ] × [0 , k ] × [0 , σ ] , X ) . (67) and if the Runge-Kutta tableau corresponds to a method of classical order ≥ , it followsthat ρ n = O ( k ) . Moreover, if f ∈ C ([0 , T ] × X, X ) , u ∈ C ([0 , T ] , X ) , (34) holds andthe Runge-Kutta tableau corresponds to a method of classical order ≥ , it follows that A − ρ n = O ( k ) . In a similar way to the proof of Theorem 16, it follows that
Theorem 21.
Under the same hypotheses of the first part of Theorem 20, and assumingalso that, for t ∈ [0 , T ] , i = 1 , . . . , s and k ∈ [0 , k ] , A l u ( t ) ∈ Z, l = 0 , , . . . , , A l f ( t, u ( t )) ∈ Z, l = 1 , , ,A l f ( t + c i k, u ( t ) + c i k ˙ u ( t )) ∈ Z, l = 1 , ,Af (cid:18) t + c i k, u ( t ) + c i kAu ( t ) + c i k A u ( t )+ k i − X j =1 a ij [ f ( t + c j k, u ( t ) + c j k ˙ u ( t )) + ( c i − c j ) kAf ( t, u ( t ))] (cid:19) ∈ Z, (68) it happens that ρ n,h = O ( k + kε h ) where the constant in Landau notation is independentof k and h . Moreover, under the additional hypotheses of the second part of Theorem20, together with condition (41), A − h, ρ n,h = O ( k + kη h ) . In a similar way to Subsection 4.2.3,
Theorem 22.
Under hypotheses of the first part of Theorem 21 and assuming that (61)holds, when considering Dirichlet boundary conditions, e n,h = O ( k + ε h + kν h + kµ k, )29 nd, with Robin/Neumann boundary conditions, e n,h = O ( k + ε h + kν h + kµ k, + k µ k, ) ,where ν h and µ k, , µ k, are the errors coming from numerical differentiation in spaceand time according to Remark 18. On the other hand, under the hypotheses of thesecond part of Theorem 21, but assuming also (22) and that f ∈ C ([0 , T ] × X, X ) , u ∈ C ([0 , T ] , D ( A )) ∩ C ([0 , T ] , X ) ,A l f ( · , u ( · )) ∈ C ([0 , T ] , X ) , l = 1 , , ,A l f ( t + c i k, u ( t ) + c i k ˙ u ( t )) ∈ C ([0 , T ] , X ) , l = 1 , , i = 1 , . . . , s, k ∈ (0 , k ] ,Af (cid:18) · + c i k, u ( · ) + c i kAu ( · ) + c i k A u ( · )+ k i − X j =1 a ij [ f ( · + c j k, u ( · ) + c j k ˙ u ( · )) + ( c i − c j ) kAf ( · , u ( · ))] (cid:19) ∈ C ([0 , T ] , X ) ,A l ˙ u ( t ) ∈ Z, l = 0 , , . . . , , ddt A l f ( t, u ( t )) ∈ Z, l = 0 , , , , t ∈ [0 , T ] ,ddt A l f ( t + c i k, u ( t ) + c i k ˙ u ( t )) ∈ Z, l = 0 , , , k ∈ [0 , k ] ,ddt A l f (cid:18) t + c i k, u ( t ) + c i kAu ( t ) + c i k A u ( t )+ k i − X j =1 a ij [ f ( t + c j k, u ( t ) + c j k ˙ u ( t )) + ( c i − c j ) kAf ( t, u ( t ))] (cid:19) ∈ Z, l = 0 , , (69) it follows that, with Dirichlet boundary conditions, e n,h = O ( k + kε h + η h + kν h + kµ k, ) and, with R/N boundary conditions, e n,h = O ( k + kε h + η h + kµ k, + k µ k, + kν h ) .Proof. As in the proof of Theorem 17, we must consider the decomposition (53) andthen study the difference U h,n +1 − ¯ U n +1 taking into account that the boundaries for U h,n +1 in (63) are just calculated approximately with an error which is given throughTable 1.More precisely, with Dirichlet boundary conditions, U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k ϕ ( kA h, ) C h O ( ν h + e n,h h γ )+ k ϕ ( kA h, ) C h O ( ν h + e n,h h γ k + µ k, h γ )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, [ f ( t n + c i k, K n,h,i ) − f ( t n + c i k, ¯ K n,h,i )]+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h O ( k ν h + k e n,h h γ )+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( ν h + e n,h h γ + kµ k, h γ )+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( ν h + e n,h h γ k + µ k, h γ ) (cid:21) , (70)30here K n,h,i − ¯ K n,h,i = e c i kA h, e n,h + c i k ϕ ( kA h, ) C h O ( ν h + e n,h h γ )+ k i − X j =1 a ij (cid:20) e ( c i − c j ) kA h, [ f ( t n + c j k, K n,h,j ) − f ( t n + c j k, ¯ K n,h,j )]+( c i − c j ) k ϕ (( c i − c j ) kA h, ) C h O ( ν h + e n,h h γ ) (cid:21) = O ( e n,h + k ν h ) , and, for the last equality, (7), (H2c) and (61) have been used. Inserting this in (70)and using again (7), (H2c) and (61), it follows that U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + O ( k ν h + ke n,h + k µ k, ) . From here, e n +1 ,h = e kA h, e n,h + O ( ke n,h ) + O ( k ν h + k µ k, ) + ρ n +1 ,h , and using a discrete Gronwall Lemma and the first part of Theorem 21, the first part ofthe theorem follows for Dirichlet boundary conditions. For the second part, the secondpart of Theorem 21 must be used, apart from (22) and the additional regularity (69).As for Robin/Neumann boundary conditions, with similar arguments, K n,h,i − ¯ K n,h,i = e c i kA h, e n,h + c i k ϕ ( kA h, ) C h O ( e n,h ) + c i k ϕ ( kA h, ) C h O ( µ k, + e n,h k + ν h + e n,h h γ )+ k i − X j =1 a ij (cid:20) e ( c i − c j ) kA h, [ f ( t n + c j k, K n,h,j ) − f ( t n + c j k, ¯ K n,h,j )]+( c i − c j ) kϕ (( c i − c j ) kA h, ) C h O ( e n,h + k ( µ k, + e n,h k ))+( c i − c j ) k ϕ (( c i − c j ) kA h, ) C h O ( µ k, + e n,h k + ν h + e n,h h γ ) (cid:21) = e c i kA h, e n,h + O ( ke n,h + k µ k, + k ν h ) = O ( e n,h + k µ k, + k ν h ) , U n +1 ,h − ¯ U n +1 ,h = e kA h, e n,h + k ϕ ( kA h, ) C h O ( e n,h ) + k ϕ ( kA h, ) C h O ( µ k, + e n,h k + ν h + e n,h h γ )+ k ϕ ( kA h, ) C h O ( µ k, + µ k, + e n,h k + ν h + e n,h kh γ )+ k s X i =1 b i (cid:20) e (1 − c i ) kA h, [ f ( t n + c i k, K n,h,i ) − f ( t n + c i k, ¯ K n,h,i )]+(1 − c i ) kϕ ((1 − c i ) kA h, ) C h O ( ke n,h + k µ k, + k µ k, + k ν h + k h γ e n,h )+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( kµ k, + kµ k, + e n,h k + ν h + e n,h h γ )+(1 − c i ) k ϕ ((1 − c i ) kA h, ) C h O ( µ k, + µ k, + e n,h k + ν h + e n,h kh γ ) (cid:21) = e kA h, e n,h + O ( ke n,h + k µ k, + k µ k, + k ν h ) . From this, e n +1 ,h = e kA h, e n,h + O ( ke n,h + k µ k, + k µ k, + k ν h ) + ρ n +1 ,h , so that, using the first part of Theorem 21 and the classical argument of convergence, e n,h = O ( k + ε h + kν h + kµ k, + k µ k, ). Again, under the second set of hypothesesin Theorem 21 and using (22) and the regularity (69), the finer result e n,h = O ( k + kε h + η h + kµ k, + k µ k, + kν h ) is achieved. In this section, we show some numerical experiments which corroborate the previousresults. For that, we have considered the following set of problems u t ( t, x ) = u xx ( t, x ) + u ( t, x ) + h ( t, x ) , x ∈ [0 , , t ∈ [0 , ,u (0 , x ) = u ( x ) , (71)where the boundary conditions are either Dirichlet u ( t,
0) = g ( t ) , u ( t,
1) = g ( t ) , (72)or mixed (Dirichlet/Neumann), u ( t,
0) = g ( t ) , u x ( t,
1) = g ( t ) , (73)and where h ( t, x ), u ( x ) and g ( t ) , g ( t ) are such that the exact solution of the problemis either u ( x, t ) = x ( x −
1) cos( x + t ) or u ( x, t ) = cos( x + t ) . (74)32n the first place, as space discretization we have considered the symmetric 2nd-orderdifference scheme for which, in the Dirichlet case (72), A h, = tridiag(1 , − , /h , C h [ g ( t ) , g ( t )] T = [ g ( t ) , , . . . , , g ( t )] T /h , (75)and, in the Dirichlet/Neumann case (73), A h, = − . . . − − − , C h (cid:20) g ( t ) g ( t ) (cid:21) = h g ( t )0...0 h g ( t ) . (76)All differential problems (71) with boundary conditions (72)-(73) satisfy hypotheses(A1)-(A5) with X = C ([0 , Z = C ([0 , ε h , η h being O ( h ) for (75) and ε h = O ( h ), η h = O ( h ) for (76). Besides, the considered solutions(74) and f are so smooth that all conditions of regularity in the paper are satisfied.Finally, although we do not provide a proof for the conditions (23) and (41) with themaximum norm, it can be numerically verified that those conditions hold uniformly on h for A h, in (75) and (76). Moreover, (41) points out that (34) is also satisfied in thecontinuous case for the regular functions u ( t ) which are considered, although its proofis not an aim of this paper either. k 1e-3 5e-4 2.5e-4 1.25e-4Local error 9.7450e-4 4.8158e-4 2.3693e-4 1.1580e-4Order 1.02 1.02 1.03Global error 1.3461e-3 6.6579e-4 3.2779e-4 1.6034e-4Order 1.02 1.02 1.03Table 2: Local and global error when integrating Dirichlet problem with vanishingboundary conditions with the classical approach (16) associated to the second-ordermethod (77).We first show the results which are obtained when integrating problem (71) associ-ated to the Dirichlet boundary conditions with the Lawson method which is constructedwith the second-order RK tableau 0 01 1 0
12 12 . (77)33 k=1e-3 k=5e-4 k=2.5e-4 k=1.25e-42e-3 1.2404e+2 6.1563e+1 3.0339e+1 1.4751e+11e-3 4.9902e+2 2.4903e+2 1.2404e+2 6.1563e+15e-4 1.9990e+3 9.9902e+2 4.9902e+2 2.4903e+2Table 3: Local error when integrating Dirichlet problem with non-vanishing boundaryconditions with the classical approach (16) associated to the second-order method (77).h k=1e-3 k=5e-4 k=2.5e-4 k=1.25e-42e-3 6.7023e+1 3.3264e+1 1.6394e+1 7.9729e+01e-3 2.6962e+2 1.3455e+2 6.7023e+1 3.3264e+15e-4 1.0801e+3 5.3977e+2 2.6962e+2 1.3455e+2Table 4: Global error when integrating Dirichlet problem with non-vanishing boundaryconditions with the classical approach (16) associated to the second-order method (77).When considering the solution in (74) which vanishes at the boundary, the classicalapproach (16) shows local and global order 1 in time, as shown in Table 2 for h =5 × − , for which the error in space is negligible. (This corroborates Theorems 2 and3.) However, when the solution does not vanish at the boundary, although the local andglobal orders are still 1, the errors are very big and even grow when h diminishes. Thiswas justified in Theorems 6 and 7 and can be observed in Tables 3 and 4. However, thatbad behaviour can be solved by using the suggested approach (30), where every termon the boundary can be calculated in terms of data, without resorting to numericaldifferentiation. In such a way, local and global order 2 is obtained in Table 5 with h = 5 × − , for which the error in space is again negligible and does not grow with h .This corroborates Theorems 11 and 12. On the other hand, with this method, it is evenpossible to achieve local order 3 with formula (58), although numerical differentiationis required to approximate the boundary of the first derivative in space of the exactsolution, as it is thoroughly explained in Remark 13. For that, we have considered the2-BDF formula, for which ν h = O ( h ) and, as already predicted by the first part ofTheorems 16 and 17, Table 6 shows local order near 3 and global order 2, but with asize of errors quite smaller than those of Table 5. In this subsection we show the results which are obtained when considering (71) asso-ciated to the mixed Dirichlet/Neumann boundary conditions in (73) and integrating itin time with the Lawson method associated to the third order Heun method0
13 1323 . (78)34 1e-03 5e-4 2.5e-4 1.25e-4Local error 1.5664e-7 3.9176e-8 9.7933e-9 2.4473e-9Order 2.00 2.00 2.00Global error 8.2929e-7 2.0714e-7 5.1712e-8 1.2903e-8Order 2.00 2.00 2.00Table 5: Local and global error when integrating Dirichlet problem with nonvanishingboundary conditions with the suggested approach (30) associated to the second-ordermethod (77), with which no numerical differentiation is required.k 8e-3 4e-3 2e-3 1e-3Local error 1.3126e-7 1.6797e-8 2.1350e-9 2.7857e-10Order 2.97 2.98 2.94Global error 5.9892e-7 1.4972e-7 3.7367e-8 9.2309e-9Order 2.00 2.00 2.02Table 6: Local and global error when integrating Dirichlet problem with nonvanishingboundary conditions with the suggested approach (58) associated to the second-ordermethod (77), for which numerical differentiation is required.We have centered on the solution of (74) which does not vanish at the boundary. Theclassical approach shows no convergence either on the local or global error where thetimestepsize diminishes, as it is justified in Theorems 4 and 5, and it is shown in Table7. (Notice the different behaviour with respect to the classical approach in Tables 3and 4. Here the errors do not diminish with k but, although not shown here for thesake of brevity, they neither grow when h diminishes as it happens in those tables. Thisis due to the fact that now every c i is different from 1). However, we can get local andthus global order 3 with our modified approach (58), by calculating the terms on theboundary following again Remark 13. For the Dirichlet boundary condition, we haveused numerical differentiation in space with the 2-BDF formula and, for the Neumannone, numerical differentiation in time with the 3-BDF scheme. In such a case, Theorem16 as well as the second part of Theorem 17 apply, with ν h = O ( h ) and µ k, = O ( k ).Therefore, when the error in space is negligible, order 3 in the timestepsize should beseen when, as Table 8 corroborates. 35 0.2 0.1 0.05 0.025Local error 9.7639e-1 9.8964e-1 9.9108e-1 9.8877e-1Order -0.02 -0.00 0.00Global error 5.3822e-1 5.3736e-1 5.3613e-1 5.3439e-1Order 0.00 0.00 0.00Table 7: Local and global error when integrating mixed D/N problem with nonvan-ishing boundary conditions with the classical approach associated to the third-ordermethod (78), h = 10 − .k 0.2 0.1 0.05 0.025Local error 1.3911e-3 1.7489e-3 2.1806e-5 2.7212e-6Order 2.99 3.00 3.00Global error 1.5136e-3 2.3369e-4 2.9913e-5 3.6533e-6Order 2.70 2.97 3.03Table 8: Local and global error when integrating mixed D/N problem with nonvanish-ing boundary conditions with the suggested approach (58) associated to the third-ordermethod (78), for which numerical differentiation is required, h = 10 − . Finally, we show that local and global order 4 can be obtained when integrating intime with the Lawson method associated to the fourth-order RK method0
13 1323 −
11 1 −
18 38 38 18 . (79)Nevertheless, we point out that it is necessary to take condition (61) into account.In our problem (71), γ = 1, as it was justified in Remarks 13 and 18. For the sakeof brevity, we will center on the Dirichlet boundary condition in (72), and we willdirectly integrate that problem with the suggested formulas (63) by inserting the neededboundaries in an exact way from the known solution. As Table 9 shows, local and globalorder 4 are observed in that way. However, when not knowing the exact solution, thoseboundaries must be calculated in terms of data following Remark 18. For that, wehave considered again the 3 (resp. 2)-BDF formula for the numerical differentiationin time (resp. in space) and, according to Theorems 21 and 22, global order 4 in thetimestepsize should be observed when the error in space is negligible and (61) holdsfor some constant C . However, as the error in space is just of second order, in orderthat the error in space is negligible with respect to that in time, h must be quite smallwith respect to k , and then the global error exploits to infinity with the parameters of36 0.2 0.1 0.05 0.025Local error 1.8356e-4 1.0396e-5 6.1679e-7 3.7509e-8Order 4.14 4.08 4.04Global error 1.9072e-4 9.3054e-6 5.4646e-7 3.5333e-8Order 4.36 4.09 3.95Table 9: Local and global error when integrating Dirichlet problem with nonvanishingboundary conditions with the suggested approach (63) associated to the fourth-ordermethod (79), when the terms on the boundary are exactly provided, h = 5 × − .k 2.5e-2 1.25e-2 6.25e-3 3.125e-3Local error 3.4537e-8 2.0441e-9 1.1954e-10 6.8247e-12Order 4.08 4.10 4.13Global error 3.3314e-8 2.0054e-9 1.1968e-10 7.0050e-12Order 4.05 4.07 4.09Table 10: Local and global error when integrating Dirichlet problem with nonvanishingboundary conditions with the suggested approach (63) associated to the fourth-ordermethod (79), when the terms on the boundary are calculated through numerical dif-ferentiation and Gauss-Lobatto collocation space discretization is used.Table 9 because condition (61) is not satisfied for a suitable constant C .In spite of all this, the problem can be solved by considering a more accurate spacediscretization. Thus we have considered a Gauss-Lobatto collocation space discretiza-tion with 17 nodes, for which the error in space is nearly of the order of rounding errorsfor this problem. Besides, the space grid is quite moderate, so condition (61) is veryweak in this case. Considering then numerical differentiation in time as before and nu-merical differentiation in space through the derivation of the corresponding collocationpolynomials, the results in Table 10 are obtained, where both local and global order 4are achieved. Acknowledgements
This research has been supported by Ministerio de Ciencia e Innovaci´on and RegionalDevelopment European Funds through project MTM2015-66837-P and by Junta deCastilla y Le´on and Feder through project VA024P17.
References [1] I. Alonso–Mallo, B. Cano and N. Reguera,
Avoiding order reduction when inte-grating linear initial boundary value problems with exponential splitting methods ,37ublished online in IMA J. Num. Anal., doi: 10.1093/imanum/drx047[2] I. Alonso–Mallo, B. Cano and N. Reguera,
Analysis of order reduction when in-tegrating linear initial boundary value problems with Lawson methods , Applied Nu-merical Mathematics, 118 (2017) , pp. 64-74.[3] I. Alonso–Mallo, B. Cano and N. Reguera,
Avoiding order reduction when inte-grating linear initial boundary value problems with Lawson methods , IMA Journal ofNumerical Analysis 37 (2017), pp. 2091-2119.[4] I. Alonso–Mallo, B. Cano and N. Reguera,
Avoiding order reduction when integrat-ing reaction-diffusion boundary value problems with exponential splitting methods ,arXiv:1705.01857, submitted for publication.[5] I. Alonso–Mallo and C. Palencia,
On the convolutions operators arising in thestudy of abstract initial boundary value problems , Proc. Royal Soc. Edinburgh. (1996), 515–539.[6] B. Cano and M. J. Moreta,
Exponential quadrature rules without order reductionfor integrating linear initial boundary value problems , to be published in SIAM J.Num. Anal.[7] B. Cano and N. Reguera,
Avoiding order reduction when integrating nonlinearSchr¨odinger equation with Strang method , J. Comp. Appl. Math. 316 (2017), 86–99.[8] L. Einkemmer and A. Ostermann,
Overcoming order reduction in diffusion-reactionsplitting. Part 1: Dirichlet boundary conditions , SIAM J. Sci. Comput. (3) (2015),A1577–A1592.[9] L. Einkemmer and A. Ostermann, Overcoming order reduction in diffusion-reactionsplitting. Part 2: Oblique boundary conditions , SIAM J. Sci. Comput. (2016)A3741–A3757.[10] E. Faou, A. Ostermann and K. Schratz, Analysis of exponential splitting methodsfor inhomogeneous parabolic equations , IMA J. Numer. Anal. (1) (2015), 161–178.[11] M. Hochbruck and A. Ostermann, Explicit exponential Runge-Kutta methods forsemilinear parabolic problems , SIAM J. Num. Anal. (2005), 1069–1090.[12] M. Hochbruck and A. Ostermann, Exponential integrators , Acta Numerica (2010)209-286.[13] J. D. Lawson,
Generalized Runge-Kutta processes for stable systems with largeLipschitz constants , SIAM J. Numer. Anal. (1967) 372–380.[14] C. Palencia and I. Alonso–Mallo, Abstract initial-boundary value problems, Pro-ceedings of the Royal Society of Edinbourgh. Section A- Mathematics. 124, (1994)879 - 908. 3815] A. Pazy, Semigroups of Linear Operators and Applications to Partial DifferentialEquations , Series: Applied Mathematical Sciences, Vol. 44, Springer, New York,Berlin, Heidelberg, Tokyo, 1983.[16] J. M. Sanz-Serna,