Efficient exponential Rosenbrock methods till order four

Journal of Computational and Applied Mathematics 453 (2025) 116158 Available online 23 July 2024 0377-0427/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam Efficient exponential Rosenbrock methods till order four B. Cano a,∗,1, M.J. Moretab aIMUVA, Departamento de Matemática Aplicada, Facultad de Ciencias, Universidad de Valladolid, Paseo de Belén 7, 47011 Valladolid, Spain bIMUVA, Departamento de Análisis Económico y Economía Cuantitativa, Facultad de Ciencias Económicas y Empresariales, Universidad Complutense de Madrid, Campus de Somosaguas, Pozuelo de Alarcón, 28223 Madrid, Spain ARTICLE INFO Keywords: Exponential Rosenbrock methods Nonlinear reaction–diffusion problems Avoiding order reduction in time Efficiency ABSTRACT In a previous paper, a technique was described to avoid order reduction with exponential Rosenbrock methods when integrating initial boundary value problems with time-dependent boundary conditions. That requires to calculate some information on the boundary from the given data. In the present paper we prove that, under some assumptions on the coefficients of the method which are mainly always satisfied, no numerical differentiation is required to approximate that information in order to achieve order 4 for parabolic problems with Dirichlet boundary conditions. With Robin/Neumann ones, just numerical differentiation in time may be necessary for order 4, but none for order ≤3. Furthermore, as with this technique it is not necessary to impose any stiff order conditions, in search of efficiency, we recommend some methods of classical orders 2, 3 and 4 and we give some comparisons with several methods in the literature, with the corresponding stiff order. 1. Introduction It is well known that stiff ordinary differential systems are typically those for which an explicit integration with standard methods is not possible due to a lack of 𝐴-stability [1]. When integrating in space partial differential equations, the space discretized system is stiff because the eigenvalues associated to the discretization of the spatial differentiation operator can be infinitely large in modulus when the space grid is refined. Exponential methods have been developed in the literature in order to get a stable integration of stiff systems in an ‘explicit’ way [2]. Although exponential-type functions of matrices applied over vectors have to be calculated with these methods, the recent improvement of Krylov techniques to approximate them has made these methods a valuable tool to integrate some initial boundary value problems [3–5]. However, the phenomenon of order reduction which already turns up when integrating this type of problems with standard methods also turns up with exponential methods. More explicitly, when the boundary conditions are not periodic or do not satisfy enough conditions of annihilation on the boundary, the order of accuracy of the time integrators when integrating the space discretized system is smaller than the order of accuracy which is observed when integrating a nonstiff ODE. Some times convergence is even lost [6]. Because of that, restrictive stiff order conditions on the coefficients of exponential Runge–Kutta methods have been firstly suggested in the literature so as to achieve the desired order of accuracy when integrating semilinear parabolic problems with vanishing boundary conditions [7]. Later, another technique has been given in order to avoid that order reduction without having to impose restrictions on the coefficients [8–11]. Moreover, the latter technique is valid for time-dependent boundary conditions. ∗Corresponding author. E-mail addresses: [email protected] (B. Cano), [email protected] (M.J. Moreta). 1This research was supported by Junta de Castilla y León and Feder through project VA169P20 https://doi.org/10.1016/j.cam.2024.116158 Received 26 July 2023; Received in revised form 12 March 2024 Journal of Computational and Applied Mathematics 453 (2025) 116158 2 B. Cano and M.J. Moreta Furthermore, by assuming that the Jacobian of the vector field can be easily calculated, that extra information can be used with Rosenbrock methods so as to achieve a desired accuracy with less stages than their Runge–Kutta counterparts. Because of that, for the integration of reaction–diffusion initial boundary value problems, stiff order conditions have also been studied in [12] for exponential Rosenbrock methods. Moreover, many particular exponential Rosenbrock methods have been constructed satisfying those conditions in order to achieve a desired accuracy while trying to be as efficient as possible [3,12–16]. In contrast, in [17], a technique is suggested to avoid order reduction with any exponential Rosenbrock method of any classical order without imposing those stiff order conditions. For that, some intermediate initial boundary problems are considered for which the boundary values have to be calculated in terms of data. If analytic expressions for the latter are known, in [17] it is stated that those boundary values can be exactly calculated in order to get local order ≥2. However, in order to get local order ≥3, numerical differentiation is in general required either in space or in time. When that in space is necessary, a weak CFL condition is needed to prove convergence [17]. Although the latter condition is much less restrictive than that required when integrating with standard explicit Runge–Kutta methods, it would be better not to require it and not to resort to numerical differentiation for the ease of implementation. Our aim in this paper is to prove that, under some simplifying assumptions (which are satisfied by mainly all already built exponential Rosenbrock methods), numerical differentiation in space is not required to achieve local order 3and 4, so that no CFL condition is necessary then. Moreover, numerical differentiation in time is just necessary to achieve local order 4with Robin/Neumann boundary conditions, but not with Dirichlet ones and that order and neither to achieve local order 3. We remark that, when integrating with standard Rosenbrock methods, a technique was also given in [18] to avoid order reduction and no numerical differentiation in space was either required to achieve local order ≤4. Therefore, the conclusions in this paper are similar to those in [18] for standard Rosenbrock methods. The advantage of using exponential Rosenbrock methods instead of standard ones correspond to problems where a good preconditioner is not available to solve linear systems in an efficient way and it is therefore more recommendable to use Krylov techniques to approximate exponentials of matrices applied over vectors [16,19–21]. In any case, coming back to the above considerations on exponential Rosenbrock methods, in this paper we also recommend some particular ones to get global orders 2, 3 and 4 (the latter for parabolic problems in which, by a summation-by-parts argument, the local order coincides with the global one). We remind that the local error corresponds to the error after just one step while the global error is the error after the required steps to get a final time. As no stiff order conditions have to be imposed, there are more parameters to play with in order to achieve a cheaper computation. The paper is structured as follows. Section 2gives some preliminaries and particularizes the obtained formulas for the full discretization in [17] in order to get local order 𝑝+ 1 for 𝑝= 1,2,3when the method has classical order at least 𝑝. Section 3justifies how those formulas greatly simplify under the precise conditions on the coefficients of the method which are stated there. (That simplification must not be seen in the length of formulas, but in the fact that the required boundary values are easier to calculate and that the linear combination of 𝜑𝑗-functions of matrices applied over vectors is given, so that Krylov subroutines can be directly applied.) Section 4gives some theorems and remarks which justify that the simplifying assumptions are satisfied by basically every constructed method. Section 5gives a recommendation for methods of order 2,3and 4and, finally, Section 6shows convergence tables and a numerical comparison in CPU time with other methods in the literature. 2. Preliminaries We will assume that the problem to integrate can be written as 𝑢(𝑡) = 𝐴𝑢(𝑡) + 𝛹(𝑢(𝑡)) + ℎ(𝑡),0≤𝑡≤𝑇 , (1) 𝑢(0) = 𝑢0, 𝜕𝑢(𝑡) = 𝑔(𝑡), where ⋅denotes differentiation with respect to time and where 𝐴∶𝐷(𝐴)⊂ 𝑋 →𝑋is a linear differential operator defined in a subset of the Banach space 𝑋=𝐿∞( 𝛺, C), where 𝛺is a bounded domain in R𝑑and the supremum norm is considered. On the other hand, 𝜕∶𝑋→𝑌corresponds to a boundary operator which arrives at another Banach space 𝑌=𝐿∞(𝜕 𝛺)with the same supremum norm. We will assume hypotheses (A1)–(A9) in [17] and, more particularly, if (1) is real and the solution stays in the interval 𝐼, that 𝛹∈𝐶2(𝐼, R)and, if (1) is complex, that 𝛹is holomorphic in a region where the solution stays. Moreover, we will assume that ℎ∈𝐶2([0, 𝑇 ], 𝑋). As stated in [17,18] through [22–24], (A1)–(A4) guarantee that problem (1) is well-posed. However, more particular assumptions on the regularity of 𝑢and 𝛹are stated in (A5)–(A9) and in Theorem 3.1 in [17] in order to get the desired local order 𝑝+1 whenever the method has non-stiff global order ≥𝑝. We do not state them here for the sake of brevity, but it is there justified that all expressions which turn up in the rest of the paper exist. Another issue would be to justify in terms of the data of the problem that the exact solution is regular enough. Although that would be very interesting, it is not an aim of this paper. In any case, we refer to [25–28] for bounds of some derivatives of the exact solution in the linear case. (The last papers are focused on singularly perturbed problems, which are not necessarily the case here, but the bounds which are obtained there are also valid for not so small parameters.) General Rosenbrock exponential methods are determined by some coefficients 𝑐1,…, 𝑐𝑠and some coefficient functions in a Butcher tableau 𝑎𝑖,𝑗 (𝑧) = 𝑟 ∑ 𝑙=1 𝜆𝑖,𝑗,𝑙𝜑𝑙(𝑐𝑖𝑧), 𝑖 = 1,…, 𝑠, 𝑗 = 1,…, 𝑖 − 1, Journal of Computational and Applied Mathematics 453 (2025) 116158 3 B. Cano and M.J. Moreta 𝑏𝑖(𝑧) = 𝑟 ∑ 𝑙=1 𝜇𝑖,𝑙𝜑𝑙(𝑧), 𝑖 = 1,…, 𝑠, (2) where 𝜑𝑙(𝑧) = ∫1 0 𝑒(1−𝜃)𝑧𝜃𝑙−1 (𝑙− 1)! 𝑑𝜃, 𝑙 ≥1, and where it is assumed that 𝑖−1 ∑ 𝑗=1 𝑎𝑖𝑗 (0) = 𝑐𝑖.(3) For an autonomous ODE differential system of the type  𝑈(𝑡) = 𝐹(𝑈(𝑡)),(4) the numerical solution 𝑈𝑛+1 at time 𝑡𝑛+1 =𝑡𝑛+𝑘is given from the numerical solution 𝑈𝑛at time 𝑡𝑛through the following formulas 𝐾𝑛,𝑖 =𝑒𝑐𝑖𝑘𝐽𝑛𝑈𝑛+𝑘 𝑖−1 ∑ 𝑗=1 𝑎𝑖𝑗 (𝑘𝐽𝑛)𝐺𝑛(𝐾𝑛,𝑗 ), 𝑖 = 1,…, 𝑠, (5) 𝑈𝑛+1 =𝑒𝑘𝐽𝑛𝑈𝑛+𝑘 𝑠 ∑ 𝑖=1 𝑏𝑖(𝑘𝐽𝑛)𝐺𝑛(𝐾𝑛,𝑖),(6) where 𝐽𝑛=𝐹′(𝑈𝑛), 𝐺𝑛(𝑈) = 𝐹(𝑈) − 𝐽𝑛𝑈. (7) When the problem is non-autonomous, by rewriting  𝑈(𝑡) = 𝐹(𝑡, 𝑈(𝑡)) as an autonomous one, the method reads 𝐾𝑛,𝑖 =𝑒𝑐𝑖𝑘𝐽𝑛𝑈𝑛+𝑐𝑖𝑘𝑡𝑛𝜑1(𝑐𝑖𝑘𝐽𝑛)𝑉𝑛+𝑘 𝑖−1 ∑ 𝑗=1 𝑟 ∑ 𝑙=1 𝜆𝑖,𝑗,𝑙[𝜑𝑙(𝑐𝑖𝑘𝐽𝑛)[𝐹(𝑡𝑛,𝑗 , 𝐾𝑛,𝑗 ) − 𝑡𝑛,𝑗 𝑉𝑛−𝐽𝑛𝐾𝑛,𝑗 ] + 𝑐𝑖𝑘𝜑𝑙+1(𝑐𝑖𝑘𝐽𝑛)𝑉𝑛], 𝑈𝑛+1 =𝑒𝑘𝐽𝑛𝑈𝑛+𝑘𝑡𝑛𝜑1(𝑘𝐽𝑛)𝑉𝑛+𝑘 𝑠 ∑ 𝑖=1 𝑟 ∑ 𝑙=1 𝜇𝑖,𝑙[𝜑𝑙(𝑘𝐽𝑛)[𝐹(𝑡𝑛,𝑖, 𝐾𝑛,𝑖) − 𝑡𝑛,𝑖𝑉𝑛−𝐽𝑛𝐾𝑛,𝑖] + 𝑘𝜑𝑙+1(𝑘𝐽𝑛)𝑉𝑛],(8) where 𝑡𝑛,𝑗 =𝑡𝑛+𝑐𝑗𝑘, and 𝑉𝑛=𝜕𝐹 𝜕𝑡 (𝑡𝑛, 𝑈𝑛), 𝐽𝑛=𝜕𝐹 𝜕𝑈 (𝑡𝑛, 𝑈𝑛). On the other hand, we will consider a general space discretization for the differential operator 𝐴, such that, when applied to the elliptic problem 𝐴𝑢 =𝐹 , 𝜕𝑢 =𝑔, the nodal values on the grid are given by the solution 𝑈ℎ∈C𝑁of this system 𝐴ℎ,0𝑈ℎ+𝐶ℎ𝑔=𝑃ℎ𝐹+𝐷ℎ𝜕𝐹 , where 𝐴ℎ,0discretizes 𝐴restricted to Ker(𝜕),𝑃ℎis the nodal projection and 𝐶ℎ, 𝐷ℎ∶𝑌→C𝑁are other linear operators over functions on the boundary. We will assume hypotheses (H1)–(H3) in [17] and will consider the matrix 𝐽𝑛,ℎ,0=𝐴ℎ,0+𝛹′(𝑈𝑛 ℎ), which discretizes the Jacobian with respect to 𝑈of the vector field which defines (1) at 𝑡=𝑡𝑛=𝑛𝑘, where 𝑘is the timestepsize. In [17], a technique is suggested to achieve local order 𝑝+ 1 with a method which has classical global order at least 𝑝. For the precise values 𝑝= 1,2,3, the modified exponential Rosenbrock method to achieve that goal when integrating the non-autonomous problem (1) reads as follows. For 𝑝= 1, 𝐾𝑛,𝑖,ℎ =𝑒𝑐𝑖𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑐𝑖𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛)] +𝑘 𝑖−1 ∑ 𝑗=1 𝑟 ∑ 𝑙=1 𝜆𝑖,𝑗,𝑙[𝜑𝑙(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝐺𝑛,𝑗,ℎ +𝑐𝑖𝑘𝜑𝑙+1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝑃ℎ ℎ(𝑡𝑛)], 𝑈𝑛+1 ℎ=𝑒𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)] +𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) + 𝑡𝑛 ℎ(𝑡𝑛)] Journal of Computational and Applied Mathematics 453 (2025) 116158 4 B. Cano and M.J. Moreta +𝑘 𝑠 ∑ 𝑖=1 𝑟 ∑ 𝑙=1 𝜇𝑖,𝑙[𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐺𝑛,𝑖,ℎ +𝑘𝜑𝑙+1(𝑘𝐽𝑛,ℎ,0)[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕  𝐺𝑛]](9) where 𝐺𝑛,𝑗,ℎ =𝛹(𝐾𝑛,𝑗,ℎ) + 𝑃ℎℎ(𝑡𝑛,𝑗 ) − 𝑡𝑛,𝑗 𝑃ℎ ℎ(𝑡𝑛) − diag(𝛹′(𝑈𝑛 ℎ))𝐾𝑛,𝑗,ℎ, with 𝑡𝑛,𝑗 =𝑡𝑛+𝑐𝑗𝑘, and   𝐺𝑛=𝛹(𝑢(𝑡𝑛)) + ℎ(𝑡𝑛) − 𝑡𝑛 ℎ(𝑡𝑛) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛).(10) As for 𝑝= 2, we get 𝐾𝑛,𝑖,ℎ =𝑒𝑐𝑖𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ +𝑐𝑖𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)] + (𝑐𝑖𝑘)2𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) + 𝑡𝑛 ℎ(𝑡𝑛)] +𝑘 𝑖−1 ∑ 𝑗=1 𝑟 ∑ 𝑙=1 𝜆𝑖,𝑗,𝑙[𝜑𝑙(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝐺𝑛,𝑗,ℎ +𝑐𝑖𝑘𝜑𝑙+1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕  𝐺𝑛]], 𝑈𝑛+1 ℎ=𝑒𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)] +𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)[𝐶ℎ𝜕[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) + 𝑡𝑛 ℎ(𝑡𝑛)] − 𝐷ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛) ℎ(𝑡𝑛)]] +𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛) ℎ(𝑡𝑛)] +𝑘 𝑠 ∑ 𝑖=1 𝑟 ∑ 𝑙=1 𝜇𝑖,𝑙[𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐺𝑛,𝑖,ℎ +𝑘𝜑𝑙+1(𝑘𝐽𝑛,ℎ,0)[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕  𝐺𝑛−𝐷ℎ𝜕 𝐽(𝑡𝑛)  𝐺𝑛] +𝑘2𝜑𝑙+2(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)  𝐺𝑛+ ℎ(𝑡𝑛)]].(11) Finally, for 𝑝= 3, we have 𝐾𝑛,𝑖,ℎ =𝑒𝑐𝑖𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ +𝑐𝑖𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)] + (𝑐𝑖𝑘)2𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝐶ℎ𝜕[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) + 𝑡𝑛 ℎ(𝑡𝑛)] − 𝐷ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛) ℎ(𝑡𝑛)]] + (𝑐𝑖𝑘)3𝜑3(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛) ℎ(𝑡𝑛)] +𝑘 𝑖−1 ∑ 𝑗=1 𝑟 ∑ 𝑙=1 𝜆𝑖,𝑗,𝑙[𝜑𝑙(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝐺𝑛,𝑗,ℎ +𝑐𝑖𝑘𝜑𝑙+1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕  𝐺𝑛−𝐷ℎ𝜕 𝐽(𝑡𝑛)  𝐺𝑛] +(𝑐𝑖𝑘)2𝜑𝑙+2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)  𝐺𝑛+ ℎ(𝑡𝑛)]]. 𝑈𝑛+1 ℎ=𝑒𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)] +𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)[𝐶ℎ𝜕[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) + 𝑡𝑛 ℎ(𝑡𝑛)] − 𝐷ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛) ℎ(𝑡𝑛)]] +𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)[𝐶ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛) ℎ(𝑡𝑛)] − 𝐷ℎ𝜕[ 𝐽(𝑡𝑛)3𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛)2 ℎ(𝑡𝑛)]] +𝑘4𝜑4(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)3𝑢(𝑡𝑛) + 𝑡𝑛 𝐽(𝑡𝑛)2 ℎ(𝑡𝑛)] +𝑘 𝑠 ∑ 𝑖=1 𝑟 ∑ 𝑙=1 𝜇𝑖,𝑙[𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐺𝑛,𝑖,ℎ +𝑘𝜑𝑙+1(𝑘𝐽𝑛,ℎ,0)[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕     𝐺𝑛,𝑖 −𝐷ℎ𝜕 𝐽(𝑡𝑛)  𝐺𝑛] +𝑘2𝜑𝑙+2(𝑘𝐽𝑛,ℎ,0)[𝐶ℎ𝜕[ 𝐽(𝑡𝑛)  𝐺𝑛+ ℎ(𝑡𝑛)] − 𝐷ℎ𝜕[ 𝐽(𝑡𝑛)2  𝐺𝑛+ 𝐽(𝑡𝑛) ℎ(𝑡𝑛)]] +𝑘3𝜑𝑙+3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[ 𝐽(𝑡𝑛)2  𝐺𝑛+ 𝐽(𝑡𝑛) ℎ(𝑡𝑛)]].(12) where     𝐺𝑛,𝑖 =𝛹(𝑢(𝑡𝑛)) + ℎ(𝑡𝑛) − 𝑡𝑛 ℎ(𝑡𝑛) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) + 𝑐2 𝑖𝑘2 2[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)]. Now we must see how to calculate the terms on the boundary taking into account that our data is just 𝜕𝑢(𝑡) = 𝑔(𝑡),ℎ(𝑡)and 𝑢0in (1). As stated in [17], for 𝑝≥2, it is in principle necessary to resort to numerical differentiation either in space or in time for both Dirichlet and Robin/Neumann boundary conditions in order to approximate those boundary values. In order to avoid that as far as possible, we will see that many times some simplifications can be performed which allow to calculate the required boundaries exactly in terms of data. Journal of Computational and Applied Mathematics 453 (2025) 116158 5 B. Cano and M.J. Moreta 3. Further simplifications and calculation of required boundaries in terms of data In this section, we will see that, under the assumptions 𝑠 ∑ 𝑖=1 𝜇𝑖,1= 1, 𝑠 ∑ 𝑖=1 𝜇𝑖,𝑙 = 0 (𝑙= 2,…, 𝑟),(13) 𝑖−1 ∑ 𝑗=1 𝜆𝑖,𝑗,1=𝑐𝑖, 𝑖−1 ∑ 𝑗=1 𝜆𝑖,𝑗,𝑙 = 0 (𝑙= 2,…, 𝑟), 𝑖 = 1,…, 𝑠, (14) some terms in the general expressions (9),(11),(12) can be simplified. In a first place, the required terms on the boundary can be calculated in a much more direct way in terms of data, i.e. without resorting to numerical differentiation either in space nor in time, except for the case of Robin/Neumann boundary conditions and 𝑝= 3, for which numerical differentiation for the first derivative in time will be required. Secondly, some other terms not concerning the boundary also simplify under those assumptions and moreover, the evaluation of 𝐺𝑛,𝑖,ℎ can be performed except for terms which lead to 𝑂(𝑘𝑝+2)-residues in 𝑈𝑛+1 ℎand therefore do not change either the local neither the global order. (This last simplification just concerns 𝑝= 1 for both the equations on the stages and the solution and 𝑝= 2 for the stages. It corresponds to the use of  𝐺𝑛,𝑖,ℎ in the formulas below.) Moreover, we will gather together all terms in the same 𝜑𝑗so that Krylov subroutines can be directly applied to calculate a linear combination of those matrix functions applied over the corresponding vectors. 3.1. 𝑝= 1 We notice that, in this case, for the stages, just the term 𝜕𝑢(𝑡𝑛) = 𝑔(𝑡𝑛)on the boundary is required (see (9)). Moreover, taking into account that 𝐺𝑛,𝑗,ℎ =𝛹(𝐾𝑛,𝑗,ℎ) + 𝑃ℎℎ(𝑡𝑛) − 𝑡𝑛𝑃ℎ ℎ(𝑡𝑛) − diag(𝛹′(𝑈𝑛 ℎ))𝐾𝑛,𝑗,ℎ +𝑂(𝑘2),(15) and the first part of (14), the terms in 𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)without boundaries can be simplified to 𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑐𝑖𝑃ℎℎ(𝑡𝑛) + ∑𝜆𝑖,𝑗,1 𝐺𝑛,𝑗,ℎ], where  𝐺𝑛,𝑗,ℎ =𝛹(𝐾𝑛,𝑗,ℎ) − diag(𝛹′(𝑈𝑛 ℎ))𝐾𝑛,𝑗,ℎ. On the other hand, the term in 𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0), considering the second part of (14) is 𝑘𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,2 𝐺𝑛,𝑗,ℎ +𝑐2 𝑖𝑘𝑃ℎ ℎ(𝑡𝑛)]. Summing up, 𝐾𝑛,𝑖,ℎ =𝑒𝑐𝑖𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[𝑐𝑖[𝑃ℎℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛)] + ∑𝜆𝑖,𝑗,1 𝐺𝑛,𝑗,ℎ] +𝑘𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,2 𝐺𝑛,𝑗,ℎ +𝑐2 𝑖𝑘𝑃ℎ ℎ(𝑡𝑛)] +𝑘 𝑟 ∑ 𝑙=3 𝜑𝑙(𝑐𝑖𝑘𝐽𝑛,ℎ,0)∑𝜆𝑖,𝑗,𝑙  𝐺𝑛,𝑗,ℎ. As for 𝑈𝑛+1 ℎin (9), multiplying 𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)𝐶ℎ, just 𝜕𝑢(𝑡𝑛) = 𝑔(𝑡𝑛)turns up again; on the other hand, the term in 𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)can be written as −𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕[𝑢(𝑡𝑛) − 𝛹(𝑢(𝑡𝑛)) + 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) − ℎ(𝑡𝑛)].(16) This term can be exactly calculated when considering Dirichlet boundary conditions and, with Robin/Neumann boundary conditions, it can be approximated through the numerical approximation at the boundary given by the space discretization of (1) itself, but without resorting to numerical differentiation. On the other hand, notice that, by using the first part of (13), some terms cancel and the term multiplying 𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐶ℎis just 𝜕 𝑢(𝑡𝑛). Therefore, that boundary can be calculated exactly in terms of data as 𝑔(𝑡𝑛). As for the terms in 𝑘2𝜑𝑙+1(𝑘𝐽𝑛,ℎ,0)𝐶ℎwith 𝑙≥2, notice that they vanish in (9) because of the second part of (13). Moreover, considering (15) and the first part of (13), the terms in 𝜑1(𝑘𝐽𝑛,ℎ,0)without boundaries can be simplified to 𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[𝑃ℎℎ(𝑡𝑛) + ∑𝜇𝑖,1 𝐺𝑛,𝑖,ℎ]. In a similar way, but using now also the second part of (13), the terms in 𝜑2(𝑘𝐽𝑛,ℎ,0)and 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)without boundaries, with 𝑙≥3, can be simplified to 𝑘𝜑2(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,2 𝐺𝑛,𝑖,ℎ +𝑘𝑃ℎ ℎ(𝑡𝑛)] 𝑘𝜑𝑙(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,𝑙  𝐺𝑛,𝑖,ℎ], 𝑙 ≥3. Summing up, under assumptions (13),𝑈𝑛+1 ℎin (9) can be simplified to 𝑈𝑛+1 ℎ=𝑒𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[𝑃ℎℎ(𝑡𝑛) + ∑𝜇𝑖,1 𝐺𝑛,𝑖,ℎ +𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) − 𝛹(𝑢(𝑡𝑛)) + 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) − ℎ(𝑡𝑛)]] Journal of Computational and Applied Mathematics 453 (2025) 116158 6 B. Cano and M.J. Moreta +𝑘𝜑2(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,2 𝐺𝑛,𝑖,ℎ +𝑘[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕 𝑢(𝑡𝑛)]]+𝑘 𝑟 ∑ 𝑙=3 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)∑𝜇𝑖,𝑙  𝐺𝑛,𝑖,ℎ. 3.2. 𝑝= 2 With similar arguments as those for 𝑝= 1, the stages for 𝑝= 2 in (11) can be simplified to the following formulas: 𝐾𝑛,𝑖,ℎ =𝑒𝑐𝑖𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ +𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,1 𝐺𝑛,𝑗,ℎ +𝑐𝑖[𝑃ℎℎ(𝑡𝑛) + 𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) − 𝛹(𝑢(𝑡𝑛)) + 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) − ℎ(𝑡𝑛)]]] +𝑘𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,2 𝐺𝑛,𝑗,ℎ +𝑐2 𝑖𝑘[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕 𝑢(𝑡𝑛)]]+𝑘 𝑟 ∑ 𝑙=3 𝜑𝑙(𝑐𝑖𝑘𝐽𝑛,ℎ,0)∑𝜆𝑖,𝑗,𝑙  𝐺𝑛,𝑗,ℎ. As for 𝑈𝑛+1 ℎin (11), now the term multiplying 𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐶ℎis again 𝜕 𝑢(𝑡𝑛). As for 𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐷ℎ, we have −𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) +  𝐽(𝑡𝑛)[𝛹(𝑢(𝑡𝑛)) + ℎ(𝑡𝑛) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛)]] = −𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)=−𝜕[𝑢(𝑡𝑛) −  ℎ(𝑡𝑛)], where, for the last equality, we have just considered the differentiation of Eq. (1). On the other hand, the term multiplying in 𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎis 𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) +  𝐽(𝑡𝑛)[𝛹(𝑢(𝑡𝑛)) + ℎ(𝑡𝑛) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛)] +  ℎ(𝑡𝑛)] =𝜕[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) +  ℎ(𝑡𝑛)] = 𝜕 𝑢(𝑡𝑛), and the term in 𝜑3(𝑘𝐽𝑛,ℎ,0)𝐷ℎvanishes because ∑𝜇𝑖,2= 0. In a similar way, the terms in 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐶ℎand 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐷ℎfor 𝑙≥4 vanish because of the second part of (13). Considering this, 𝑈𝑛+1 ℎin (11) can be calculated as 𝑈𝑛+1 ℎ=𝑒𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,1  𝐺𝑛,𝑖,ℎ +𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) − 𝛹(𝑢(𝑡𝑛)) + 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) − ℎ(𝑡𝑛)] ] +𝑘𝜑2(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,2  𝐺𝑛,𝑖,ℎ +𝑘[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕 𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) −  ℎ(𝑡𝑛)]]] +𝑘𝜑3(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,3  𝐺𝑛,𝑖,ℎ +𝑘2𝐶ℎ𝜕 𝑢(𝑡𝑛)] +𝑘 𝑟 ∑ 𝑙=4 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)∑𝜇𝑖,𝑙   𝐺𝑛,𝑖,ℎ, where   𝐺𝑛,𝑗,ℎ =𝛹(𝐾𝑛,𝑗,ℎ) − diag(𝛹′(𝑈𝑛 ℎ))𝐾𝑛,𝑗,ℎ +𝑃ℎℎ(𝑡𝑛,𝑗 ) − 𝑐𝑗𝑘𝑃ℎ ℎ(𝑡𝑛). We notice that again a term like (16) turns up, which can either be calculated exactly or approximated without resorting to numerical differentiation. As for the other terms on the boundary, they can be calculated in terms of data with both Dirichlet and Robin/Neumann boundary conditions since 𝜕 𝑢(𝑡) = 𝑔(𝑡)and 𝜕 𝑢(𝑡) = 𝑔(𝑡). 3.3. 𝑝= 3 Similarly to the calculation of 𝑈𝑛+1 ℎwith 𝑝= 2, but using (14), the stages in (12) can be simplified to 𝐾𝑛,𝑖,ℎ =𝑒𝑐𝑖𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,1  𝐺𝑛,𝑗,ℎ +𝑐𝑖[𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) − 𝛹(𝑢(𝑡𝑛)) + 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) − ℎ(𝑡𝑛)]]] +𝑘𝜑2(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,2  𝐺𝑛,𝑗,ℎ +𝑐2 𝑖𝑘[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕 𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) −  ℎ(𝑡𝑛)]]] +𝑘𝜑3(𝑐𝑖𝑘𝐽𝑛,ℎ,0)[∑𝜆𝑖,𝑗,3  𝐺𝑛,𝑗,ℎ +𝑐3 𝑖𝑘2𝐶ℎ𝜕 𝑢(𝑡𝑛)] +𝑘 𝑟 ∑ 𝑙=4 𝜑𝑙(𝑐𝑖𝑘𝐽𝑛,ℎ,0)∑𝜆𝑖,𝑗,𝑙   𝐺𝑛,𝑗,ℎ.(17) As for the terms concerning boundaries to calculate 𝑈𝑛+1 ℎin (12), using the left part of (13) and (14), the term in 𝜑2(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕can be simplified to 𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[𝑢(𝑡𝑛) + 𝑘2 2(∑𝜇𝑖,1𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)]], which can again be calculated exactly in terms of data with Dirichlet boundary conditions and can be approximated with Robin/Neumann ones considering the error from the approximation itself at the boundary and resorting to numerical differentiation just for the first time derivative 𝑢. On the other hand, the term in 𝜑2(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕can be simplified to −𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) +  𝐽(𝑡𝑛)[𝛹(𝑢(𝑡𝑛)) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) + ℎ(𝑡𝑛)]] Journal of Computational and Applied Mathematics 453 (2025) 116158 7 B. Cano and M.J. Moreta = −𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛)=−𝑘2𝜑2(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕[𝑢(𝑡𝑛) −  ℎ(𝑡𝑛)], which is exactly calculable in terms of data for both Dirichlet and Robin/Neumann boundary conditions. As for the terms in 𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕, they can be written as 𝑘2𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[𝑘2 2(∑𝜇𝑖,2𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)] +𝑘[ 𝐽(𝑡𝑛)2𝑢(𝑡𝑛) +  𝐽(𝑡𝑛)[𝛹(𝑢(𝑡𝑛)) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) + ℎ(𝑡𝑛)] +  ℎ(𝑡𝑛)]] =𝑘2𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[𝑘2 2(∑𝜇𝑖,2𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)] + 𝑘[ 𝐽(𝑡𝑛)𝑢(𝑡𝑛) +  ℎ(𝑡𝑛)]] =𝑘2𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[𝑘2 2(∑𝜇𝑖,2𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)] + 𝑘𝑢(𝑡𝑛)]. Similarly, the term in 𝜑3(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕can be written as −𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕[ 𝐽(𝑡𝑛)3𝑢(𝑡𝑛) +  𝐽(𝑡𝑛)2[𝛹(𝑢(𝑡𝑛)) + ℎ(𝑡𝑛) − 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛)] +  𝐽(𝑡𝑛) ℎ(𝑡𝑛)] = −𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕[ 𝐽(𝑡𝑛)[  𝐽(𝑡𝑛)𝑢(𝑡𝑛) +  ℎ(𝑡𝑛)]] = −𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕 𝐽(𝑡𝑛)𝑢(𝑡𝑛) = −𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕[… 𝑢(𝑡𝑛) −  ℎ(𝑡𝑛) − 𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2], where the last equality comes from differentiating (1) three times. In a similar way, it can be deduced that the term in 𝜑4(𝑘𝐽𝑛,ℎ,0)𝐶ℎis 𝑘4𝜑4(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕[1 2(∑𝜇𝑖,3𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)]+ … 𝑢(𝑡𝑛) −  ℎ(𝑡𝑛) − 𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2] , but that in 𝜑4(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕vanishes, as well as the possible terms 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕and 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)𝐷ℎ𝜕for 𝑙≥5. Summing up, under the assumptions in (14),𝑈𝑛+1 ℎin (12) can be written as 𝑈𝑛+1 ℎ=𝑒𝑘𝐽𝑛,ℎ,0𝑈𝑛 ℎ+𝑘𝜑1(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,1  𝐺𝑛,𝑖,ℎ +𝐶ℎ𝜕𝑢(𝑡𝑛) − 𝐷ℎ𝜕[𝑢(𝑡𝑛) − 𝛹(𝑢(𝑡𝑛)) + 𝛹′(𝑢(𝑡𝑛))𝑢(𝑡𝑛) − ℎ(𝑡𝑛)]] +𝑘𝜑2(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,2  𝐺𝑛,𝑖,ℎ +𝑘[𝑃ℎ ℎ(𝑡𝑛) + 𝐶ℎ𝜕[𝑢(𝑡𝑛) + 𝑘2 2(∑𝜇𝑖,1𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)]]−𝐷ℎ𝜕[𝑢(𝑡𝑛) −  ℎ(𝑡𝑛)]]] +𝑘𝜑3(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,3  𝐺𝑛,𝑖,ℎ +𝑘2[𝐶ℎ𝜕[𝑘 2(∑𝜇𝑖,2𝑐2 𝑖)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)] + 𝑢(𝑡𝑛)] −𝐷ℎ𝜕[… 𝑢(𝑡𝑛) −  ℎ(𝑡𝑛) − 𝛹′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2]]] +𝑘𝜑4(𝑘𝐽𝑛,ℎ,0)[∑𝜇𝑖,4  𝐺𝑛,𝑖,ℎ +𝑘3𝐶ℎ𝜕[… 𝑢(𝑡𝑛)+(1 2∑𝜇𝑖,3𝑐2 𝑖− 1)[𝛹′′(𝑢(𝑡𝑛)) 𝑢(𝑡𝑛)2+ ℎ(𝑡𝑛)]]] +𝑘 𝑟 ∑ 𝑙=5 𝜑𝑙(𝑘𝐽𝑛,ℎ,0)∑𝜇𝑖,𝑙   𝐺𝑛,𝑖,ℎ (18) Again all the terms at the boundary in this expression can be exactly calculated in terms of data with Dirichlet boundary conditions and in an approximated way with Robin/Neumann boundary ones taking into account the approximated values of the space discretization of (1) at the boundary and the approximation of 𝑢 through numerical differentiation in time. 3.3.1. Concluding remarks Remark 3.1. We notice that, in any case, no numerical differentiation in space is required to approximate the required boundary values, so that no weak CFL condition is required to prove the classical order of the method, as it was in principle necessary in the more general case [17]. Remark 3.2. Although, for the sake of brevity, we do not show the calculations here, for 𝑝= 4, numerical differentiation in space would be required to approximate boundary values with both Dirichlet and Robin/Neumann boundary conditions, even under assumptions (13)–(14). 4. Conditions under which the simplifying assumptions are satisfied In this section, we will see why assumptions (13) and (14) are nearly always satisfied for methods of classical order ≤4. For that, we will consider the classical order conditions on the coefficients of the Rosenbrock method (2) and, from them, we will justify when (13) and (14) are guaranteed. Classical order conditions till order four were derived in [13], although just assuming that 𝑠= 2 and that ∑𝑏𝑖(𝑧) = 𝜑1(𝑧). In the general case, just considering Taylor expansions of (5)–(6) on the timestepsize, the fact that 𝐺′ 𝑛(𝑈𝑛)=0 because of (7), and comparing with the Taylor expansion of the exact solution of (4), the order conditions of Table 1 turn up. Then, we have the following result: Journal of Computational and Applied Mathematics 453 (2025) 116158 8 B. Cano and M.J. Moreta Table 1 Classical order conditions for exponential Rosenbrock methods. Order Conditions 1∑𝑏𝑖(0) = 1 2∑𝑏′ 𝑖(0) = 1 2 3∑𝑏′′ 𝑖(0) = 1 3∑𝑏𝑖(0)𝑐2 𝑖=1 3 4∑𝑏′′′ 𝑖(0) = 1 4∑𝑏′ 𝑖(0)𝑐2 𝑖=1 12 ∑𝑏𝑖(0)𝑐𝑖𝑎′ 𝑖𝑗 (0) = 1 8∑𝑏𝑖(0)𝑐3 𝑖=1 4 Theorem 4.1. If 𝑞denotes the classical order of an exponential Rosenbrock method (1≤𝑞≤4) and 𝑟in (2) satisfies 𝑟≤𝑞, then (13) is satisfied. Proof. We firstly notice that, considering (2), the first column of conditions in Table 1 is equivalent to 𝑟 ∑ 𝑙=1 1 𝑙! 𝑠 ∑ 𝑖=1 𝜇𝑖,𝑙 = 1, 𝑟 ∑ 𝑙=1 1 (𝑙+ 1)! 𝑠 ∑ 𝑖=1 𝜇𝑖,𝑙 =1 2, 𝑟 ∑ 𝑙=1 1 (𝑙+ 2)! 𝑠 ∑ 𝑖=1 𝜇𝑖,𝑙 =1 6, 𝑟 ∑ 𝑙=1 1 (𝑙+ 3)! 𝑠 ∑ 𝑖=1 𝜇𝑖,𝑙 =1 24 .(19) Therefore, when 𝑞= 1, just the first equation must hold and, when 𝑟= 1,(13) follows directly. When 𝑞= 2, the first two equations must hold. When 𝑟= 1, both equations are the same and again (13) follows immediately. When 𝑟= 2, we have a linear system of two equations in the two unknowns ∑𝜇𝑖,1and ∑𝜇𝑖,2. The matrix associated to that system is clearly nonsingular. Because of that, there is a unique solution of that system, which obviously corresponds to ∑𝜇𝑖,1= 1 and ∑𝜇𝑖,2= 0. When 𝑞= 3, the first three equations must hold. When 𝑟≤2, the first two of them lead to (13) with the same arguments than before and the corresponding solution happen to also satisfy the third equation. When 𝑟= 3, we have a linear system of three equations and three unknowns, which matrix is again non-singular and the unique solution of the system is ∑𝜇𝑖,1= 1,∑𝜇𝑖,2= 0 and ∑𝜇𝑖,3= 0. For 𝑞= 4, a similar argument leads to the result. □ Let us see now under which conditions (14) is guaranteed. Theorem 4.2. (i) If 𝑟= 1 or 𝜆𝑖,𝑗,𝑙 = 0 for 𝑙≥2,(14) is always satisfied. (ii) If 𝑟= 2 or 𝜆𝑖,𝑗,𝑙 = 0 for 𝑙≥3,𝑠= 2 and 𝑞= 4,(14) is satisfied. Proof. (i) comes directly from (3), which can be written like this considering (2) 𝑟 ∑ 𝑙=1 1 𝑙! 𝑖−1 ∑ 𝑗=1 𝜆𝑖,𝑗,𝑙 =𝑐𝑖, 𝑖 = 1,…, 𝑠. (20) In order to prove (ii), we notice that, for 𝑞= 4 and 𝑠= 2, the last two conditions in the last row of Table 1 read 𝑏2(0)𝑐2𝑎′ 21(0) = 1 8, 𝑏2(0)𝑐3 2=1 4. Considering (2) again, the latter can also be written as 𝑐2 2( 𝑟 ∑ 𝑙=1 𝜇2,𝑙 1 𝑙!)( 𝑟 ∑ 𝑙=1 𝜆2,1,𝑙 1 (𝑙+ 1)! )=1 8, 𝑐3 2 𝑟 ∑ 𝑙=1 𝜇2,𝑙 1 𝑙!=1 4, which imply that 𝑟 ∑ 𝑙=1 𝜆2,1,𝑙 1 (𝑙+ 1)! =𝑐2 2. Taking here 𝑟= 2 as well as in (20) with 𝑖= 2, a uniquely solvable linear systems of two equations and two unknowns turn up, which lead to 𝜆2,1,1=𝑐2and 𝜆2,1,2= 0.□ Journal of Computational and Applied Mathematics 453 (2025) 116158 9 B. Cano and M.J. Moreta Remark 4.3. We notice that the conditions which guarantee that the simplifying assumptions are satisfied mainly concern the maximum value 𝑟for the index 𝑙in the 𝜑𝑙-functions, which must be small enough with respect to the classical order which wants to be achieved. This is not a serious drawback since, for the sake of simplicity, in the construction of methods, 𝑟is taken as small as possible. Only in case (ii) of Theorem 4.2 there is also a restriction for order 𝑞= 4 on the number of stages if 𝜑2is turning up in the coefficients 𝑎𝑖𝑗 . However, this is not restrictive either since, as stated in the introduction, with Rosenbrock methods, very few stages are required to get a desired accuracy. In particular, classical order 4can be obtained with just 2stages. Remark 4.4. We also remark that the simplifying assumptions (13)–(14) are equivalent to the simplifying assumptions in [12] 𝑠 ∑ 𝑖=1 𝑏𝑖(𝑧) = 𝜑1(𝑧), 𝑖−1 ∑ 𝑗=1 𝑎𝑖𝑗 (𝑧) = 𝑐𝑖𝜑1(𝑐𝑖𝑧),1≤𝑖≤𝑠, under which it was assured that equilibria of autonomous problems were preserved and which allowed to simplify stiff order conditions in that paper. More particularly, stiff order 2was assured under those assumptions when integrating that type of problems when considering vanishing boundary conditions. Because of that, assumptions (13)–(14) are satisfied by mainly all already constructed methods. As distinct, in this paper, we justify through Theorems 4.1 and 4.2 that those assumptions are assured to be satisfied in many cases, without the need to resort to stiff order conditions of any kind or the preservation of equilibria. 5. Recommended methods depending on the desired accuracy In this section, considering the results on the previous ones, we will suggest what we think is the best choice up to the moment of exponential Rosenbrock methods to achieve the particular orders of accuracy 𝑞= 2,3and 4. 5.1. 𝑞= 2 Rosenbrock-Euler method, which just has one stage and corresponds to 𝑏1(𝑧) = 𝜑1(𝑧), is well-known to have classical order 2. By looking at Table 1, we can see that not only the conditions for classical order 2are satisfied, but also one of the conditions for classical order 3. Moreover, the other condition to achieve the latter accuracy cannot be satisfied with any method which just has one stage since, in such a case, because of (3),𝑐1= 0. Therefore, this seems to be an unbeatable second-order method. What is more, it happens to satisfy (13) and (14). (In fact, that could also be deduced directly from Theorems 4.1 and 4.2.) Because of that, the technique in [17] with 𝑝= 2 can be applied to achieve local order 3without resorting to numerical differentiation to calculate the required boundary values. (We notice that Euler-Rosenbrock method has stiff order 2according to [12], but just shows local order 2when implemented through the standard method of lines when the boundary condition 𝑔(𝑡)is time-dependent). Numerical results in [17] show the big advantage in computational time of using the modified Rosenbrock-Euler method against Rosenbrock-Euler with the standard method of lines. The simplifications for that particular simple modified Rosenbrock-Euler method were already done in [17], where it was observed that the difference with the standard method of lines just consisted on adding a term of the form 𝑘3𝜑3(𝑘𝐽𝑛,ℎ,0)𝐶ℎ𝜕 𝑢(𝑡𝑛), when calculating 𝑈𝑛+1 ℎfrom 𝑈𝑛 ℎ. 5.2. 𝑞= 3 As stated before, it is impossible to get an exponential Rosenbrock method of classical order 3with just one stage. Because of that, we look for one with two stages. Trying to be as efficient as possible, we take 𝑐2= 1 so that possible evaluations at 𝑡=𝑡𝑛+𝑐2𝑘 can also be used at the next step. The simpler function 𝑎21(𝑧)of the form (2) satisfying (3) is then 𝑎2,1(𝑧) = 𝜑1(𝑧). If we now look for functions 𝑏1(𝑧)and 𝑏2(𝑧)satisfying the four necessary conditions in Table 1, we can see that we can achieve that just by considering 𝑟= 1. Although a linear system of four equations with two unknowns is obtained, three of them are equivalent and altogether lead to 𝑏1(𝑧) = 2 3𝜑1(𝑧), 𝑏2(𝑧) = 1 3𝜑1(𝑧). This method again satisfies conditions (13) and (14), as it was also assured through Theorems 4.1 and 4.2. By considering 𝑟= 2, a one-parameter family of methods turn up, which correspond to 𝜇1,1=2 3+𝜇2,2 2, 𝜇1,2= −𝜇2,2, 𝜇2,1=1 3−𝜇2,2 2.(21) We notice that with all these methods, the first equation in the last row of Table 1 is satisfied. As for the third and fourth equation in the same row, they are never satisfied. However, the second equation in that row is just satisfied for 𝜇2,2= 1, and that leads to 𝑏1(𝑧) = 7 6𝜑1(𝑧) − 𝜑2(𝑧), 𝑏2(𝑧)=−1 6𝜑1(𝑧) + 𝜑2(𝑧). We may therefore expect that this leads to the smallest local errors inside the family (21). All previous methods have stiff order 2but not stiff order 3according to [12] since 𝑏1(𝑧) + 𝑏2(𝑧) = 𝜑1(𝑧), 𝑎21(𝑧) = 𝑐2𝜑1(𝑐2𝑧), 𝑏2(𝑧)𝑐2 2≠2𝜑3(𝑧). However, the technique in [17] can be applied to avoid order reduction. Again the simplifying assumptions (13) and (14) are satisfied and therefore, the advantages of the simplified formulas in Section 3can be used.