# Implementing Adams Methods with Preassigned Stepsize Ratios

Authors: David J. López; José G. Romay
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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.


---

## Body

## 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)y′(x)=f(x,y(x)),y(x0)=y0.
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(1.2)y(x+h)=y(x)+∫xx+hf(t,y(t))dt≃y(x)+∫xx+hp(t)dt.When the stepsize is constant, the explicitk-step Adams-Bashforth method can be written in terms of the backward differences of f,(1.3)yn+1=yn+h∑j=0k-1γj*∇jfn,
where the coefficients γj* are obtained by the recurrence(1.4)γ0*=1,γj*=1-∑i=1jγj-i*i+1.The formulation of the implicitk-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 hn=xn+1-xn 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(1.5)pn+1=yn+hn∑j=0k-1gj(n)βj(n)ϕj(n),
and the k-step Adams-Moulton corrector can be expressed as(1.6)yn+1=pn+1+hngk(n)ϕk(n+1).The coefficients can be calculated easily by the recurrences(1.7)β0(n)=1,βj(n)=βj-1(n)xn+1-xn-j+1xn-xn-j,(1.8)ϕ0(n)=fn,ϕj(n)=ϕj-1(n)-βj-1(n-1)ϕj-1(n-1).
The calculation of gj(n) requires a double loop in terms of j and q:(1.9)c0,q(xn+1)=1q,cj,q(xn+1)=cj-1,q(xn+1)-cj-1,q+1(xn+1)hnxn+1-xn-j+1,gj(n)=cj,1(xn+1).
The computational cost of the recalculation of gj(n) 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 “β(n)” and “g(n)” are scalar and they are independent of the function f. Let us denote the ratios(2.1)ri=hi+1hi.
Then it is well known that dependence on the points xi in (1.7) and (1.9) is actually a dependence only on the ratios ri.


### 2.1. The Coefficients in Terms of the Stepsize Ratios

Let us study the number of ratios needed for each coefficient.Proposition 2.1.
Forj≥1 the coefficient βj(n) only depends on the j ratios rn-1,rn-2,…,rn-j, and for j≥2 the coefficient gj(n) only depends on the j-1 ratios rn-1,rn-2,…,rn-j+1.Proof.
We only need to apply the relation (withp>m)
(2.2)xp-xm=∑j=0p-m-1hm+j=∑j=0p-m-1hm∏i=0j-1rm+i=hm·φp-m(rp-2,rp-3,…,rm)
to deduce that the fraction in (1.7) can be rewritten as
(2.3)xn+1-xn-j+1xn-xn-j=hn-j+1·φj(rn-1,rn-2,…,rn-j+1)hn-j·φj(rn-2,rn-3,…,rn-j)=hn-j·rn-j·φj(rn-1,rn-2,…,rn-j+1)hn-j·φj(rn-2,rn-3,…,rn-j)=ψj(rn-1,rn-2,…,rn-j)
and only depends on the j ratios rn-1,rn-2,…,rn-j. Using an inductive reasoning, starting with the constant value β0(n)=1, it follows that if j≥1, then
(2.4)βj(n)=Ψj(rn-1,rn-2,…,rn-j).
A similar proof can be developed for gj(n). Applying (2.2) to the fraction in (1.9),
(2.5)hnxn+1-xn-j+1=hn-j+1·∏i=0j-2rn-j+i+1hn-j+1·φj(rn-1,rn-2,…,rn-j+1)=∏i=0j-2rn-j+i+1φj(rn-1,rn-2,…,rn-j+1)=θj(rn-1,rn-2,…,rn-j+1).
From this equation and the starting values
(2.6)c0,q(xn+1)=1q,c1,q(xn+1)=1q(q+1),
the inductive reasoning shows that, if j≥2,
(2.7)cj,q(xn+1)=Θj,q(rn-1,rn-2,…,rn-j+1),
and so
(2.8)gj(n)=cj,1(xn+1)=Θj,1(rn-1,rn-2,…,rn-j+1),
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,(2.9)ri∈{ω1,ω2,…,ωλ},∀i≥0.As the ratiosri belong to a finite subset, then the coefficients βj(n) and gj(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 Proposition2.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, g0(n)=1, and g1(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(2.10)λ+λ2+⋯+λk-1=λk-λλ-1.Again from Proposition2.1, if j≥2, the coefficient gj(n) can only take λj-1 different values. But in this case gk(n) is used in the implicit corrector (1.6), so it is necessary to compute g2(n),…,gk(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 thek-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)LE0(n)=hn(gk(n)-gk-1(n))ϕk(n+1).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(2.12)LE-1(n)=hn(gk-1(n)-gk-2(n))ϕk-1(n+1),LE+1(n)=hn(gk+1(n)-gk(n))ϕk+1(n+1),
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 gk+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
(2.13)λ+λ2+⋯+λk-1+λk=λk+1-λλ-1.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 usuallyk+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 gMO(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
(2.14)λ+λ2+⋯+λMO-2=λMO-1-λλ-1.The amount of RAM needed to store and read the coefficients can be an important disadvantage of afixed ratios implementation. In the following section we search some useful choices of λ and ωm in (2.9).


## 2.1. The Coefficients in Terms of the Stepsize Ratios

Let us study the number of ratios needed for each coefficient.Proposition 2.1.
Forj≥1 the coefficient βj(n) only depends on the j ratios rn-1,rn-2,…,rn-j, and for j≥2 the coefficient gj(n) only depends on the j-1 ratios rn-1,rn-2,…,rn-j+1.Proof.
We only need to apply the relation (withp>m)
(2.2)xp-xm=∑j=0p-m-1hm+j=∑j=0p-m-1hm∏i=0j-1rm+i=hm·φp-m(rp-2,rp-3,…,rm)
to deduce that the fraction in (1.7) can be rewritten as
(2.3)xn+1-xn-j+1xn-xn-j=hn-j+1·φj(rn-1,rn-2,…,rn-j+1)hn-j·φj(rn-2,rn-3,…,rn-j)=hn-j·rn-j·φj(rn-1,rn-2,…,rn-j+1)hn-j·φj(rn-2,rn-3,…,rn-j)=ψj(rn-1,rn-2,…,rn-j)
and only depends on the j ratios rn-1,rn-2,…,rn-j. Using an inductive reasoning, starting with the constant value β0(n)=1, it follows that if j≥1, then
(2.4)βj(n)=Ψj(rn-1,rn-2,…,rn-j).
A similar proof can be developed for gj(n). Applying (2.2) to the fraction in (1.9),
(2.5)hnxn+1-xn-j+1=hn-j+1·∏i=0j-2rn-j+i+1hn-j+1·φj(rn-1,rn-2,…,rn-j+1)=∏i=0j-2rn-j+i+1φj(rn-1,rn-2,…,rn-j+1)=θj(rn-1,rn-2,…,rn-j+1).
From this equation and the starting values
(2.6)c0,q(xn+1)=1q,c1,q(xn+1)=1q(q+1),
the inductive reasoning shows that, if j≥2,
(2.7)cj,q(xn+1)=Θj,q(rn-1,rn-2,…,rn-j+1),
and so
(2.8)gj(n)=cj,1(xn+1)=Θj,1(rn-1,rn-2,…,rn-j+1),
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,(2.9)ri∈{ω1,ω2,…,ωλ},∀i≥0.As the ratiosri belong to a finite subset, then the coefficients βj(n) and gj(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 Proposition2.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, g0(n)=1, and g1(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(2.10)λ+λ2+⋯+λk-1=λk-λλ-1.Again from Proposition2.1, if j≥2, the coefficient gj(n) can only take λj-1 different values. But in this case gk(n) is used in the implicit corrector (1.6), so it is necessary to compute g2(n),…,gk(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 thek-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)LE0(n)=hn(gk(n)-gk-1(n))ϕk(n+1).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(2.12)LE-1(n)=hn(gk-1(n)-gk-2(n))ϕk-1(n+1),LE+1(n)=hn(gk+1(n)-gk(n))ϕk+1(n+1),
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 gk+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
(2.13)λ+λ2+⋯+λk-1+λk=λk+1-λλ-1.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 usuallyk+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 gMO(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
(2.14)λ+λ2+⋯+λMO-2=λMO-1-λλ-1.The amount of RAM needed to store and read the coefficients can be an important disadvantage of afixed ratios implementation. In the following section we search some useful choices of λ and ωm 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)ω≤ri≤Ω
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)r=max(ω,min(Ω,0.9(tolest)1/(k+1))).The fixed ratios implementation that we propose selects the ratio as the largestωm in (2.9) lower to r in (3.2):(3.3)ri=max{ωm≤r}.
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 ri 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 classicalTwo-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.(2)
AnArenstorf 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)
TheLorenz problem in the interval [0,16], extremely sensitive to errors in the first steps [11, pages 120, 245].(4)
A fictitious seven-body problem, denotedPleiades [11, pages 245-246]. In this problem the function f(x,y) is expensive to evaluate.In the figures we compare the decimal logarithm of theCPU 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).(iii)
FRABM methods withλ=5 ratios, ri∈{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 orbite=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 Figure1 diminish. In any case, the structure of Figure 2 shows little change.In Table1 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.Table 1
Percentiles 25, 50, and 75 of the ratios for the highly eccentrice=0.9 two-body problem.
ToleranceP25P50P7510-30.941.041.1710-40.961.111.1410-50.951.071.1110-60.931.041.0910-70.931.031.0710-80.941.021.0610-90.951.021.0410-100.951.021.0410-110.961.011.0310-120.971.011.0310-130.971.011.02Table 2
Relevant percentiles of the ratios of accepted steps in a global sample.
P10P20P25P30P33P40P50P60P66P70P75P80P900.940.960.970.980.981.001.011.021.031.031.041.041.09We 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 "β(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 gj(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 gj(n).We present different FRABM integrations for the Arenstorf problem, all of them withλ=4 prefixed ratios, ri∈{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 ri 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, bothreading “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.Figure 3
FRABM, computing versus reading “β(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 ω1=ω<0.9 is required for the rejected steps, and the maximum ωλ=Ω will be used specially in the first steps, as the first stepsize h0 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 ri∉{ω1,ω2,…,ωλ} forces the code to compute the next k “g(n)” coefficients from (1.9) to eliminate ri 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, ri∈{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, ri∈{0.5,0.9,1.1,2} (from now on denoted FRABM4, diamonds in Figure 4), performed significantly better than the one with ri∈{0.5,0.98,1.03,2} obtained with percentiles 33 and 66 in Table 2 (stars in Figure 4).Figure 4
Two differentλ=4, 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 ω3=1 to the ratios in FRABM4. Then for FRABM5 we selected ri∈{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 ri∈{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 MO=11.
λNumber of coefficientsRAM cost3295230.23 MB43495242.67 MB5244140518.63 MB61209323492.26 MBFigure 5
FRABMλ with 3≤λ≤6, MO=11 in the Lorenz problem.In Figure5, 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, ri∈{0.5,0.9,1,1.1,2} performed better than λ=6, ri∈{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 Table4), 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.
λMORAMCategory3132.03 MBLow cost4112.67 MBLow cost5103.73 MBLow cost41210.67 MBMiddle cost51118.63 MBMiddle cost41342.67 MBHigh cost51293.13 MBHigh costWe rejected the caseMO=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 ordersMO=12 and MO=13 in the Pleiades problem.The behaviour of the methods shown in Figure6 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.Figures7 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 withMO=13, FRABM4 with MO=11 and FRABM5 with MO=10 in a middle eccentric Two-Body problem.Figure 8
High cost, FRABM3 withMO=13, FRABM4 with MO=13, and FRABM5 with MO=12 in the Arenstorf problem.


## 3.1. The Ratios in the Classical Implementation

We present some results of five revolutions of a highly eccentric orbite=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 Figure1 diminish. In any case, the structure of Figure 2 shows little change.In Table1 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.Table 1
Percentiles 25, 50, and 75 of the ratios for the highly eccentrice=0.9 two-body problem.
ToleranceP25P50P7510-30.941.041.1710-40.961.111.1410-50.951.071.1110-60.931.041.0910-70.931.031.0710-80.941.021.0610-90.951.021.0410-100.951.021.0410-110.961.011.0310-120.971.011.0310-130.971.011.02Table 2
Relevant percentiles of the ratios of accepted steps in a global sample.
P10P20P25P30P33P40P50P60P66P70P75P80P900.940.960.970.980.981.001.011.021.031.031.041.041.09We 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 "β(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 gj(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 gj(n).We present different FRABM integrations for the Arenstorf problem, all of them withλ=4 prefixed ratios, ri∈{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 ri 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, bothreading “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.Figure 3
FRABM, computing versus reading “β(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 ω1=ω<0.9 is required for the rejected steps, and the maximum ωλ=Ω will be used specially in the first steps, as the first stepsize h0 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 ri∉{ω1,ω2,…,ωλ} forces the code to compute the next k “g(n)” coefficients from (1.9) to eliminate ri 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, ri∈{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, ri∈{0.5,0.9,1.1,2} (from now on denoted FRABM4, diamonds in Figure 4), performed significantly better than the one with ri∈{0.5,0.98,1.03,2} obtained with percentiles 33 and 66 in Table 2 (stars in Figure 4).Figure 4
Two differentλ=4, 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 ω3=1 to the ratios in FRABM4. Then for FRABM5 we selected ri∈{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 ri∈{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 MO=11.
λNumber of coefficientsRAM cost3295230.23 MB43495242.67 MB5244140518.63 MB61209323492.26 MBFigure 5
FRABMλ with 3≤λ≤6, MO=11 in the Lorenz problem.In Figure5, 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, ri∈{0.5,0.9,1,1.1,2} performed better than λ=6, ri∈{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 Table4), 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.
λMORAMCategory3132.03 MBLow cost4112.67 MBLow cost5103.73 MBLow cost41210.67 MBMiddle cost51118.63 MBMiddle cost41342.67 MBHigh cost51293.13 MBHigh costWe rejected the caseMO=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 ordersMO=12 and MO=13 in the Pleiades problem.The behaviour of the methods shown in Figure6 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.Figures7 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 withMO=13, FRABM4 with MO=11 and FRABM5 with MO=10 in a middle eccentric Two-Body problem.Figure 8
High cost, FRABM3 withMO=13, FRABM4 with MO=13, and FRABM5 with 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 forMO=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.Error versus CPU time. FRABM5MO=10: green line with stars, FRABM5 MO=12: red line with diamonds, VSVOABM MO=13: blue line with squares, DOPRI853: black line with triangles.Error versus evaluations. FRABM5MO=10: green line with stars, FRABM5 MO=12: red line with diamonds, VSVOABM 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 samey-coordinate (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 andCPU 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 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 thenumber 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 functionf(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 functionf(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 gains about 20% CPU time with respect to DOPRI853.The Pleiades problem is expensive to evaluate because it is formed by 28long 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 discretizedBrusselator with diffusion [11, pages 248-249], constituted by 882 short first-order differential equations.(6)
A discretized hangingRope problem. In this problem we let the dimension n in equation (10.8) of [11, pages 247-248] as a free parameter. Then this problem consists of 2nlong 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 withMO=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 beatsMO=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 whenf(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.(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 f(x,y) 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.



---

*Source: 765620-2010-04-11.xml*