Implementing Adams Methods with Preassigned Stepsize Ratios

Abstract

Runge-Kutta and Adams methods are the most popular codes to solve numerically nonstiff ODEs. The Adams methods are useful to reduce the number of function calls, but they usually require more CPU time than the Runge-Kutta methods. In this work we develop a numerical study of a variable step length Adams implementation, which can only take preassigned step-size ratios. Our aim is the reduction of the CPU time of the code by means of the precalculation of some coefficients. We present several numerical tests that show the behaviour of the proposed implementation.

1. The Variable-Stepsize Adams Method

The Adams methods (Bashforth and Adams [1], Moulton [2])are the most popular multistep formulas to solve a non-stiff initial valueproblem

𝑦^′(𝑥)=𝑓(𝑥,𝑦(𝑥)),𝑦(𝑥₀)=𝑦₀. (1.1) These codes, and other interesting multistep methods, are based on the replacement of 𝑓(𝑥,𝑦(𝑥)) by an interpolating polynomial 𝑝(𝑥), and the approximation

𝑦(𝑥+ℎ)=𝑦(𝑥)+_𝑥^𝑥+ℎ𝑓(𝑡,𝑦(𝑡))𝑑𝑡≃𝑦(𝑥)+_𝑥^𝑥+ℎ𝑝(𝑡)𝑑𝑡. (1.2)

When the stepsize is constant, the explicit 𝑘-step Adams-Bashforth method can be written in terms of the backward differences of 𝑓,

𝑦_𝑛+1=𝑦_𝑛+ℎ^𝑘−1_𝑗=0𝛾^∗_𝑗∇^𝑗𝑓_𝑛, (1.3) where the coefficients 𝛾^∗_𝑗 are obtained by the recurrence

𝛾^∗₀=1,𝛾^∗_𝑗=1−^𝑗_𝑖=1𝛾^∗_𝑗−𝑖𝑖+1. (1.4)

The formulation of the implicit 𝑘-step Adams-Moulton is straightforward, as can be seen, for example, by Lambert in [3, Section 3.9]. The implicit equation is usually solved using a 𝑃𝐸𝐶𝐸 fixed-point implementation [3, page 144], starting from an Adams-Bashforth formula. The order of the 𝑘-step ABM pair is 𝑘+1 (an exhaustive study of the order of a predictor-corrector pair can be seen in [3, Section 4.2] ), and it is equivalent to the 𝑘-order pair based on 𝑘-step Adams-Bashforth and (𝑘−1)-step Adams-Moulton with local extrapolation [3, pages 114-115].

The 𝛾^∗_𝑗 coefficients do not depend on the step length, so they remain constant in the whole integration and need only be computed once. However, they do vary when the stepsize ℎ_𝑛=𝑥_𝑛+1−𝑥_𝑛 is not constant, so they must be recalculated in each step of the integration.

There are lots of different ways to implement a variable-stepsize variable-order Adams code (actually one different for each developer). For example, codes as DEABM [4], VODE [5] and LSODE [6, 7] include some controls to diminish the computational cost in some situations by keeping the stepsize constant, or by avoiding the growth of the order. Nowadays, symbolic programs such as Matlab and Mathematica as well as their numerical solvers of ODEs ode113 [8] and NDSolve [9, pages 961–972,1068,1215-1216], are very popular. We do not want to restrict our study to a particular code, so we will use in the tests a purer implementation of the variable-stepsize Adams code, developed by ourselves, based on the Krogh divided differences formulation [10].

Krogh developed a way to compute the Adams formulas in terms of modified divided differences. Following the notation of Hairer et al. [11, pages 397–400], the variable-stepsize explicit 𝑘-step Adams-Bashforth predictor is

𝑝_𝑛+1=𝑦_𝑛+ℎ_𝑛^𝑘−1_𝑗=0𝑔_𝑗(𝑛)𝛽_𝑗(𝑛)𝜙_𝑗(𝑛), (1.5) and the 𝑘-step Adams-Moulton corrector can be expressed as

𝑦_𝑛+1=𝑝_𝑛+1+ℎ_𝑛𝑔_𝑘(𝑛)𝜙_𝑘(𝑛+1). (1.6)

The coefficients can be calculated easily by the recurrences

𝛽₀(𝑛)=1,𝛽_𝑗(𝑛)=𝛽_𝑗−1𝑥(𝑛)^𝑛+1−𝑥_{𝑛−𝑗+1}𝑥_𝑛−𝑥_𝑛−𝑗,𝜙₀(𝑛)=𝑓_𝑛,𝜙_𝑗(𝑛)=𝜙_𝑗−1(𝑛)−𝛽_𝑗−1(𝑛−1)𝜙_𝑗−1(𝑛−1). (1.7)(1.8) The calculation of 𝑔_𝑗(𝑛) requires a double loop in terms of 𝑗 and 𝑞:

𝑐_0,𝑞(𝑥_𝑛+11)=𝑞,𝑐_𝑗,𝑞(𝑥_𝑛+1)=𝑐_{𝑗−1,𝑞}(𝑥_𝑛+1)−𝑐_{𝑗−1,𝑞+1}(𝑥_𝑛+1)ℎ_𝑛𝑥_𝑛+1−𝑥_{𝑛−𝑗+1},𝑔_𝑗(𝑛)=𝑐_𝑗,1(𝑥_𝑛+1). (1.9) The computational cost of the recalculation of 𝑔_𝑗(𝑛) in each step is the major disadvantage of the variable-stepsize Adams code when the cost is measured in CPU time.

2. The Stepsize Ratios

The coefficients “𝛽(𝑛)” and “𝑔(𝑛)” are scalar and they are independent of the function 𝑓. Let us denote the ratios

𝑟_𝑖=ℎ_𝑖+1ℎ_𝑖. (2.1) Then it is well known that dependence on the points 𝑥_𝑖 in (1.7) and (1.9) is actually a dependence only on the ratios 𝑟_𝑖.

2.1. The Coefficients in Terms of the Stepsize Ratios

Let us study the number of ratios needed for each coefficient.

Proposition 2.1. For 𝑗≥1 the coefficient 𝛽_𝑗(𝑛) only depends on the 𝑗 ratios 𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_𝑛−𝑗, and for 𝑗≥2 the coefficient 𝑔_𝑗(𝑛) only depends on the 𝑗−1 ratios 𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1}.

Proof. We only need to apply the relation (with 𝑝>𝑚) 𝑥_𝑝−𝑥_𝑚=^{𝑝−𝑚−1}_𝑗=0ℎ_𝑚+𝑗=^{𝑝−𝑚−1}_𝑗=0ℎ_𝑚^𝑗−1_𝑖=0𝑟_𝑚+𝑖=ℎ_𝑚⋅𝜑_𝑝−𝑚(𝑟_𝑝−2,𝑟_𝑝−3,…,𝑟_𝑚) (2.2) to deduce that the fraction in (1.7) can be rewritten as 𝑥_𝑛+1−𝑥_{𝑛−𝑗+1}𝑥_𝑛−𝑥_𝑛−𝑗=ℎ_{𝑛−𝑗+1}⋅𝜑_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1})ℎ_𝑛−𝑗⋅𝜑_𝑗(𝑟_𝑛−2,𝑟_𝑛−3,…,𝑟_𝑛−𝑗)=ℎ_𝑛−𝑗⋅𝑟_𝑛−𝑗⋅𝜑_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1})ℎ_𝑛−𝑗⋅𝜑_𝑗(𝑟_𝑛−2,𝑟_𝑛−3,…,𝑟_𝑛−𝑗)=𝜓_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_𝑛−𝑗) (2.3) and only depends on the 𝑗 ratios 𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_𝑛−𝑗. Using an inductive reasoning, starting with the constant value 𝛽₀(𝑛)=1, it follows that if 𝑗≥1, then 𝛽_𝑗(𝑛)=Ψ_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_𝑛−𝑗). (2.4) A similar proof can be developed for 𝑔_𝑗(𝑛). Applying (2.2) to the fraction in (1.9), ℎ_𝑛𝑥_𝑛+1−𝑥_{𝑛−𝑗+1}=ℎ_{𝑛−𝑗+1}⋅∏^𝑗−2_𝑖=0𝑟_{𝑛−𝑗+𝑖+1}ℎ_{𝑛−𝑗+1}⋅𝜑_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1})=∏^𝑗−2_𝑖=0𝑟_{𝑛−𝑗+𝑖+1}𝜑_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1})=𝜃_𝑗(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1}). (2.5) From this equation and the starting values 𝑐_0,𝑞(𝑥_𝑛+11)=𝑞,𝑐_1,𝑞(𝑥_𝑛+11)=𝑞(𝑞+1), (2.6) the inductive reasoning shows that, if 𝑗≥2, 𝑐_𝑗,𝑞(𝑥_𝑛+1)=Θ_𝑗,𝑞(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1}), (2.7) and so 𝑔_𝑗(𝑛)=𝑐_𝑗,1(𝑥_𝑛+1)=Θ_𝑗,1(𝑟_𝑛−1,𝑟_𝑛−2,…,𝑟_{𝑛−𝑗+1}), (2.8) and the proof is complete.

2.2. The Number of Coefficients in Terms of the Number of Ratios

Let us fix only 𝜆≥2 options to choose the ratios, that is,

𝑟_𝑖∈{𝜔₁,𝜔₂,…,𝜔_𝜆},∀𝑖≥0. (2.9)

As the ratios 𝑟_𝑖 belong to a finite subset, then the coefficients 𝛽_𝑗(𝑛) and 𝑔_𝑗(𝑛) belong to finite spaces of possible events. As there are only a finite number of “𝛽(𝑛)” and “𝑔(𝑛)” coefficients, they can be precalculated, and read at the start of the numerical integration. In particular, the major computational disadvantage of the variable-stepsize Adams code (the double loop in the computation of (1.9) in each integration step) is removed.

The importance of Proposition 2.1 is not the calculation of the functions Ψ_𝑗 and Θ_𝑗,1, but the existence of those functions and the number of the ratios for each one. Let us count the number of different possible coefficients in terms of 𝜆 and the order 𝑘+1. For simplicity we ignore the constants 𝛽₀(𝑛)=1, 𝑔₀(𝑛)=1, and 𝑔₁(𝑛)=1/2.

For the predictor (1.5) we require from 𝛽₁(𝑛) to 𝛽_𝑘−1(𝑛). According to (1.8), the value 𝛽_𝑘(𝑛) is not necessary for 𝜙_𝑘(𝑛+1) in the corrector (1.6). From Proposition 2.1 it follows that there are only 𝜆^𝑗 different values for 𝛽_𝑗(𝑛), so the total length is

𝜆+𝜆²+⋯+𝜆^𝑘−1=𝜆^𝑘−𝜆.𝜆−1 (2.10)

Again from Proposition 2.1, if 𝑗≥2, the coefficient 𝑔_𝑗(𝑛) can only take 𝜆^𝑗−1 different values. But in this case 𝑔_𝑘(𝑛) is used in the implicit corrector (1.6), so it is necessary to compute 𝑔₂(𝑛),…,𝑔_𝑘(𝑛) and (2.10) is also valid for the “𝑔(𝑛)” coefficients. We can formulate the first lemma.

Lemma 2.2. The “𝛽(𝑛)” and “𝑔(𝑛)” coefficients needed in the 𝑘-step pair (1.5) and (1.6) can be stored in two arrays of length (2.10).

A variable-stepsize algorithm needs an estimator of the error to change the stepsize. This estimator can be computed in different ways. We follow Lambert's [3, section 4.4] and extend them to variable-stepsizes.

As the 𝑘-step pair is actually the 𝑘-order pair with local extrapolation, we can choose the same local error estimator for both pairs, which leads to

LE₀(𝑛)=ℎ_𝑛(𝑔_𝑘(𝑛)−𝑔_𝑘−1(𝑛))𝜙_𝑘(𝑛+1). (2.11)

In practice the variable-stepsize Adams method is also endowed with the capability of changing its order, so it is necessary to obtain some local error estimators for decreasing or increasing the order. By modifying the previous equation we obtain the estimators

LE₋₁(𝑛)=ℎ_𝑛(𝑔_𝑘−1(𝑛)−𝑔_𝑘−2(𝑛))𝜙_𝑘−1(𝑛+1),LE₊₁(𝑛)=ℎ_𝑛(𝑔_𝑘+1(𝑛)−𝑔_𝑘(𝑛))𝜙_𝑘+1(𝑛+1), (2.12) for decreasing or increasing the order, respectively. The value of 𝛽_𝑘(𝑛) is necessary for 𝜙_𝑘+1(𝑛+1) in LE₊₁(𝑛), and also 𝑔_𝑘+1(𝑛), so Lemma 2.2 can be updated to the following lemma.

Lemma 2.3. The “𝛽(𝑛)” and “𝑔(𝑛)” coefficients needed for the predictor (1.5), the corrector (1.6), and the estimators (2.11) and (2.12) can be stored in two arrays of length 𝜆+𝜆²+⋯+𝜆^𝑘−1+𝜆^𝑘=𝜆^𝑘+1−𝜆.𝜆−1 (2.13)

The CPU time used by the code will increase drastically if we let the order grow without control. For this reason it is necessary to impose a maximum order. In double precision (8 bytes per data, 15-16 decimal digits) the maximum order is usually 𝑘+1=13 (Lambert [3, page 144]). For our purposes we let this maximum as a parameter MO (Maximum Order). We must emphasize that, when the maximum order 𝑘+1=MO is reached, the estimator LE₊₁(𝑛) is not necessary, and then we do not need to calculate (and store) 𝛽_MO−1(𝑛) and 𝑔_MO(𝑛). This fact produces the final lemma.

Lemma 2.4. The “𝛽(𝑛)” and “𝑔(𝑛)” coefficients needed for a Variable-Step Variable-Order Adams 𝑘-step code, with order up to 𝑘+1=MO, can be stored in two arrays of length 𝜆+𝜆²+⋯+𝜆^MO−2=𝜆^MO−1−𝜆.𝜆−1 (2.14)

The amount of RAM needed to store and read the coefficients can be an important disadvantage of a fixed ratios implementation. In the following section we search some useful choices of 𝜆 and 𝜔_𝑚 in (2.9).

3. The Fixed Ratios Strategy

Usually, the implementation of the variable-stepsize Adams method includes a restriction in the ratios of the kind of

𝜔≤𝑟_𝑖≤Ω (3.1) that appears in the stability of generic multistep variable-stepsize methods (Hairer et al. [11, page 403]) and in the convergence of the variable-stepsize variable-order Adams method (Shampine [12]). Then, the stepsizes ratio is obtained by adapting equation (391a) by Butcher in [13, page 293]:

𝑟=max(𝜔,min(Ω,0.9(tol)est^1/(𝑘+1))). (3.2)

The fixed ratios implementation that we propose selects the ratio as the largest 𝜔_𝑚 in (2.9) lower to 𝑟 in (3.2):

𝑟_𝑖=max{𝜔_𝑚≤𝑟}. (3.3) The pathologic case 𝜔_𝑚>𝑟 for 1≤𝑚≤𝜆 can be easily avoided by including 𝜔₁=𝜔 as the minimum of the prefixed ratios.

For the numerical study we programmed a code, denoted as VSVOABM, with a classical implementation. VSVOABM uses (3.2) for the stepsize selection, and lets the order vary freely up to 𝑘+1=13. In addition, we denote by FRABM (fixed ratios Adams-Bashforth-Moulton) a code based on VSVOABM, changing (3.2) by (3.3) to select the new stepsize. As the prefixed ratio 𝑟_𝑖 in (3.3) is lower than 𝑟 in (3.2), the stepsize selected in FRABM will be smaller than the one that VSVOABM would choose from the same estimator, so the error of FRABM is expected to be lower (of course using more steps). As the coefficients “𝛽(𝑛)” and “𝑔(𝑛)” do not depend on the problem, they were precalculated only once, stored in a file and put in RAM at the start of the integration in FRABM. Then, in each step, they are read from memory instead of being calculated. As FRABM was programmed from VSVOABM, the comparison of the two will show the behaviour of the classical strategy let the ratios free and calculate the coefficients in each step versus the new fix a finite number of ratios and read the coefficients from RAM in each step.

The authors presented in [14] a large number of experiments. In this paper we summarize them. First of all we tested the following problems. In all of them the variation of the stepsize is crucial.

(1)Five revolutions of the classical Two-Body problem used in many relevant books, such as Butcher [13, pages 4–7], Dormand [15, pages 86–90], or Shampine [16, pages 90–95]. The variation of the stepsize grows with the eccentricity of the ellipsis. We test a pair of different problems, one with middle eccentricity 𝑒=0.6 and another one with high eccentricity 𝑒=0.9.(2)An Arenstorf problem [17] with the popular orbit printed by Butcher in [13, page 9]. of Hairer et al. [11, page 130, Figure 0.1] shows that a variable-stepsize strategy is very important in this problem.(3)The Lorenz problem in the interval [0,16], extremely sensitive to errors in the first steps [11, pages 120, 245].(4)A fictitious seven-body problem, denoted Pleiades [11, pages 245-246]. In this problem the function 𝑓(𝑥,𝑦) is expensive to evaluate.

In the figures we compare the decimal logarithm of the CPU time versus the decimal logarithm of the error, with tolerances from 10⁻³ to 10⁻¹³, except in the Lorenz problem in which tolerances were restricted from 10⁻⁷ to 10⁻¹³ (see Hairer et al. [11, page 245] for an explanation). Nowadays most of these problems can be solved quickly in a standard computer, so CPU time was measured with extreme precision in “uclocks” (about 840 nanoseconds, less than 1 microsecond).

The main conclusions were the following.

(i)Ratios 𝜔_𝑚 must be close to 1, even if the stepsize may vary a lot in the integration (Section 3.1).(ii)Coefficients “𝛽(𝑛)” must be calculated from (1.7) in each step and coefficients “𝑔(𝑛)” must be precalculated once and read in each step (Section 3.2).(iii)FRABM methods with 𝜆=5 ratios, 𝑟_𝑖∈{0.5,0.9,1,1.1,2}, and maximum orders MO=10 and MO=12 are two good methods with very different RAM requirements (Section 3.3).

3.1. The Ratios in the Classical Implementation

We present some results of five revolutions of a highly eccentric orbit 𝑒=0.9 of the Two-Body problem with VSVOABM. In Figure 1 we see significant changes in the stepsize. In the beginning of a revolution (perigee) it is about one hundred times less than in the middle of a revolution (apogee). However, the ratios in Figure 2 remain near to 1 everywhere, even those less than 0.9 from the rejected steps.

Figure 1 

Decimal logarithm of the stepsizes in a highly eccentric two-body problem.

Figure 2 

Ratios in a highly eccentric two-body problem.

We have tested some other eccentricities in the Two-Body problem. As the eccentricity reduces to 0, the changes in the stepsizes in Figure 1 diminish. In any case, the structure of Figure 2 shows little change.

In Table 1 we order the ratios and select the percentiles 25, 50 (the median), and 75 of the accepted steps for different tolerances. The percentiles approach to unity as the tolerance reduces, probably because the selected orders increase with small tolerances (see in Shampine [16, page 262, Figure 5.3]). In Table 2 we show the behaviour of the ratios in a global sample obtained with tolerances from 10⁻³ to 10⁻¹³ for all the proposed problems.

Table 1 

Percentiles 25, 50, and 75 of the ratios for the highly eccentric 𝑒=0.9 two-body problem.

Table 2 

Relevant percentiles of the ratios of accepted steps in a global sample.

We must conclude that a good choice of a finite number of prefixed ratios must select them near to 1, even for problems with big changes in the stepsize.

3.2. The Difference between the "𝛽(𝑛)” and the “𝑔(𝑛)” Coefficients

The coefficients “𝛽(𝑛)” and “𝑔(𝑛)” are archived in two arrays of the same size, so the CPU time required when we search for a particular 𝛽_𝑗(𝑛) or 𝑔_𝑗(𝑛) must be the same. However, there is an important difference when we compute them from their algorithms, because while 𝛽_𝑗(𝑛) is calculated by means of the single loop (1.7), a double loop is required in (1.9) for 𝑔_𝑗(𝑛).

We present different FRABM integrations for the Arenstorf problem, all of them with 𝜆=4 prefixed ratios, 𝑟_𝑖∈{0.5,0.9,1.1,2}, and maximum order MO=13. We change the way to evaluate (computing or reading) the (“𝛽(𝑛)”/“𝑔(𝑛)”) coefficients, which produces four different implementation modes. In all of them the coefficients must be the same, and also the numerical results and errors (round-off errors can produce in (3.2) slightly different estimators est and different selections of the ratio 𝑟_𝑖 in (3.3), but it does not happen often). As the CPU time can be different for the four FRABM methods, the comparison must provide parallel graphs.

In all the tests we made, both reading “𝑔(𝑛)” strategies (diamonds and squares in Figure 3) won to both computing “𝑔(𝑛)” strategies (stars and triangles). In Figure 3 the best implementation mode resulted to be reading “𝛽(𝑛)” and “𝑔(𝑛)”. In other problems we tested computing “𝛽(𝑛)” and reading “𝑔(𝑛)” was slightly superior. Taking into account all the tests, both strategies were virtually tied, and they were doubtlessly more efficient than both computing “𝑔(𝑛)” modes. But as storing “𝛽(𝑛)” requires an important amount of RAM, we discard the reading mode for “𝛽(𝑛)”, so from now on the FRABM implementation will compute “𝛽(𝑛)” from (1.7) in each step, and only precalculate, store, and read in each step “𝑔(𝑛)”. In this way RAM requirements are halved.

Figure 3 

FRABM, computing versus reading “$\beta (n)$” and “$g(n)$” in the Arenstorf problem.

3.3. Proposing Some FRABM Algorithms

Now we will develop some concrete FRABM algorithms. It is impossible to reproduce all of the experiments we made in [14] comparing sets of ratios obtained by statistical and intuitive processes, so in this paper we only present the final conclusions.

First of all, we looked for the minimum and maximum ratios 𝜔 and Ω. These ratios will not be used often, but they are indispensable. The minimum 𝜔₁=𝜔<0.9 is required for the rejected steps, and the maximum 𝜔_𝜆=Ω will be used specially in the first steps, as the first stepsize ℎ₀ must be very small because the code starts with the minimum order 𝑘+1=2. Computing (and not reading) the “𝑔(𝑛)” coefficients directly from (1.9) only in these special cases looks like an interesting idea, because it sets two preassigned ratios free. However, we must note that using one nonpreassigned ratio 𝑟_𝑖∉{𝜔₁,𝜔₂,…,𝜔_𝜆} forces the code to compute the next 𝑘 “𝑔(𝑛)” coefficients from (1.9) to eliminate 𝑟_𝑖 from (2.8), and the computational cost will be increased.

In our experiments there were no clear winners, but maybe the pair 𝜔=0.5 and Ω=2 was more robust than other compared ratios, so we assumed the criterion of Dormand [15, page 73] and Butcher [13, page 293] and selected this pair.

With the addition of the ratio 𝜔₂=1 we obtained the basic doubling and halving FRABM3 method with 𝜆=3 prefixed ratios, 𝑟_𝑖∈{0.5,1,2}. Lemma 2.4 shows that for maximum order MO=13 this method only requires 265719 “𝑔(𝑛)” coefficients and 2.03 MB of RAM to store them in double precision.

In the vast majority of occasions the stepsize is accepted, not rejected. That is why we decided to keep 𝜔₁=0.5 as the unique prefixed ratio for rejections, and looked for 𝜔_𝑚≥0.9 if 𝑚≥2. In fact, taking into account the numerical results of our tests, for 𝜆>3 we decided to fix 𝜔₂=0.9, and 𝜔_𝜆−1=1.1 as its companion. In particular, the method with 𝜆=4 ratios, 𝑟_𝑖∈{0.5,0.9,1.1,2} (from now on denoted FRABM4, diamonds in Figure 4), performed significantly better than the one with 𝑟_𝑖∈{0.5,0.98,1.03,2} obtained with percentiles 33 and 66 in Table 2 (stars in Figure 4).

Figure 4 

Two different $\lambda =4\text{,}$$\text{MO}=13$ FRABM methods in the Lorenz problem.

We tested several different strategies for 𝜆=5 ratios, including the one with percentiles 25, 50, and 75 from Table 2. However, the best strategy turned out to be the simplest: adding 𝜔₃=1 to the ratios in FRABM4. Then for FRABM5 we selected 𝑟_𝑖∈{0.5,0.9,1,1.1,2}. The methods with ratios obtained by means of a statistical criterion from Table 2 performed clearly worse than others selected intuitively in all the tests we made.

For no specific reason, we add the ratio 1.05 to develop a 𝜆=6 method, with 𝑟_𝑖∈{0.5,0.9,1,1.05,1.1,2}. In Figure 5 we compared the behaviour of FRABM methods with 3≤𝜆≤6 ratios, all of them with maximum order MO=11 (the RAM requirements of the methods are in Table 3).

Table 3 

RAM cost in terms of the number of prefixed ratios 𝜆 for maximum order 𝑀𝑂=11.

Figure 5 

FRABM$\lambda$ with $3\le \lambda \le 6\text{,}$$\text{MO}=11$ in the Lorenz problem.

In Figure 5, and in all the problems we tested, FRABM5 was the best of the compared methods. We cannot assure that 𝜆=5 ratios are a better choice than 𝜆=6; we assure that in all of our tests the combination 𝜆=5, 𝑟_𝑖∈{0.5,0.9,1,1.1,2} performed better than 𝜆=6, 𝑟_𝑖∈{0.5,0.9,1,1.05,1.1,2}. This situation did not change when we raised 𝜆. We compared methods with 3≤𝜆≤10 ratios, all of them with the same MO=9 (reduced through the fault of the RAM cost for 𝜆=10), and FRABM5 was always the winner. It is surprising because FRABM5, MO=9 only requires an amount of 0.75 MB of RAM.

We cannot conclude yet that the best method is FRABM5, because when 𝜆<5 bigger maximum orders MO are attainable with lower RAM cost. However, it appears that when 𝜆>5 the improvement of the error is less than the growth of the CPU time (and of course the RAM requirements), so we restrict our study to 3≤𝜆≤5.

We divided the methods into three categories in terms of their RAM requirements (see Table 4), with the usual restriction MO≤13 for 𝜆=3,4, and MO=12 for FRABM5. As the case 𝜆=3 reaches the maximum order MO=13 with low-cost, FRABM3 was considered in the three categories.

Table 4 

RAM cost in terms of the number of prefixed ratios 𝜆 and the maximum order MO.

We rejected the case MO=13 for FRABM5 because it needs the enormous amount of 465.66 MB of RAM. Besides, the RAM requirements worked against the efficacy of the code. Figure 6 shows the comparative of FRABM5 methods with maximum orders MO=12 and MO=13 in the Pleiades problem. In all the problems under consideration the case MO=13 was not competitive with MO=12, specially for large tolerances.

Figure 6 

FRABM5, maximum orders $\text{MO}=12$ and $\text{MO}=13$ in the Pleiades problem.

The behaviour of the methods shown in Figure 6 is an exception. The efficacy of FRABM𝜆 methods in Table 4 was consistent with its maximum order, that is, it was better with higher MO. However the middle cost methods were between low-cost and high-cost methods for all the problems in almost every tolerance, so we discarded them because they were not the winners in any circumstance.

Figures 7 and 8 show representative examples of the performance of the methods in all the problems under consideration. FRABM3 with MO=13 was clearly the poorest of the compared methods in both categories. For low- and high-cost FRABM5, MO=10,12 won FRABM4, MO=11,13, respectively. Thus, combinations of 𝜆=5 prefixed ratios {0.5,0.9,1,1.1,2} with maximum orders MO=10 and MO=12 are officially chosen as our two FRABM candidates to solve general non-stiff ODEs.

Figure 7 

Low cost, FRABM3 with $\text{MO}=13$, FRABM4 with $\text{MO}=11$ and FRABM5 with $\text{MO}=10$ in a middle eccentric Two-Body problem.

Figure 8 

High cost, FRABM3 with $\text{MO}=13$, FRABM4 with $\text{MO}=13$, and FRABM5 with $\text{MO}=12$ in the Arenstorf problem.

4. The Behaviour of FRABM5

In this section we compare the pair of selected FRABM5 methods (green line with stars for MO=10, red line with diamonds for MO=12) and the classical implementation VSVOABM with maximum order MO=13 (blue line with squares). As the most popular representative of Runge-Kutta methods, we include in the tests the version of the Dormand and Prince DOPRI853 code [18] downloaded from E. Hairer's website (black line with triangles). Figure 9 shows the efficacy of these four methods when the computational cost is measured in CPU time. We also present in Figure 10 the efficacy graphics with the computational cost measured in number of function calls.

Figure 9 

Error versus CPU time. FRABM5 $\text{MO}=10$: green line with stars, FRABM5 $\text{MO}=12$: red line with diamonds, VSVOABM $\text{MO}=13$: blue line with squares, DOPRI853: black line with triangles.

Figure 10 

Error versus evaluations. FRABM5 $\text{MO}=10$: green line with stars, FRABM5 $\text{MO}=12$: red line with diamonds, VSVOABM $\text{MO}=13$: blue line with squares, and DOPRI853: black line with triangles.

Let us interpret the figures. At a point where a pair of methods has the same 𝑦-coordinate (error) the method with graph 𝛼 units to the left in the 𝑥-coordinate gains 100⋅(1−10^−𝛼)% in the cost.

Taking into account only the Adams methods and CPU time, both FRABM5 codes clearly beat VSVOABM in both Two-Body, Arenstorf, Lorenz and Pleiades problems. The difference fluctuates between 0.1–0.15 units in the 𝑥-coordinate (a gain of 20%–30% in CPU time). Let us remark that FRABM5, MO=10 requires the insignificant amount of 3.73 MB of RAM and is the most efficient multistep method for large tolerances. The case MO=12 is more efficient for small tolerances, but its RAM cost grows to 93.13 MB.

When comparing the number of function calls, VSVOABM does not start very well in large tolerances, but it can take advantage of its superior maximum order MO=13 to win FRABM5, MO=10 in small tolerances. In the Two-Body problem with 𝑒=0.9 the case MO=12 is slightly better than VSVOABM even in these small tolerances. In the rest of the tests VSVOABM is the winner. We made some experiments (not included in this paper) to compare the efficacy in number of function calls of VSVOABM with FRABM4 and FRABM5, all of them with MO=13. There was not a clear winner in our tests, which allows us to assert that the main variable that produces differences among the Adams methods is not the implementation mode, but the maximum order MO. This can be seen in the figures of Two-Body with 𝑒=0.6 and Pleiades.

It is well known that Runge-Kutta methods, like DOPRI853, can hardly compete with multistep methods when comparing the number of evaluations. This can clearly be seen in all of the experiments. However, as Runge-Kutta methods do not need to recalculate their coefficients in each step, they are very superior in CPU time unless the function 𝑓(𝑥,𝑦) is expensive to evaluate. DOPRI853 shows its potential in both Two-Body, Arenstorf,and Lorenz problems, all of them with a cheap function 𝑓(𝑥,𝑦). The difference between DOPRI853 and VSVOABM is about 0.4–0.6 units in the 𝑥-axis, so the gain of DOPRI853 oscillates between 60%–75%, which means that it is about 2.5–4 times faster than VSVOABM. It is true that the fixed ratios strategy improves the efficacy of the Adams method, but the difference between DOPRI853 and VSVOABM was so great that both FRABM5 methods are clearly worse than DOPRI853 in these problems.

The function 𝑓(𝑥,𝑦) in the Pleiades is more complicated, and its evaluation requires a significant amount of CPU time. DOPRI853 beats ABMVSVO in all tolerances, but the difference is small, which allows both FRABM5 methods to win in this problem. For small tolerances FRABM5, MO=12 gains about 20% CPU time with respect to DOPRI853.

The Pleiades problem is expensive to evaluate because it is formed by 28 long first-order differential equations. Hairer et al. [11, page 427] indicated that the efficacy of Adams methods gets worse as the number of differential equations grows. To study this question in our fixed ratios implementation we add two new expensive to evaluate ODEs extracted from [11].

(5) A discretized Brusselator with diffusion [11, pages 248-249], constituted by 882 short first-order differential equations.(6) A discretized hanging Rope problem. In this problem we let the dimension 𝑛 in equation (10.8) of [11, pages 247-248] as a free parameter. Then this problem consists of 2𝑛long first-order differential equations.

The comparison in number of function calls in the Brusselator problem does not show any new significant conclusion. However, in the Rope examples FRAMB5, MO=12 beats ABMVSVO even for small tolerances.

In the Brusselator problem, DOPRI853 clearly beats all FRABM5 methods in CPU time, which corresponds to the indication of Hairer, Nørsett, and Wanner. In the Rope problem we can also see the negative influence of the dimension in the behaviour of the multistep methods. With 8 first-order equations FRABM5 methods are more efficient than DOPRI853, specially with MO=12, which gains up to 40% in small tolerances. For 24 equations FRABM5, MO=12 still remains in the first position in these tolerances, but now all of the methods have similar results. However, DOPRI853 takes the lead for 80 equations, and the difference grows for 320 equations. We must conclude that the statement in Hairer et al. [11, page 427] about dimension and efficacy of the Adams code is also valid for the fixed ratios implementation.

VSVOABM is not far from FRABM5 methods in the Brusselator and Rope problems, but it only beats MO=10 in small tolerances. FRABM5, MO=12 is superior in all tolerances.

We sum up the conclusions of this efficacy section.

(1)When the CPU time is used to measure the efficacy, DOPRI853 is the best choice when 𝑓(𝑥,𝑦) is short or it has a high number of components. However, the Pleiades and Rope problems show that FRABM5 can beat DOPRI853 in CPU time when 𝑓(𝑥,𝑦) has a moderate number of long components. For these kinds of ODEs we propose our fixed ratios implementation, particularly FRABM5 with maximum order MO=10 for large tolerances and MO=12 for small ones. VSVOABM was not the best method in any test.(2)When comparing the number of evaluations, fixed ratios and VSVO strategies are comparable; only the maximum order MO produces differences between large and small tolerances. In this chapter, specially relevant if the evaluation of 𝑓(𝑥,𝑦) has an associated monetary cost, all the multistep methods won DOPRI853. In particular, both FRABM5 were slightly better than VSVOABM in large tolerances. In small tolerances FRABM5, MO=12 was the leader in some problems and VSVOABM in others, while the case MO=10 did not win in any test. FRABM4, MO=13 is a choice with an acceptable RAM cost that performed slightly better than FRABM5, MO=12 in these small tolerances, and was virtually tied to ABMVSVO.