Avoiding order reduction phenomenon for general linear methods when integrating linear problems with time dependent boundary values

Journal of Computational and Applied Mathematics 439 (2024) 115629 Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam Avoiding order reduction phenomenon for general linear methods when integrating linear problems with time dependent boundary values I. Alonso-Mallo a,∗,1, N. Reguera b,1 aIMUVA, Departamento de Matemática Aplicada, Universidad de Valladolid, Paseo de Belén, 7, 47011 Valladolid, Spain bIMUVA, Departamento de Matemáticas y Computación, Escuela Politécnica Superior, Universidad de Burgos, Avda. Cantabria, 09006 Burgos, Spain article info Article history: Received 7 June 2023 Received in revised form 12 September 2023 Keywords: Order reduction General linear methods Initial boundary value problems abstract When applied to stiff problems, the effective order of convergence of general linear methods is governed by their stage order, which is less than or equal to the classical order of the method. This produces an order reduction phenomenon, present in all general linear methods except those with high stage order, in a manner similar to that observed in other time integrators with internal stages. In this paper, we investigate the order reduction which arises when general linear methods are used as time integrators when using the method of lines for solving numerically initial boundary value problems with time dependent boundary values. We propose a technique, based on making an appropriate choice of the boundary values for the internal stages, with which it is possible to recover one unit of order, as we prove in this work. As expected, this implies a considerable improvement for the general linear methods suffering order reduction. Moreover, numerical experiments show that the improvement is not only in these cases, but that, even when the order reduction is not expected, the size of the errors is drastically reduced by using the technique proposed in this paper. ©2023 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/). 1. Introduction General linear methods (GLMs), introduced by Burrage and Butcher in [1], are a class of multistep-multistage methods used for the time integration of systems of ordinary differential equations [2–7]. Many of the time integrators most commonly used in practice for numerically solving ordinary differential equations are particular cases of GLM, including Runge–Kutta methods and linear multistep methods (see e.g. [8,9]). When a GLM is used for the time discretization of a stiff system of ordinary differential equations, the order observed in practice can be lower than the classical order of the GLM. This order reduction phenomenon is ubiquitous, and it has been observed for many schemes with internal stages [10–13], and also for exponential schemes [14–16]. ∗Corresponding author. E-mail addresses: [email protected] (I. Alonso-Mallo), [email protected] (N. Reguera). 1This research has been supported by project PGC2018-101443-B-100, Ministerio de Ciencia, Innovación y Universidades (Spain) and by Junta de Castilla y León (Grant numbers VA169P20 and VA193P20). https://doi.org/10.1016/j.cam.2023.115629 0377-0427/©2023 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC license (http://creativecommons.org/ licenses/by-nc/4.0/). I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 The order reduction suffered by GLMs has been studied in [9] for the particular case when they are applied to Prothero–Robinson test problems [17] u′(t)=λ(u(t)−φ(t)) +φ′(t),t∈ [t0,T], u(t0)=φ(t0). More precisely, in [9], it is shown that the effective global order of a GLM is equal to its stage order qand, for GLMs with Runge–Kutta stability which are A(0)-stable, it is equal to q+1 when the stability function R(z) of the underlying Runge–Kutta method is such that R(+∞)= 1. Therefore, for GLMs with these stability properties, the order reduction phenomenon is observed when q+1 is less than the classical order p, and, for general GLMs, order reduction is observed if q<p. It follows that only with GLMs with high stage order, i.e., q=por q=p−1, the order reduction phenomenon will not appear. Consequently, recent literature on general linear methods applied to stiff problems has dealt with methods that have high stage order, such as Diagonally Implicit Multistage Integration methods (DIMSIMs) [18–20] and also implicit-explicit methods (IMEX) [21]. An important case of very stiff high dimensional systems of ordinary differential equations arises when a time evolution partial differential equation is discretized in space with the idea of using the method of lines to obtain a full discretization (see e.g. [22,23], where GLMs are used for the time integration). Although the order reduction is not observed in the special case when the solution and its derivatives vanish sufficiently at the boundary of the spatial domain where the solutions are defined (see [24] for the case in which Runge–Kutta methods are used as integrators in time), order reduction arises in the more general and interesting case of time dependent boundary conditions. Due to the importance of this phenomenon, we can find in the literature the proposal of different techniques that allow to avoid the order reduction, when a initial boundary value problem (IBVP) is full discretized in time and space, for the case of some schemes with internal stages, such as Runge–Kutta [10,12], Rosenbrock methods [11], Spectral-fractional step Runge–Kutta [12], as well as for exponential splitting schemes [14–16]. The aim of this paper is, on the one hand, to study the order reduction that GLMs suffer when they are used to solve well posed linear abstract IBVP (see [24,25] for a detailed description of this type of problems). Since it is assumed that the differential operator of the problem is the infinitesimal generator of a C0−semigroup when the boundary values vanish, both hyperbolic and parabolic cases are included. We prove that the effective local order is q+1; then the effective global order is also q+1 when the stability function R(z) of the underlying Runge–Kutta method is such that R(+∞)= 1, similarly to the conclusions in [9] for Prothero–Robinson test problems. The numerical experiments that we show in this paper confirm this results. On the other hand, we also prove that, in the general case of time dependent boundary values, a unit of order can be recovered by an appropriate choice of the boundary values of the GLM stages. In this way, with the technique we propose, we prove that the local order is min(q+2,p+1) instead of order q+1 that would be obtained with a standard implementation. Then, for GLMs that present order reduction, one unit of order is gained, with the improvement that this entails. But, as we will shown in the numerical experiments, even in the particular case of GLMs for which no order reduction is expected, the size of the errors is much smaller with our technique, also providing an improvement to the standard implementation. Furthermore, we note that the computational cost of the technique we propose is negligible compared to that of the method itself. The organization of the paper is as follows. Section 2is devoted to present in an abstract way the IBVP which we want to solve numerically, and the GLMs that we will use for this. In Section 3, we prove the order reduction that GLMs suffer when using a standard implementation We also propose the boundary values for the stages that must be used to avoid the order reduction at least in one unit. The numerical experiments, which confirm the previous results, are included in Section 4. 2. Discretization of an IBVP by means of a GLM Let Xand Ybe two Banach spaces, D(A)⊂Xa dense subspace of X, and let A:D(A)⊂X→X,B:D(A)⊂X→Y, be two linear operators. We consider an abstract non homogeneous linear IBVP u′(t)=Au(t)+f(t),t∈ [0,T], u(0) =u0, Bu(t)=g(t)∈Y. (1) where Aand Bare operators for which we assume, as in [24,25], the following hypothesis: (A1) The operator Bis onto. (A2) Ker(B)=D(A0) is a dense subspace of Xand A0=A|D(A0)is the infinitesimal generator of a C0semigroup {S(t)}t≥0 in Xwith type ω. (A3) Let If z> ω, the steady state problem Ax =zx, 2 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 Bx =v∈Y, possesses a unique solution denoted by x=K(z)vand there exists a constant Csuch that the linear operator K(z):Y→D(A) satisfies, for any z≥ω0> ω, ∥K(z)v∥X≤C∥v∥Y. (A4) The solution uin (1) satisfies u(t)∈D(A) for t∈ [0,T]and it is smooth enough in time. We will consider a GLM for the numerical time integration of the IBVP (1). For this, let Nbe a natural number and k=T/Nwill denote the time step and tn=nk,n=0,...,N. A GLM of sstages and rsteps applied to (1) is given by the equations U[n] i=k s ∑ j=1 aij[AU[n] j+f(tn+cjk)] + r ∑ j=1 uiju[n] j,i=1,2,...,s,(2) u[n+1] i=k s ∑ j=1 bij[AU[n] j+f(tn+cjk)] + r ∑ j=1 viju[n] j,i=1,2,...,r,(3) where the stages U[n] iare approximations to the values u(tn+cik) of order qof the solution, U[n] i=u(tn+cik)+O(kq+1),i=1,2,...,s and u[n] iare approximations of order pto certain linear combinations of derivatives of the solution at the point tn, u[n] i= p ∑ l=0 qilklu(l)(tn)+O(kp+1),i=1,2,...,r.(4) The values pand q, which are the classical order and the stage order (see Definition 2.1) of the GLM, will be important throughout this work. These equations of a GLM can be rewritten by means of the abscissa vector c= [c1,...,cs]T, the vectors q0= [q1,0,...,qr,0]T,q1= [q1,1,...,qr,1]T, ..., qp= [q1,p,...,qr,p]T, and the matrices A= [aij] ∈ Rs×s,U= [uij] ∈ Rs×r, B= [bij] ∈ Rr×s,V= [vij] ∈ Rr×rand C=diag(c1,...,cs)∈Rs×s. Next, we are going to mention some definitions and properties of GLMs that will be of interest for the study we will carry out next. The GLM is strictly zero-stable when all the eigenvalues of the coefficient matrix Vare inside of the unit circle an it has a unique eigenvalue on the unit circle. As a consequence of the strict zero-stability, we can assure that there exists lim n→∞ Vn=˜ V=uvT where u,v∈Rrare vectors satisfying ˜ Vu =u,vT˜ V=vT, and vTu=1 (cf. [8]). We use the notation e= [1,...,1]Tand cl= [cl 1,...,cl s]Tand we consider the following values: γ0=e−Uq0 γk=ck k!−Ack−1 (k−1)!−Uqk,k=1,...,p.⎫ ⎬ ⎭ Then, a way to define the stage order is as follows. Definition 2.1. The maximum number q≤psuch that γk=0, k=0,1,...,q, is the stage order of the GLM method. In the rest of this work, we assume the stage consistency conditions γ0=0⇔Uq0=e, γ1=0⇔Ae +Uq1=c,}(5) or, equivalently, we assume that q≥1. Also, we will make the hypothesis that, for the GLM that is being used, the stage order q<p, so that the order reduction phenomenon will appear for the local error (see Theorem 3.1). On the other hand, we also consider the values ˆγ0=q0−Vq0 ˆγk= k ∑ l=0 ql (k−l)!−Bck−1 (k−1)!−Vqk,k=1,...,p.⎫ ⎪ ⎬ ⎪ ⎭ 3 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 Table 1 Order conditions for GLMs of order p≤4. Order porder condition 0ˆγ0=0, γ0=0 1ˆγ1=0, γ1=0 2ˆγ2=0 3ˆγ3=0 ˜ VBγ2=0 or Bγ2=0 4ˆγ4=0 ˜ VBγ3=0 ˜ VBAγ2=0 ˜ VBCγ2=0 Then, the conditions ˆγk=0,k=0,1,...,p,(6) are necessary (but not sufficient) for the GLM to have order p. In the rest of this work, we also assume the consistency conditions ˆγ0=0⇔Vq0=q0, ˆγ1=0⇔Be +Vq1=q0+q1, which imply, along with (5), that the GLM has classical order p≥1. Order conditions can be deduced in order to obtain a GLM of any desired order. See for example [8], where they are displayed up to order p=6. We have listed in Table 1 these conditions up to the order p=4. To approximate the solution of (1) by means of a GLM, the first step is to solve the equation for the internal stages (2), which can be rewritten as (I⊗I−A⊗kA)U[n]=(U⊗I)u[n]+k(A⊗I)F[n+c](7) where U[n]= [U[n] 1,...,U[n] s]Tis the vector of internal stages, u[n]= [u[n] 1,...,u[n] r]Tis the approximation at t=tnand F[n+c]= [f(tn+c1k),...,f(tn+csk)]T. When the internal stages have been computed, the approximation u[n+1]is obtained with (3), which can be rewritten as u[n+1]=(B⊗kA)U[n]+k(B⊗I)F[n+c]+(V⊗I)u[n].(8) To obtain a unique solution of the Eqs. (7), it is necessary to incorporate the boundary values of the internal stages. If we choose the boundary values G[n]= [G[n] 1,...,G[n] s]T, we obtain the following equations for the stages, (I⊗I−A⊗kA)U[n]=(U⊗I)u[n]+k(A⊗I)F[n+c], BU[n]=G[n].}(9) In order to see that (9) has a unique solution, it suffices to consider the problems (I⊗I−A⊗kA0)U[n,0]=(U⊗I)u[n]+k(A⊗I)F[n+c], BU[n,0]=0,}(10) and (I⊗I−A⊗kA)U[n,b]=0, BU[n,b]=G[n].}(11) We endow the space Xswith the usual norm product. Then, the solvability of (10) and (11) comes from the following two lemmas (see [10]). Lemma 2.2. There exists a constant C >0such that the operator (I⊗I−A⊗kA0)−1:Xs→Xsis boundedly invertible and ∥(I⊗I−A⊗kA0):D(As 0)⊂Xs→Xs∥ ≤ C for k >0small enough. Lemma 2.3. For k >0small enough and W = [W1,...,Ws]T∈Ys, the problem, (I⊗I−A⊗kA)V=0, BV =W.}4 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 possesses a unique solution V := K((kA)−1)W∈Xs. Moreover, there exists a constant C >0such that, ∥K((kA)−1)W∥ ≤ C∥W∥. 3. Local error analysis of the time semidiscrete problem In order to define the time semidiscrete local error, we consider U[n+c]= [u(tn+c1k),...,u(tn+csk)]T, and we use the following notation G[n+c]=BU[n+c]= [Bu(tn+c1k),...,Bu(tn+csk)]T = [g(tn+c1k),...,g(tn+csk)]T.(12) These boundary values are the ones used for the stages in the standard method of lines. Using the boundary values (12), we define U[n+c]= [U[n+c] 1,...,U[n+c] s]Tas the solution of (I⊗I−A⊗kA)U[n+c]=(U⊗I)˜ u(tn)+k(A⊗I)F[n+c], BU[n+c]=G[n+c],}(13) where ˜ u(tn)= [˜ u1(tn),...,˜ ur(tn)]T, is given taking into account (4) with exact values of the solution and its derivatives, that is ˜ ui(tn)= p ∑ l=0 qilklu(l)(tn),i=1,...,r.(14) Notice that we can consider ˜ u(tn) as the optimal approximation of the exact solution at t=tnthat we expect to obtain with the GLM. Then, let u[n+1]be given by u[n+1]=(B⊗kA)U[n+c]+k(B⊗I)F[n+c]+(V⊗I)˜ u(tn),(15) Then, we define the time semidiscrete local error at tn+1as follows. ρn+1=˜ u(tn+1)−u[n+1]. Theorem 3.1. Let u be the solution of (1) and let us assume the hypotheses (A1)–(A4). We consider the GLM given by the stages Eq. (13), with the boundary values G[n+c]defined in (12), and the final step (15). We assume that q, the stage order of the GLM, is smaller that p, the classical order. Then the time semidiscrete local errors ρn=˜ u(tn)−u[n],1≤n≤N, satisfy ∥ρn∥ ≤ Ckq+1,for k >0,(16) where the constant C depends only on the derivatives of u, the GLM method and the differential operator A. Proof. Let εn,kbe the value satisfying (I⊗I−A⊗kA)U[n+c]=(U⊗I)˜ u(tn)+k(A⊗I)F[n+c]+εn,k.(17) Then, εn,k=U[n+c]−(A⊗kI)(AU[n+c]+F[n+c])−(U⊗I)˜ u(tn) =U[n+c]−(A⊗kI)U′[n+c]−(U⊗I)˜ u(tn) where we have used the notation U′[n+c]= [u′(tn+c1k),...,u′(tn+csk)]T. Expanding U[n+c]and U′[n+c]into Taylor’s series around tn, we obtain εn,k=γ0+ p ∑ l=0 γlklu(l)(tn)+O(kp+1) = p ∑ l=q+1 γlklu(l)(tn)+O(kp+1),(18) where qis the stage order of the GLM and we have assumed that q<p. Denoting En,k=U[n+c]−U[n+c], we subtract (17) from (13) and, taking into account that U[n+c]and U[n+c]share the same boundary values, we deduce that (I⊗I−A⊗kA)En,k=εn,k BEn,k=0}5 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 and, using Lemma 2.2, En,k=(I⊗I−A⊗kA0)−1εn,k. On the other hand, u[n+1]=(V⊗I)˜ u(tn)+(B⊗kA)U[n+c]+k(B⊗I)F[n+c] =(V⊗I)˜ u(tn)+(B⊗kA)U′[n+c] =˜ u(tn+1)−˜ u(tn+1)+(V⊗I)˜ u(tn)+(B⊗kA)U′[n+c] =˜ u(tn+1)−ˆγ0− p ∑ l=0ˆγlklu(l)(tn)+O(kp+1) =˜ u(tn+1)+O(kp+1),(19) where, as for the stages equation, we have expanded ˜ u(tn+1) and U′[n+c]into Taylor’s series around tnand we have used the order conditions (6). Finally, we subtract (15) from (19), ρn+1=˜ u(tn+1)−u[n+1] =(B⊗kA0)(I⊗I−A⊗kA0)−1εn,k and, from (18) and using that (B⊗kA0)(I⊗I−A⊗kA0)−1are well-defined and bounded for k>0, we deduce that ρn+1=O(kq+1). We will now propose a technique to increase the order by one unit in cases where order reduction is present. The key is to choose properly the boundary values for the stages. For this, we consider the following vector U[n+c,1]=(U⊗I)˜ u(tn)+(A⊗kA)U[n+c]+k(A⊗I)F[n+c] =(U⊗I)˜ u(tn)+(A⊗kA)U′[n+c](20) and we denote G[n+c,1]=BU[n+c,1]=(U⊗I)B˜ u(tn)+(A⊗kA)BU′[n+c] =(U⊗I)˜ G(tn)+(A⊗kA)G′[n+c],(21) where G′[n+c]= [g′(tn+c1k),...,g′(tn+csk)]T, and ˜ G(tn)= [B˜ u1(tn),...,B˜ ur(tn)]T, with B˜ u1(tn)= p ∑ l=0 qilklg(l)(tn),i=1,...,r. Using these boundary values, we define the s-dimensional vector U[n+c,1]as the solution of (I⊗I−A⊗kA)U[n+c,1]=(U⊗I)˜ u(tn)+k(A⊗I)F[n+c] BU[n+c,1]=G[n+c,1]}(22) and the second step of the GLM is given by u[n+1,1]=(V⊗I)˜ u(tn)+(B⊗kA)U[n+c,1]+k(B⊗I)F[n+c].(23) Then, we define the time semidiscrete local error at tn+1 ρn+1,1=˜ u(tn+1)−u[n+1,1]. Theorem 3.2. Let u be the solution of (1) and let us assume the hypotheses (A1)–(A4). We use (13) with the boundary values G[n+c,1]defined in (21) and the final step (23). Then the time semidiscrete local error ρn,1=˜ u(tn)−u[n,1],1≤n≤N, satisfy ∥ρn,1∥ ≤ Ckmin(q+2,p+1),for k >0,(24) where q and p are, respectively, the stage order and the classical order of the GLM, and the constant C depends only on the derivatives of u, the GLM method and the differential operator A. 6 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 Proof. We subtract (20) from (17) and we deduce that εn,k=U[n+c]−U[n+c,1].(25) Let εn,k,1be the value satisfying U[n+c,1]=(U⊗I)˜ u(tn)+(A⊗kA)U[n+c,1]+k(A⊗I)F[n+c]+εn,k,1.(26) Now, we subtract (20) from (26) and we have 0=(A⊗kA)(U[n+c,1]−U[n+c])+εn,k,1(27) and, using (18), εn,k,1=(A⊗kA)(U[n+c]−U[n+c,1]) =(A⊗kA)εn,k = p ∑ l=q+1 Aγlkl+1Au(l)(tn)+O(kp+1) =O(kq+2)+O(kp+1). On the order hand, we define u[n+1,1]=(V⊗I)˜ u(tn)+(B⊗kA)U[n+c,1]+k(B⊗I)F[n+c] =(V⊗I)˜ u(tn)+(B⊗kA)U[n+c]+k(B⊗I)F[n+c] +(B⊗kA)(U[n+c,1]−U[n+c]) =u[n+1]+(B⊗kA)εn,k=˜ u(tn+1)+O(kp+1)+O(kq+2).(28) where we have used (19). Now, we use the notation En,k,1=U[n+c,1]−U[n+c,1]. Then, from (20),(25) and (26), En,k,1=(I⊗I−A⊗kA0)−1εn,k,1. Finally, from (28) ρn+1,1=˜ u(tn+1)−u[n+1,1]+u[n+1,1]−u[n+1,1] =O(kp+1)+O(kq+2)+(B⊗kA)(U[n+c,1]−U[n+c,1]) =O(kp+1)+O(kq+2)+(B⊗kA0)(I⊗I−A⊗kA0)−1εn,k,1 =O(kp+1)+O(kq+2). 4. Numerical experiments In the previous section, Theorem 3.1 and 3.2 show that the technique we propose, based on using appropriate boundary conditions for the stages (20), allows us to increase the local order by one unit in the general case in which there is order reduction. In this section we show numerical experiments that support these theoretical results. Moreover, we show that also in the case of GLMs for which order reduction is not expected, our technique provides an improvement because, also in this case, it allows to obtain much smaller errors than with a standard implementation. In this section we are going to solve numerically the initial boundary value problem ut(x,t)=uxx(x,t)+f(x,t),x∈ [0,1] u(x,0) =u0(x) u(0,t)=g0(t),u(1,t)=g1(t) (29) with f(x,t)=et(x2−1), u0(x)=x2+1, g0(t)=etand g1(t)=2et, whose exact solution is u(x,t)=et(x2+1). This problem fits into the abstract IBVPs theory developed in [25], taking the spaces X=L2(0,1),Y=R2and defining, for u∈D(A)=H2(0,1),Au =uxx and Bu = [g0(0),g1(0)]T. For the spatial discretization of (29) we will consider the second-order symmetric difference scheme. We take J∈N, J≥3, and h=1/Jand we consider the nodes xj=jh,j=0,...,J−1. We denote Ah,0=1 h2tridiag(1,−2,1),Chg(t)=1 h2[g0(t),0,...,0,g1(t)]T. Notice that, since the solution of (29) is a second degree polynomial in space, the error in space vanishes since we are using second-order symmetric differences. In this way, we can better study the error in time, which is what we are interested in. 7 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 Table 2 Local and global error when integrating problem (29) with GLM31 with a standard implementation. For the discretization in space h=4×10−3has been considered. k=4×10−2k=2×10−2k=10−2k=5×10−3k=2.5×10−3 L2-local error 1.2931e−01 2.7007e−02 5.6031e−3 1.1553e−3 2.3611e−4 Order 2.26 2.27 2.28 2.29 L-global error 5.6650e−02 1.2217e−02 2.5743e−3 5.3406e−4 1.0924e−4 Order 2.21 2.25 2.27 2.29 Table 3 Local and global error when integrating problem (29) with GLM31 avoiding order reduction. For the discretization in space h=4×10−3has been considered. k=4×10−2k=2×10−2k=10−2k=5×10−3k=2.5×10−3 L2-local error 1.5281e−03 1.5434e−04 1.6616e−5 1.9643e−6 2.5371e−7 Order 3.31 3.22 3.10 3.00 L2-global error 9.8820e−04 1.2439e−04 1.5806e−5 2.0039e−6 2.5335e−7 Order 2.99 2.98 2.98 2.98 After applying this spatial discretization, we obtain the semidiscrete problem U′ h(t)=Ah,0Uh(t)+fh(t)+Chg(t) (30) where Uh(t)= [U1(t),...,UJ−1(t)]Tand fh(t)= [f(x1,t),...,f(xJ−1,t)]T. We are going to carry out the numerical integration of (30) with two different GLM studying in each case the order in time that is observed. First, we use the method GLM31, deduced in [9] that has classical order p=3 and stage order q=1. For this method, Vis a rank one matrix such that V=e1vT, where e1= [1,0,...,0]Tis the first vector of the canonical basis of Rp+1 and vTe1=1. In this way, Vn=e1vT=Vand we deduce that ˜ V=V. Then, the condition for order three is given by ˜ VBγ2=VBγ2=0. We include here the coefficients of the GLM because there is a small error in [9]: A=[0.5 0 0 0.6114715765267838 0.5 0 34.95200172589030 −1.545721169870557 0.5], U=⎡ ⎣ 1−0.5−0.2509325024749301 0 1−0.6114715765267838 −1 4 1 3 1−32.90628055601974 1 2−1 8 ⎤ ⎦, B=⎡ ⎢ ⎣ 36.20771429329408 −2.724800143643149 1.340965575511875 001 1−4 3 4−8 4 ⎤ ⎥ ⎦, v=[1−33.82387972516280 0.5214344963096989 −0.1632161031338774 ]T We can observe in Table 2 that, when we solve numerically the problem (29) with the spatial discretization previously explained and the scheme GLM31 for the time integration, with its standard implementation, local and global orders in time q+1=2 are obtained. This local order is expected from Theorem 3.1. Although in the general case the global order is expected to be 1, we get again the value 2 due to the summation by parts. The numerical results observed in Table 2 show that order reduction is taking place when we use the standard implementation of GLM31. In order to avoid this order reduction effect, now we are going to carry out the same numerical experiment, but using the implementation (22) that we propose in this work for avoiding the order reduction phenomenon. In Table 3, we display the local and global orders in time, observing that now, both of them are q+1+1=3, so the order reduction is completely avoided. These numerical results are in agreement with Theorem 3.2. Note that, when we use the boundary values proposed in Theorem 3.2 to avoid order reduction, we get the added advantage that the errors are much smaller (compare the errors in Table 3 with those in Table 2). This is clearly seen in Fig. 1. Now, we are going to consider the method GLM32 deduced in [9] with p=3 and q=2. The coefficients of such a method, are as follows A=[0.5 0 0 19.03844038247158 0.5 0 −3.961809715877235 −0.1911525347850310 0.5], U=[1−0.5 0 0 1−19.03844038247158 −0.125 0.333333333333333 1 4.652962250662266 0.09557626739251550 −0.125 ], 8 I. Alonso-Mallo and N. Reguera Journal of Computational and Applied Mathematics 439 (2024) 115629 Fig. 1. Global errors when integrating problem (29) with GLM31. Table 4 Local and global error when integrating problem (29) with GLM32 with a standard implementation. For the discretization in space h=4×10−3has been considered. k=4×10−2k=2×10−2k=10−2k=5×10−3k=2.5×10−3 L2-local error 2.3028e−04 2.43442e−05 2.5611e−06 2.6760e−07 2.7773e−08 Order 3.24 3.25 3.26 3.27 L2-global error 7.0327e−04 7.6149e−05 8.0456e−06 8.3339e−07 8.4674e−08 Order 3.21 3.24 3.27 3.30 Table 5 Local and global error when integrating problem (29) with GLM32 avoiding order reduction. For the discretization in space h=4×10−3has been considered. k=4×10−2k=2×10−2k=10−2k=5×10−3k=2.5×10−3 L2-local error 1.5405e−05 1.0631e−06 7.0663e−08 4.5948e−09 2.9492e−10 Order 3.86 3.91 3.94 3.96 L2-global error 1.4321e−04 7.9134e−06 4.3307e−07 2.3743e−08 1.3305e−09 Order 4.18 4.19 4.19 4.16 B=⎡ ⎢ ⎣ −1.624407520530026 −0.1715153809360207 0.5228985152929018 0 0 1 1−4 3 4−8 4 ⎤ ⎥ ⎦, v=[1 2.273024386173144 0.06285917517510861 −0.07334316836278161 ]T We can observe in Table 4 that when we integrate problem (29) considering the method GLM32 for the time integration, with its standard implementation, local and global order q+1=3 are obtained. The global order is 3, and it seems that there is no order reduction. However, if there were no order reduction, local order p+1=4 would be expected and this is not the case in Table 4. Therefore, we have used the implementation of the GLM32 that we propose in order to avoid order reduction and that is based in the boundary values of Theorem 3.2. The local and global errors displayed in Table 5 show that the local order predicted by Theorem 3.2,q+2=4, is reached. In addition, the global error sizes are much smaller than when using the standard implementation, which justifies the use of the new boundary values even when q+1=p, in which case no order reduction is expected (see Fig. 2). 9