Two-step peer methods with equation-dependent coefficients

We introduce a new class of explicit two-step peer methods with the aim of improving the stability properties of already existing peer methods, by making use of coefficients depending on the Jacobian of the Ordinary Differential Equations (ODEs) system to solve. Numerical tests highlight the best stability and accuracy properties of the new methods compared to the classical and equation-dependent ones proposed in Conte et al. (Lect Notes Comput Sci 12949:309–324, 2021).


Introduction
It is the purpose of this work to derive new explicit two-step s-stage peer methods with improved stability properties for the numerical solution of first-order ODEs of the type (1) 2017), implicit-explicit (Schneider et al. 2021;Weiner 2017, 2018), and parallelizable (Kulikov and Weiner 2010;Schmitt and Wiener 2010;Schmitt et al. 2009;Weiner et al. 2012) peer methods have been derived. There are also specific techniques that can allow to build peer methods adapted to the problem to solve. For example, if it is a priori known that the analyzed problem has an oscillating solution (Budroni et al. 2021a, b), the Exponential Fitting (EF) technique (Ixaru 1997;Ixaru and Vanden 2004) leads to peer methods with coefficients depending on the oscillation frequency (Conte et al. 2018(Conte et al. , 2019(Conte et al. , 2020b, that are much more accurate than the classical ones. Peer methods are characterized by several stages like Runge-Kutta schemes, but unlike the latter all of them show the same accuracy and stability properties. They have been introduced with the aim of combining the advantages of Runge-Kutta and multistep methods and are very convenient. In fact, peer methods do not suffer from order reduction if applied to ODEs systems with a decidedly high stiffness. Furthermore, they are suitable for parallel implementation, as the actual stages rely only on the previous ones. In this paper, we focus on explicit peer methods with fixed step-size h. Explicit methods are less expensive than implicit ones, but they usually have worse stability properties. For this, the main aim of our work lies in the derivation of new peer methods that are still explicit, but with improved stability properties than the classical ones. Considering a time discretization {t n , n = 1, ..., N } of the integration interval [t 0 , T ] related to the ODE (1), classic explicit peer methods can be expressed as where A = (a i, j ) s i, j=1 , B = (b i, j ) s i, j=1 and R = (r i, j ) s i, j=1 are matrices containing the coefficients of the method. I d is the identity matrix of order d. In our case, since we analyze explicit methods, R is strictly lower triangular, i.e., r i, j = 0 ∀i ≤ j. The notation used to represent these methods in vectorial form (2) is the following: , Y n,i ≈ y(t n,i ), t n,i = t n + hc i , where the nodes c i are assumed to be distinct with c s = 1. The stages Y [n] , and therefore F(Y [n] ), are column vectors.
To determine the methods proposed in this work, we apply a technique that leads to new coefficients with respect to classical schemes, depending on the Jacobian of the problem. This technique briefly consists in modifying the classical case by considering a different expression of the stages of the method. By imposing the order conditions using the new form of the stages, the new coefficients of the method are obtained. We were inspired by Ixaru (2012), in which a similar procedure was applied to explicit Runge-Kutta methods improving their accuracy and stability properties, as evidenced also by the numerous numerical tests conducted in Conte et al. (2020a). An extension of this methodology has also been applied to peer methods in Conte et al. (2021), in which the authors have managed to derive Jacobiandependent coefficients, slightly improving the stability properties of classical schemes. In this work, we show that by changing the approach proposed in Conte et al. (2021), it is possible to obtain New Explicit Jacobian-Dependent Peer (NEJDP) methods with much better stability and accuracy properties than those obtained in Conte et al. (2021), which will be called Old Explicit Jacobian-Dependent Peer (OEJDP) methods.
Moreover, this paper also completes the work done in Conte et al. (2021), where the coefficients of the proposed methods have been derived only in the scalar case and when the accuracy order and the number of stages are equal to two. Now, we derive the old methods coefficients (which become matrices) also for non-scalar problems. We also provide the order conditions of these methods even when the accuracy order p and number of stages s satisfy p = s ≥ 2. Therefore, Sect. 2 of this paper is devoted to the extension of OEJDP methods of Conte et al. (2021), while the following ones are devoted to the derivation of NEJDP methods.
Specifically, this paper is organized as follows: in Sect. 2, we extend the work done in Conte et al. (2021), formulating the OEJDP methods in the non-scalar case and deriving the related order conditions for p = s ≥ 2; in Sect. 3, we derive two-stage NEJDP methods of order two, showing the relative linear stability properties; in Sect. 4, we carry out numerical tests highlighting the advantages of the new methods over the old Jacobian-dependent and classical ones; in Sect. 5, we discuss the obtained results and the possible future research.

Old equation-dependent peer methods
In this section, after recalling the formulation of the two-stage OEJDP methods for scalar ODEs (Sect. 2.1), we first of all extend these methods to the case of s-stages (Sect. 2.2), and afterwards, we derive two-stage OEJDP methods able to solve differential problems of any dimension d (Sect. 2.3).

Two-stage OEJDP methods for scalar ODEs
Consider explicit two-stage peer methods: + hr 21 f (t n,1 , Y n,1 ). (3) By defining the error operators L 1 and L 2 associated, respectively, with the first and second stages as where Y 1 and Y 2 represent the continuous expression of Y n,1 and Y n,2 it has been shown that annihilating L 1 (t k ) and L 2 (t k ), k = 0, 1, 2, at t = 0, i.e., annihilating the moments L i,k , i = 1, 2, k = 0, 1, 2, leads to the following order conditions for classic peer methods of accuracy order equal to two, respectively Instead, in the derivation of equation-dependent methods, it was assumed that the first stage computed at t n,1 , i.e. Y n,1 , is affected by an error err 1 , which has the form To prove this, the following relationship between error operators and moments has been used: By doing so, the expression of Y 2 (and therefore that of L 2 ) has changed The function j 1 is the Jacobian related to the ODE (1) at Y 1 , i.e., j 1 (t) = f y (t + hc 1 , y) |y=Y 1 (t) . By annihilating L 1,k , k = 0, 1, 2, and the new L 2,k , k = 0, 1, 2, 3, the same conditions (6) and (7) related to the classical case are obtained, plus an additional one The resulting two-stage peer methods are still explicit and of order two, but they exhibit better accuracy and stability properties than the classical ones. The Jacobian-dependent scalar coefficients of these methods, which can be derived by imposing the order conditions (6), (7), and (11), are a 11 = − b 11 (−1 + c 1 ) 2 + c 2 1 / 2(−1 + c 1 ) , Note that the moments associated with L 2 (t k ) have been annihilated also for k = 3, i.e., the second stage is calculated with accuracy order equal to three. That is why there is an additional order condition (11) than the classical case and why this technique can also lead to better accuracy properties.

s-stage OEJDP methods for scalar ODEs
In this subsection, we derive the order conditions of the OEJDP methods in the general sstage case where the order is p = s. To do this, recall the definition of accuracy order of an explicit peer method ).

Definition 1 The peer method
The values of i are called residuals. Imposing that the method (13) has accuracy order p = s, leads to express its coefficients in the following form ): B1 = 1, 1 = (1, ..., 1) T (pre-consistency condition), What we want to do now is to express the coefficients of OEJDP methods with accuracy order p = s in a similar form. Let us consider the conditions (6), (7), and (11) that come out by canceling the moments L 1,k , k = 0, 1, 2, and L 2,k , k = 0, 1, 2, 3, in the simple case of two-stage schemes of order two. All the equations of (6) and (7) correspond to the classical order conditions. Equation (11) returns an additional order condition with respect to the classical case, which arises by assuming that in general, for s-stage methods, the first s − 1 stages are affected by error and the last stage is more accurate than the previous ones.
Define the continuous expression of the ith stage The related error operator is .., s assume the following form: Note that err i depends on the product To understand why the error associated with Y i takes this form, just extend the procedure shown in Conte et al. (2021) for err 1 (8), observing that it depends on the moment L 1,3 . Still referring to the simple two-stage case, look at Eq. (11), in which at the end, the product between j 1 and L 1,3 (deprived of h 3 which has been simplified) just appears. In fact, here, the moment L 2,3 is annihilated, whose expression contains L 1,3 , i.e., it depends on t err 1 . In the general s-stage case, therefore, the moment L s,s+1 must be annihilated, which depends on all the errors made at the previous stages.
Define the following column vector operr ∈ R s×1 to derive a compact notation for the moments associated with the error operators: The vector operr essentially contains in the ith component the moment (deprived of the powers of h) associated with L i (t s+1 ). Therefore, once we have determined the first s − 1 (in consecutive order) components of operr, we can calculate its last component operr s . Annihilating operr s corresponds to the additional order condition we have for the s-stage OEJDP methods. Summarizing, these methods have order p = s if their coefficients satisfy

Two-stage OEJDP methods for systems of ODEs
In this subsection, we show how the coefficients of the OEJDP methods can be obtained in the multi-dimensional case (d ≥ 1), practically deriving them for two-stage methods of order two.
As seen above, explicit peer methods can be expressed in the vectorial form (2), in which the coefficients matrices A, B, and R are tensorially multiplied by the identity matrix of order d (which from now on we simply denote with I instead of I d ). However, since the Jacobian related to the ODE (1) has dimension d, the Jacobian-dependent coefficients must be re-formulated appropriately, and there is no need to multiply them by the identity matrix.
For example, in the case of two-stage methods, the coefficients take the form (12). Specifically, a 21 , a 22 and r 21 depend on the Jacobian, and must therefore be re-formulated as matrices. For this reason, let us indicate them with capital letters A 21 , A 22 , and R 21 , respectively. To understand how to derive such matrices, let us accurately show the computation of R 21 .
The coefficient r 21 of (12) can be expressed as r num 21 /r den 21 , where Now, isolate the coefficients of r num 21 and r den 21 that do not multiply the Jacobian j 1 . Considering r den 21 , the only coefficient that respects this property corresponds to −6c 1 and appears in the last term, while r num 21 does not contain j 1 . Therefore, define the coefficient r c 21 = r num 21 /(−6c 1 ). The next step involves highlighting the term h j 1 in r den 21 , dividing everything by the coefficient −6c 1 previously determined. This operation must be done, in general, also for r num 21 (dividing it by the numerator of r c 21 ). However, in this case, we can avoid that as r num 21 does not depend on the Jacobian. This leads to By indicating the Jacobian with the capital letter J in the multi-dimensional case, we obtain the following form for the matrix R 21 : Similarly, the matrices A 21 and A 22 can be derived Here, a num 21 , a num 22 , and a den 21 correspond to the numerators and denominators of a 21 and a 22 of (12), deprived of the coefficients that do not multiply h j 1 .
Note that the matrices A den 21 and A den 22 match, as the denominators of a 21 and a 22 of (12) are equal. This is a good thing, as these matrices need to be inverted at each time-step. Thus, the computational cost of this two-stage method is halved.
When we deal with numerical tests in the next sections, we will show that these Jacobiandependent peer methods work well even in the multi-dimensional case, resulting once again better than the classical ones.

New equation-dependent peer methods
As mentioned in the Introduction, in the paper Conte et al. (2021) the authors derived the OEJDP methods shown in the previous section by applying a technique similar to that adopted in Ixaru (2012) for Runge-Kutta methods. However, the procedure applied for Runge-Kutta and peer methods could not be the same, as for the latter the stages in the current interval depend on those in the previous one. In the work Conte et al. (2021) this was ignored. That is, the OEJDP methods have been derived by assuming that only the stages Y n, j , j < i, in the computation of Y n,i , are affected by error. Now, we derive NEJDP methods by assuming that also the stages Y n−1, j , j < i, are affected by error in the computation of Y n,i .

Two-stage NEJDP methods of order two
In this subsection, we determine two-stage NEJDP methods by following the observations just made. The continuous expressions of the stages Y n,1 and Y n,2 in this case are Assuming that Y n−1,1 is affected by error (which we call err −1 , while with t err −1 we indicate the local truncation error which is obtained by imposing that the first stage has order two) means that the following writings hold: Once again, the Jacobian j 1 (t −h) = f y (t +h(c 1 −1), y) |y=Y 1 (t−h) arises from the application of Taylor series expansion. Substituting the expressions (25) in the continuous stages (24) leads to By defining the error operator related to Y n−1,1 as it's possible to determine err −1 and therefore t err −1 . Following a similar approach to the one detailed in Conte et al. (2021), to have an order-two first stage we have to annihilate the moments related to L −1 (t k ), k = 0, 1, 2 (it is easy to see that annihilating them leads to conditions equivalent to (6)). The resulting error err −1 is given by the following expression, that involves the moment L −1,3 : In order to determine err −1 we have used the property (9) that binds error operators to the related moments. Note that by defining L −1 as done in (27) we have implicitly assumed that, for simplicity, there are no further errors at the earlier stages Y n−k,1 , k ≤ 2. Therefore, the additional error we consider compared to the OEJDP methods concerns Y n−1,1 .
It is necessary to define the error operators associated with the stages Y n,1 and Y n,2 to determine the order conditions of the NEJDP methods. Such operators have the same form as (4), but obviously the continuous expression of the stages is now given by (26).
To calculate the first stage with accuracy order equal to two, the moments L 1,k , k = 0, 1, 2, must be annihilated. Note that, by (26), L 1 is completely known, as the expression of t err −1 is known (28). This operator, evaluated at y(t) = t k , assumes the following form: Canceling the moments L 1,k , k = 0, 1, 2, leads to order conditions equivalent to that obtained in the classical case (6) for the first stage. Leveraging property (9), it is possible to conclude that, for NEJDP methods, the error associated with the first stage Y n,1 is The final step now consists in deriving the conditions coming from the imposition that the second stage has accuracy order equal to three. Therefore, we consider the continuous expression of the second stage obtained previously (26), assuming that Y n,1 is affected by the error (30): As usual, the Jacobian j 1 (t) = f y (t + hc 1 , y) |y=Y 1 (t) appears by applying Taylor expansion of y (t + hc 1 ) = f t + hc 1 , Y 1 (t) + err 1 (t) at the second component.
By defining the error operator associated with Y n,2 as before (4), considering this time the new continuous expression of the second stage (31), we can determine the moments L 2,k , k ≥ 0, which are obtained by canceling L 2 (t k ) at t = 0. Report the expression of L 2 (t k ): Now, the moments L 2,0 , L 2,1 , L 2,2 and L 2,3 can be easily calculated. While canceling the first three, we obtain exactly the same order conditions of the classical case (7), annihilating the last moment leads to the following additional order condition, different from that derived for OEJDP methods: To conclude, therefore, the NEJDP methods that we propose in this paper are obtained by imposing the order conditions (6), (7) and (33). This leads to coefficients which, compared to the OEJDP methods, depend on the Jacobian calculated at the current and the previous grid points: r cden 21 =12(−1 + c 1 )c 1 . The coefficients a 11 and a 12 are scalar and correspond to that of OEJDP methods (12). Furthermore, as in the previous case, the denominators of A 21 and A 22 correspond. For convenience, we have reported the coefficients of the NEJDP methods directly in non-scalar case.
The coefficients dependence on the Jacobian also at the previous point does not cause an additional computational cost (with respect to the OEJDP methods) in terms of function evaluations. In fact, the Jacobian determined at the previous grid points can be stored in an array, without the need to re-evaluate it at each step of the method. We note that, as with the OEJDP methods, we still have two-stage methods of order two, and the second stage is always computed by imposing that it has order three. Since we have also considered the errors on the stages in the previous interval, it is reasonable to expect that the new methods are more accurate than the classic and old equation-dependent ones. This observation will be confirmed in numerical tests.

Stability analysis of NEJDP methods
In the current subsection, we show the absolute stability region of the NEJDP methods obtained in this paper. Through a careful selection of the parameters left free, the NEJDP methods result more stable than the OEJDP methods (12) obtained in Conte et al. (2021), and therefore also than the classical ones. Furthermore, in numerical tests, we will see that they also have better accuracy.
To perform the linear stability analysis of the considered peer methods, we need to solve with them the test equation where λ is a complex number with negative real part, by evaluating for which values of hλ the numerical solution does not explode. The application of the peer method (2) to the test Eq. (35) leads to By defining the matrix M(z) = (I − z R) −1 (B + z A), it is possible to conclude that the numerical solution does not explode if its spectral radius is in modulus less than one. Obviously the NEJDP, the OEJDP, and the classical method have three different stability matrices M, respectively, as they have different coefficients. Consequently, in order to analyze the stability properties of these three methods, the eigenvalues of the respective matrices M must be determined, evaluating for which values of z they assume a value less than one in modulus.
Regarding the OEJDP and classical peer methods, the values of the free parameters for which the relative stability region has real part as large as possible have already been determined in Conte et al. (2021). For classical methods, these parameters are while for the OEJDP schemes they assume the following values: Conducting the same analysis for the NEJDP methods derived in this work, using the Matlab fmincon function (considering random initial values for the coefficients), we have ascertained that the parameters that maximize the absolute stability region real part are b 11 = −0.24, b 21 = −0.31, c 1 = 0.2.
(39) Figure 1 shows the absolute stability regions of the three methods with the values of the parameters just reported. Note that the considered OEJDP method has stability region that in some places contains that of the NEJDP method. However, the latter has a larger area and includes a much larger part of the real axis.

Numerical tests
In this section we perform numerical tests to evaluate the advantages of NEJDP methods derived in this paper, comparing them with the OEJDP and the classic ones. Since in Conte et al. (2021) numerical tests were conducted exclusively on scalar ODEs, we now also consider non-scalar problems. Specifically, we will focus on the methods versions (39), (38) and (37), which guarantee optimal stability properties as demonstrated in the previous subsection. We evaluate for each method the absolute error committed at the final grid point T , comparing the computed numerical solution with the exact one (if known), or with that determined by the Matlab function ode15s, requiring maximum accuracy. In addition, we report tables with the estimate of the order of each method, calculated using the formula where cd(h) = −log 10 (absolute error) represents the number of correct digits obtained with step-size h.
In this case, we do not report tables with errors due to the application of the considered methods. In fact, we will evaluate the greater accuracy of the new methods obtained in this paper through the relative application on non-scalar problems. Now, we want to evaluate the better stability of the NEJDP methods, both with respect to OEJDP and classical schemes. On the other hand, for the same values of the parameter and of the step-size, both the OEJDP method and the classical one produce numerical solutions that either explode, or in any case that have a completely different trend with respect the exact one.

Euler problem
Euler problem (Euler 1758) is represented by the following non-stiff ODEs system: In this case, from Table 1, it is evident the greater accuracy of the NEJDP methods compared to the OEJDP and the classic ones. In fact, the error of the new methods is smaller than the others by at least two orders of magnitude. Specifically, there are two orders of difference between the new methods and the old ones, and even three orders of difference between the new methods and the classic peer schemes.
As we expected, however, the order of NEJDP methods remains equal to two, despite the fact that the second stage is calculated with order three. In fact, as mentioned in the Introduction, peer methods are characterized by having all stages with the same properties in terms of accuracy and stability. Therefore, if one stage is determined with order two, and the other with order three, the overall order of the method is still equal to two. We can ascertain this by looking at Table 2. However, note that, although the order tends to two, it is initially close to three. Thus, the higher accuracy in the second-stage calculation is responsible for the higher overall accuracy of the NEJDP methods.

Brusselator model
The last ODEs system we consider is the Brusselator (Nicolis and Prigogine 1977) problem y 2 (t) =3y 1 (t) − y 1 (t) 2 y 2 (t),  In this case, Tables 3 and 4 confirm the observations already made previously. Therefore, even for non-scalar problems with non-constant and non-linear Jacobian, the NEJDP methods work well.

Non-linear Burgers equation
To show the trend of the new methods also on ODEs systems of higher dimensions than those discussed so far, let us consider a Partial Differential Equation (PDE) and discretize it in space using the method of lines. In particular, we consider the following formulation of the non-linear Burgers PDE (Burgers 1948): We set periodic boundary conditions and as initial conditions, then employing an order four finite differences spatial semidiscretization, both for the second-and first-order spatial derivatives, fixing a grid of M = 2 5 spatial intervals. The semi-discretized is space problem, which is an ODEs system of dimension M, takes the following form: Here, L 1 and L 2 are constant matrices having the discretization coefficients of the secondand first-order spatial derivatives, respectively. The Jacobian is the non-constant matrix J y(t) = L 1 − L 2 y(t).
We solve the problem (44) with the NEJDP methods, also comparing the CPU time they take with that of the Matlab ode45 routine, with parity of errors. The errors are evaluated, as before, with respect to the reference solution determined by the ode15s function applied to the semi-discretized ODEs system (44), requiring maximum accuracy. The corresponding results are reported in Table 5.
Since ode45 has an adaptive step-size control, the comparison is not very fair, as NEJDP methods are with fixed step-size. In fact, from Table 5, note that the number of grid points taken by ode45 is significantly lower than that of the NEJDP methods. However, the latter still manage to be competitive (they take on average a CPU time four times bigger than ode45, but still low). To show that thanks to the step control ode45 saves a lot of time, in Table 6, we have reported the minimum step-size h min used by it during the integration. Then, note that even though ode45 takes far fewer grid points than NEJDP methods, it is sometimes forced to choose an even smaller step than them to obtain a comparable error.
Furthermore, note that in the far right column of Table 5, we have reported the results obtained by applying the NEJDP methods by freezing the Jacobian as J = L 1 . This choice is usual when the system of ODEs to be solved is of the form (44), and it is also known that its stiffness is mainly concentrated in the first term. In this case, NEJDP methods are advantageous with respect to ode45 in terms of computational time, with parity of errors. This happens despite the choice J = L 1 does not look at accuracy, and therefore could be done in a better way, leading to much lower errors. This suggests that, with adequate Jacobian freezing strategies and adaptive step-size control, NEJDP methods could definitely perform even better than they already do. Indeed, in this test, we have shown that, only with frozen (non-optimized) Jacobian and without adaptive step, the NEJDP methods are already competitive with ode45. Surely, with adaptive step-size control, the computing times would be lowered even more drastically.

Conclusions and future perspectives
In this paper, we have proposed a new class of explicit peer methods with improved stability and accuracy properties. Furthermore, we have proved that such methods are very advantageous both with respect to the classic ones, and to those derived in Conte et al. (2021), which have exactly the same computational cost. Furthermore, we have definitely improved the work Conte et al. (2021), deriving the OEJDP methods also in non-scalar case and generalizing their order conditions. Moreover, we have shown how it is possible to generalize in an appropriate way the technique proposed in Ixaru (2012) for Runge-Kutta methods also on peer methods, which unlike the former are characterized by stages that in the current interval depend on those determined in the previous one. In fact, we considered the stages in the intervals [t n , t n+1 ] and [t n−1 , t n ] to be affected by error. Due to the shape of peer methods, we think that the Jacobian dependency must be applied in some way to all points of the grid in [0, t n ], to obtain even better stability and accuracy properties.
Finally, the results reported in the last numerical test lead us to think about investigating a variable step-size formulation of the NEJDP methods derived in this paper, also with a Jacobian freezing strategy. article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.