Implementing Adams Methods with Preassigned Stepsize Ratios

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.


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 value problem y x f x, y x , y x 0 y 0 . 1.1 These codes, and other interesting multistep methods, are based on the replacement of f x, y x by an interpolating polynomial p x , and the approximation y x h y x x h x f t, y t dt y x x h x p t dt.

Mathematical Problems in Engineering
When the stepsize is constant, the explicit k-step Adams-Bashforth method can be written in terms of the backward differences of f, where the coefficients γ * j are obtained by the recurrence The formulation of the implicit k-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 PECE fixed-point implementation 3, page 144 , starting from an Adams-Bashforth formula.The order of the k-step ABM pair is k 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 k-order pair based on k-step Adams-Bashforth and k − 1 -step Adams-Moulton with local extrapolation 3, pages 114-115 .The γ * j 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 h n x n 1 − x n 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 k-step Adams-Bashforth predictor is and the k-step Adams-Moulton corrector can be expressed as The coefficients can be calculated easily by the recurrences The calculation of g j n requires a double loop in terms of j and q: 1.9 The computational cost of the recalculation of g j n in each step is the major disadvantage of the variable-stepsize Adams code when the cost is measured in CPU time.

The Stepsize Ratios
The coefficients "β n " and "g n " are scalar and they are independent of the function f.Let us denote the ratios Then it is well known that dependence on the points x i in 1.7 and 1.9 is actually a dependence only on the ratios r i .

The Coefficients in Terms of the Stepsize Ratios
Let us study the number of ratios needed for each coefficient.
Proposition 2.1.For j ≥ 1 the coefficient β j n only depends on the j ratios r n−1 , r n−2 , . . ., r n−j , and for j ≥ 2 the coefficient g j n only depends on the j − 1 ratios r n−1 , r n−2 , . . ., r n−j 1 .
Proof.We only need to apply the relation with p > m Mathematical Problems in Engineering to deduce that the fraction in 1.7 can be rewritten as and only depends on the j ratios r n−1 , r n−2 , . . ., r n−j .Using an inductive reasoning, starting with the constant value β 0 n 1, it follows that if j ≥ 1, then A similar proof can be developed for g j n .Applying 2.2 to the fraction in 1.9 ,

2.5
From this equation and the starting values the inductive reasoning shows that, if j ≥ 2, c j,q x n 1 Θ j,q r n−1 , r n−2 , . . ., r n−j 1 , 2.7 and so and the proof is complete.

The Number of Coefficients in Terms of the Number of Ratios
Let us fix only λ ≥ 2 options to choose the ratios, that is, As the ratios r i belong to a finite subset, then the coefficients β j n and g j n belong to finite spaces of possible events.As there are only a finite number of "β n " and "g n " 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 Ψ j and Θ j,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 k 1.For simplicity we ignore the constants β 0 n 1, g 0 n 1, and g 1 n 1/2.For the predictor 1.5 we require from β 1 n to β k−1 n .According to 1.8 , the value β k n is not necessary for φ k n 1 in the corrector 1.6 .From Proposition 2.1 it follows that there are only λ j different values for β j n , so the total length is Again from Proposition 2.1, if j ≥ 2, the coefficient g j n can only take λ j−1 different values.But in this case g k n is used in the implicit corrector 1.6 , so it is necessary to compute g 2 n , . . ., g k n and 2.10 is also valid for the "g n " coefficients.We can formulate the first lemma.
Lemma 2.2.The "β n " and "g n " coefficients needed in the k-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 k-step pair is actually the k-order pair with local extrapolation, we can choose the same local error estimator for both pairs, which leads to

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 for decreasing or increasing the order, respectively.The value of β k n is necessary for φ k 1 n 1 in LE 1 n , and also g k 1 n , so Lemma 2.2 can be updated to the following lemma.
Lemma 2.3.The "β n " and "g n " 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 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 k 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 k 1 MO is reached, the estimator LE 1 n is not necessary, and then we do not need to calculate and store β MO−1 n and g MO n .This fact produces the final lemma.
Lemma 2.4.The "β n " and "g n " coefficients needed for a Variable-Step Variable-Order Adams k-step code, with order up to k 1 MO, can be stored in two arrays of length 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 ω m in 2.9 .

The Fixed Ratios Strategy
Usually, the implementation of the variable-stepsize Adams method includes a restriction in the ratios of the kind of 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 :

3.2
The fixed ratios implementation that we propose selects the ratio as the largest ω m in 2.9 lower to r in 3.2 : The pathologic case ω m > r for 1 ≤ m ≤ λ can be easily avoided by including ω 1 ω 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 k 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 r i in 3.3 is lower than r 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 "β n " and "g n " 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 e 0.6 and another one with high eccentricity e 0.9.
In the figures we compare the decimal logarithm of the CPU time versus the decimal logarithm of the error, with tolerances from 10 −3 to 10 −13 , except in the Lorenz problem in which tolerances were restricted from 10 −7 to 10 −13 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 ω m must be close to 1, even if the stepsize may vary a lot in the integration Section 3.1 .
ii Coefficients "β n " must be calculated from 1.7 in each step and coefficients "g n " must be precalculated once and read in each step Section 3.2 .

The Ratios in the Classical Implementation
We present some results of five revolutions of a highly eccentric orbit e 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.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.

Mathematical Problems in Engineering
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 −3 to 10 −13 for all the proposed problems.
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.

The Difference between the "β n " and the "g n " Coefficients
The coefficients "β n " and "g n " are archived in two arrays of the same size, so the CPU time required when we search for a particular β j n or g j n must be the same.However, there is an important difference when we compute them from their algorithms, because while β j n is calculated by means of the single loop 1.7 , a double loop is required in 1.9 for g j n .
We present different FRABM integrations for the Arenstorf problem, all of them with λ 4 prefixed ratios, r i ∈ {0.5, 0.9, 1.1, 2}, and maximum order MO 13.We change the way to evaluate computing or reading the "β n "/"g n " 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 r i 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 "g n " strategies diamonds and squares in Figure 3 won to both computing "g n " strategies stars and triangles .In Figure 3 the best implementation mode resulted to be reading "β n " and "g n ".In other problems we tested computing "β n " and reading "g n " was slightly superior.Taking into account all the tests, both strategies were virtually tied, and they were doubtlessly more efficient than both computing "g n " modes.But as storing "β n " requires an important amount of RAM, we discard the reading mode for "β n ", so from now on the FRABM implementation will compute "β n " from 1.7 in each step, and only precalculate, store, and read in each step "g n ".In this way RAM requirements are halved.

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 ω 1 ω < 0.9 is required for the rejected steps, and the maximum ω λ Ω will be used specially in the first steps, Mathematical Problems in Engineering  as the first stepsize h 0 must be very small because the code starts with the minimum order k 1 2. Computing and not reading the "g n " 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 r i / ∈ {ω 1 , ω 2 , . . ., ω λ } forces the code to compute the next k "g n " coefficients from 1.9 to eliminate r i 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 ω 2 1 we obtained the basic doubling and halving FRABM3 method with λ 3 prefixed ratios, r i ∈ {0.5, 1, 2}.Lemma 2.4 shows that for  maximum order MO 13 this method only requires 265719 "g n " 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 ω 1 0.5 as the unique prefixed ratio for rejections, and looked for ω m ≥ 0.9 if m ≥ 2. In fact, taking into account the numerical results of our tests, for λ > 3 we decided to fix ω 2 0.9, and ω λ−1 1.1 as its companion.In particular, the method with λ 4 ratios, r i ∈ {0.5, 0.9, 1.1, 2} from now on denoted FRABM4, diamonds in Figure 4 , performed significantly better than the one with r i ∈ {0.5, 0.98, 1.03, 2} obtained with percentiles 33 and 66 in Table 2 stars in Figure 4 .
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 ω 3 1 to the ratios in FRABM4.Then for FRABM5 we selected r i ∈ {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 r i ∈ {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 .
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, r i ∈ {0.5, 0.9, 1, 1.1, 2} performed better than λ 6, r i ∈ {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 Mathematical Problems in Engineering  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.
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.
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.

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.
Let us interpret the figures.At a point where a pair of methods has the same ycoordinate error the method with graph α units to the left in the x-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.15units in the x-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 e 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 e 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 f x, y is expensive to evaluate.DOPRI853 shows its potential in both Two-Body, Arenstorf,and Lorenz problems, all of them with a cheap function f x, y .The difference between DOPRI853 and VSVOABM is about 0.4-0.6 units in the x-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 f x, y 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   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 f x, y 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 f x, y 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.

9 Figure 1 : 9 Figure 2 :
Figure 1: Decimal logarithm of the stepsizes in a highly eccentric two-body problem.

Figure 6 :
Figure 6: FRABM5, maximum orders MO 12 and MO 13 in the Pleiades problem.

6 Figure 7 :Figure 8 :
Figure 7: Low cost, FRABM3 with MO 13, FRABM4 with MO 11 and FRABM5 with MO 10 in a middle eccentric Two-Body problem.

Table 1 :
Percentiles 25, 50, and 75 of the ratios for the highly eccentric e 0.9 two-body problem.

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

Table 3 :
RAM cost in terms of the number of prefixed ratios λ for maximum orderMO 11.

Table 4 :
RAM cost in terms of the number of prefixed ratios λ and the maximum order MO.
gains about 20% CPU time with respect to DOPRI853.The Pleiades problem is expensive to evaluate because it is formed by 28 long firstorder 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 .