1 Introduction

1.1 State of the art on ELDIRK methods

Most physical problems involve a time-dependent process and consequently require a reliable numerical integration in time of the underlying ordinary differential equations. 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 to tackle these complex tasks is provided by so-called \(S\)-stage Runge-Kutta (RK) methods, where explicit and implicit schemes are distinguished, [6, 9, 17]. All the individual RK schemes require some strategy on step size selection. For this purpose, Fehlberg suggested an efficient strategy on the basis of a posteriori local error estimates, by means of so-called embedded Runge–Kutta methods, [10,11,12]. 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 Runge–Kutta schemes have been developed, see e.g. Bogacki & Shampine [4, 13].

Based on the classification of RK-methods for Diagonal Implicit Runge–Kutta (DIRK)-methods, in Kværnø, Nørsett and Owren[16], a new class of ELDIRK methods is introduced in Mahnken [20, 21]. Here, the single explicit first stage introduced in Kennedy and Carpenter [15] in the general Butcher array is replaced by an Explicit Last stage. As a consequence to the numerical implementation, only the implicit leading stages have to be solved iteratively, e.g. by a Newton method, whereas the additional explicit last stage requires only one extra function calculation which is in accordance to Fehlberg’s idea. In this way, a general procedure has been proposed in Mahnken [21] to obtain second-order ELDIRK schemes for the embedding of the explicit and implicit Euler method and to obtain third-order ELDIRK schemes for the embedding of the implicit Euler-method, the trapezoidal-rule, the midpoint -rule and the Ellsiepen method, respectively. Moreover, in Westermann and Mahnken [23] a fourth-order ELDIRK scheme for embedding of an implicit third-order RK method is developed.

In Mahnken [20, 21]; Westermann and Mahnken [23] various numerical examples of physical problems in dynamics, continuum mechanics and material sciences verify the expected higher convergence rates. However, regarding stability the investigations in Westermann and Mahnken [23] reveal severe numerical problems of several ELDIRK methods. As a remedy, an accompanying eigenvalue determination allows a step size selection for the critical time step thus leading to stability for each of the related ELDIRK schemes.

1.2 Key objectives and outline of the present work

The present paper is directed to the extensive investigation of stability properties of the general ELDIRK methods with second-order and third-order properties. As a point of departure, we are confronted with the above-mentioned numerical instabilities revealed in Westermann and Mahnken [23], which at first sight could be explained by the semi-explicit nature of the ELDIRK Butcher tableaus. Moreover, there exists a general statement on the stability of RK-methods, see e.g. Alexander [3]: ’To be A-stable, and possibly useful for stiff systems, a Runge–Kutta formula must be implicit’. On this basis, the question remains, whether it is possible to construct ELDIRK methods, which are A-stable, although explicit in the last stage. To the best knowledge of the authors, no semi-explicit Runge–Kutta schemes with the desired A-stability property, typical for implicit methods, has been proposed so far. Thus, the key contributions are as follows:

  • We perform stability investigations based on the general formulations in Mahnken [21] for second- and third-order ELDIRK methods.

  • We construct novel second- and third-order A-stable ELDIRK methods in agreement with order and stage conditions.

  • Regarding the third-order methods, A-stable and not A-stable methods will be distinguished. For the latter class, the regions of absolute stability are obtained, which allow a step size selection by means of an accompanying eigenvalue determination.

  • Two numerical examples are concerned with the curing of a thermosetting material and the representative volume element (RVE)-modeling based on a Kobayashi-Warren-Carter phase field formulation for crystallinity and orientation.

The outline of the paper reads as follows. Section 2 presents some preliminaries on embedded Runge–Kutta methods applied for adaptive step size refinement. Section 3 recalls some basics of Mahnken [20] for a general approach to embed a Runge-Kutta method of order 1 (respectively order 2) by an ELDIRK method of order 2 (respectively order 3). Then, as a main contribution of this paper, in Sect. 4 the stability of ELDIRK methods is investigated. In particular, we will construct novel A-stable ELDIRK methods for the general approach proposed in Mahnken [20]. In Sect. 5 the expected results on convergence order and A-stability are verified for two examples on the curing of a thermo-setting material and a Kobayashi-Warren-Carter phase-field model for grain boundary migration based on orientation and crystallinity in polycrystals. Section 6 briefly summarizes the findings of the present paper and discusses potential improvements to be considered in the future.

1.3 Notation

Throughout the manuscript, several mathematical formulations occur, particularly for the presentation of a coupled phase-field approach. For simplicity and ease of reading the essential formulations are declared in advance. The notation of the present paper includes derivative with respect to time \({}^{.}\left( \cdot \right) \), higher order solution \(\hat{(\cdot )}\), scalar a, first order tensor \(\varvec{a}\), second order tensor \(\varvec{A}\), fourth-order tensor \({\mathbb {A}}\), gradient \(\nabla \left( \cdot \right) \), divergence \(\nabla \cdot \left( \cdot \right) \), and Laplacian \(\nabla ^2\left( \cdot \right) \). The inner product on a domain \(\Omega \) of two functions \(a[\varvec{x}]\) and \(b[\varvec{x}]\), \(\varvec{x}\in \Omega \), which are elements of an appropriate Sobolov space \({\mathcal {V}}\) is defined as a space integral

$$\begin{aligned} \int _\Omega (a\cdot b) dv =: <a, b>_\Omega , \end{aligned}$$
(1)

and where the chevron notation \(<,>\) is used for brief notation. The Gâteaux differentiation of a non-linear form \(F[\cdot ]:\) \(\underline{y} \rightarrow F[\underline{y}]\) at the point \(\underline{y}\) in the direction of a variation \(\Delta \underline{y}\) is defined as

$$\begin{aligned}&\text {D}_{\underline{y}} F\big [\underline{y}; \Delta \underline{y}\big ] = \lim _{\varepsilon \rightarrow 0} \frac{1}{\varepsilon } \left( F\big [\underline{y} + \varepsilon \Delta \underline{y}\big ] - F\big [\underline{y}\big ]\right) . \end{aligned}$$
(2)

Here and in subsequent presentations, a linear dependency of a functional on a variable \(\Delta \underline{y}\) is indicated by a semicolon.

2 Preliminaries on embedded Runge–Kutta methods

2.1 Generalized Runge–Kutta scheme

Consider the initial value problem (IVP)

(3)

where the solution typically represents a physical process of n variables dependent on time t. Moreover, is a constant scaling matrix, is a non-linear vector-valued function. The source (or respectively loading) term \( {f}\) and the initial value are given data. We note, that the representation (3) can be generalized to parabolic equations, see e.g. Schmich and Vexler [22]. An example for phase-field simulations is given in Sect. 5.2.

For brevity of presentation, in the remainder of this section, we restrict ourselves to temporal discretization. To this end, we perform a partition of the time interval as and use the notations

(4)

and where N is the number of time steps. Additionally, denotes the numerical approximations of a q-order RK-method for the initial value problem (3) at times . Then, starting with at the time-stepping algorithm of the RK method is generally written in terms of a slope function (or increment function, see Ernst Hairer [9]), of the method as

(5)

For an \(S\)-stage RK-method the algorithm (5) specifies as

(6a)
(6b)
(6c)

and where the weights \(b_j\), the knots \(c_j\) and the RK-coefficients \(a_{jl}\), respectively, satisfy the restrictions

$$\begin{aligned}{} & {} 1.~~ \sum \limits _{j=1}^{S} b_j = 1, \qquad 2.~~ 0\le c_j \le 1, \nonumber \\{} & {} 3.~~ \sum \limits _{l=1}^{S} a_{jl} = c_j,~~j = 1,\ldots ,S. \end{aligned}$$
(7)

Local and global errors are distinguished in the literature. The convergence studies at a later stage in this work are performed for the global error of RKq defined as, see e.g. Lambert [17]

(8)

2.2 Embedded Runge–Kutta methods

Two kinds of embedded RK methods are known in the literature, standard embedded and reversed embedded RK-methods.

The former, standard embedded RK(r)q as introduced in Fehlberg [10,11,12], can be represented by an extended Butcher tableau, Butcher [6]

$$\begin{aligned} ~~{{\underline{B}}}(r)q= \begin{array}{ll|llllllll} c_1 &{} &{} &{} a_{11} &{} \ldots &{} \ldots &{} \ldots &{} ~~a_{1S} \\ c_2 &{} &{} &{} a_{21} &{} &{} &{} &{} ~~a_{2S} \\ \vdots &{}&{} &{} \vdots &{} \vdots &{} \ddots \\ c_{S} &{} &{} &{} a_{S1} &{} a_{S2} &{} \ldots &{} ~~a_{SS-1} &{} ~~a_{SS} \\ \hline &{}&{}&{} \\ &{} &{} &{} {\hat{b}}_{1} &{} {\hat{b}}_{2} &{} \ldots &{}~~ {\hat{b}}_{S-1} &{}~~~ {\hat{b}}_{S} \\ &{} &{} &{} b_{1} &{} b_{2} &{} \ldots &{} ~~b_{S-1} &{} ~~~ b_{S} \\ \end{array}. \end{aligned}$$
(9)

Here, \(a_{jl}, b_j, c_j,\) with \(j, l \in \{1,\ldots , S\}\) are the coefficients in the Equations (7). Additionally the first row of coefficients \({\hat{b}}_{j}\) below the horizontal line in (9) are the weights according to Fehlberg, [10,11,12]. In this way, the applied integration of order q is embedded by a second method with higher order \(r > q\), however, identical Runge-Kutta coefficients \(a_{jl}\). For \(a_{jl} = 0, l\ge j\) the slopes \( {k}_j\) in Eq.(6b) 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 Mahnken [20].

The second class of embedded RK-methods is the reversed embedded RK(q)r also known as local extrapolation, see e.g. Lambert [17]. It is given as a pair of time-stepping algorithms according to (5),

(10)

Note, that the lower order method RKq,r, \(q\le r-1\) and stage \(Q\) uses the initial value of the higher order method RKr at each time step, which explains the upper index \(q,r\).

Extensive investigations for both types of algorithms with a focus on ELDIRK methods have been done in Mahnken [20]. In this paper, only reversed embedded RK(q)r methods will be considered.

From the above pair of algorithms (10)

(11)

becomes a posteriori estimate for the local error of RK(q)r. These are used to control the step-size in Eq.(10) according to the formula

(12)

where \(||\bullet ||\) denotes some norm of the related quantity, e.g. the L2-norm or the maximum-norm. Also, \({\varepsilon }\) is a tolerance for the estimated error, 0.9 is a safety factor and 2 and 0.2 are values for the upper and lower bounds of the step length. Since the error estimate in Eq.(11) is dependent on the (unknown) final time-step , in practice is obtained in a trial and error procedure. We point out, that the approach for local error estimation by means of the extended Butcher tableau in Eq.(9) is applicable to both, explicit and implicit embedded RK-methods.

3 Review of some low-order ELDIRK methods

3.1 The class of ELDIRK methods

A classification of RK methods is introduced in Kværnø et al. [16]. Here, in particular Diagonal Implicit Runge-Kutta (DIRK)-type methods are of particular interest, where in the Butcher tableau in Eq.(9) one has lower triangular coefficient matrix with elements \(a_{ij} = 0\) for all \(j > i\) and \(a_{ii} \ne 0\) for at least one i.

An embedding of implicit DIRK methods has been considered for some prominent low-order methods in Mahnken [20, 21]. Let \(Q\) and \(R> Q\) be the stages for the pair of a \(q\)-order method embedded by an \(r\)-order RK-method, where \( r> q\). Then, according to Mahnken [20, 21] 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.(9) are zero. For a systematic generalization of this convenient property the class of RK-methods is defined as

$$\begin{aligned}{} & {} \text{ ELDIRK: } \text{ Explicit } \text{ Last } \text{ stage } \text{ DIRK },\nonumber \\ {}{} & {} \quad \qquad \textrm{where}\quad ~ a_{jj}= 0, \quad Q< j \le R. \end{aligned}$$
(13)

This is in contrast to the class of EDIRK methods, which according to Kværnø et al. [16] have an explicit first stage in the general Butcher array in Eq.(9), that is (at least) \(a_{11} = 0\). The extended Butcher tableau for ELDIRK in Eq.(9) reads

$$\begin{aligned} {{\underline{B}}}(r)q= \begin{array}{lc|lll|llll} &{} &{} &{} &{} &{} 0 &{} \ldots &{} &{} 0 \\ {{{\underline{c}}}_{Q}}&{} &{} &{} \underline{{{\underline{A}}}}_{Q}&{} &{} 0 &{} \ldots &{} &{} 0 \\ {} &{}&{} &{} &{} &{} 0 &{} \ldots &{} &{} 0 \\ \hline \\ c_{Q+1} &{}&{} a_{Q+1,1}&{} \ldots &{} a_{Q+1\ Q} &{} 0 &{} \ldots &{} &{} 0 \\ c_{Q+2} &{}&{} a_{Q+2,1}&{} \ldots &{} a_{Q+2\ Q} &{} a_{Q+2 \ Q+1} &{} 0 &{} &{} 0 \\ &{}&{} \vdots &{} \vdots &{}&{} \ldots &{} \ddots &{} \ddots &{} \vdots \\ c_{R} &{}&{} a_{R1} &{} \ldots &{} a_{RQ} &{} a_{RQ+1} &{} \ldots &{} a_{R\ R-1} &{}0 \\ \hline &{} &{} &{} &{} \\ &{} &{} {{\hat{b}}}_1 &{} \ldots &{} {{\hat{b}}}_{Q} &{} {{\hat{b}}}_{Q+1} &{}\ldots &{} &{} {{\hat{b}}}_{R} \\ &{} &{} b_1 &{} \ldots &{} b_{Q} &{} 0 &{}\ldots &{} &{} 0 \\ \end{array} \nonumber \\ \end{aligned}$$
(14)

where \({{{\underline{c}}}_{Q}}= \left( c_{j}\right) _{j=1}^{Q}\). The lower triangular coefficient submatrix \(\underline{{{\underline{A}}}}_{Q}\) with elements \(a_{ij} = 0\) for all \(j > i\) and \(a_{ii} \ne 0\) for at least one i, and where \(i,j = 1,\ldots ,Q\), represents the underlying \(Q\)-stage method. In compact form, the full matrix of RK-coefficients is

$$\begin{aligned} \underline{{{\underline{A}}}}= \left[ \begin{array}{lclll} \underline{{{\underline{A}}}}_{Q}&{} {\underline{0}} \\ \underline{{{\underline{A}}}}_{RQ}&{} \underline{{{\underline{A}}}}_{RR}\\ \end{array} \right] , \end{aligned}$$
(15)

where \(\underline{{{\underline{A}}}}_{RQ}= \left( a_{ij}\right) _{i=\tiny {Q}+1,j=1}^{R,Q}\) is a submatrix which can be fully structured, and \(\underline{{{\underline{A}}}}_{RR}\) is a strictly lower triangular coefficient submatrix with elements \(a_{ij} = 0\) for all \(j \ge i\), and where \(i,j = Q+1,\ldots ,S\).

The practical implementation of ELDIRK can be split into two steps, Mahnken [20, 21]: Firstly, the unknowns \( {k}_j\), \(j=1,\ldots ,Q< R\) of the related implicit \(Q\)-stage method are determined by a nonlinear Newton algorithm as explained e.g. in Mahnken [20]. Then, in a second (post-processing) step, the additional slopes are obtained in an explicit manner by exploiting definition (13) from Eq.(6b) as

(16)

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 (3.1).

3.2 Second-order ELDIRK for embedding of RK1

This section recalls some basics for a general procedure to embed a Runge-Kutta method of order 1, subsequently denoted as RK1 with Butcher array \({{\underline{B}}}1\), by an ELDIRK-method of order 2, subsequently denoted as RK2. The resulting standard embedded RK-method RK(2)1 has an extended Butcher array \({{\underline{B}}}(2)1\):

$$\begin{aligned} ~~ {{\underline{B}}}1= \begin{array}{ll|llllllll} c_1 &{} &{} &{} a_{11} &{} \\ \hline &{} &{} &{} &{} \\ &{} &{} &{} b_1 &{} \\ \end{array}, \qquad ~~ {{\underline{B}}}(2)1= \begin{array}{ll|llllllll} c_1 &{} &{} &{} a_{11} &{} \\ {{\hat{c}}}_2 &{} &{} &{} {{\hat{a}}}_{21} &{} &{} \\ \hline &{} &{} &{} &{} \\ &{} &{} &{} {{\hat{b}}}_1 &{} {{\hat{b}}}_2 \\ &{} &{} &{} b_1 &{} \\ \end{array}. \end{aligned}$$
(17)

In these tableaus, RK-coefficients without a hat are assumed to satisfy the stage condition, Butcher [5] and Dormand and Prince [8]

$$\begin{aligned} \begin{array}{lllllll} j=1:&s^{(1)}= & {} \sum _l a_{1l} - c_1= & {} 0 ~~~\Longrightarrow ~~ a_{11} = c_1. \end{array} \end{aligned}$$
(18)

The additional three conditions for the four RK-coefficients with a hat (\( {{\hat{c}}}_2, {{\hat{a}}}_{21},\) \({{\hat{b}}}_1,\) \({{\hat{b}}}_2\)), to obtain order \(r=2\) with stage \(R=2\), are the 2 order conditions and 1 stage condition taken from Butcher [5]; Dormand and Prince [8] as:

$$\begin{aligned} a_1^{(1)}&= \sum \limits _{j} b_{j} - 1 = {{\hat{b}}}_1 + {{\hat{b}}}_2 - 1 = 0, \end{aligned}$$
(19a)
$$\begin{aligned} a_1^{(2)}&= \sum \limits _{j} b_{j} c_{j} - \frac{1}{2} = {{\hat{b}}}_1 c_1 + {{\hat{b}}}_2 {{\hat{c}}}_2 - \frac{1}{2} = 0, \end{aligned}$$
(19b)
$$\begin{aligned} s^{(2)}&= \sum \limits _l a_{2l} - {{\hat{c}}}_2 = {{\hat{a}}}_{21} - {{\hat{c}}}_2 = 0. \end{aligned}$$
(19c)

Here, additionally, the ELDIRK-condition Eq.(13) has been used leading to \({{\hat{a}}}_{22} = 0\). This results in the relations, Mahnken [21]

$$\begin{aligned} {{\hat{a}}}_{21}&= {{\hat{c}}}_2, \end{aligned}$$
(20a)
$$\begin{aligned} {{\hat{b}}}_2&= \displaystyle \frac{\frac{ 1}{ 2}- a_{11}}{{{\hat{a}}}_{21}-a_{11}}, \end{aligned}$$
(20b)
$$\begin{aligned} {{\hat{b}}}_1&= 1- {{\hat{b}}}_2. \end{aligned}$$
(20c)

The unknowns \({{\hat{a}}}_{21}, {{\hat{b}}}_1, {{\hat{b}}}_2\) in Eq.(20a) to Eq.(20c) can be evaluated, where the variable \({{\hat{c}}}_2\) can be chosen arbitrarily, however, subject to the condition, Mahnken [21]:

$$\begin{aligned} {{\hat{c}}}_2 \ne c_1. \end{aligned}$$
(21)

Table  1 summarizes two Butcher-Tableaus of second-order ELDIRK methods, namely the embedded explicit Euler-method RK(2)1-Eul-exp and the embedded implicit Euler-method RK(2)1-Eul-imp, developed in Mahnken [21].

3.3 Third-order ELDIRK for embedding of RK2

This section recalls some basics for a general procedure to embed a Runge-Kutta method of order 2, subsequently denoted as RK2 with Butcher array \({{\underline{B}}}2\), by an ELDIRK-method of order 3, subsequently denoted as RK3. The resulting embedded RK-method RK(3)2 has an extended Butcher array \({{\underline{B}}}(3)2\):

$$\begin{aligned} ~~ {{\underline{B}}}2= \begin{array}{ll|llllllll} c_1 &{} &{} &{} a_{11} &{} \\ c_2 &{} &{} &{} a_{21} &{} a_{22} &{} \\ \hline &{} &{} &{} &{} \\ &{} &{} &{} b_1 &{} b_2 \\ \end{array}, \qquad ~~ {{\underline{B}}}(3)2= \begin{array}{ll|llllllll} c_1 &{} &{} &{} a_{11} &{} \\ c_2 &{} &{} &{} a_{21} &{} a_{22} &{} \\ {{\hat{c}}}_3 &{} &{} &{} {{\hat{a}}}_{31} &{} {{\hat{a}}}_{32} &{} \\ \hline &{} &{} &{} &{} \\ &{} &{} &{} {{\hat{b}}}_1 &{} {{\hat{b}}}_2 &{} {{\hat{b}}}_3 \\ &{} &{} &{} b_1 &{} b_2 \\ \end{array}. ~~~~~ \end{aligned}$$
(22)

As before, in these tableaus RK-coefficients for RK2 without a hat are assumed to satisfy the stage conditions taken from Butcher [5]; Dormand and Prince [8] as

$$\begin{aligned} j=1: ~~&s^{(1)} = \sum _l a_{1l} - c_1 = 0 \end{aligned}$$
(23a)
$$\begin{aligned} j=2: ~~&s^{(2)} = \sum _l a_{2l} - c_2 = 0. \end{aligned}$$
(23b)

The additional conditions for the six unknowns with a hat (\( {\hat{c}}_3, {{\hat{a}}}_{31},\) \({{\hat{a}}}_{32},\) \({\hat{b}}_1, {\hat{b}}_2, {\hat{b}}_3\)) to obtain order \(r=3\) with stage \(R=3\), are 4 order conditions and 1 stage condition taken from Butcher [5]; Dormand and Prince [8] as:

$$\begin{aligned} a_1^{(1)}&= \sum \limits _{j} {\hat{b}}_{j} - 1 = {\hat{b}}_1 + {\hat{b}}_2 + {\hat{b}}_3 - 1 = 0 \end{aligned}$$
(24a)
$$\begin{aligned} a_1^{(2)}&= \sum \limits _{j} {\hat{b}}_{j} c_{j} - \frac{1}{2} = {\hat{b}}_1 c_1 + {\hat{b}}_2 c_2 + {\hat{b}}_3 {{\hat{c}}}_3 - \frac{1}{2} = 0 \end{aligned}$$
(24b)
$$\begin{aligned} a_1^{(3)}&= \frac{1}{2}\sum \limits _{j} {\hat{b}}_{j} c_{j}^2 - \frac{1}{6} = \frac{1}{2} \left[ {\hat{b}}_1 c_1^2 + {\hat{b}}_2 c_2^2 + {\hat{b}}_3 {{\hat{c}}}_3^2 \right] - \frac{1}{6} = 0 \end{aligned}$$
(24c)
$$\begin{aligned} a_2^{(3)}&= \sum \limits _{jl} {\hat{b}}_{j} a_{jl} c_{l} - \frac{1}{6} = {\hat{b}}_1 a_{11} c_1 + {\hat{b}}_2( a_{21} c_1 + a_{22} c_2 ) \\ \nonumber&\quad + {\hat{b}}_3 \left( {{\hat{a}}}_{31} c_1 + {{\hat{a}}}_{32} c_2 \right) - \frac{1}{6} = 0 \end{aligned}$$
(24d)
$$\begin{aligned} s^{(3)}&= \sum _l a_{3l} - {{\hat{c}}}_3 = {{\hat{a}}}_{31} + {{\hat{a}}}_{32} - {{\hat{c}}}_3 = 0. \end{aligned}$$
(24e)

This results in the relations, Mahnken [21]

$$\begin{aligned} {\hat{b}}_1&= - \displaystyle \frac{3c_2 - 2 + 3 {{\hat{c}}}_3 (1- 2 c_2)}{6 (c_1 - {{\hat{c}}}_3 ) ( c_1- c_2) } \end{aligned}$$
(25a)
$$\begin{aligned} {\hat{b}}_3&= -\displaystyle \frac{1}{{{\hat{c}}}_3 ({{\hat{c}}}_3 - c_2)} \left( {\hat{b}}_1 (c_1^2 - c_1 c_2) + \frac{1}{2} c_2 - \frac{1}{3} \right) \end{aligned}$$
(25b)
$$\begin{aligned} {\hat{b}}_2&= \displaystyle \frac{1}{c_2} \left( \frac{1}{2} - {\hat{b}}_3 {{\hat{c}}}_3 - {\hat{b}}_1 c_1 \right) \end{aligned}$$
(25c)
$$\begin{aligned} {{\hat{a}}}_{31}&= \displaystyle \frac{-1}{{\hat{b}}_3 ( c_1 - c_2 ) } \left( {\hat{b}}_1 a_{11} c_1 + {\hat{b}}_2( a_{21} c_1 + a_{22} c_2 ) + {\hat{b}}_3 c_2 {{\hat{c}}}_3 - \frac{1}{6} \right) \end{aligned}$$
(25d)
$$\begin{aligned} {{\hat{a}}}_{32}&= {{\hat{c}}}_3 - {{\hat{a}}}_{31}. \end{aligned}$$
(25e)

For given value \({{\hat{c}}}_3\) the unknowns \( {{\hat{a}}}_{31},\) \({{\hat{a}}}_{32},\) \({\hat{b}}_1, {\hat{b}}_2, {\hat{b}}_3\) in Eq.(22) can be evaluated by means of Eq.(25a) to Eq.(25e), however, subject to the conditions

$$\begin{aligned} 1.~~ {{\hat{c}}}_3 \ne c_1, ~~~~~~~~ 2.~~ {{\hat{c}}}_3 \ne c_2. \end{aligned}$$
(26)

According to (19) the weights for the low order method RK2 can be obtained as

$$\begin{aligned} b_1 = \displaystyle \frac{(1/2-c_2)}{(c_1-c_2)}, \qquad b_2 = 1- b_1, \qquad b_3 = 0. \end{aligned}$$
(27)

Table  2 summarizes some Butcher-Tableaus of third-order ELDIRK methods developed in Mahnken [20, 21]. RK3-Eul being a three-stage ELDIRK-method requires two extra function evaluations to increase the order from \(q= 1\) for the classical Euler methods (backward and forward) to \(r=3\). 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\), RK(3)2-Mid is a three-stage method, however, it requires two extra function evaluations to increase the order from \(q= 2\) to \(r=3\). Recall, that these features hold without solving another nonlinear system at all valid for nonlinear ordinary differential equations of the form (3.1).

Table 1 Summary of Butcher-tableaus for second order explicit and implicit Euler embedded Runge–Kutta-methods
Full size table
Table 2 Summary of Butcher-tableaus for third order implicit embedded Runge–Kutta-methods
Full size table

4 Stability of ELDIRK methods

4.1 Preliminaries

An important feature of any RK method is its stability property. The classical way for the related investigation considers the scalar test equation of Dahlquist, see e.g. Dahlquist [7]

$$\begin{aligned} 1. ~~\dot{y} = \lambda y(t), ~~~~ 2.~~y(0) = 1, ~~~ 3.~~ \text{ Re }(\lambda ) < 0. \end{aligned}$$
(28)

The \(S\)-stage RK-method (6) applied to the test problem (28) yields, Hairer and Wanner [14]

(29a)
(29b)
(29c)
(29d)

Here \({\bar{\tau }} = \tau \lambda \) is a complex variable, is a vector with values one, \({{\underline{b}}}= \left( b_{j}\right) _{j=1}^{S} \) is the vector of weighting factors and \(\underline{{{\underline{A}}}}= \left( a_{ij}\right) _{i,j}^{S}\) is the matrix of RK-coefficients in the Butcher array in Eq.(9). \( U({\bar{\tau }} )\) in Eq.(29b) is a rational function (one polynomial divided by another) of degree \(\le S\) with numerator \( {V({\bar{\tau }} )}\) and denominator \( {W({\bar{\tau }} )}\).

Two important stability criteria are distinguished in the literature, see e.g. Hairer and Wanner [14]

$$\begin{aligned}&\text{ A-Stability }: |U({\bar{\tau }})| \le 1,~~ \forall ~ {\bar{\tau }} \in {\mathbb {C}}^{-}. \end{aligned}$$
(30a)
$$\begin{aligned}&\text{ L-Stability }: \text{ A-Stability } \text{ and }~ |U({\bar{\tau }} )|\rightarrow 0 \text{ as } \text{ Re } ({\bar{\tau }}) \rightarrow -\infty . \end{aligned}$$
(30b)

To be A-stable, the degrees of nominator and denominator in Eq.(29b) must satisfy at least

$$\begin{aligned} \hbox {deg} V\le \hbox {deg} W. \quad \end{aligned}$$
(31)

As noted in Hairer and Wanner [14] condition (30a) alone means stability on the imaginary axis and there it is called I-stability. It is equivalent to the condition for the E-polynomial defined as

$$\begin{aligned}{} & {} E(y):= |W(iy)|^2 -|V(iy)|^2 = W(iy)W(-iy)\nonumber \\{} & {} \quad - V(iy)V(-iy) \ge 0 ~~\forall y \in {\mathbb {R}}. \end{aligned}$$
(32)

For RK methods, which are not A-stable, the definition of the minimal weighted time step is essential

$$\begin{aligned} \begin{array}{llll} {\bar{\tau }}_{\min }\in {\mathbb {R}}^{-} \text{ such } \text{ that } ~~~~~&|U({\bar{\tau }}_{\min })| = 1, \end{array} \end{aligned}$$
(33)

and which renders the critical time step for the Dahlquist problem (28) as

$$\begin{aligned} {\bar{\tau }}_{\min }= \text{ Re }(\lambda ) \tau _{crit}\quad \Longleftrightarrow \quad \tau _{crit}= {\bar{\tau }}_{\min }/ \text{ Re }(\lambda ). \end{aligned}$$
(34)

For general systems of nonlinear differential equations of the form (3) in the requirements (30) and (34), respectively, \(\lambda \) is replaced by the smallest negative eigenvalue \(\lambda _{\min }\) of the generalized eigenvalue problem related to (3) as

$$\begin{aligned} \tilde{\underline{J}} \cdot \underline{x} = \underline{\lambda } \underline{M} \cdot \underline{x}, \qquad \textrm{where} \qquad \tilde{\underline{J}} = \frac{\partial \underline{A}_j[t,\underline{y}\left[ t\right] ]}{\partial \underline{y}}\Big |_{\underline{y}_j}. \end{aligned}$$
(35)
Table 3 Stability functions \(U({\bar{\tau }})\) Eq.(29b) for some implicit and explicit RK-methods: From literature: 1. RK2-Mid/-Trap, 2. SDIRK3, 3. RK4-exp, 4. RK2-Ell. ELDIRKS from Mahnken [20, 21]): 5. RK2-Eul-impl, 6. RK3-Eul-impl, from Table  1; 7. RK3-Trap, 8. RK3-Ell from Table 2. New methods: 9. RK2-stab (\({{\hat{c}}}_2 = 1\)), 10. RK2-stab (\({{\hat{c}}}_2 = 0.8\)), 11. RK3-nstab (\(a_{22}=1/6\)), 12. RK3-stab (\(a_{22}=1\)) from Table 6
Full size table
Table 4 Polynomials \(E(y)\) Eq.(32) for all RK-methods listed in Table 3
Full size table

As noted e.g. in Hairer and Wanner [14], for the implicit midpoint rule RK2-Mid and the trapezoidal rule given in Table 2 we have \( {V(y)}\) \(=\) \( {W(-y)}\) and therefore \(E(y) =0\), meaning that both methods are A-stable. Elaboration of the polynomials \( {V({\bar{\tau }} )}\) and \( {W({\bar{\tau }} )}\) in (29) as well as for \(E(y)\) in Eq.(32) for the additional RK-schemes in Table 2 is somewhat laborious. Therefore this task has been done with the symbolic feature of the software MATLAB.

A summary of the results on stability characteristics is given in Table  3 for the stability functions \(U({\bar{\tau }})\) according to Eq.(29) as well as for the polynomials \(E(y)\) according to Eq.(32) in Table 4, respectively. Corresponding illustrations are provided in Fig. 1. All results for RK2-Mid, RK2-Trap, are identical to those given e.g. in Hairer and Wanner [14]. For comparison (and moreover, for validation of our MATLAB routine) we present results for the implicit third-order SDIRK3 and the classical fourth-order explicit RK-method RK4-exp. The result for the E-polynomial \(E(y)\) for SDIRK3 in Table 4 is not given in Hairer and Wanner [14]. To the authors’ knowledge, also results for RK2-Ell have not been published before. Additionally, Table 5 summarizes the minimal weighted time steps \({\bar{\tau }}_{\min }\) according to Definition (33) for the above-mentioned RK methods.

4.2 Stability of ELDIRK-methods

According to Eq.(15) ELDIRK methods have a strictly lower triangular coefficient matrix \(\underline{{{\underline{A}}}}_{R}\) such that the denominator \(W({\bar{\tau }}) \) in Eq.(29d) becomes

$$\begin{aligned} \begin{array}{llllllll} &{}&{}W^{R}({\bar{\tau }}) = \det ( {I} - {\bar{\tau }} \underline{{{\underline{A}}}}_{R}) = \det \left( {I} - {\bar{\tau }} \left[ \begin{array}{lclll} \underline{{{\underline{A}}}}_{Q}&{} {\underline{0}} \\ \underline{{{\underline{A}}}}_{RQ}&{} \underline{{{\underline{A}}}}_{RR}\\ \end{array} \right] \right) \\ &{}&{}\quad = \det \left[ \begin{array}{ccc|cccc} &{} &{} &{} 0 &{} 0 &{} \ldots &{} 0 \\ &{} {I}^Q -{\bar{\tau }} \underline{{{\underline{A}}}}_{Q}&{} &{} 0 &{} &{} \ldots &{} 0 \\ &{} &{} &{} 0 &{} &{} \ldots &{} 0 \\ \hline \\ -{\bar{\tau }} a_{S+1, 1} &{} \ldots &{} \ldots &{} 1 &{} 0 &{}\ldots &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} &{} \ddots &{} \vdots &{} \\ -{\bar{\tau }} a_{R~1} &{} \ldots &{} \ldots &{} -{\bar{\tau }} a_{R~Q+1} &{} &{} \ldots &{} 1 \end{array}\right] \\ &{}&{}\quad = \det ( {I}^Q - {\bar{\tau }} \underline{{{\underline{A}}}}_{Q}) = W^{Q}({\bar{\tau }}). \end{array} \end{aligned}$$
(36)

This means that the embedded scheme RKq and its ELDIRK extension RKr have the same denominator polynomial, such that its degrees have the property

$$\begin{aligned} \hbox {deg} W^R= \hbox {deg} W^Q. \end{aligned}$$
(37)

The nominator \(V^{R}({\bar{\tau }}) \) in Eq.(29c) of the ELDIRK extension

(38)

also is a polynomial. In general, it is different to \(V^{Q}({\bar{\tau }})\) of the embedded scheme RKq. Let us assume, that the embedded scheme RKq is A-stable, the unfavorable case of equality \(\hbox {deg} V^Q= \hbox {deg} W^Q\) in Eq.(31) holds, and the case of practical interest \( V^{R}({\bar{\tau }}) \ge V^{Q}({\bar{\tau }})\) holds, then in view of Eq.(37) it follows

$$\begin{aligned} \hbox {deg} V^R\ge \hbox {deg} V^Q= \hbox {deg} W^Q= \hbox {deg} W^R. \quad \quad \end{aligned}$$
(39)

Consequently, their stability function \(U^{R}\) is a polynomial of degree \(\hbox {deg} U^R\ge 1\), which proves

$$\begin{aligned} \text{ Proposition: } \text{ In } \text{ general } \text{ an } \text{ ELDIRK } \text{ method } \text{ is } \text{ not } \text{ A-stable }. \end{aligned}$$
(40)

A summary of results on stability characteristics of ELDIRK methods RK2-Eul-impl, RK3-Eul-impl, RK3-Trap and RK3-Ell according to Table 1 and Table 2 is listed in Table  3 and Table 4, respectively. As before, we provide results for the stability function \(U({\bar{\tau }})\) and the polynomial \(E(y)\), respectively, which are illustrated in Fig. 1. These results verify Proposition (40), namely, that ELDIRK methods in general are not A-stable. Additionally, Table 5 summarizes the minimal weighted time steps \({\bar{\tau }}_{\min }\) according to Definition (33) for the RK-methods.

Please, note, that the above proposition (40) holds for general ELDIRK methods. To be A-stable, the embedded scheme RKq and its ELDIRK extension RKr must have nominator polynomials with the same degree

$$\begin{aligned} \hbox {deg} V^R= \hbox {deg} V^Q. \end{aligned}$$
(41)

In the ensuing subsections, it will be shown, that indeed for certain cases A-stable ELDIRK-methods can be constructed.

4.3 Construction of second order A-stable ELDIRKs (RK2-stab)

We consider the two-stage ELDIRK method RK2 according to Eq.(17). Its denominator polynomial \(W^{R}({\bar{\tau }})\) in Eq.(36) is

$$\begin{aligned} \begin{array}{lllllll} W^{R}({\bar{\tau }}) = W^{Q}({\bar{\tau }})= & {} 1 - a_{11} {\bar{\tau }}. \end{array} \end{aligned}$$
(42)

The nominator polynomial \(V^{R}({\bar{\tau }})\) in Eq.(38) (obtained from our MATLAB symbolic routine) is

$$\begin{aligned} \begin{array}{lllllll} V^{R}({\bar{\tau }}) &{}=&{} 1+ {\bar{\tau }}( {{\hat{b}}}_{1} + {{\hat{b}}}_{2} - a_{11} ) + {\bar{\tau }}^2{{\hat{b}}}_{2} ( {{\hat{a}}}_{21} - a_{11} ) \\ &{}=&{} V^{R}({\bar{\tau }}, {{\hat{b}}}_i, a_{ij}, {{\hat{a}}}_{ij}). \end{array} \end{aligned}$$
(43)

Next, we exploit the stage condition (20a) and the order condition (20c) as well as Eq.(20b) to obtain

$$\begin{aligned} \begin{array}{lllllll} V^{R}({\bar{\tau }}) \ {}= & {} 1+ {\bar{\tau }}( 1 - a_{11} ) + {\bar{\tau }}^2 \left( \frac{ 1}{ 2}- a_{11}\right) = V^{R}({\bar{\tau }},a_{11}). \end{array}\nonumber \\ \end{aligned}$$
(44)

Interestingly, this representation shows no dependence of \( V^{R}\) on the additional coefficients \( {{\hat{b}}}_i, {{\hat{a}}}_{ij}\) of the ELDIRK method. Therefore, by recalling Eq.(29b) and Eq.(42) we conclude

$$\begin{aligned} \begin{array}{lllll} \text{ Proposition: } \text{ The } \text{ stability } \text{ function } U^{R}= \displaystyle {V^{R}}/{W^{R}}\hbox { of }\\ \hspace{-1mm}{{RK2}} \hbox { can be formulated solely in terms of the} \\ \hspace{-1mm} \text{ RK-coefficient } a_{11} \text{ of } \text{ the } {{RK1}} \text{ scheme. } \end{array} \end{aligned}$$
(45)
Fig. 1
figure 1

Stability characteristics of RK-methods: Left, RK2-Mid/-Trap, SDIRK3, RK4-exp, RK2-Ell, taken from literature. Right, second order ELDIRK RK2-Eul-impl from Table 1 and third order ELDIRKs RK3-Eul-impl, RK3-Trap, RK3-Ell from Table 2

Full size image
Table 5 Minimal weighted time step \({\bar{\tau }}_{\min }\) according to Definition (33) for all RK-methods listed in Table 3
Full size table

In order to construct an A-stable method we need according to Eq.(41) \(\hbox {deg} V^R= \hbox {deg} V^Q= 1\) and therefore from Eq.(44) we can conclude \(a_{11} = 1/2\). Moreover, setting \(b_1 = 1\) for B(2)1 in Eq.(17) would render the second-order Midpoint rule RK2-Mid in Table 2. However, we are interested in a first-order RK1 for an embedded method with its second-order ELDIRK extension. Therefore, in contrast to Eq.(17) we choose two values \(b_1 = b_{2} = 1/2\) as well as \({{\hat{c}}}_{2} = 1\), \({{\hat{b}}}_{2} = 1\) to obtain the extended Butcher array in Eq.(17) labeled as RK2-stab(\({{\hat{c}}}_2=1\)) in Table 6. Analogously, we choose (arbitrarily) \({{\hat{c}}}_2=0.8\), to obtain a second scheme labeled as RK2-stab(\({{\hat{c}}}_2=0.8\)) in Table 6.

A summary of results on stability characteristics of the new ELDIRK methods namely, RK2-stab(\({{\hat{c}}}_2 = 1\)) and RK2-stab(\({{\hat{c}}}_2 = 0.8\)) according to Table 6, is listed in Table  3 and Table 4, respectively. As before, we provide results for the stability function \(U({\bar{\tau }})\) and the polynomial \(E(y)\), respectively, which are illustrated in Fig. 2. Additionally, Table 5 summarizes the minimal weighted time steps \({\bar{\tau }}_{\min }\) according to Definition (33) for the new RK methods.

Interestingly, analogously to Eq.(44), the stability functions \(U({\bar{\tau }})\) and the polynomials \(E(y)\), for both, RK2-stab\(({\hat{c}}_2=0.8)\) and RK2-stab\(({\hat{c}}_2=1)\), are the same as for the midpoint rule and the trapezoidal rule, respectively. The results verify, that the stability characteristics of RK2-stab are indeed independent of the choice for \({{\hat{c}}}_2\), which is a statement of Proposition (45). In summary, all results verify the properties (30a) or (32), respectively, for A-stability, of the new ELDIRK methods RK2-stab although they are explicit in the second stage. Accordingly, Table 5 verifies, that the minimal weighted time steps defined in (33) obtain \({\bar{\tau }}_{\min }\rightarrow -\infty \) for both methods.

Table 6 Summary of Butcher-tableaus for new ELDIRK methods RK(2)1-stab(\({{\hat{c}}}_2=1\)), RK(2)1-stab(\({{\hat{c}}}_2=0.8\)), RK(3)2-nstab(\(a_{22}=1/6\)), RK(3)2-stab(\(a_{22}=1\))
Full size table
Fig. 2
figure 2

Stability characteristics of new ELDIRK-methods in Table 6: Left, second order methods RK2-stab(\({{\hat{c}}}_2\)=1) and RK2-stab(\({{\hat{c}}}_2\)=0.8). Right, third order methods, RK3-nstab(\(a_{22}\) =1/6) and RK3-stab(\(a_{22}\) = 1). All methods show the A-stability characteristic, execpt RK3-nstab(\(a_{22}\) = 1/6), since it violates condition (52)

Full size image

4.4 Construction of third order A-stable ELDIRKs (RK3-stab)

We consider the three-stage ELDIRK method RK3 according to Eq.(22). Its denominator polynomial \(W^{R}({\bar{\tau }})\) in Eq.(36) is

$$\begin{aligned} W^{R}({\bar{\tau }})= & {} (a_{11} {\bar{\tau }}- 1) (a_{22} {\bar{\tau }}- 1)\nonumber \\ {}= & {} 1 - {\bar{\tau }}(a_{11} + a_{22}) + {\bar{\tau }}^2{a_{11} a_{22}}. \end{aligned}$$
(46)

Its nominator polynomial \(V^{R}\) in Eq.(38) is

$$\begin{aligned} V^{R}({\bar{\tau }})= & {} {{\hat{b}}}_{1} {\bar{\tau }}- a_{22} {\bar{\tau }}- a_{11} {\bar{\tau }}+ {{\hat{b}}}_{2} {\bar{\tau }}+ {{\hat{b}}}_{3} {\bar{\tau }}+ a_{11} a_{22} {\bar{\tau }}^2\nonumber \\{} & {} - a_{11} {{\hat{b}}}_{2} {\bar{\tau }}^2 - a_{11} {{\hat{b}}}_{3} {\bar{\tau }}^2 + a_{21} {{\hat{b}}}_{2} {\bar{\tau }}^2 - a_{22} {{\hat{b}}}_{1} {\bar{\tau }}^2 \nonumber \\{} & {} - a_{22} {{\hat{b}}}_{3} {\bar{\tau }}^2 + a_{31}l {{\hat{b}}}_{3} {\bar{\tau }}^2 + {{\hat{a}}}_{32} {{\hat{b}}}_{3} {\bar{\tau }}^2 + a_{11} a_{22} {{\hat{b}}}_{3} {\bar{\tau }}^3\nonumber \\{} & {} - a_{11} {{\hat{a}}}_{32} {{\hat{b}}}_{3} {\bar{\tau }}^3 + a_{21} {{\hat{a}}}_{32} {{\hat{b}}}_{3} {\bar{\tau }}^3 - a_{22} a_{31}l {{\hat{b}}}_{3} {\bar{\tau }}^3 + 1\nonumber \\= & {} V^{R}({\bar{\tau }}, {{\hat{b}}}_i, a_{ij}, {{\hat{a}}}_{ij}). \end{aligned}$$
(47)

Next, we exploit the order conditions (25c) to (25d) and obtain after some algebra:

$$\begin{aligned} \begin{array}{lllllll} V^{R}({\bar{\tau }}) &{}=&{} 1 + {\bar{\tau }}(1 - a_{11} - a_{22} ) + {\bar{\tau }}^2\left( \frac{ 1}{ 2}- a_{11} + a_{22} + a_{11} a_{22} \right) \\ &{}&{} + {\bar{\tau }}^3\left( \displaystyle \frac{1}{6} - \frac{ 1}{ 2}a_{11} - \frac{ 1}{ 2}a_{22} + a_{11} a_{22} \right) = V^{R}({\bar{\tau }}, a_{11},a_{22}). \end{array} \end{aligned}$$
(48)

Interestingly, the result shows no dependence of \( V^{R}({\bar{\tau }}) \) on the additional coefficients \( {{\hat{b}}}_i, {{\hat{a}}}_{ij}\) of the ELDIRK method. Therefore, by recalling Eq.(29b) and Eq. (46) we conclude

$$\begin{aligned} \begin{array}{lllll} \text{ Proposition: } \text{ The } \text{ stability } \text{ function } U^{R}= \displaystyle {V^{R}}/{W^{R}} \text{ of } {{RK3}}\\ \text{ can } \text{ be } \text{ formulated } \text{ solely } \text{ in } \text{ terms } \text{ of } \text{ the } \\ \text{ diagonal } \text{ RK-coefficients } a_{11}, a_{22} \text{ of } \text{ the } {{RK2}}\text{-scheme }. \end{array} \end{aligned}$$
(49)

In order to construct an A-stable method we need according to Eq.(41) \(\hbox {deg} V^R= \hbox {deg} V^Q= 2\). From the third-order term in Eq.(48) we obtain:

$$\begin{aligned} \displaystyle \frac{1}{6} - \frac{ 1}{ 2}a_{11} - \frac{ 1}{ 2}a_{22} + a_{11} a_{22} = 0 ~~\Longleftrightarrow ~~ a_{11} = \displaystyle \frac{1/3 - a_{22}}{1-2 a_{22}}.\nonumber \\ \end{aligned}$$
(50)

After some algebra, this leaves Eq.(46) and Eq. (48) as

$$\begin{aligned} W^{R}({\bar{\tau }})&= \displaystyle \frac{-1}{{6(2 a_{22} - 1)}} ( \underbrace{6 - 12 a_{22}}_{f_{W0}} + {\bar{\tau }}\underbrace{(12 a_{22}^2 -2 )}_{f_{W1}} \nonumber \\&\quad + {\bar{\tau }}^2 \underbrace{(2a_{22} - 6 a_{22}^2)}_{f_{W2}} ) \end{aligned}$$
(51a)
$$\begin{aligned} V^{R}({\bar{\tau }})&= \displaystyle \frac{-1}{6(2 a_{22} - 1)} (\underbrace{6 - 12 a_{22}}_{f_{V0}} + {\bar{\tau }}\underbrace{( 4 + 12 a_{22}^2 - 12 a_{22})}_{f_{V1}} \nonumber \\&\quad + {\bar{\tau }}^2 \underbrace{(1 + 6 a_{22}^2 - 4 a_{22} )}_{f_{V2}} ). \end{aligned}$$
(51b)

The condition for Eq.(29) on A-Stability in Eq.(30) renders a requirement for the second order \({\bar{\tau }}^2\) terms \({f_{W2}}\) and \({f_{V2}}\) in Eq.(51). Then, after some further algebra finally we obtain the result:

$$\begin{aligned} \begin{array}{llllllll} \lim \limits _{\bar{\tau }\rightarrow -\infty }|U({\bar{\tau }})| = \lim \limits _{\bar{\tau }\rightarrow -\infty } \left| \displaystyle \frac{V({\bar{\tau }} )}{W({\bar{\tau }} )}\right| = \displaystyle \frac{|f_{V2}|}{|f_{W2}|} \le 1 \quad \Longrightarrow \quad a_{22} \ge \frac{ 1}{ 2}. \end{array} \nonumber \\ \end{aligned}$$
(52)

From the equations (51), we see, that the case \( a_{22} = 1/2\) is not defined. The extended Butcher array is summarized as RK3-stab in Table 6, where in view of the result (52) \(a_{22}=1 \) has been chosen. For later comparison, also the (arbitrarily chosen) case \(a_{22}= {1}/{6}< {1}/{2}\) has been introduced.

A summary of results on stability characteristics of the new third-order ELDIRK methods RK3-stab(\(a_{22} = 1.0\)), and RK3-nstab(\(a_{22} = 1/6\)) according to Table 6 is listed in Table 3 and Table 4, respectively. As before, we provide results for the stability function \(U({\bar{\tau }})\) and the polynomial \(E(y)\), respectively, which are illustrated in Fig. 2. The results verify the properties (30a) or (32), respectively, for A-stability, of the new ELDIRK method RK3-stab(\(a_{22} = 1.0\)) although it is explicit in the third stage. Contrary RK3-nstab(\(a_{22} = 1/6\)) is not A-stable, which is due to a violation of the condition \(a_{22} \ge 1/2\) in Eq.(52). Accordingly, Table 5 shows, that the minimal weighted time steps \({\bar{\tau }}_{\min }\) according to Definition (33) is limited only for RK3-nstab(\(a_{22} = 1/6\)).

5 Representative numerical examples

In the following, numerical results of the embedded Ellsiepen method RK(3)2-Ell in Table  2, and all four new embedded schemes of the classes RK(2)1 and RK(3)2 summarized in Table  6, are presented for two test problems: Firstly, curing for a thermosetting material and, secondly, phase field RVE-modeling for crystallinity and orientation. Results are presented for equidistant step lengths as well as for non-equidistant step lengths based on adaptive time-step control, where as mentioned before, only reversed embedded RK(q)r methods according to Eq.(10) will be considered. Regarding the following nomenclature, occationally also results for equidistant step lengths (without adaptive time-step control) are labeled as RK(q)r, instead of RKr.

As an error measure, we use the quantity of interest for the global error in Eq.(8)

(53)

As an alternative to the above time point goal function, results for a time-integrated goal function are presented in Mahnken [20]. Evaluation of Eq. (53) requires an analytical solution which is not available for both test problems. Therefore, in both cases a highly accurate solution, subsequently denoted as,,reference solution”, is generated with the classical 4th-order RK-method having Butcher-Tableau

$$\begin{aligned} {{\underline{B}}}= \begin{array}{ll|llllllll} 0 &{} &{} &{} &{}&{} &{}&{} &{} \\ \frac{1}{2} &{} &{} &{} \frac{1}{2} &{}&{} &{}&{} \\ \frac{1}{2} &{} &{} &{} 0 &{}&{} \frac{1}{2} &{}&{} \\ 1 &{} &{} &{} 0 &{}&{} 0 &{}&{} 1\\ \hline &{} &{} &{} &{}&{} &{}&{} \\ &{} &{} &{} \frac{1}{6} &{}&{} \frac{1}{3} &{}&{} \frac{1}{3} &{}&{} \frac{1}{6} \end{array}. \end{aligned}$$
(54)

5.1 Curing for a thermosetting material

The first numerical example is concerned with the evolution of curing for a thermosetting material with IVP as, see Lion and Höfer [18]

$$\begin{aligned} \begin{array}{lllllllllll} &{}1.\quad {{\dot{z}}} = A(z), \quad \qquad \textrm{where}\quad \quad 2.\quad A(z) = K z^a (1-z)^n, \\ {} &{}3.\quad K = k_A \exp \left[ \displaystyle \frac{-E_A}{R\theta }\right] , \quad 4.\quad z(0) = z_0. \end{array} \end{aligned}$$
(55)

Here \(k_A,E_A,R\) are respectively the Arrhenius constant, the activation energy and the universal gas constant. Furthermore, a and n are power constants. In order to restrict the investigations in this paper purely to discretization errors rather than model errors, only numerical results are presented. A value \(z_0 = 10^{-3}\) in Eq.(55.4) is chosen as an initial condition in order to trigger curing. A discussion on the influence of the initial value is not in the scope of this paper. The material constants in (55) are chosen as: \(k_A\) \(=\) 3.0e9 [-], \(E_A\) \(=\) 89110 J/mol, R \(=\) 8.3144621, J/(mol K), \(\theta \) \(=\) 410 K, a = 1.2 [-], n = 3.01 [-]. The initial time is \(t_0 = \alpha = 0\) and the total time is \(T = \beta -\alpha =12000\) s. The iterative solution of the resulting nonlinear residual by means of Newtons method is described in Mahnken [20]. On solution procedures for an IVBP with related boundary conditions, the interested reader is referred to Mahnken [19]. Since an analytical solution is not available, a,,reference solution” is generated by means of the classical 4th-order RK-method with Butcher tableau in Eq. (54) in \(N=10^5\) equidistant time steps with length \(\tau _i = 0.12\) s.

Figure 3 shows results for all five reversed embedded RK-methods RK(q)r for \(N= 200\) equidistant time steps. For the curing variable versus time z(t) in Fig. 3.a) a difference to the exact solution is hardly visible for all RK-methods.

Table 7 Curing for a thermosetting material for time point goal function in Eq. (53): Results for five low-order methods RKq and its extensions RKr with equidistant step sizes. Convergence orders for the global error in Eq.(8)
Full size table

Table 7 summarizes the convergence orders according to Eq. (8) for all low-order methods RKq and its extensions RKr with equidistant step sizes. A visualization of the results for RKr in provided in Fig. 3.b). It can be seen, that all RK methods reach their theoretical convergence orders.

Fig. 3
figure 3

Curing for a thermosetting material: Results for five RK(q)r methods with equidistant step sizes. a) Curing versus time, b) global error in Eq.(53)

Full size image

5.2 A Kobayashi-Warren-Carter phase-field model capturing grain boundary migration

The Kobayashi-Warren-Carter model for grain boundary migration in Abrivard, Busso, Forest and Appolaire [1] and Abrivard, Busso, Forest and Appolaire [2] is used to demonstrate the stability properties and accuracy in time of the novel A-stable ELDIRK time integration schemes with the finite-element method.

5.2.1 Evolutions equations

The evolution equations for the order parameters of the crystallinity \(\eta \) and the orientation \(\theta \) are

$$\begin{aligned}{} & {} \tau _\eta \frac{\partial \eta }{\partial t} = \left( \alpha ^2\nabla ^2\eta -\frac{\partial f\left[ \eta \right] }{\partial \eta } -\frac{\partial g\left[ \eta \right] }{\partial \eta }|\nabla \theta | -\frac{\partial h\left[ \eta \right] }{\partial \eta }|\nabla \theta |^2\right) , \nonumber \\{} & {} \tau _\theta \frac{\partial \theta }{\partial t} = \nabla \cdot \left( 2h\left[ \eta \right] \nabla \theta +g\left[ \eta \right] \frac{\nabla \theta }{|\nabla \theta |}c \right) , \end{aligned}$$
(56)

where the functions \(f\left[ \eta \right] \), \(g\left[ \eta \right] \) and \(h\left[ \eta \right] \) are given as

$$\begin{aligned}{} & {} f\left[ \eta \right] = \frac{\omega ^2}{2}\left( 1-\eta \right) ^2, \qquad g\left[ \eta \right] = a_1\eta +a_2\eta ^2+a_3\eta ^3, \nonumber \\{} & {} h\left[ \eta \right] = d \eta ^2 + e, \end{aligned}$$
(57)

and \(\tau _\eta \) as well as \(\tau _\theta \), are positive inverse mobilities, related to the grain boundary mobility in a non-trivial manner

$$\begin{aligned} \tau _\eta = \beta _\eta \left[ T\right] \omega ^2, \qquad \tau _\theta [\eta ] = \beta _\theta \left[ T\right] \omega ^2. \end{aligned}$$
(58)

Here, \(\omega \) is associated with the grain boundary energy. The temperature-dependent parameters \(\beta _\eta \left[ T\right] \) and \(\beta _\theta \left[ T\right] \) define the time scales for the evolutions of crystallinity and orientation and read

$$\begin{aligned} \beta _\eta \left[ T\right] = \beta ^0_\eta \exp \left[ \frac{Q}{kT}\right] , \qquad \beta _\theta \left[ T\right] = \beta ^0_\theta \exp \left[ \frac{Q}{kT}\right] , \end{aligned}$$
(59)

where an Arrhenius law captures the temperature influence with the activation energy Q and the universal gas constant k.

The singularity of the orientation gradient \(\nabla \theta \) leads to a discontinuous partial differential equation. The discontinuity is resolved by introducing the scaling function

$$\begin{aligned} c = \tanh \left[ \gamma |\nabla \theta |\right] . \end{aligned}$$
(60)

5.2.2 Variational formulation, temporal discretization and Gâteaux derivatives

The derivation of a variational formulation of the evolution equations (56) is explained in Westermann and Mahnken [23]. Here, the non-linear vector-valued function in Eq. (3) is filled with the order parameters in Eq. (56) and yields a pair of primal variables

$$\begin{aligned} \underline{y} = \left[ \eta ,\theta \right] . \end{aligned}$$
(61)

The variational formulations for the crystallinity \(\eta \) and orientation \(\theta \) in Eq. (56) with the corresponding linear and semilinear forms are written in chevron notation of Eq.(1) by means of test functions \(\delta \underline{y}\)

(62)

and in compact form reads

$$\begin{aligned}<\delta \underline{y}\cdot \underline{M}\cdot \dot{\underline{y}}>+<\delta \underline{y}\cdot \underline{A}\big [\underline{y}\big ]> = <\delta \underline{y}\cdot \underline{f}>, \quad \forall \delta \underline{y}.\nonumber \\ \end{aligned}$$
(63)

For a discretization in time the S-stage Runge–Kutta formulation in Eq. (6) applied to Eq. (63) becomes

$$\begin{aligned}{} & {} <\delta \underline{y}\cdot \underline{M}\cdot \underline{k}_j>+<\delta \underline{y}\cdot \underline{A}\big [\underline{y}_j\big [\underline{k}_j\big ]\big ]> \nonumber \\{} & {} \quad = <\delta \underline{y}\cdot \underline{f}_{ij}>, \quad \underline{y}_j= \underline{y}_{i-1}+\tau _i\sum _{l=1}^S a_{jl}\underline{k}_l, \quad j,l = 1,\dots ,S.\nonumber \\ \end{aligned}$$
(64)

Eq. (64) may be used to define a residual where the coefficients \(\underline{k}_{j}\) of the generalized S-stage RK-method are unknown:

$$\begin{aligned}{} & {} \underline{\rho }_{ij}\big [\underline{k}_{j};\delta \underline{y}\big ]=<\delta \underline{y}\cdot \underline{M}\cdot \underline{k}_j>+<\delta \underline{y}\cdot \underline{A}\big [\underline{y}_j\big [\underline{k}_j\big ]\big ]>\nonumber \\{} & {} \quad -<\delta \underline{y}\cdot \underline{f}_{ij}>=\underline{0},\qquad j,l=1,\dots ,S. \end{aligned}$$
(65)

This non-linear residual equation is solved with a Newton method by means of the Gâteaux derivative in Eq.(2) applied to the residual \(\underline{\rho }_{ij} \big [\underline{k}_{j};\delta \underline{y}\big ]\) in Eq.(65)

$$\begin{aligned}{} & {} D_{\underline{k}_j} \underline{\rho }_{ij}\big [\underline{k}_j;\Delta \underline{k}_j,\delta \underline{y}\big ] = D_{\underline{k}_j} \underline{\rho }_{ij}\big [\underline{k}_j;\delta \underline{y}\big ]\cdot \Delta \underline{k}_j \nonumber \\{} & {} \quad =<\delta \underline{y}\cdot \underline{M}\cdot \Delta \underline{k}_j> +<\delta \underline{y}\cdot \frac{\partial \underline{A}}{\partial \underline{k}_j}\cdot \Delta \underline{k}_j> \nonumber \\{} & {} \quad =<\delta \underline{y}\cdot \underline{M}\cdot \Delta \underline{k}>+<\delta \underline{y}\cdot \frac{\partial \underline{A}}{\partial \underline{y}_j} \cdot \frac{\partial \underline{y}_j}{\partial \underline{k}_j}\cdot \Delta \underline{k}_l>\nonumber \\{} & {} \quad =<\delta \underline{y}\cdot \underline{M}\cdot \Delta \underline{k}>+<\delta \underline{y}\cdot \frac{\partial \underline{A}}{\partial \underline{y}_j} \tau _ia_{jl}\cdot \Delta \underline{k}_l>. \end{aligned}$$
(66)

The discretization in space as well as the derivatives of the Gâteaux concept are found in Appendix A and Westermann and Mahnken [23], respetively.

5.2.3 Convergence investigations for equidistant step sizes

The variational form in Eq. (62) is applied to meaningful examples henceforth. The Kobayashi-Warren-Carter model is chosen based on its severe need for fine discretizations in time. A two-dimensional squared representative volume element (RVE) with a side length a of 1 mm replaces the microscopic domain \(\Omega \) introduced in the IVP in Eq. (3). A straightforward bicrystal example is used to illustrate the Kobayashi-Warren-Carter approach’s outcomes. In this case, periodic boundary conditions are used to account for the periodicity of crystalline microstructures. These constraints are implemented as primary/secondary set Dirichlet boundary conditions. The 256x256 quadrilateral elements of the spatial discretization have linear shape functions. The source term is eliminated as the microstructure changes across the entire time \(T=1\,\hbox {s}\), driven exclusively by temperature activation.

Fig. 4
figure 4

Bicrystal example - Constitutive variables for four time increments: a–d Orientation \(\theta \), e–h crystallinity \(\eta \)

Full size image

The example takes into account a small grain (\(\theta _1=-0.8\)) that is enclosed by a larger grain (\(\theta _1=0.8\)), which is distinguished by the orientation (\(\theta \)). The grains’ orientation is decided upon at random. Grain boundaries arise between two grains, which are identified by the crystallinity (\(\eta \gg 0\)), according to the model equations in Eq. (56). The initial values are displayed in Fig. 4.a and Fig. 4.e for the orientation \(\theta \) and crystallinity \(\eta \), respectively. The orientation gradient (\(\nabla \theta \)) between both grains causes the inner grain to rotate until eventual consumption by the outer grain. Thus, the orientation becomes the same across the domain \(\omega \) and the orientation gradient (\(\nabla \theta \)) vanishes. Since the crystallinity displays grain boundaries (\(\eta \gg 0\)) only for \(\nabla \theta \gg 0\), the grain boundary vanishes as well. After the consumption of the inner grain, the microstructure reaches a steady state.

The convergence orders for various Butcher tableaus are obtained for the equidistant step sizes

$$\begin{aligned} \underline{\tau } = \left[ {}{1e-4}{}, {}{5e-5}{}, {}{2e-5}{}, {}{1e-5}{}, {}{1e-6}{}\right] .\nonumber \\ \end{aligned}$$
(67)

Only 0.01 s of the simulation in Fig. 4 is taken into consideration for determining the convergence orders due to the significant numerical effort. The reference solution is produced via the traditional fourth-order explicit Runge–Kutta technique RK4-Class in Eq. (54) with N=100000 equidistant time increments in the current situation because analytical solutions are not attainable. On a logarithmic scale, the global temporal error \(\underline{E}\) is shown in comparison to the equidistant step sizes \(\underline{\tau }\) in the Fig. 5. In the legend, the convergence orders are shown next to the corresponding procedure.

Fig. 5
figure 5

Equidistant step size - Global error \(\underline{E}\) over equidistant time steps \(\tau \)

Full size image

Even for the non-linear partial differential equations of the Kobayashi-Warren-Carter problem, all four new methods summarized in Table  6, RK(1)2-stab(\({{\hat{c}}}_2=1\)), RK(1)2-stab(\({{\hat{c}}}_2=0.8\)), RK(2)3-nstab(\(a_{22}=1/6\)), and RK(2)3-stab(\(a_{22}=1\)) reveal the respective increased orders, determined by linear regression, with only one additional explicit evaluation in the last RK-stage.

5.2.4 Stability investigations for equidistant step sizes

The stability of the new ELDIRK methods is of interest, to determine A-stability. In Westermann and Mahnken [23], the explicit nature of the last stage requires adherence to sufficiently small step sizes. Instability problems result from the growth of the Kobayashi-Warren-Carter model for insufficiently small step sizes. ELDIRK methods may result in oscillation build-up, as shown in Fig. 6. In this case, RK(2)3-nstab(\({{\hat{a}}}_{22}=1/6\)) is used with N=250 . evenly spaced steps over \(T=1\,\hbox {s}\). Both unknowns start to oscillate in time and space, until simulation termination.

Fig. 6
figure 6

Stability investigation - RK(2)3-nstab(\({{\hat{a}}}_{22}=1/6\)): Order parameter for four time increments: a–d Orientation \(\theta \), e–h crystallinity \(\eta \)

Full size image

The potential stability problems of the ELDIRK methods are quantified in Fig. 7, Fig. 8, and Fig. 9, which show the results for RK(2)3-Ell, RK(2)3-nstab(\(a_{22}=1/6\)), and RK(2)3-stab(\(a_{22}=1\)), respectively. Each contains the order parameters \(\eta \) and \(\theta \) in a, the minimum eigenvalue \(\lambda _{\textrm{min}}\) in b, the step size \(\overline{\tau }\) in c, and the stability function \(U(\overline{\tau })\) in d for two different step sizes. Some step sizes are chosen to invoke stability problems. In Fig. 7 the results for the RK(2)3-Ell method proposed in Mahnken [21] are shown. Here, stability problems are observed for \(\tau =5\times 10^{-4}\,\hbox {s}\) but not for \(\tau =2\times 10^{-4}\,\hbox {s}\).

Fig. 7
figure 7

Stability investigation - RK3-Ell: Crystallinity \(\eta \) and orientation \(\theta \), minimum eigenvalue \(\lambda _{\textrm{min}}\), step size \(\overline{\tau }\), stability function \(U(\overline{\tau })\): a–d \(\tau =5\times 10^{-4}\,\hbox {s}\), e–h \(\tau =2\times 10^{-4}\,\hbox {s}\)

Full size image

Both order parameters \(\eta \) and \(\theta \) show oscillatory behavior. By solving the eigenvalue problem in Eq. (35), the minimal eigenvalues \(\lambda _{\textrm{min}}\) are determined. Numerical stability is evaluated based on the stability function \(U(\overline{\tau })\) following Eq. (29). In Fig. 7.d the function \(U(\overline{\tau })>1\) explains the occurrence of oscillatory behavior and in Fig. 7.h \(U(\overline{\tau })<1\) explains the absence of oscillatory behavior.

Fig. 8
figure 8

Stability investigation - RK3-nstab(\(a_{22}=1/6\)): Crystallinity \(\eta \) and orientation \(\theta \), minimum eigenvalue \(\lambda _{\textrm{min}}\), step size \(\overline{\tau }\), stability function \(U(\overline{\tau })\): a–d \(\tau =1\times 10^{-3}\,\hbox {s}\), e–h \(\tau =2.5\times 10^{-4}\,\hbox {s}\)

Full size image

Similar to RK(2)3-Ell, numerical problems arise for RK(2)3-nstab(\(a_{22}=1/6\)). Figure 8 shows the results for step sizes \(\tau =1\times 10^{-3}\,\hbox {s}\) and \(\tau =2.5\times 10^{-3}\,\hbox {s}\). Since the critical step size \(\tau _{crit}\) for RK(2)3-nstab(\(a_{22}=1/6\)) shows greater stability than the critical step size for RK(2)3-Ell, larger step sizes are possible.

In contrast to RK(2)3-nstab(\(a_{22}=1/6\)), RK(2)3-stab(\(a_{22}=1\)), shows A-stable behavior for the same step sizes. The choices of \(\tau =1\times 10^{-3}\,\hbox {s}\) in Fig. 9.a-d as well as \(\tau =2.5\times 10^{-4}\,\hbox {s}\) in Fig. 9.e-h both yield stable results.

Fig. 9
figure 9

Stability investigation - RK3-stab(\(a_{22}=1\)): Crystallinity \(\eta \) and orientation \(\theta \), minimum eigenvalue \(\lambda _{\textrm{min}}\), step size \(\overline{\tau }\), stability function \(U(\overline{\tau })\): a–d \(\tau =1\times 10^{-3}\,\hbox {s}\), e–h \(\tau =2.5\times 10^{-4}\,\hbox {s}\)

Full size image

5.2.5 Stability-controlled adaptively refined step sizes

The bicrystal example in Fig. 4 is applied to adaptive step size refinement with stability-controlled criteria. The increment size for the next time increment is determined by the critical step size \(\tau _{crit}\) in Table 5 for the respective integration scheme as well as the minimum eigenvalue of the Jacobian.

Fig. 10
figure 10

Stability-controlled adaptivity - RK(2)3-Ell: a Order parameters y, b minimum eigenvalue \(\lambda _{\textrm{min}}\), c step size \(\tau \), d stability function \(U(\overline{\tau })\)

Full size image

In Fig. 10 the adaptively refined step size \(\tau \) for RK(2)3-Ell is displayed. As expected, the stability function \(U(\overline{\tau })\) in Fig. 10.d is constant to 1, leading to numerical stability, which results in an adaptively refined step size \(\tau \) in Fig. 10.c. In Fig. 11 the adaptively refined step size \(\tau \) for RK(2)3-nstab(\(a_{22}=1/6\)) is displayed. As expected, the stability function \(U(\overline{\tau })\) in Fig. 11.d is constant to 1, leading to numerical stability, which results in an adaptively refined step size \(\tau \) in Fig. 11.c.

Fig. 11
figure 11

Stability-controlled adaptivity - RK(2)3-nstab(\(a_{22}=1/6\)): a Order parameters y, b minimum eigenvalue \(\lambda _{\textrm{min}}\), c step size \(\tau \), d stability function \(U(\overline{\tau })\)

Full size image
Fig. 12
figure 12

Adaptively refined step size - Second-order methods (RK(1)2-stab(\({{\hat{c}}}_{2}=1\)), RK(1)2-stab(\({{\hat{c}}}_{2}=0.8\))): Norm of order parameters \(||\underline{y}||_2\), step size \(\tau \), estimated error \(\underline{{\hat{e}}}\): a–c \(\varepsilon ={}{1e-2}{}\), d–f \(\varepsilon ={}{1e-4}{}\), g–i \(\varepsilon ={}{1e-6}{}\)

Full size image
Fig. 13
figure 13

Adaptively refined step size - Third-order methods (RK(2)3-Ell, RK(2)3-nstab(\(a_{22}=1/6\)), RK(2)3-stab(\(a_{22}=1\))): Norm of order parameters \(||\underline{y}||_2\), step size \(\tau \), estimated error \(\underline{{\hat{e}}}\): a–c \(\varepsilon ={}{1e-6}{}\), d–f \(\varepsilon ={}{1e-8}{}\), g–i \(\varepsilon ={}{1e-10}{}\)

Full size image

5.2.6 Error-controlled adaptively refined step sizes

The ELDIRK schemes are now used for adaptive discretization in time with error-controlled methods. The results of adaptive step size selection based on a posteriori error estimation with embedded RK methods are presented for different maximum allowable errors in Fig. 12 for second-order methods RK(1)2-stab(\({{\hat{c}}}_2=1\)) and RK(1)2-stab(\({{\hat{c}}}_2=0.8\)). Here, the L2-norm of the order parameters \(||\underline{y}||_2\), the adaptively refined step size \(\tau \), and the estimated error \(\hat{\underline{e}}\) are displayed for the allowable errors chosen are \(\underline{\varepsilon }=[{1e2\,\textrm{,}} {1e4\,\textrm{,}} {1e6\,\mathrm{]}}\) in Fig. 12.a-c, Fig. 12.d-f, and Fig. 12.g-i, respectively. Recall, the step size is refined based on the local error estimate \({{\hat{e}}}\) according to Eq. (12), which is determined by the lower-order and higher-order solutions of the embedded integration scheme.

Figure 13 exhibits results for the three third-order embedded ELDIRK methods RK(2)3-Ell, RK(2)3-nstab(\(a_{22}=1/6\)), and RK(2)3-stab(\(a_{22}=1\)). Figure 13.a displays the L2-norm of the order parameters \(||\underline{y}||_2\), Fig. 13.b displays the adaptively refined step size \(\tau \) and Fig. 13.c displays the estimated error for the allowable error \(\varepsilon ={}{1e-6}{}\). In the same way, the allowable errors \(\varepsilon ={}{1e-8}{}\) and \(\varepsilon ={}{1e-10}{}\) are considered in Fig. 13.d-f and Fig. 13.g-i, respectively. The step size \(\tau \) varies heavily during the evolution of the Kobayashi-Warren-Carter problem for all methods. This confirms the importance of adaptivity in time for this problem. Recall, that RK(2)3-Ell and RK(2)3-nstab(\(a_{22}=1/6\)) are not A-stable in contrast to RK(2)3-stab(\(a_{22}=1\)). This fact is demonstrated e.g. in Fig. 13.e, where the refined step size \(\tau \) of RK(2)3-Ell is limited by the stability requirements, whereas the refined step size \(\tau \) of RK(2)3-stab(\(a_{22}=1\)) is only limited by error requirements. Therefore, much larger step sizes can be used with RK(2)3-stab(\(a_{22}=1\)) and the discretization of the total time T requires fewer time steps, resulting in less computational effort. The difference in step size between the instable methods RK(2)3-Ell and RK(2)3-nstab(\(a_{22}=1/6\)) is explained by the critical step sizes \(\tau _{crit}\) obtained from Table 5.

6 Conclusion and outlook

In this paper, we perform stability investigations for the ELDIRK methods, proposed in Mahnken [20, 21]. In particular, we investigate stability properties for the general approach of ELDIRK methods developed in Mahnken [21] for second- and third-order ELDRIK methods. As a first main result we show for both cases, that the stability properties of the higher-order ELDIRK methods can be written solely in terms of the lower-order embedded scheme. This allows, as a second main result, to construct A-stable second- and third-order ELDRIK methods. Moreover, regarding the third-order methods, A-stable and not A-stable methods are distinguished. For the latter case, regions of absolute stability are obtained, which allow a step size selection by means of an accompanying eigenvalue determination. Two numerical examples are concerned with the curing for a thermosetting material and phase-field RVE-modeling for crystallinity and orientation. The numerical results confirm the theoretical results on convergence order and stability.

The investigations regarding convergence order and stabilty in Mahnken [20, 21]; Westermann and Mahnken [23] and in this work are restricted to second and third-order ELDIRK methods. Therefore, an extension of the proposed strategy to higher order ELDIRK methods would be of high interest in future work.