Singly TASE Operators for the Numerical Solution of Stiff Differential Equations by Explicit Runge–Kutta Schemes

In this paper new explicit integrators for numerical solution of stiff evolution equations are proposed. As shown by Bassenne, Fu and Mani in (J Comput Phys 424:109847, 2021), the action on the original vector field of the stiff equations of an appropriate time-accurate and highly-stable explicit (TASE) linear operator, allows us to use explicit Runge–Kutta (RK) schemes with these modified equations so that the resulting algorithm becomes stable for the original stiff equations. Here a new family of TASE operators is considered. The new operators, called Singly TASE, have the advantage over the TASE operators of Bassenne et al. that the action on the vector field depends on the powers of the inverse of only one matrix, which can be computationally more simple, without loosing stability properties. A complete study of the linear stability properties of k–stage, kth–order explicit RK schemes under the action of Singly TASE operators of the same order is carried out for k≤4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \le 4$$\end{document}. For orders two, three and four, particular schemes that are nearly strongly A–stable and therefore suitable for stiff problems are devised. Further, explicit RK schemes with orders three and four that can be implemented with only two storage locations under the action of Singly TASE operators of the same order are discussed. A particular implementation of the classical four–stage fourth–order RK scheme with two Singly TASE operators is presented. A set of numerical experiments has been conducted to demonstrate the performance of the new schemes by comparing with previous RKTASE and other established methods. The main conclusion is that the new integrators provide a very simple solver for stiff systems with good stability properties and avoids the difficulties of using implicit algorithms.


Introduction
Numerical methods for the solution of stiff Initial Value Problems must have adequate stability properties. Unfortunately, for general problems such properties imply that the method must be implicit, which makes the implementation of the integrators complicated and the computational cost can be large for high dimensional problems. However, for particular classes of stiff problems, explicit Runge-Kutta (RK) schemes with suitable stability properties can be constructed. For example, in the semidiscretization of parabolic PDEs, that are known to have large negative real spectrum, the so called (damped) Chebyshev Runge-Kutta methods [1][2][3] can be adequate. These methods are explicit, and with s stages and order less or equal than two they have a stability interval [−α, 0] with α ≤ 2 s 2 (see e.g. [1]), and a thin stability region around the real axis and they have given raise to an integration code [4]. This class of schemes can solve stiff problems with eigenvalues located near the real axis, but at the price of having probably a large number of (explicit) stages.
In [5,6], a family of explicit Runge-Kutta schemes, referred to as Paired Explicit Runge-Kutta schemes, that are suitable for the solution of stiff systems of equations, is proposed. The P-ERK approach allows Runge-Kutta schemes to have a (bounded) large stability region at the price of having a large number of stages, and this make the schemes suitable for midly stiff problems. Highly stiff problems would require a very large number of (explicit) stages.
For stiff problems for which the eigenvalues of the Jacobian matrix are separated into two clusters, one containing the "stiff" or fast components, and the other containing the slow components, Bocher et al. [7] proposed a family of explicit Runge-Kutta schemes that use exponential fitting techniques. The stability region of these methods contains two regions, one close to the origin and the other one fitting the large eigenvalues.
In this paper we are considering a different approach to solve stiff problems by means of explicit RK methods. This approach was proposed by Bassenne, Fu and Mani in a recent paper [8], where, together the differential system associated to the step from (t n , Y n ) to (t n + t, Y n+1 ) the authors consider a second IVP that here will be called the stabilized IVP where T = T ( t L) is a linear operator, referred to as the TASE operator, that depends on the product of the time step size t and a d-dimensional linear operator L that approximates the Jacobian F u (t n , U n ) at the initial values at each time-step. The TASE operator is choosen so that the numerical solution of Eq. (2) U RK (t n + t) obtained with an explicit RK of order p approximates the exact solution of (1), Y (t n + t) at t n + t with order p i.e. and also U RK (t n + t) satisfies some stability requirements, such as A-or L-stability, that are necessary for solving stiff systems. Thus the introduction of the TASE operator into the original governing equation allows us to overcome the numerical stability restrictions of explicit RK time-advancing schemes for solving stiff systems.
In the simplest case of a linear system, i.e., F(t, Y ) = L Y , by defining the TASE operator as T = (I − t L) −1 and solving Eq. (2) with the explicit Euler with the time step size t, we get the approximation U n+1 = U RK (t n + t) ≈ U (t n+1 ) given by which is equivalent to apply the first order BDF method to the original problem (1) and as shown in [9] is L-stable and therefore does not have stability restrictions on the step size t for diffusion problems and other stiff systems with strongly damped components in the solution. On the other hand by choosing T = (I − ( t/2)L) −1 and solving again the stabilized IVP by explicit Euler's scheme we get that is the implicit midpoint rule (see [10], p. 204) that is suitable for convective problems.
Observe that according to the nonlinear variation-of-constants formula [10, p. 96], assuming that F is sufficiently smooth the solutions of (1) and (2) satisfy is the solution of (1). Hence if the TASE operator T = T ( t L) satisfies for some integer q ≥ 1 (the order of approximation of the operator to the identity), it follows from (3) that Moreover if the explicit RK scheme applied to the stabilized equation (2) has order p, then and from Eqs. (5) and (6) it follows that That is, the numerical solution of the stabilized system Eq. (2), U RK (t n ) has at least order min( p, q). A crucial point in the above approach is the construction of TASE operators that satisfy (5) together with the A-stability and/or L-stability requirements. In [8], the authors propose the TASE operator T = T p (α t) defined recursively by Here T p (α t) has an accuracy of order p and α > 0 is a real parameter to be chosen by the stability requirements. By choosing α ≥ α min = (2 p − 1) C −1 p , where C p is the stability abscissae of the explicit RK scheme of order p, some specific schemes that combine explicit p-stage RK schemes with orders p ≤ 4 and TASE operators with the same orders have been derived in [8], and their linear stability properties have been studied as well. Moreover applications to some benchmark equations showed that this approach leads to time marching algorithms that can be easily implemented and do not suffer from the stability restrictions on the step size.
More recently, Calvo et al. [11] have considered a more general family of TASE operators with the form where β j and α j > 0, j = 1, . . . , p, are free real parameters with α j distinct between them, and β j uniquely determined by the condition T p = I + O( t p ). Here the free parameters α j have been selected to improve the linear stability properties of an explicit RK scheme with order p for the scaled equation Eq. (2). Even though TASE operators are very recent, they have already been considered in a number of related papers ( [12][13][14]) which encourages us to advance in this research.
Note that, in the TASE operators of types described by Eqs. (8) and (9), each T p evaluation involves the inversion of the p linear operators (I − t α j L), j = 1, . . . , p, for given positive α j . In fact, for an s-stage explicit RK scheme defined by the Butcher tableau [15] 0 c 2 a 21 . . . (2) is advanced from (t n , U n ) → (t n+1 = t n + t, U n+1 ) by the formulas where K j , j = 1, . . . , s are computed recursively from the formulas The explicit scheme defined by the Eqs. (10), (11), will be called a RKTASE scheme with s stages and order p. Now the computational cost is the s evaluations of the vector field in Eq. (11) and the computation of the p × s matrix operations Note that these operations do not require the computation of the inverse of the matrices, because as usual they can be done by solving the associated linear systems.
The aim of this paper is to study a new family of RK TASE schemes associated to TASE operators given by where α > 0 and β p, j , j = 1, . . . , p are real parameters. It will be seen that β pj , j = 1, . . . , p are uniquely determined as functions of α by the condition that ST p ( t L) = I + O( t p ). Note that this family is not a subclass of the family (9) because that family dependes on the operator (I − t α j L) −1 with the power −1 whereas in the proposed family the operator The reason of the choice (12) is that the calculation of the stabilized stages ST p F k requires a matrix operation with only the powers of the inverse of one matrix (I − t α L). Moreover it will be seen that with a suitable choice of the free parameter α > 0 the linear stability properties are similar to the RK TASE schemes with the TASE operator (9). Clearly, these operators require the solution of s linear systems all of them with the same matrix (I − t α L), similarly to the case of Singly Implicit RK methods [16,17] and the Singly Diagonally RK methods [18,19] when simplified Newton is used to solve the implicit stage equations. For this reason, the TASE operators (12) will be called Singly-TASE operators hereafter. Moreover, for high dimensional ODEs that arise in the spatial semidiscretization of some PDEs, the techniques used in their numerical solution must take into account not only the accuracy and stability properties and it may be convenient to use low-storage RK schemes [20,21]. Then we will consider the use of the schemes that combine 2N -storage RK schemes where N is the dimension of the system of ODEs, with Singly TASE operators.
The outline of this paper is as follows. In Sect. 2, we propose a new family of TASE operators of order p that requires the computation of p linear systems with the same matrix of coefficients. They depend on one free parameter d = α −1 that can be optimized to adjust the stability properties of the resulting Singly-RKTASE scheme. In Sect. 3, the stability of twostage Singly-RKTASE schemes with order two is studied, and the necessary and sufficient conditions on the free parameter d for A-stability and strong A-stability are obtained; In Sect. 4, the stability of three-stage Singly-RKTASE schemes with order three is considered. It is observed that there are no A-stable schemes, but there exist A(θ )-stable schemes with θ = 89.05 degrees and also strongly A(θ )-stable schemes with θ = 88.99 degrees; Sect. 5 deals with the case of four-stage schemes with order four. It is noticed that no scheme of this class can be A-stable, but there are A(θ )-stable and strongly A(θ )-stable schemes with θ = 87.18 and θ = 87.17 degrees, respectively. Finally, in Sect. 5, we present the results of some numerical experiments with selected stiff benchmark problems, comparing the new schemes with previous RKTASE schemes and also with Rosenbrock and Singly Diagonally Implicit RK (SDIRK) methods.

Singly-RKTASE Schemes
To simplify some expressions, Eq. (12) will be written in the equivalent form First of all, we determine the coefficients b pj of Eq. (13) so that ST p has order ( t) p , i.e., This is equivalent to say that the scalar function ST p (z), defined as should be of order z p , i.e., Then such a requirement leads to the non-singular linear system of p equations with regard It has the unique solution Furthermore, the error coefficient Q p of the Singly-TASE operator Eq. (14) with the coefficients Eq. (15) satisfies Also in view of Eq. (15), the function Eq. (14) can be written in the compact form In particular for p = 2, 3, 4, the coefficients b pj have the values Next, we study the linear stability properties of an explicit s-stage RK scheme of order p combined with a singly TASE operator of the same order. When applied to the scalar test equation y = λy with λ ∈ C, the modified equation becomes is the stability function of the explicit s-stage RK scheme with some real coefficients a j , the stability function of the stabilized scheme RST (z) will be Recall the following stability properties: its stability region contains the left hand region of the complex plane with angle θ . If in addition We first derive a necessary condition for A-stability. Since and therefore for given values of s ≥ p, a necessary condition for A-stability is that the available positive parameter d satisfies A sufficient condition for A-stability can be derived from the application of the modulus maximum theorem to RST s (z) along the imaginary axis z = i y. In fact, it follows from Eq. (16) that RST s (z) has the form where the numerator π 2s (z) is a polynomial in z of degree ≤ 2 s with coefficients depending polynomially on d. Now, is an even function of y and taking into account that the scheme has order p with μ j polynomial functions of d.
Clearly, a necessary and sufficient condition for A-stability is that the numerator of Eq. (18) should be ≤ 0 for all y 2 > 0, or else that the polynomial for all u > 0. In general, this condition becomes too complicated to be checked, but there is a necessary condition that can be easily checked, for example that μ 2 p ≤ 0.

Second Order Singly-RKTASE Schemes
We need at least a two-stage RK scheme with order two combined with a second-order Singly-TASE operator. Next we study the linear stability properties of this scheme. The stability function of a two-stage RK scheme with order two is and it can be seen that Since μ 4 < 0, μ 6 < 0 and μ 8 ≤ 0 for all d ∈ (0, 1], for all d ∈ (0, 1] all second order Singly RK schemes with s = 2 are A-stable. Moreover for all d ∈ (0, 1), |RST 2 (∞)| < 1 and the schemes are strongly A-stable. In particular min |RST 2 (∞)| = 1/2 is attained for d = 1/2 and therefore there are no L-stable second order Singly RK schemes with two stages.
If we want to have L-stable second order Singly RK schemes we must consider second order RK schemes with s = 3 three stages. Now the stability polynomial will be and the condition is satisfied if and only if A numerical study of the conditions (18) and (19)  In the right plot of Fig. 1, we display (solid line) the boundary of the stability region of the Singly-RKTASE RT 2 for d = 1 as well as (dashed line) that of the Singly-RKTASE scheme with d = 1/2. A remarkable fact is that the stability regions are not connected. In fact, apart from the component that includes the negative complex plane, they have a small component  Fig. 1 the boundary of the stability regions of the RKTASE schemes proposed in [8,11].
In addition to the plots, we include in the Figure two tables with the error coefficient Q 2 and the stability angle θ that, since all the schemes are A-stable, is 90 degrees.
In view of the Tables displayed in Fig. 1 it can be seen that the RKTASE schemes, in spite of having been selected within a family with two free parameters α 1 and α 2 , have larger error constants than those of the Singly-RKTASE schemes, for which there is only one free parameter α = 1/d. Finally, it is worth to remark that (see [10]) there exist a oneparameter family of two-stage RK schemes with order two, and each scheme of this family supplemented with Singly TASE ST 2 leads to Singly-RKTASE with the same linear stability properties.

Third Order Singly-RKTASE Schemes
We will consider three-stage third order RK schemes. The third-order Singly TASE operator is defined by and taking into account that the stability function of a three-stage RK scheme with order three is the stability function of the third-order RKTASE is Recalling that R 3 (z) given by Eq. (22) Moreover, since and for all d that satisfies Eq. (23) the coefficient of y 4 in Eq. (24) is positive, all the schemes that satisfy Eq. (23) cannot be A-stable. In fact, the stability region does not contain the imaginary axis around the origin and is nearly A-stable.
On the other hand, since R 3 (u L ) = 0 for and consequently the corresponding scheme is L(θ )-stable. Hence we conclude: Theorem 2 (1) A Singly-RKTASE scheme of order three with three stages cannot be Astable.
In the right plot of Fig. 2, we display (solid line) the boundary of the stability region of the third-order Singly-RKTASE scheme with d = d A together with (dashed line) that of the Singly-RKTASE scheme with d = d L that is nearly L-stable. In the left plot of the Fig. 2 we also display the boundary of the stability regions of the RKTASE schemes of order three proposed in [8,11]. In addition to the plots, we include in the Fig. 2 two tables with the error coefficient Q 3 and the stability angle θ in degrees.
For the third-order schemes, the error coefficients of the Singly-RKTASE schemes are smaller than those of the RKTASE schemes. However, the stability angle (none of them is A-stable) is slightly smaller for the Singly-RKTASE schemes.
To show the behavior of the boundaries around the origin, we present here a zoomed-in view around this point of the boundaries in Fig. 3.

Remark 1
We must remark that (see [10]) there exist a two-parameter family of three-stage RK schemes with order three. All of these schemes possess the same linear stability function R(z). Therefore all the associated Singly RKTASE scheme of order three have the same stability properties that depend only on the parameter d in the ST 3 operator. Hence we may select the RK scheme taking into account some additional criteria as accuracy or low storage requirements.
As an example of 2N -storage implementation of a third order RK scheme with three stages supplemented with a Singly TASE ST 3 , we may consider the RK scheme proposed by Williamson [20] and Carpenter and Kennedy [21] with the Butcher array where the coefficients have been selected so that they minimize the local truncation error among the three stage RK schemes with order three that can be cast in the 2N -storage implementation of Williamson. Recall that in this implementation the two 2N -storage are U j and dU j that are computed by (see [20]) with p = 3, U 0 = Y 0 , U 3 = Y 1 and A j and B j given by Each successive stage is written onto the same register without zeroing the previous value.

Fourth-Order Singly RKTASE
First of all we consider explicit four-stage fourth-order RK schemes supplemented with Singly TASE operators of order four. Following the notations of the previous sections, we have Taking into account that where γ is the real stability abscissae of R 4 (z), given by On the other hand, since that is attained at The optimal strong A-stability scheme will be for Concerning the second necessary condition for A-stability, we have and taking into account that the coefficient of y 6 in Eq. (27) is non negative for all d ∈ (0, −γ 4 /4], all schemes that satisfy this condition are not A-stable. In fact, small values of z = iy near the origin in the imaginary axis are not contained in the stability region of the scheme, and we conclude:  In the right plot of Fig. 4, we display (solid line) the boundary of the stability region of the fourth-order Singly-RKTASE for d = d A = −γ 4 /4 together with (dashed line) that of the Singly-RKTASE scheme with d = d L that is nearly strongly A-stable. In the left plot of the Fig. 4, we display the boundary of the stability regions of the schemes of order four proposed in [8,11]. In addition to the plots, we include in the Fig. 4 two Tables with the error coefficients Q 4 and the stability angles θ in degrees. It can be seen in the Tables that the properties of Singly-RKTASE schemes are very similar to those of RKTASE schemes.
To show the behavior of the boundaries around the origin, we present here a zoomed-in view around this point of the boundaries in Fig. 5.
One may wonder whether or not the consideration of RK schemes with five stages and order four would make possible to have L-stable schemes. For a five-stage fourth-order RK scheme the stability function is with some constant ν = 0. Then a necessary condition for L-stability is lim z→∞ R 5 z ST 4 (z) = 1 120 (120 − 480d + 960d 2 − 1280d 3 + 1280d 4 − 1024d 5 ν) = 0 and a numerical study shows that this condition is not compatible with Remark 2 It can be seen (see [10], page 135) that there exist several families of four-stage fourth-order RK schemes depending on one or two parameters. All these schemes have the same linear stability function R(z). Clearly each of these schemes combined with the Singly TASE ST 4 gives a Singly RKTASE with order four whose stability properties depend only on the parameter d. Hence we may select the RK scheme taking into account some additional criteria as accuracy or low storage requirements. We can for example take a four stages fourth order low-storage Runge Kutta scheme (see e.g. [22]) or else the classical RK of order four and combine it with the Singly TASE operator ST 4 and apply it to stiff problems.

The Implementation of Singly RKTASE Schemes
In this section we examine the main issues in the implementation of a Singly RKTASE scheme for the numerical solution of stiff IVPs (1). We will consider the case of four stages, s = 4, but the generalization to any other s is straightforward. For the Singly TASE operator A first point is the choice of the linear operator L in ST 4 . If the original IVP (1) is linear F(t, Y ) = M(t) Y +g(t), the safest option is to take at the step from t n to t n + t, L = M(t n ).
In the non linear case we could take L = ∂ Y F(t n , Y n ). However, the order of the scheme is not reduced if L is not the exact Jacobian, so we can take an approximation instead whenever the stability of the numerical solution is maintained. Something similar happens with implicit Runge-Kutta methods, where the Jacobian matrix is often frizzed for several steps unless the iterative scheme used has som problems to converge. Thus, to reduce the computational effort, we can freeze the Jacobian for several steps (even the whole integration), or even take an approximation to it that makes the solution of the linear systems cheaper.
A second point is the computation of the scaled stages (11) Taking into account that we compute (28) by using a Horner type algorithm because due to (29), K k can be written in the form And therefore given F k computed by (11) and W = d I − t L, the scaled stage K k of (28) can be computed by the recurrence 1: K k ← 0 2: for j from 4 to 1 step −1 do 3: The implementation of the complete RKTASE integrator is very simple. We can take the implementation of the explicit Runge-Kutta scheme (of order four in this case) and add the above algorithm just after each evaluation of the vector field. With this, the integrator can be used to solve stiff problems.

Numerical Experiments
In this section, some numerical experiments are presented to confirm the theoretical analyses in the above sections and to compare the performance of our proposed Singly-RKTASE methods with that of the methods in [8] and [11]. Specifically, we will consider the following fourth-order RKTASE and Singly-RKTASE methods based on the classical RK scheme of fourth-order given by the Butcher array TASE-S Strongly A(θ )-stable, with minimum value |RT 4 (∞)| = 0.270395 [11] and the parameters are given as We want also to compare the proposed methods with some Rosenbrock methods [23], which have a structure similar to RKTASE schemes, and with the SDIRK methods [18,19]. More specifically, the following methods will be used.
SDIRK-A Fourth-order, three-stage, SDIRK A-stable scheme (see Hairer et al. [9]) defined by the Butcher array SDIRK-L Fourth-order, five-stage, SDIRK L-stable scheme (see Hairer et al. [9]) defined by the Butcher array With these methods, we will integrate the following four stiff differential problems, with which two of them are linear while the other two are non-linear.
• Problem 1 (linear) The 1D diffusion of a scalar function y = y(x, t), with a time dependent source term (taken from [8] and [11]) The solution y = y(x, t) is assumed to be 2π-periodic in x. As the initial condition, we have taken For the spatial discretization, we have used fourth-order centered difference schemes with a grid resolution of N = 512 ( x = 2π/512) and periodic conditions are assumed at the two domain boundaries. The real part of the eigenvalues of the Jacobian matrix of the semi-discrete problem ranges from −3.5 × 10 4 to 2.79 × 10 −5 . • Problem 2 (linear) The two-species advection-diffusion-reaction problem in the bounded domain of 0 ≤ x ≤ 1 is taken from [8] as with K = 10 4 , D = 10 2 , and U = 10 2 . Two different initial and boundary conditions have been considered as (i) (ii) Dirichlet condition is assumed at the two boundaries, and for the spatial discretization, we used fourth-order centered difference schemes in x ∈ [0, 1], with a grid resolution of N = 256 ( x = 1/256). The real part of the eigenvalues of the Jacobian matrix of the semi-discrete problem ranges from −3.5 × 10 7 to 1.4 × 10 −9 in both cases, suggesting that the problem is highly stiff. • Problem 3 (non-linear) The 1D Burgers' equation can be written in the conservative form as where ε > 0 denotes the constant viscosity coefficient. The initial condition y(x, 0) is given as The solution is 2π-periodic, and we have solved this problem on a spatial grid of 0 ≤ x ≤ 2π using fourth-order centered difference schemes with a grid resolution of N = 512 ( x = 2π/512). Periodic conditions will be used for both boundaries. The real part of the eigenvalues of the Jacobian matrix of the semi-discrete problem ranges from −3.5 × 10 3 to 2.6 × 10 −6 , indicating that the problem is mildly stiff. • Problem 4 (non-linear) This non-linear problem is taken from [8] with the governing equation as where β is a positive constant. For the spatial discretization, we have used fourth-order centered difference schemes with a grid resolution of N = 512 ( x = 2π/512). Neumann conditions, i.e., ∂ y/∂ x| x=−5 = 2.5e −25/4 and ∂ y/∂ x| x=5 = −2.5e −25/4 , are assumed at the left and right domain boundaries, respectively. The real part of the eigenvalues of the Jacobian matrix of the semi-discrete problem ranges from −1.38 × 10 4 to 6.41 × 10 −6 .
With respect to the implementation of TASE operators, for the linear problems 1 and 2 we have taken as operator L the Jacobian matrix of the vector field F(t, Y ). Since in both cases this matrix is constant it was computed just once at the first time-step. For the nonlinear problem 4 we have taken as operator L the Jacobian matrix of the vector field F(t, Y ), computed by standard finite difference schemes and it has been re-evaluated at each time-step.
For problem 3 we have taken as operator L the Jacobian matrix of only the diffusion term, which is not the Jacobian matrix of the full vector field at any point (t n , Y n ). In this way, the matrix L is constant and need to be computed only at the first step of the integration. This has been also used in the implementation of Rosenbrock and SDIRK methods to save computational effort.
In all the cases, the corresponding linear systems were solved by factorizing the matrix in the LU form.
For the solution of the non-linear systems that appear in SDIRK methods, we used a simplified Newton method, iterated until the norm of the difference between two consecutive approximations was smaller than the threshold of 10 −9 .

Results
Regarding the comparison between Singly-RKTASE and previous RKTASE schemes, in all the cases, all the considered methods behave very similarly with respect to accuracy as well as stability. The only relevant difference is that the new Singly-RKTASE schemes required, as expected, a lower computational cost than the standard RKTASE methods.
Concerning the computational efficiency, the Singly-RKTASE schemes are clearly more efficient than the standard RKTASE methods, and, moreover, the strongly A(θ )-stable methods are more efficient for highly stiff problems.
Regarding the comparison with the classical SDIRK and Rosenbrock methods, for the highly stiff problems the strongly stable Singly-RKTASE and RKTASE schemes present efficiency plots (see Fig. 10) comparable to the other methods.
For problem 3, the efficiency of Singly-RKTASE schemes is similar to that of SDIRK methods for large time step sizes. Rosenbrock scheme (see Fig. 12) looses efficiency because it does not use the exact Jacobian matrix and this makes the method to lose order of accuracy.  (Fig. 13).
Rosenbrock and SDIRK methods present better efficiency plots than Singly-RKTASE methods for problems 1 and 4 (see Figs. 7, and 14). The reason is not the lower computational cost of these methods (indeed, the cost is similar for all the methods), but the lower global error given by the integrators. Note that the selected SDIRK and Rosenbrock methods are among the most efficient ones in their families. This encourages us to undertake a further research to find the most efficient combination of the RK schemes and the TASE operators.
Also note that SDIRK schemes require, at each time step, the solution of several non-linear systems by means of some iterative method. It runs the risk that, for a specific problem, the iterative method fails to converge. In general, any implicit method requires at each step the solution of one (or more) non-linear implicit system of equations, and therefore the overall efficiency of the integrator depends strongly on the non-linear solver used. On the other hand, An advantage of the Singly-RKTASE methods with respect to the Rosenbrock schemes is that the order of accuracy does not change if we take an approximation to the Jacobian matrix [24], which further allows to save Jacobian matrix evaluations and LU matrix factorizations, as was the case for example with problem 3.

Conclusions
In this work, a new family of STASE operators for the solution of stiff differential equations has been proposed. The advantage over previous TASE operators [8,11] is that they require the solution of linear systems with the same coefficient matrix, whereas in the previous cases there are different coefficient matrices. As a consequence, the new STASE method only relies on the inverse of a single matrix. The new schemes are computationally comparable to the classical Rosenbrock methods. A detailed study of the accuracy and the linear stability properties up to fourth-order methods has been carried out, resulting in nearly A-stable and strongly A-stable Singly-RKTASE schemes. The results of several numerical experiments have been presented by comparing two new Singly-RKTASE fourth-order schemes with the standard RKTASE methods of the same order in [8,11]. The new methods are also compared with the classical SDIRK and Rosenbrock methods.
The main conclusions derived from the numerical experiments are as follows.
• Singly-RKTASE schemes are more efficient than their corresponding RKTASE methods.
• For highly stiff problems, strong stability is an essential feature.
• In most problems, Singly-RKTASE schemes introduce larger global errors than the SDIRK and Rosenbrock methods. Further research needs to be done to find an efficient and optimal combination of the RK schemes and the STASE operators.