Derivation of third order Runge–Kutta methods (ELDIRK) by embedding of lower order implicit time integration schemes for local and global error estimation

Three prominent low order implicit time integration schemes are the first order implicit Euler-method, the second order trapezoidal rule and the second order Ellsiepen method. Its advantages are stability and comparatively low computational cost, however, they require the solution of a nonlinear system of equations. This paper presents a general approach for the construction of third order Runge–Kutta methods by embedding the above mentioned implicit schemes into the class of ELDIRK-methods. These will be defined to have an Explicit Last stage in the general Butcher array of Diagonal Implicit Runge–Kutta (DIRK) methods, with the consequence, that no additional system of equations must be solved. The main results—valid also for non-linear ordinary differential equations—are as follows: Two extra function calculations are required in order to embed the implicit Euler-method and one extra function calculation is required for the trapezoidal-rule and the Ellsiepen method, in order to obtain the third order properties, respectively. Two numerical examples are concerned with a parachute with viscous damping and a two-dimensional laser beam simulation. Here, we verify the higher order convergence behaviours of the proposed new ELDIRK-methods, and its successful performances for asymptotically exact global error estimation of so-called reversed embedded RK-method are shown.


Introduction
The mathematical modeling of several engineering problems requires the numerical solution of ordinary differential equations.Typical examples for a first order ordinary differential equation are the resulting heat equation discretized in space or inelastic constitutive equations, whilst structural dynamic problems lead to examples of second order ordinary differential equations, see e.g.[16,33].Choosing the appropriate integration method requires deciding on aspects such as the order of the integration scheme, stability properties and computational efficiency.A general class of algorithms is provided by so-called S-stage Runge-Kutta methods (RK-methods), where explicit and implicit schemes are distinguished, [5,10,19].Due to the properties of A-stability, where the region of absolute stability is unbounded, see e.g.[19], and relatively low computational cost, the implicit Euler-method and the generalized trapezoidal-rule are most prominent integration schemes in structural mechanics, see e.g.[16,33].Also the second order Ellsiepen method [7], with the additional property of L-stability, proved to be very efficient for the laser beam simulation in additive manufacturing processes, see e.g.[28,29].
A main focus of the paper is the derivation of reliable a posteriori local error estimates, which are essential for the two aims, 1. to bound the local error by adaptive timestep control, 2. to contribute to an estimate of the global error.
Regarding the above mentioned first aim, to bound the local error, two concepts are distinguished for the calculation of a posteriori local error estimates: firstly, by means of higher order methods and secondly by means of extrapolation methods, see e.g.[19] on general information and [23] on applications for creep analysis with the finite element method.In line with the first concept, Fehlberg suggested a strategy for stepsize selection by so-called embedded Runge-Kutta methods, [11][12][13].The novelty of Fehlberg's method is that the applied integration scheme is embedded by a second method with higher order, however, identical Runge-Kutta coefficients.As a consequence only one extra function calculation is required to estimate the local error of the embedded method.Thereupon, several embedded embedded Runge-Kutta schemes have been developed.For example Bogacki and Shampine [3] present a Runge-Kutta method of order three with four stages in order to embed a second-order method.A recent paper [14] investigates strong stability properties of embedded pairs of RK-methods, however, restricted to explicit methods.
A slightly different embedding compared to Fehlbergs explicit schemes has been given by Ellsiepen [7] with Butcher array ( This array defines an embedded RK-method, subsequently denoted as RK(2)1-Ell, combining a method RK2-Ell (using α below the horizontal line) of order 2 with embedded method RK1-Ell of order 1 (using β) [7].Both methods have the same stage and, different to Fehlbergs explicit schemes, both methods are implicit.Exploiting both choices below the horizontal line has the consequence, that even for nonlinear problems no additional functional evaluation becomes necessary for local error estimation.
This paper intends to combine the advantage of Fehlsbergs embedding approach with only one additional functional evaluation to obtain a local error estimate, with the advantages of A-stability and relatively low computational cost of low order implicit methods.In order to avoid an additional solution of a system of equations we introduce a certain class of RK-methods labeled as ELDIRK.These methods are defined to have an explicit last stage in the general Butcher array of Diagonal Implicit Runge-Kutta (DIRK) methods, contrary to the class of EDIRK, where according to Kvaernø et al. [18] the RK-methods have an explicit first stage in the general Butcher array In addition to Mahnken [22], a general approach for the construction of third order ELDIRK-methods is presented by embedding some well known low order implicit time integration schemes.The proposed methodology is applied to embed the first order implicit Euler-method, the second order trapezoidal rule and the second order Ellsiepen method by a third order scheme.As main results, two extra function calculations are required in order to embed the implicit Eulermethod and one extra function calculation is required for both, the trapezoidal-rule and Ellsiepen's method, in order to obtain the third order properties.In this way, we will extend the Butcher array in Eq. ( 1) to obtain a new method RK3-Ell of order 3 in order to embed RK2-Ell of order 2 with only one additional function evaluation.As verified in the numerical examples in Sect.7, the higher order properties for all new ELDIRK-methods with only one (and respectively two) additional function evaluation(s) are obtained for both linear as well as for nonlinear ordinary differential equations.
Two different types of embedded RK-methods are distinguished in the literature.The original embedded RK-method explained above will subsequently be termed standard embedded RK-method.For its variation, the roles of embedded and improved solutions are simply reversed, see [19, p. 108], and therefore this method will subsequently be denoted as reversed embedded RK-method.
Regarding the above mentioned second aim of local error estimates, to contribute to an estimate of the global error, we follow closely the procedure recently presented in [22].Here, we consider an estimate for the global error measure obtained from a quantity of interest, for both, the standard embedded and the reversed embedded RK-method.This strategy avoids the solution of a dual problem as performed e.g. in [25].The dual problem runs backwards in time [2,17,21,31], which due to additional storage requirements might become a serious drawback when applied e.g. to FEM or respectively to FDM simulations.Moreover, in [22], we perform a linearisation of the quantity of interest, related to both embedded RK-methods, in order to investigate the corresponding effectivity indices with the following main results: The error estimate for the standard RK-method approaches I eff = κ in the limit.Here κ is a constant, which can be above or below one, meaning, the error can be over-or underestimated.The effectivity index for the reversed RK-method approaches I eff = 1 in the limit, meaning that this error estimate is asymptotically exact.In this work, these results are verified numerically, for all new RK-schemes developed in this paper.
The structure of the paper is as follows: Sect. 2 summarises some preliminaries on numerical time integration of systems of ordinary equations with RK-methods including a posteriori error estimation.Here also, the class of ELDIRK-methods is introduced.Moreover, we will recall the order conditions from [4,6] to guarantee a certain convergence order of an RK-method.In Sects. 3 and 4 we present a general approach for the construction of second order embedding ELDIRK-methods and third order embedding ELDIRK-methods, respectively.The methodology is applied to three prominent low order implicit integration schemes, that is the implicit Euler-method, the generalized trapezoidal-rule and Ellsiepen's method.In this way, for each of the methods we obtain an extended Butcher tableau.Section 5 provides a brief revision of goal oriented error measures introduced in [22], and the asymptotic behaviour of error estimates for the standard embedded and the reversed embedded RK-methods.Two numerical examples in Sect.6 are concerned with a parachute with viscous damping and a two-dimensional laser beam simulation.Here, we illustrate the higher order convergence behaviour for equidistant step sizes of the proposed new Runge-Kutta methods.Finally, a summary of the main results of this paper and an outlook for future work are given in Sect.7.

Generalized Runge-Kutta scheme
Consider the initial value problem (IVP) where the solution y(t) ∈ R n , n > 0 typically represents a physical process dependent on time t.Moreover, A is a nonlinear operator defined in D(A) ⊂ V , where V is a suitable Hilbert space with corresponding norm induced by an inner product (omitting mathematical details for brevity of presentation).The source (or respectively loading) term f and the initial value y α are given data.We note, that the representation (2) can be generalized to parabolic equations, see e.g.[31].An example for the heat equation taking into account an additional spatial dependence will be given in Subsection 6.2 of this work.For brevity of presentation, in the remainder of this section we restrict ourselves to temporal discretization.
We begin with a partition of the time interval and y q i , i = 1, . . ., N denotes the numerical approximations of a q-order RK-method for the initial value problem (2) at times t i .Then, starting with y 0 = y α at t 0 = α the time stepping algorithm of the RK-method is generally written in terms of a slope function (or increment function, see [10]), of the method φ q as y q i = y For an S-stage RK-method the algorithm (4) specifies as 1. y 123 where the weights b j , the knots c j and the RK-coefficients a jl , respectively, satisfy the restrictions 1.

S l=1
a jl =c j , j=1, . . ., S. ( Then, the slope function in Eq. ( 4) becomes a weighted sum of the slopes in Eq. (5.2) as φ q (t i , y The coefficients a jl , b j , c j , with j, l ∈ {1, . . ., S} are the ones corresponding to the Butcher tableau [5] The first row of coefficients b j below the horizontal line in (7) will be explained later.For a jl = 0, l ≥ j the slopes k j in Eq. (5.2) are obtained explicitly such that the number of stages S is equal to the number of function evaluations per time-step.Otherwise the method becomes implicit, where for non-linear ordinary differential equations generally the number of function evaluations per time-step exceeds the number of stages S due to some iterative procedure.An iterative Newton scheme for determination of the slopes k j valid for implicit RK-methods is described e.g. in [22].

The class of ELDIRK-methods
A classification of RK-methods is introduced in [18].Here, in particular Diagonal Implicit Runge-Kutta (DIRK)-type methods are of particular interest, where in the Butcher tableau in Eq. ( 7) one has a i j = 0 for all j > i and a ii = 0 for at least one i.We denote by Q and R > Q the stages for the pair of a q-order method embedded by an r −order RK-method, where r > q.Such an embedding of implicit DIRK methods has been considered for some prominent low order methods in [22].There, it is shown, that an additional solution of a system of equations for the R-stage method can be avoided, if the diagonal elements at the last stage(s) in the general Butcher array in Eq. ( 7) are zero.For a systematic generalization of this convenient property we introduce in this paper a class of RK-methods defined as ELDIRK: Explicit Last stage DIRK, where a j j = 0, Q < j ≤ R. (8) This is in contrast to the class of EDIRK methods, which according to Kvaernø et al. [18] have an explicit first stage in the general Butcher array in Eq. (7), that is (at least) a 11 = 0.The practical implementation of ELDIRK can be split into two steps: Firstly, the unknowns k j , j = 1, . . ., Q < R of the related implicit Q-stage method are determined by a nonlinear Newton algorithm as explained e.g. in [22].Then, in a second (post-processing) step, the additional slopes are obtained in an explicit manner by exploiting definition (8) from Eq. (5.2) as In this way, there is no need to solve another nonlinear system at all, which holds for linear as well as for nonlinear ordinary differential equations of the form (2.1).

Order and stage conditions
Following standard text books on odes, see e.g.[19], a local error is introduced, As shown in several classical textbooks, see e.g.[19], the order of accuracy can be obtained analytically, by inserting the slope φ q (t i−1 , y(t i−1 ), y(t i ), τ i ) = Q j=1 b j k j for the Runge-Kutta method under investigation into the definition for the local error in Eq. (10) followed by a Taylor-series expansion.
For more complex RK-methods this procedure may ensue a tedious analysis, which can completely be avoided by means of the general results on RK-methods in [6].Here it is extensively shown, that if the RK-method is of order q, then the RK-coefficients c k , b k , a kl of the Butcher array in Eq. ( 5) have to satisfy certain order conditions for the RKq formula.Butcher [4] has listed the related conditions up to order 8. Table 1 contains the complete set of required equations for orders up to q = 4. Note, that the particular case q = 1 accounts for consistency, see e.g.[19].Additionally, in this table we introduce the number of elementary differentials n r thus resulting into the total number of order conditions m q = q r =1 n r , see [6].Also, the stage conditions up to stage Q are recalled from Eq. (6.3) in Table 1.
Example.Application of the order and stage conditions of Table 1 to the RK(2)1-Ell method of Ellsiepen [7] with Butcher tableau in Eq. (1) reveals: Table 1 Order conditions (up to q = 4) and stage conditions • Numbers of elementary differentials n r and order conditions m q Orders r , q 1 2 3 4 n r 1 1 2 4 n r 1 2 4 8 • Order conditions for orders q = 1, . . ., 4 which confirms, that RK1-Ell (using β) is indeed only a first order method whilst RK2-Ell (using α) is a second order method.Both RK-schemes satisfy the consistency condition.

Embedded Runge-Kutta methods
An efficient computation of error estimates is achieved by so-called embedded Runge-Kutta methods, where two conceptions are distinguished: The standard embedded RK-method in its original form introduced in [11][12][13], and subsequently denoted as RK(r)q, is given as a pair of time stepping algorithms according to (4) as, see e.g.[19], Here RKq has order q and stage Q.The higher order method RKr,q, r ≥ q + 1 and stage R > Q, uses the initial value y r ,q i−1 = y q i−1 at each time step, which explains the upper index r , q.Alternatively, a reversed embedded RK-method RK(q)r see e.g.[6], also known as local extrapolation, see e.g.[19], is given as a pair of time stepping algorithms according to (4), Note, that the lower order method RKq,r, q ≤ r −1 and stage Q uses the initial value y q,r i−1 = y r i−1 at each time step, which explains the upper index q, r .
According to Fehlberg's suggestion the RK-formulae of orders q and r > q both in (12) and in (13) (usually r = q + 1) share the same function evaluations, i.e. they have the same RK-coefficients a jl and the same knots c j , j = 1, ..., Q.In addition to b j , j = 1, ..., Q weighting factors b j , j = 1, ..., R are introduced in the Butcher array in Eq. ( 7), where the symbol ˆrefers to the higher r -order method with numerical solutions y r i according to Eq. ( 5).A reliable procedure for the numerical solution of initial value problems requires adaptive time-step selection.Let us assume the following conditions [19]: 1.The value y q i−1 at the beginning of the time interval is exact, i.e. y ), where r ≥ q + 1, 3. Terms of order O(τ q+2 ) are neglected.(14) Then, from the above embedding 1. ηr, become a posteriori estimates for the local error of e i for both RK(r)q and RK(q)r.These are used to control the step-size for both, RK(r)q in (12) and RK(q)r in (13), respectively, according to the formulae where a, b = r , q or a, b = q, r , respectively, and where • denotes some norm of the related quantity, e.g. the l2norm or the maximum-norm.Also, ε is a tolerance for the estimated error, 0.9 is a safety factor and 2 and 0.2 are values for upper and lower bounds of the step-length.Since the error estimate ηa,b i is dependent on the (unknown) final time-step τ i , in practice τ i is obtained in a trial and error procedure.We point out, that the procedure for local error estimation by 123 means of the extended Butcher tableau in Eq. ( 7) is applicable to both, explicit and implicit embedded RK-methods.

General setting
This section intends to develop a general procedure to embed a Runge-Kutta method of order 1, subsequently denoted as RK1 with Butcher array B1, by an ELDIRK-method of order 2, subsequently denoted as RK2, thus resulting into an embedded RK-method RK(2)1 with extended Butcher array B(2)1: In these tableaus RK-coefficients without a hat are assumed to be known, and RK-coefficients with a hat are assumed to be unknown.Note, that in order to become an ELDIRK-method according to definition (8), we require â22 = 0. We also assume, that the stage condition is a priori satisfied for RK1.Two well known examples for B1 are the explicit and the implicit Euler-schemes In the tableau for B(2)1 in Eq. ( 17) the higher order ELDIRK-method has four unknowns ĉ2 , â21 , b1 , b2 .To obtain order r = 2 with stage R = 2 reference to Table 1 shows n 1 = 1, n 2 = 1, such that according to Table 1 m2 = 2 order conditions and R = 2 stage conditions have to be satisfied.Since the stage condition ( 18) is a priori satisfied, with 3 conditions for 4 unknowns, there is one degree of freedom in the selection of the constants.The resulting 2 order conditions and 1 stage condition are taken from Table 1 as: Taking into account one degree of freedom in the above conditions a reformulation is performed next, with the aim to obtain explicit expressions for â21 , b1 , b2 in terms of a chosen value for ĉ2 .To this end, from (20.1) we deduce Inserting Eq. ( 21) into Eq.(20.2) renders Finally, from Eq. ( 20.3) we obtain For given value of ĉ2 the unknowns â21 , b1 , b2 in Eqs. ( 21)-( 23) can be evaluated, subject to the condition

Embedding of the explicit and implicit Euler method
As a first example, an embedding of the explicit Euler-scheme in Eq. ( 19) can be obtained as which is the classical modified Euler-method, see e.g.[19].An embedding of the implicit Euler-scheme in Eq. ( 19) can be obtained with the same choice for ĉ2 as For convenience, the resulting Butcher-Tableaus for the second order embedded Euler methods RK(2)1-Eul-exp and RK(2)1-Eul-imp are given in Table 2.

Remark 1
The choices ĉ2 = 1 2 in Eq. ( 25) and Eq. ( 26) are indeed rather arbitrary.A different choice would only result into a different Butcher-Tableau for the embedded RK(2)1-Eul-exp and RK(2)1-Eul-imp in Table 2.The procedure can be compared to the classical choice of constants for derivation of the modified Euler-method, the Heun-method and the Ralston-method, see e.g.[19].All of these are second order explicit RK-methods where none of them has an advantage over the other two.

Evaluation of the new second order Butcher-Tableaus
Comparing the new Butcher-Tableaus in Table 2 reveals the following stages R, Q for the embedded RK (r ) p schemes: According to Sect.2.2 on ELDIRK-methods, the practical implementation of both embedded methods in Table 2 can be split into two steps: Firstly, the unknown slope k 1 of the related 1-stage Euler-method is calculated, for RK1-Eul-exp as an explicit update, for RK1-Eul-imp by a nonlinear Newton algorithm as explained e.g. in [22].Then, according to Eq. ( 9) the additional slope is obtained in a second (post-processing) step, in an explicit manner as (observe In this way, both two-stage methods, the forward scheme RK2-Eul-exp and the backward scheme RK2-Eul-imp, require only one extra function evaluation to increase the order from q = 1 to r = 2, without solving another nonlinear system at all.Note, that this feature holds for linear as well as for nonlinear ordinary differential equations of the form (2.1).

General setting
This section intends to develop a general procedure to embed a Runge-Kutta method of order 2, subsequently denoted as RK2 with Butcher array B2, by an ELDIRK-method of order 3, subsequently denoted as RK3, thus resulting into an embedded RK-method RK(3)2 with extended Butcher array As before, RK-coefficients without a hat are assumed to be known, and RK-coefficients with a hat are assumed to be unknown.Note, that in order to become an ELDIRK-method according to definition (8), we require â33 = 0. We also assume that the stage conditions are a priori satisfied for RK2.Some examples for B2 are the embedded implicit modified Euler-method and the implicit Euler-scheme of Table 2 as well as the trapezoidal method or RK2-Ell of Ellsiepen [7] according to Eq. ( 1) In the tableau for B(3)2 in Eq. ( 29) the higher order ELDIRK-method has six unknowns ĉ3 , â31 , â32 , b1 , b2 , b3 .To obtain order r = 3 with stage R = 3 reference to Table 1 shows n 1 = 1, n 2 = 1, n 3 = 2, such that according to Table 1 m3 = 4 order conditions and R = 3 stage conditions have to be satisfied.Since the stage conditions (30) are a priori satisfied, there is one degree of freedom in the selection of the constants.The resulting 4 order conditions and 1 stage condition are taken from Table 1 as: Taking into account one degree of freedom in the above conditions a reformulation is performed next, with the aim to obtain explicit expressions for â31 , â32 , b1 , b2 , b3 in terms of a chosen value for ĉ3 .From (33.2) we deduce Inserting (34) into Eq.(33.3) multiplied with 2 renders With the results ( 34) and ( 35), Eq. (33.1) By means of Eq. (33.5) and we exploit Eq. (33.4) to obtain For given value ĉ3 the unknowns â31 , â32 , b1 , b2 , b3 in Eq. ( 29) can be evaluated by means of Eqs. ( 34)- (38) with the conditions

Embedding of the modified implicit Euler method
For the second order modified implicit Euler method the corresponding RK-coefficients are obtained from Table 2 as Consequently, with the choice Equations ( 36), ( 35) and (33.1) render Next, inserting (42) into Eq.( 38) and using Eq. ( 37) gives For convenience, the resulting new Butcher-Tableau for the third order embedded implicit Euler method RK(3)2-Eul is given in Table 3.

Evaluation of the new third order Butcher-Tableaus
The embedding schemes for the trapezoidal-rule RK(3)2-Trap and the Ellsiepen method RK(3)2-Ell are derived in a similar fashion as to the implicit Euler method RK(3)2-Eul in Appendix A. Comparing the new Butcher-Tableaus in Table 3 reveals the following stages R, Q for the embedded RK (r )q schemes: According to Sect.2.2 on ELDIRK-methods, the practical implementation of both embedded methods in Table 3 can be split into two steps: Firstly, the unknowns k 1 , k 2 of the related 2-stage method are calculated: For RK2-Eul as explained in Sect.3.3 including Eq. ( 28), for RK2-Trap and RK2-Ell, respectively, by a nonlinear Newton algorithm as explained e.g. in [22].Then, according to Eq. ( 9) the additional slopes are obtained in a second (post-processing) step, in an explicit manner as (observe a 33 = 0) In this way, RK3-Eul being a three-stage ELDIRK-method requires one extra function evaluation to increase the order from q = 2 for RK2-Eul to r = 3.In other words, RK3-Eul requires two extra function evaluations to increase the order from q = 1 for the classical Euler methods (backward and forward) to r = 3. RK3-Trap has already been derived in [22], however from a somewhat different viewpoint.Both, RK3-Trap and RK3-Ell are three stage ELDIRK-methods and require one extra function evaluation to increase the order from q = 2 to r = 3, without solving another nonlinear system at all.As mentioned before, these features hold for linear as well as for nonlinear ordinary differential equations of the form (2.1).

Remark 2
In [18] stability is investigated for certain subclasses of DIRK, such as EDIRK.However, due to the different structure of ELDIRK compared to EDIRK, the results of [18] cannot directly be applied.Therefore, additional stability investigations are required for the new ELDIRK-methods in Tables 2 and 3, respectively.This non-trivial task, however, is not in the scope of the current paper and therefore will be a task of future work.

Revision of goal oriented global error measures applied to RK(r)q and RK(q)r
This section intends to provide a brief summary of goal oriented global temporal error measures introduced in [22], and its asymptotic behaviour for standard embedded RK(r)q in Eq. ( 12) and for reversed embedded RK(q)r in Eq. ( 13) Global temporal error estimation can be achieved within a mathematically sound way by means of two functionals as originally proposed for partial differential equations, see e.g.[2,8,27,30,36]: Firstly, a weak form for temporal discretization of the underlying initial value problem (IVP), and, secondly, a quantity of interest (QoI) or goal function, respectively, which enables to measure the related global error.For the temporal discretization, the former can be obtained as a weak Petrov-Galerkin formulation.Here the concept of the finite element method for space discretizations is adapted to the time domain, originally introduced in [1,15,26].To this end the underlying differential equation ( 2) is reformulated into a corresponding weak form by a premultiplication with a test function and an integration over the time domain I .The solution functions as well as the test functions are then approximated by means of appropriate interpolations.Conventionally, discontinuous Galerkin (dG) and continuous Galerkin (cG) methods are distinguished, see e.g.[8].Eventually, the weak form enables to derive RK-schemes which result into the Butcher array (7).A simple example is e.g.given in [35] for the backward (implicit) Euler scheme.Regarding the mathematical details for derivation of a general explicit S-stage RK-method on the basis of a weak form the interested reader is referred e.g. to Muñoz-Matute et al. [25].The QoI approach offers a great flexibility on global error estimation, since the goal function is a user specified quantity, see e.g.[17,21].
In this section, a global error estimate is obtained directly from a quantity of interest, for both, the standard embedded and the reversed embedded RK-method by means of the functional representation for the weak form.This strategy avoids the solution of a dual problem as performed e.g. in [25].The dual problem runs backwards in time [2,17,21,31], which due to additional storage requirements might become a serious drawback when applied e.g. to FEM or respectively to FDM simulations.

Global error and effectivity index for standard embedded RK(r)q
We let y ∈ V be the solution of a related weak form of the IVP (2) where V is a suitable functional space.Also, y q τ and y r ,q τ are two solutions related to standard embedded RK(r)q of Eq. ( 12).Note, that these functions are continuous for the cG-method, whilst for the dG-method they are allowed to be discontinuous at the nodal points, but are taken to be continuous within the time intervals I i , see e.g.[8,35].As mentioned above, the resulting weak forms can be exploited to derive RK-schemes in the form of the Butcher array (7), or to formulate dual problems, which however is not in the scope of this paper, see e.g.[25].
In the sequel, the above functions are used to define two global semi-discretization (temporal) errors and a global error estimate as 1.E r ,q (t) := y(t) − y r ,q τ (t), 2. E q (t) := y(t) − y q τ (t), 3. η r ,q (t) := y r ,q τ (t) − y q τ (t), t ∈ I .( Next, following [9] and [30], we introduce a goal quantity of interest such that two global error measures E (., .): V × V → R can be distinguished, an error of the quantity and respectively a quantity of the error where for the second case the definition for E q in (46.2) has been used.Clearly, for linear goal functions, both error measures coincide.For simplicity, if there is no danger of confusion, the upper index eq or respectively qe is dropped from now on.Note, also that both functionals E (•, •) share the property A Taylor series of E (y, y q τ ) around y where Eq. ( 49) has been used.E T (y q τ , y q τ ; E q ) := D y q τ E (y q τ , y q τ ) • E q is the Gateaux-derivative (or tangent, respectively) of E (•, y q τ ) with respect to the first argument.Note, that the tangent is linear in the third argument, which is indicated by a semicolon.
As a concrete example for an QoI, in this work we will consider the time point error function For an alternative, a time integrated goal function is considered in [22].
The following relations are introduced in [22] for general goal functions in Eq. (51.1) with related error measures in Eq. (52.2) valid for "small" time steps τ defined in Eq. (3.4), 1. Q(E q )= E (y, y q τ ) ≈ E T (y q τ , y q τ ; E q ) ≈ C r ,q τ P q 2. Q(E r ,q )=, E (y, y r ,q τ ) ≈ E T (y r ,q τ , y r ,q τ ; E r ,q ) ≈ C r ,q τ P r ,q 3. Q(η r ,q )=, E (y r ,q τ , y q τ ) ≈ E T (y q τ , y q τ ; η r ,q ) ≈ C η r ,q τ Z r ,q . (52) Here, C r ,q , C r ,q and C η r ,q are some constants, and P q , P r ,q , Z r ,q are three convergence orders which satisfy, [19,22], P q = P r ,q = Z r ,q = q. (53) Remark 3 1.In Lambert, p. 58, it has been pointed out, while the local truncation error e i for an RKq-method in Eq. (52.1) is O(τ q+1 i ), the global truncation error Q(E q ) in Eq. ( 53) is O(τ q ) with τ according to Eq. (3.4).That is, one power of τ has been lost owing to the process of accumulation and which leads to P q = q in Eq. ( 53).This property is also verified for the classical methods RK1-Eul, RK2-Trap and RK2-Ell by the numerical results in the forthcoming Sect.6, see e.g. ), provided the conditions Table 4 Parachute with viscous damping for goal function in Eq. ( 60): results for four standard embedded RK-methods and its reversed methods with equidistant step sizes RK-Method RK(r)q Z r ,q P r ,q P q = q RK(2)1-Ell 1.00e+00 Global convergence orders, top, for standard embedded RK(r)q Z r ,q , P r ,q , P q = q according to Eq. ( 52) and, bottom, for reversed embedded RK(q)r Z q,r , P q,r , P r = r according to Eq. (57) (14) hold.Regarding the global error estimate Q(η r ,q ) in Eq. (52.2), the order condition (14.1) for the initial value does not hold, such that the global error estimate Q(η r ,q ) in Eq. (52.2) is only order O(τ q i ), and which leads to Z r ,q = q in Eq. ( 53), see also e.g.[32] and references therein.
The above relations allow to evaluate the quality of an error estimate.To this end, we borrow from Szabo and Babuska [34] the definition for the effectivity index originally formulated for finite element error estimation as the ratio of the estimated global error E (y r ,q τ , y r τ ) and the global error E (y, y q τ ).Then, the effectivity index becomes in the limit τ → 0 lim Consequently, the error estimate Q(η r ,q ) is a κ-multiple of the error Q(E q ).The effectivity index I q eff approaches some constant κ, which depends on the specific RK-algorithm and the individual problem to be solved.In general, the effectivity index cannot be one (even in the limit), which is a main limitation of standard embedded RK-methods.

Global error and effectivity index for reversed embedded RK(q)r
Analogously to the relations Eq. ( 46) two global errors and an error estimate are defined 1. E q,r (t) := y(t) − y q,r τ (t), 2. E r (t) := y(t) − y r τ (t), 3. η q,r (t) := y q,r τ (t) − y r τ (t), t ∈ I .( A relation between these quantities is Then, analogously to (52) the following relations are introduced in [22] for reversed embedded RK(q)r in in ( 13), valid for "small" time steps τ defined in Eq. (3.4), as where C r , C r and C η q,r are some constants.For the global convergence order constants it holds, [22], 1. P r = r = q + 1, 2. P q,r = Z q,r = q = r − 1. (58) These properties are verified for the new ELDIRK methods RK3-Eul, RK3-Trap and RK3-Ell by the numerical results in the forthcoming Sect.6, see e.g.Table 4.Then, by means of Eq. ( 56) and exploiting the linearity of the tangent form in the third argument a related effectivity index becomes in the limit τ → 0 lim The interpretation of Eq. ( 59) is, that the error estimate Q(η q,r ) is asymptotically exact for the global error Q(E q,r ) of RKq,r.This is a main advantage of the reversed embedded method RK(q)r compared to standard embedded RK(r)q, [22].

Representative numerical examples
In the following, numerical results of the algorithms summarized in Table 3, the embedded implicit Euler method RK(3)2-Eul, the embedded trapezoidal rule RK(3)2-Trap and the embedded Ellsiepen method RK(3)2-Ell are presented for two test problems: Firstly, a parachute with viscous damping and, secondly, a laser beam simulation with the finite difference method.Results are presented for equidistant step lengths as well as for non-equidistant step lengths based on adaptive time-step control.We will also verify the expected order properties of the new third order ELDIRK methods RK3-Eul, RK3-Trap and RK3-Ell, and we will present results for the reversed embedded Runge-Kutta methods, denoted as RK(2)3-Eul, RK(2)3-Trap and RK(2)3-Ell.Thereby, the nonlinear equations resulting from the implicit RK-schemes are solved with a Newton method, as explained in [22].
For comparison, we present also results for the embedded RK-method of Ellsiepen according to Eq. ( 1) [7], denoted as RK(2)1-Ell.As explained before, both RK1-Ell and RK2-Ell are implicit, however, the Newton method iterates on the same slopes k j , which satisfy a residual obtained from Eq. (5.2).In analogy to previous notations, RK (1)2-Ell denotes the corresponding reversed RK-method.
For both, standard embedded RK(r)q and reversed embedded RK(q)r the error measure in Eq. (51.2) is investigated.Thus, by means of the notation • N := •(t N ) and a 2-norm, for standard embedded RK(r)q the following three global error functions are introduced: with global error estimate η r ,q N according to Eq. (46.3).We will verify the result lim τ →0 I q eff = κ in Eq. (54).For reversed embedded RK(q)r the following three global error functions are introduced: with global error estimate η q,r N according to Eq. (55.3).We will verify the result lim τ →0 I q,r eff = 1 in Eq. ( 59).As an alternative to the above time point goal function, results for a time integrated goal function are presented in [22].

Initial value problem (IVP)
The first numerical example is concerned with a singledegree-of-freedom system in dynamics.Here the second order ODE typically attains the form m ẍ + d ẋ + cx = F(t), where m is the mass, x is the position, d is a damping constant, c is an elastic constant and F(t) is some applied force, see e.g.[20].As a particular problem we consider the motion of a parachute with viscous damping d ẋ (alternatively a quadratic ansatz could be applied for the damping term, see e.g [20].)Moreover we set c = 0 and F = mg, where g is the gravitational acceleration.Using the definition v := ẋ for the velocity, and assuming zero velocity at time t = 0, renders a first order IVP of the form (2) as: In [20] its analytical solution is derived as Here v st is the stationary velocity, obtained from the condition v(t = ∞) = 0 = g − dv st /m.In the subsequent calculations the constants are chosen as: m = 70 kg, d = 20.5 kg/s, g = 9.81 m/s 2 .The initial time is t 0 = α = 0 and the total time is T = β − α = 10 s.

Global convergence investigations for equidistant step sizes
Figure 1 shows results for all four standard embedded RKmethods on the left hand side and its reversed embedded RK-methods on the right hand side for N = 10 equidistant time steps.For the velocity versus time v(t) in Fig. 1a) a difference to the exact solution is hardly visible for all RKmethods.The same holds for all reversed methods in Fig. 1b).Table 4 summarizes the global convergence orders Z r ,q , P r ,q , P q defined in Eq. ( 52) for all four embedded RKmethods and Z q,r , P q,r , P r defined in Eq. (57) for its reversed methods with equidistant step sizes.For standard embedded RK(r)q the approximation P q ≈ P r ,q ≈ Z r ,q ≈ q according to Eq. ( 53) is obtained, whereas for reversed embedded RK(q)r P r ≈ r = q + 1 > P q,r ≈ Z q,r ≈ q according to Eq. ( 58) is obtained.In this way, all RK-methods reach (and occasionally even surpass) its theoretical convergence orders.
Regarding the new ELDIRK methods in Table 3 the following is pointed out from the results in Table 4: The order for RK3-Ell and RK3-Trap is increased by one compared to the implicit methods RK2-Ell and RK2-Trap, respectively, with only one additional function evaluation according to Eq. (45), whilst the order of RK3-Eul is increased by two compared to the implicit method RK1-Eul with two additional function evaluations according to Eq. ( 28) and Eq.(45), respectively.There is no need to solve another nonlinear system, although the IVP (62) is nonlinear.1 Parachute with viscous damping: Results for four standard embedded RK-methods and its reversed methods with equidistant step sizes.Left for standard embedded RK(r)q: a velocity versus time, c exact global error E q in (46.2) and estimated global error η r ,q in Eq. (46.3), e effectivity index I q eff in Eq. (54).Right for reversed embedded RK(q)r: b Velocity versus time, d exact global error E q,r in (55.1) and estimated global error η q,r in Eq. (55.3), f effectivity index I q,r eff in Eq. ( 59) 123 Table 5 Parachute with viscous damping for goal function in Eq. ( 60): results for four standard embedded RK-methods and its reversed methods with equidistant step sizes RK-Method  Top standard embedded RK(r)q: exact global errors Q(E q ), Q(E r ,q ), estimated global error Q(η r ,q ) in Eq. ( 60), and effectivity index I q eff in Eq. ( 54).Bottom reversed embedded RK(q)r: exact global errors Q(E r ), Q(E q,r ), estimated global error Q(η q,r ) in Eq. ( 61) and effectivity index I q,r eff in Eq. ( 59) A detailed documentation of further results on the the global convergence analysis for all eight methods with equidistant step sizes is provided in Table 5. Regarding standard embedded RK(r)q we present values for the exact global errors Q(E q ), Q(E r ,q ), the estimated global error Q(η r ,q ) in Eq. ( 60) and the effectivity index I q eff in Eq. ( 54).Regarding reversed embedded RK(q)r we present the exact global errors Q(E r ), Q(E q,r ), the estimated global error Q(η q,r ) Eq. ( 61) and the effectivity index I q,r eff in Eq. (59). Figure 1 and Figure 2 provide some corresponding graphical illustrations for the goal function in Eq. (60). Figure 1.c) shows the exact global error E q in (46.2) and the estimated global error η r ,q in Eq. (46.3).The corresponding effectivity index I q eff in Eq. ( 54) is displayed in Fig. 1.e).Figure 1.d) shows the exact global error E q,r in (55.1) and the estimated global error η q,r in Eq. (55.3).The corresponding effectivity index I q,r eff in Eq. ( 59) is displayed in Fig. 1.f).
The results in Table 5 and in Fig. 1.e) reveal, that the effectivity index for standard embedded RK(r)q tends to some constant, that is I q eff → κ according to Eq. (54).However, this constant is quite a distance away from 1.0, where I q eff may be as bad as κ ≈ 2.51, that is, only overestimation occurs in this example.This however is not a general trend of standard RK-methods, see [22], where both under-and overestimation occurs even for the same example.Consequently, the error estimates for standard embedded RK(r)q are not very accurate or robust.Contrary, the effectivity index I q,r eff for reversed embedded RK(q)r in Table 5 and in Fig. 1.f) approaches one, meaning it is asymptotically exact.We point out, that this result holds for all reversed embedded RK-methods, RK(1)2-Ell, RK(2)3-Trap, RK(2)3-Ell and RK(2)3-Eul.Table 5 further reveals, that the error estimates Q(η r ,q ) and Q(η q,r ) are some distance away from each other for τ i = 1.0 s, whilst for τ i = 1.0 × 10 −3 s the relation Q(η r ,q ) ≈ Q(η q,r ) is achieved.
Figure 2 depicts global error measures vs. time step of standard embedded RK-methods, and its reversed methods with equidistant step sizes.The subfigures a, c, e on the left refer to standard embedded RK(r)q with exact global errors Q(E r ,q ), Q(E q ) and estimated global error Q(η r ,q ).The subfigures b, d, f on the right refer to reversed embedded RK-method RK(q)r with exact global errors Q(E q,r ), Q(E r ) and estimated global error Q(η q,r ). Figure 2d compared to Fig. 2c demonstrates the superior convergence results of the new methods RK(1)2-Ell, RK(2)3-Eul, RK(2)3-Trap and RK(2)3-Ell which is documented by its higher slopes.

Adaptively refined step sizes
Next, results of adaptive step size selection according to Eq. ( 16) based on a posteriori error estimation with embedded RK-methods are presented.The tolerance for the estimated error in Eq. ( 16) is ε = 0.1.Figure 3 exhibits results for all four standard embedded RK-methods on the left hand side and its reversed embedded RK-methods on the right hand side.For the velocity versus time v(t), both in Fig. 3a and in Fig. 3b, a difference to the exact solution is hardly visible for all RK-methods. Figure 3c shows the exact global error E q in (46.2) and the estimated global error η r ,q in Eq. (46.3). Figure 3d shows the exact global error E q,r in (55.1) and the estimated global error η q,r in Eq. (55.3).Moreover, the Fig. 3e, f depict the time step lengths τ i resulting from the adaptive selection according to Eq. ( 16).
Table 6 summarizes the values for the exact global errors Q(E r ,q ), Q(E q ) and the estimated global error Q(η r ,q ) in Eq. ( 60) with corresponding effectivity index I q eff in Eq. (54) for RK(r)q; as well as the exact global errors Q(E r ), Q(E q,r ) Fig. 2 Parachute with viscous damping: results for four standard embedded RK-methods and its reversed methods with equidistant step sizes.Left for standard embedded RK(r)q: a, c exact global errors Q(E r ,q ), Q(E r ) and e estimated global error Q(η r ,q ) in Eq. (60).Right for reversed embedded RK(q)r: b, d exact global errors Q(E q,r ), Q(E r ) and f estimated global error Q(η q,r ) in Eq. ( 61) Fig. 3 Parachute with viscous damping: Results for four standard embedded RK-methods and its reversed methods with adaptively refined step sizes: left for standard embedded RK(r)q: a velocity versus time, c exact global error E q in (46.2) and estimated global error η r ,q in Eq. (46.3).e time steps τ i .Right for reversed embedded RK(q)r: b velocity versus time, d exact global error E q,r in (55.1) and estimated global error η q,r in Eq. (46.3).f Time steps τ i Table 6 Parachute with viscous damping: results for four standard embedded RK-methods and its reversed methods with adaptively refined step sizes  Top, standard embedded RK(r)q: exact global errors Q(E q ), Q(E r ,q ), estimated global error Q(η r ,q ) in Eq. ( 60) and effectivity index I q eff in Eq. ( 54).Bottom, reversed embedded RK(q)r: exact global errors Q(E r ), Q(E q,r ) and estimated global error Q(η q,r ) in Eq. ( 61), effectivity index I q,r eff in Eq. ( 59) and estimated global error Q(η q,r ) in Eq. ( 61), with corresponding effectivity index I q,r eff for RK(q)r in Eq. (59) for RK(q)r.The first order standard RK-method RK1-Ell has a comparatively large global error Q(E q ) = 1.67 for standard RK-method.The third order reversed RK-method RK3-Ell achieves the lowest values for the exact global error The effectivity indices I q eff of the standard embedded methods RK(r)q range from 1.74 to 2.49, meaning that in this example all errors are overestimated with a range that is not very robust.Contrary, the range of I q,r eff is from 0.938 to 1.14, that is, these values have a close range around one, which compared to the range of I q eff appears to be much more robust.

Remark 4
Of course, the results in this subsection on adaptively refined step sizes depend strongly on the maximum allowable error ε in Eq. (16).Thereby, the approximation for "small" time step sizes as e.g. in (52) or (57), respectively, is not satisfied anymore compared to the previous subsection on equidistant time steps in the limit τ → 0. As discussed in [22], adaptive time step control by means of the scheme (16) only controls the local error.Bounding the global error by means of adaptive strategies requires further techniques, which are not within the scope of this paper, see e.g.[21].

Initial boundary value problem
The second example is concerned with the temperature distribution on a plate due to a laser beam movement typically occurring in selective laser melting (SLM) or respectively selective electron beam melting (SEBM) in additive manufacturing processes, [24,29,38].Since our focus lies on the performance investigation of different Runge-Kutta methods the numerical simulations are restricted to a solid phase of a two-dimensional area as illustrated schematically in Fig. 4. Within a rectangular sheet, the laser beam follows a sequence of straight lines representing the letters LTM.The underlying initial boundary value problem (IBVP) for the temperature θ reads In Eq. (64.1) ρ is the density of the material, c s is the specific heat capacity, λ is the conductivity and r is the heat source.Moreover, = [a × b], a = 400 mm, b = 200 mm, and ∂ describe the spatial domain and its boundary, respectively, of the rectangular sheet, and is the Laplace operator, which in Cartesian coordinates x, y reads The Dirichlet boundary condition Eq. (64.2) prescribes a constant temperature θ R at the boundary.In Eq. (64.3) initial conditions θ 0 are prescribed.In our simulations the temperatures for boundary and initial conditions in (64) are θ R = 172 o C and θ 0 = 172 o C, respectively.The initial time is t 0 = α = 0 and the total time is T = β − α = 25 s.The material parameters for the simulation are given in Table 7.In order to model the energy input of the electron beam in the SEBM process, we let b represent the total beam path as illustrated in Fig. 4 and represents the actual position of the laser beam along the path.Here x 0 ∈ is the initial position of the laser beam, related to the point P 0 in Fig. 4. Then the energy input of the electron beam in the SEBM process is modeled as a time dependent volumetric heat source r (t, x) moving along b with power density Here we have used: the absorption factor A b , the absorption length β p , the laser beam power P b and the beam radius r b , see e.g.[28,29].For the two-dimensional simplification in our numerical example we neglect the penetration depth z of the beam.Moreover, in Eq. (66) the δ-function is being used to describe a unit point source positioned at the point x, see e.g.[37].The total arc-length of the beam scan path illustrated   Table 8 Laser beam simulation: convergence orders Z r ,q , P r ,q , P q = q according to Eq. ( 52) and Z q,r , P q,r , P r = r according to Eq. (57) for four standard embedded RK-methods and its reversed methods with equidistant step sizes RK-Method RK(r)q Z r ,q P r ,q P q = q RK(2)1-Ell 1.05e+00The spatial discretization and the time integration for the IVBP (64) are extensively described in [22].Thereby, we also account for the effect of order reduction by means of an effective temporal "node-to-node" discretization technique.
Figure 5 provides the resulting temperature distribution at four time steps, here generated with the classical 4th-order Runge-Kutta method with Butcher-Tableau All other RK-methods introduced above did not show any visible difference for this illustrative simulation.The laser beam starts at the upper point P 0 of the letter L, introduced in Fig. 4, such that the current position of the beam spot can be seen from the current status of the letters LTM.The final phase from laser scanning at t b = 15.986s to the total time T = 25 s shows some redistribution of the temperature field as a result of stationary effects.

Global convergence investigations for equidistant step sizes
This subsection provides convergence investigations for equidistant step sizes analogously to the previous example.Since an analytical solution is not available and in order to reduce the numerical effort, only the initial phase of the laser beam movement is considered.The total time is set to T R = 0.21 s, such that the distance from the initial node P 0 to the final point P 1 , illustrated in Fig. 4, becomes P 0 P 1 = v b T R = 10.5 mm.A highly accurate solution, subsequently denoted as "reference solution", is generated using the classical 4th-order Runge-Kutta method in Eq. ( 67) in N = 200 equidistant time steps.Table 8 shows that all RK-methods attain (and in most cases even surpass) its theoretical convergence orders according to Eq. (53) for standard RK(r)q and according to Eq. (58) for reversed RK(q)r.As in the previous example, the new ELDIRK methods RK3-Eul, RK3-Trap and RK3-Ell attain the third order properties, with no need to solve another nonlinear system.
Table 9 provides the exact global errors Q(E q ), Q(E r ,q ), the estimated global error Q(η r ,q ) in Eq. ( 60), and the corresponding effectivity index I q eff in Eq. (54).Regarding reversed embedded RK(q)r we present the exact global errors Q(E r ), Q(E q,r ), the estimated global error Q(η q,r ) in Eq. ( 61) and the effectivity index I q,r eff in Eq. ( 59).The results in Table 9 reveal, that the effectivity index for standard embedded RK(r)q tends to some constant, that is I q eff → κ = 1 according to Eq. ( 54).This effect, different to the previous example, has not been investigated so far.
It might be attributed to the fact, that the matrix A in heat transfer equation ( 64) is constant.As expected, the effectivity index I q,r eff for reversed embedded RK(q)r in Table 9 approaches one, meaning it is asymptotically exact.We point out, that this result holds for all reversed embedded RK-methods, RK(1)2-Ell, RK(2)3-Trap, RK(2)3-Ell and RK(2)3-Eul.
Figure 6 depicts errors vs. time step of standard embedded RK-methods, and its reversed methods with equidistant step sizes.The subfigures a, c, e on the left refer to standard embedded RK(r)q with exact global errors Q(E r ,q ), Q(E q ) and estimated global error Q(η r ,q ).The subfigures b, d, f on the right refer to reversed embedded RK-method RK(q)r with exact global errors Q(E q,r ), Q(E r ) and estimated global error Q(η q,r ). Figure 6d compared to Fig. 6c demonstrates the superior convergence results of the new methods RK(2)3-Eul, RK(2)3-Trap and RK(2)3-Ell documented by its higher slopes.

Summary and outlook
In this paper we have proposed a general approach to embed low order implicit time integration methods into third order Runge-Kutta schemes.In order to obtain a systematic representation for the higher order methods, we introduce the notion of ELDIRK-methods.These are defined to have an explicit last stage in the general Butcher tableau, contrary to the class of EDIRK methods, which according to Kvaernø et al. [18] have an explicit first stage in the general Butcher array.The consequence is, that there is no need to solve an additional nonlinear system of equations for the ELDIRKmethod compared to the embedded schemes.This feature Table 9 Laser beam simulation: results for four standard embedded RK-methods and its reversed methods with equidistant step sizes RK-Method τ i = 5.250e−02 s τ i = 2.100e−03 s RK(r)q Q(E q ) Q(E r ,q ) Q(η r ,q ) I q eff → κ Q(E q ) Q(E r ,q ) Q(η r ,q ) I q eff → κ  Top, standard embedded RK(r)q: exact global errors Q(E q ), Q(E r ,q ), estimated global error Q(η r ,q ) in Eq. ( 60), and effectivity index I q eff in Eq. (54).Bottom, reversed embedded RK(q)r: exact global errors Q(E r ), Q(E q,r ), estimated global error Q(η q,r ) in Eq. (61) and effectivity index I q,r eff in Eq. ( 59) holds for linear as well as for nonlinear ordinary differential equations.
The approach has been applied to three prominent low order implicit methods, that is, the implicit Euler-method, the trapezoidal-rule and the Ellsiepen method.In this way, for each of the methods we obtain an extended Butcher tableau.
The new embedded standard RK-methods, RK(3)2-Eul, RK(3)2-Trap and RK(3)2-Ell and the reversed methods RK(2)3-Eul, RK(2)3-Trap and RK(2)3-Ell have been tested extensively in two numerical examples: Firstly, a parachute with viscous damping and, secondly, two-dimensional temperature evolution due to laser beam movement, and they have been compared to the RK(2)1 method of Ellsiepen [7], Eq. (1) in this work.For both examples it has been verified, that with only two additional function evaluations for RK1-Eul, the ELDIRK-method RK3-Eul attains the third order properties.Similarly, third order properties are obtained for the new ELDIRK-methods RK3-Trap and respectively for RK3-Ell with only one additional function evaluation for RK2-Trap and respectively for RK2-Ell.We point out, that this behaviour holds also for nonlinear IVPs.
Concerning further extensions, being only half implicit, a stability analysis of all new ELDIRK-methods is of interest.Also, strategies to control the global error, see e.g.[21], should be considered for embedded RK-methods, based on the proposed global error measures of [22] and recalled in this work.6 Laser beam simulation: errors vs. time step for four standard embedded RK-methods and its reversed methods with equidistant step sizes.Left standard embedded RK(r)q: a, c exact global errors Q(E r ,q ), Q(E r ) and e estimated global error Q(η r ,q ) in Eq. (60).Right reversed embedded RK(q)r: b, d exact global errors Q(E q,r ), Q(E r ) and f estimated global error Q(η q,r ) in Eq. (61) RK(3)2-Trap: (R, S) = (3, 2), (r , p) = (3, 2)

Fig.
Fig.1Parachute with viscous damping: Results for four standard embedded RK-methods and its reversed methods with equidistant step sizes.Left for standard embedded RK(r)q: a velocity versus time, c exact global error E q in (46.2) and estimated global error η r ,q in Eq. (46.3), e

in Fig. 4
is l b = | b | = 799.3mm, such that with the scan velocity v b in Table 7 the resulting heating time for laser scan-ning becomes t b = l b /v b = 15.986s.The total number for the spatial discretization visualized in Fig. 4.b is n h = 1891, see[22] on more details.Details on a reduced total beam path length and the corresponding step size selections for a convergence analysis are given in Sect.6.2.2.Contrary to the previous example, results are presented only for equidistant step lengths.Results for non-equidistant step lengths based on adaptive time-step control would give no new insight for this example compared to those presented recently in[22].

Fig.
Fig.6 Laser beam simulation: errors vs. time step for four standard embedded RK-methods and its reversed methods with equidistant step sizes.Left standard embedded RK(r)q: a, c exact global errors Q(E r ,q ),

Table 7
Laser beam simulation: parameters for the simulation