Frequency evaluation for adapted peer methods

In this paper, we consider exponentially fitted peer methods for the numerical solution of first order differential equations and we investigate how the frequencies can be tuned in order to obtain the maximal benefit. We will show that the key is analyzing the error’s behavior. Formulae for optimal frequencies are computed. Numerical experiments show the properties of the proposed algorithm.

the same accuracy with significantly larger step sizes, permitting to develop efficient and accurate numerical methods.
Despite the enormous efforts on Ordinary Differential Equations (ODEs) problems, the critical question is how to choose frequencies to maximize the benefits of EF methods. In Ixaru et al. (2001Ixaru et al. ( , 2002Ixaru et al. ( , 2003 for EF multistep methods, an algorithm was presented, and the answer was provided to numerically tune frequencies as optimally as possible for systems of ODEs. The optimal frequencies were derived by exploiting the expression of the error, depending on some combination of higher order derivative of the solution. The present paper considers a general class of EF two-step peer methods (Conte et al. 2018(Conte et al. , 2020 of the form (2) for numerical integration of ODEs (1) with oscillatory solution. We investigate Ixaru's frequency evaluation algorithm for adapted EF peer methods. We take the behavior of the leading term of the error as a starting point of the whole investigation and we develop frequency formulae.
The paper is organized as follows: in Sect. 2, we give a short review of EF peer methods. Section 3 investigates some fundamental concepts of frequency evaluation and explores the development of the EF peer method frequency assessment, while Sect. 4 describes the estimation of derivatives. Some numerical examples in Sect. 5 illustrate the efficiency of the new frequency assessment algorithm, and, lastly, some conclusions are drawn in Sect. 6. Appendix A contains the expressions of the leading term of the error.

EF peer methods
We consider systems of ODEs of the form where f : R × R d → R d is sufficiently smooth to ensure that the solution exists and it is unique. In addition, we assume that the problem (1) possesses an oscillatory solution. Let {t n := t 0 + nh, n = 0, . . . , N } be a discretization of the interval [t 0 , T ] with fixed stepsize h > 0, 0 ≤ c 1 ≤ · · · ≤ c s = 1 be fixed distinct nodes and define t ni = t n + c i h, i = 1, . . . , s. Two-step peer methods can be formulated as follows: where h > 0 is the stepsize, I is the identity matrix of dimension d, A = a i j As c s = 1, Y ns is the approximation of the solution at grid point t n+1 .
The order of EF peer methods has been analyzed by Schmitt and Weiner (2004) by introducing the condition and providing the following Theorem.
Theorem 1 If AB( p + 1) is verified, the implicit s-stage peer method (2) has order of consistency p.

Corollary 1
The peer method (2) has order p ≥ s if We briefly review the procedure of construction of EF peer methods introduced in Conte et al. (2018Conte et al. ( , 2020 in the following Algorithm.

Algorithm 1 : Construction of EF peer methods
Start: Choose the fitting space: with μ = i ω where ω ∈ R is problem's oscillating frequency. Use: Define the linear difference operator where the vector w contains all the coefficients of the method (2). The coefficient matrices A, B and R of peer methods (2) are obtained by imposing that the operator (6) annihilates when the function y belongs to the fitting space (5).
Step 1: Construct the classic moments and the dimensionless classic moments: which have the form: Step 2: We now look for the maximum value M that ensures the compatibility of the system which is equivalent to annihilating the difference operator (6) on polynomials with a degree less or equal to M − 1.
Therefore, we may construct an s-order peer method if M = s + 1, due to the Theorem 1 and Corollary 1.
Step 3: Construct the formal expressions of Step 4: Construct the possible expressions for the fitting space (5) taking into account that M = s + 1 and the self-consistency condition K + 2P = M − 3 has to be verified. The number of stages s and the dimension M of the system (9) are of different parities, so the number K + 1 = s − 1 − 2P of classic functions in the fitting space is odd or even if s is even or odd.
Step 5: Solve formally the linear systems with Z dependent coefficients. The numeric values of a i j and b i are computed either for real or imaginary μ-values.
Step 6: Compute the leading term of the local truncation error for CL and EF peer methods (Choices for K and P summarized in Table 1): where we denote D the derivative with respect to time.

Basic elements
The question of how frequencies need to be tuned to achieve maximum benefit from EF methods has not been answered for a long time. Ixaru et al. (2001Ixaru et al. ( , 2002Ixaru et al. ( , 2003, and Vanden Berghe et al. (2001a, b) have proposed a frequency evaluation algorithm for EF multistep methods and EF Runge-Kutta methods (respectively) which enable to tune of the frequency μ in the way that the principal local truncation error vanishes. Therefore, the analysis of the error behavior is a required step. We refer to the articles Ixaru et al. (2001Ixaru et al. ( , 2002Ixaru et al. ( , 2003 and Vanden Berghe et al. (2001a, b) for technical information and even some practical points and we restrict the discussion to only three relevant cases: A0, A1, and A2 algorithms that exactly incorporate all linear combinations from the reference set of functions.
where K is specified by the considered method. The choice A0 covers the purely algebraic classical method; Algorithm A1 is especially of importance whenever the solution exhibits a purely exponential behavior, while Algorithm A2 describes oscillatory solutions if μ is strictly imaginary. As said, for the investigation of frequency, analyzing the behavior of the error is a necessary stage. We compute the expression of the leading term of the error (lte) for these algorithms by using the general procedure described in Sect. 2 and appropriate options for K and P in the fitting space based on the above algorithms.

Frequency evaluation for adapted EF peer methods
To start, we consider EF peer methods derived in Algorithm 1 to the scalar equation y (t) = f (t, y(t)), then the leading term of their error (lte) assumes the form When we apply three types of algorithms A0, A1 and A2 to the scalar equation (1) and assume appropriate options for K and P in the fitting space summarized in Table 1 and by attention to Eq. (13) the lte is derived as follows.
• Algorithm A0: in this case K = s and P = −1, therefore lte is described as follows • Algorithm A1: in this case K = s − 1 and P = 0 and lte is given by the following expressions assuming s even or odd.
-if s = 2, then K = 1, P = 0 • Algorithm A2: for this case with K = s − 2 and P = 0, assuming that s is even or odd, the following expressions provide lte.
s = 2, then K = 0, P = 0 Our principal purpose is to determine the μ value that guarantees maximum accuracy when the classical A0 case is replaced by one of the two EF Algorithms. For the calculation of parameter μ, we restrict discussion to the expressions for the lte of A1 and A2 cases especially for s = 2 and s = 3: In all cases, we see that the lte consists of a product of three factors, i.e.
• a general h 3 (in the case s = 2) or general h 4 factor (in the case s = 3), • a function depending on Z which tends to the classical value when Z tends to zero.
• a factor that involves two derivatives of the solution.
The important thing for our studies is the different behavior of the third factor. This differential factor can make a real difference in accuracy. Let us then introduce the following functionals: If a μ exists such that D 1 identically vanishes on the quoted interval then the version A1 corresponding to that μ will be exact. The reason is that identically vanishing D 1 is equivalent to looking at the differential equation y (4) (t) − μ 2 y (t) = 0 and y (5) (t) − μ 2 y (3) (t) = 0 for s = 2 and s = 3, respectively.
In general, no constant μ can be found such that D 1 identically vanishes but it makes sense to address the problem of finding that value of μ which ensures that the values of D 1 are kept as close to zero as possible for t n in the considered interval. When D 1 is held close to zero, the optimal μ are given by In A1 case the frequency must be real.
For algorithm A2, the same considerations can be repeated for D 2 . The reason is that identically vanishing D 2 are equivalent to looking at the differential equation y (3) (t) − μ 2 y (t) = 0 and y (4) (t) − μ 2 y (t) = 0 for s = 2 and s = 3, respectively. When D 2 is held close to zero, the optimal μ are given by The frequencies are either real or imaginary if Algorithm A2 is chosen. We have succeeded in proposing formulae for the optimal μ = i ω (ω is frequency) value in (20)-(21). Having found the optimal value for μ, we have the optimal frequency ω.

Estimation of the derivatives
The evaluation of the optimal μ value and then optimal frequency ω requires knowing the total derivatives appearing in the expressions of the formulae (20)-(21). At first, it seems like a very simple task: the first-order derivative is equivalent to the right-hand sides of the Eq. (1) i.e. f (t, y(t)), after that, it is very straightforward to calculate the higher-order derivatives. This technique works well on many problems, but in Ixaru et al. (2002), authors demonstrate this should be avoided on stiff problems. They also demonstrate that it is sufficient to use finite difference approximations of the derivatives. Since the expressions of the formulae (20)-(21) contain derivatives of orders three and four (for s = 2) and four and five (for s = 3), we will estimate the derivatives in each integration point t n with five points finite difference formulae for s = 2 and six points finite difference formulae for s = 3.
Since in Van de Vyver (2005), the author has shown approximation of the derivatives by five points finite difference formulae, in this work we just mention the formulae of the six points finite difference with data at t n−4 , t n−3 , t n−2 , t n−1 , t n and t n+1 for the input.
The data y n−4 , y n−3 , y n−2 , y n−1 , y n for the input points are the value of the numerical solution at t n−4 , t n−3 , t n−2 , t n−1 , t n . The estimation for the y n+1 at t n+1 is determined by Milne-Simpson two-step formulae (Ixaru et al. 2002) as follows: It is necessary that we make this approximation with a method that has a higher or the same order as the EF peer method. Therefore it is sufficient to choose Milne-Simpson's twostep formulae. We only use the result for the calculation of the derivatives and not for the propagation of the solution.

Case s = 2
If t n is a root of y (t), the algorithm A1 is not defined, while A2 is not defined when t n is a root of y (t). In a very special case, when t n is a root of both y (t) and y (t) EF algorithms are not suitable and so the classic A0 Algorithm must be activated. In general, a logical way to choose between A1 and A2 involves comparing |y (t)| and |y (t)|. If |y (t n )| < |y (t n )| then A1 is selected, otherwise A2.

Case s = 3
In this case, if t n is a root of y (3) (t), the algorithm A1 is not defined, whereas A2 when t n is a root of y (t), not defined. It happens that both y (t) and y (3) (t) change the sign with the same time interval, neither the two EF algorithms are not suitable, so the classic A0 should be enabled. In general, as in the previous case, there is a reasonable way to choose between A1 and A2 involves comparing |y (t)| and |y (3)

Numerical experiments
In this section we present numerical experiments showing the behaviour of the new "optimal EF peer methods", with Z -dependent coefficients whose are computed for optimal μ-values in (20)-(21). In the following examples, we apply the described optimal EF peer methods and the classic and EF peer methods to solve the test cases. We compare errors of the optimal EF implicit peer methods with errors of classic and EF implicit peer methods of Conte et al. (2020) in Examples 1 and 2, and we also compare achieved results from Example 3 with reported results by Ixaru et al. (2002). The error will be estimated as the infinite norm of the difference between the numerical solution and the exact solution at the endpoint and reported in the tables. In addition, we will use the following notation to represent the used numerical methods: • CL = classic, • EF = exponentially fitted, • IM P2 = implicit peer method of order 2.
In the following examples, we use the notation reported in Conte et al. (2020). We consider s = 2. In this case K = 0 and P = 0. We fix c 1 = 0, c 2 = 1. According to Z = μ 2 h 2 = −ω 2 h 2 , the numerical values of a i j and b i are computed either for real or imaginary μ-values. The corresponding optimal EF IM peer method and EF IM peer method are: Example 1 Let us consider the Prothero-Robinson problem y (t) = λ ( y(t) − sin(51 t) ) + (51) cos(51 t), t ∈ 0, π 2 ,  Optimal EF IM P2 ω op 9.49e−09 3.08e−10 9.62e−12 whose exact solution is The oscillating behavior of the exact solution leads us to utilize the EF methods with the parameter μ characterizing the functions belonging to the fitting space equal to μ = i ω. So the problem is integrated by the EF peer methods, where the parameter ω is chosen equal to the frequency of the exact solution, i. e. ω = 50.
We also consider the case in which the oscillatory frequency ω is not known exactly. Therefore by finding the frequency from the formulae (20)-(21) and denoting with ω op , we employ the EF peer methods whose coefficients are computed in correspondence of a frequency ω op and μ = i ω op value. We used the initial conditions and carried out with CL, EF IM peer methods, and optimal EF IM peer methods, whose algorithms are constructed by the procedure described in Algorithm 1. We consider interval [0, π 2 ] with different grid points N = 320, 640, 1280. Tables 2, 3 represent the absolute errors from the considered methods, for λ = −1 (non stiff case) and λ = −10 6 (stiff case). We see that the optimal EF IM peer method works much better than CL and is close to EF IM peer methods, irrespective of whether the problem is stiff or non-stiff.
For additional confirmation, we present some graphs for both cases λ = −1 and λ = −10 6 . In Fig. 1, we depict the variation of ω op at each integration point for the problem when the A2 algorithm is chosen. It is seen from Fig. 2, (as expected) the obtained ω op is close to ω = 50. It is instructive to mention that Fig. 2 shows the efficiency curve for this problem obtained by the CL IM, EF IM, and optimal EF IM peer methods.

Example 2 Consider the Prothero-Robinson problem
where exact solution is y(t) = sin(101 t) = sin(100 t) cos(t) + cos(100 t) sin(t).  The oscillating behavior of the exact solution leads us to utilize the EF methods with the parameter μ = i ω. This system is integrated by the EF peer methods, where the parameter ω is chosen equal to the frequency of the exact solution, i.e. ω = 100. In this example, we also consider the case in which the oscillatory frequency ω is not known exactly and we integrate the system by the optimal EF peer methods with ω op . In detail, we utilized the initial conditions and performed the experiments with CL and EF IM peer methods, as well as optimal EF IM peer methods, whose algorithms are constructed in the process outlined in Algorithm 1. We examine the interval [0, π 2 ] with various grid points N = 320, 640, 1280. The absolute errors from the considered approaches are listed in the Tables 4 and 5 for λ = −1 (non-stiff case) and λ = −10 6 (stiff case).
Whether the problem is stiff or non-stiff, the optimal EF IM peer methods perform much better than CL IM and are close to EF IM peer methods. We offer some graphs for both scenarios λ = −1 and λ = −10 6 for extra validation. It is instructive to mention that Fig. 3 shows the efficiency curve for this problem obtained by the CL, EF, and optimal EF IM peer methods. In Fig. 4, we depict the variation of ω op at each integration point for the problem when the A2 algorithm is chosen. As predicted, it is obvious from Fig. 4, the obtained ω op is close to ω = 100.  Example 3 Consider the following test case with the exact solution This is a simple differential equation, and it helps to illustrate some exciting aspects when (24) is approached by optimal EF IM peer methods (with Z -dependent coefficients whose are computed for optimal μ-values) introduced in this work. We can construct the optimal μ value for both A1 and A2 algorithms by using formulae (20) and (21). For instance, when s = 2, by using formulae (20) the optimal μ value for A1 is the constant representation  while by using formulae (21) the optimal μ 2 value for A2 is It follows that by using the A1 algorithm when s = 2, we don't have all of the data to construct the optimal μ value. For this reason, in this test case, we construct the solution by the A2 algorithm along with the whole interval. In Fig. 5, we depict the variation of optimal μ 2 at each integration point for the problem. As predicted, it is evident from Fig. 5, μ 2 0.024 when t = 10, s = 2, h = 0.0125, as the theoretically expected value 0.024 is also given from ϕ 2 (10).
This test case has been employed by Ixaru et al. (2002). They used EF multistep methods for Eq. (24). In Table 6, we report the absolute errors at t = 1, t = 5 and t = 10, with the fixed stepsize h = 0.0500, 0.0250 and h = 0.0125 by optimal EF IM peer methods and compare them with reported results in Ixaru et al. (2002). From this table, we observe that for s = 2, the optimal EF IM peer methods have the same accuracy behavior concerning with EF multistep methods (Ixaru et al. 2002).

Conclusion
In this paper, we applied EF peer methods for the numerical solution of first-order ODEs and examined the problem of how the frequencies should be tuned to obtain the maximal benefit from the exponential fitting versions. To answer this question, we analyzed the error behavior of EF peer methods. We have succeeded in proposing formulae for optimal μ values. Under this condition and with the determination of optimal μ values, we achieved the "optimal EF peer methods". The introduced methods were tested on some examples, and the efficiency of optimal EF peer methods was shown.