Local error estimation and step size control in adaptive linear multistep methods

In a k-step adaptive linear multistep methods the coefficients depend on the k − 1 most recent step size ratios. In a similar way, both the actual and the estimated local error will depend on these step ratios. The classical error model has been the asymptotic model, chp+ 1y(p+ 1)(t), based on the constant step size analysis, where all past step sizes simultaneously go to zero. This does not reflect actual computations with multistep methods, where the step size control selects the next step, based on error information from previously accepted steps and the recent step size history. In variable step size implementations the error model must therefore be dynamic and include past step ratios, even in the asymptotic regime. In this paper we derive dynamic asymptotic models of the local error and its estimator, and show how to use dynamically compensated step size controllers that keep the asymptotic local error near a prescribed tolerance tol. The new error models enable the use of controllers with enhanced stability, producing more regular step size sequences. Numerical examples illustrate the impact of dynamically compensated control, and that the proper choice of error estimator affects efficiency.


Introduction
We shall consider a standard initial value probleṁ y = f (t, y), y(0) = y 0 , and let y n denote the numerical approximation to y(t n ). Further, we let y n denote f (t n , y n ). Thus derivatives of functions are denoted by a dot, while samples of the vector field are denoted by prime. The methodology developed here does not depend on whether the problem is scalar or a system of equations. We will study linear multistep methods on nonuniform grids for this initial value problem, of the normalized form where h n = t n+1 − t n . Normalization will also be discussed in connection with error constants. The subscript n of the coefficients indicates that these are rational functions of the step ratios, ρ n−1 = h n /h n−1 . The method is implicit if b k,n = 0. We shall consider the grid-independent representation, [1], where every linear multistep method is defined in terms of a fixed set of parameters and interpolation conditions. Given any grid point distribution {t n } ∞ 0 , this representation generates the coefficients a k−j,n and b k−j,n .
Time-step adaptivity is of key importance to make initial value problem solvers efficient. In particular, in stiff problems, where time scales may vary by several orders of magnitude, it is necessary to vary the step size accordingly. If possible, an implementation should offer tolerance proportionality, i.e., the global error should be proportional to a preset accuracy requirement, TOL. The classical approach has been to assume that the local error can be represented asymptotically by r n = ϕ n h q n , (2) where ϕ n is the norm of the principal error function, h n is the current time step, and r n is the norm of the error estimate. The time step is then typically changed according to the elementary control law For a method of order p, we take q = p + 1 or q = p depending on whether the local error per step (EPS) or the local error per unit step (EPUS) is controlled. The scheme is usually accompanied by safety heuristics, and detailed descriptions can e.g. be found in monographs such as [3,6,12,13,17]. Starting in the late 80s, adaptivity based on control theory was developed, [9][10][11], leading to [18] and a framework for using digital filters and signal processing, [19,21]. Computational stability was further studied in [20]. Special controllers for explicit and implicit geometric integration were also developed in [14] and [23]. Combinations of these techniques have been applied to near-Hamiltonian problems with weak Rayleigh damping, [16].
This theory has focused on one-step methods, for which the asymptotic error model (2) applies whenever the tolerance is small enough. Taking logarithms, (2) becomes log r n = q · log h n + log ϕ n , suggesting that a step size change affects the error (approximately) by a constant "gain" factor of q, without any dependence on step size history. While this static error model is correct for one-step methods in the asymptotic regime, it is no longer correct outside asymptotic operation, [10]. For a k-step linear multistep method it is not even correct in the asymptotic regime. As the method coefficients depend on k − 1 past step ratios, errors are similarly affected, and the error's step size dependence becomes dynamic. Nevertheless, the elementary controller (3) has been applied to multistep methods, with varying success. Unfortunately, its stability may deteriorate in combination with e.g. the Adams-Bashforth methods. Thus, depending on the choice of error estimator, the dynamic interaction of method and controller may become resonant and oscillatory. The problems can be overcome by a careful combination of method, estimator and controller.
In order to improve adaptive multistep methods, we shall investigate new dynamic error models of multiplicative form, Here the first factor is identical to the static model (2). The second factor accounts for step size history in terms of step ratios ρ n−j = h n−j +1 /h n−j . The number of past step ratios is s = k − 1 for a k-step method if r n represents the actual error, while it is s = k for error estimators that use additional data from one previous step. The parameters δ j are characteristic of each multistep method and its error estimator. In multistep methods, there is a concern over zero stability when the step size varies, focusing on how large step ratios can be allowed without causing instability, see e.g. [5,7,8,22]. While this issue is of significance, it does not reflect that feedback control interacts with the method and manages overall stability. Thus, in a strongly zero stable method a stable computational process can be maintained by generating a smooth step size sequence (i.e., having step ratios ρ n ≈ 1). This can be achieved by a well designed controller, working with incremental step size changes, avoiding traditional step size halving/doubling.
The grid-independent approach to multistep methods has also been implemented in a proof-of-concept software platform, [2], where any multistep method can be evaluated under ceteris paribus conditions. While the theory makes it possible to find method-specific dynamic error estimators in terms of current and past step sizes, the implementation has so far only offered controllers for the static model (2). Our main objectives in this paper are therefore 1) to analyze dynamic error estimators, based on models of the form (5); and 2) to devise suitable controllers that manage the local error as well as the stability of the interaction of method and controller. Our approach is based on using discrete control theory as outlined in [19].
The main result is that we construct a dynamic compensator that extracts the static part of the asymptotic error model from the estimator, allowing us to employ standard digital filters and other controllers to achieve results comparable to those of onestep methods. We further show that without this dynamic compensation, conventional techniques are not completely reliable. Thus, for example, the elementary controller (3) is unsuitable for the control of adaptive linear multistep methods, as the stability margin of the process deteriorates with increasing step size history dependence.

Dynamic error models and control objectives
In the classical constant step size theory, the local error of a multistep method has an asymptotic expansion with leading term as h → 0, where c * is referred to as the method's (normalized) error constant. But this assumes that all steps shrink simultaneously. In adaptive methods, however, having reached the point t n+1 = t n + h n , we consider changing the next step size, h n+1 , but none of the previously accepted steps h n , . . . , h n−k+1 . As a consequence, the asymptotic error will depend on previous step ratios. Local error control is justified by the possibility of linking a local error bound to the global error via tolerance proportionality. Global error accumulation can be modeled by the variational equation,u = J (t)u + w, where J (t) = ∂f/∂y along the exact solution y(t). Here u(t) and w(t) represent global and local errors, respectively. Let μ[J (t)] denote the logarithmic norm of J (t), and assume that μ[J (t)] ≤ M for all t ≥ 0. Note that M may be negative. Further, let w ∞ = max t≥0 w . Then a global error bound can be obtained by integrating the differential inequality d dt u ≤ M · u + w ∞ over a single step to obtain In nonstiff computation, we have |hM| 1. Apart from the error propagation, the global error is then incremented by approximately h w ∞ in a step of length h. By using the control objective w ∞ ≤ TOL we keep the global error growth rate in check and achieve tolerance proportionality. This corresponds to local error per unit step (EPUS) control.
If on the other hand −hM 1, modeling stiff computation and strong damping, the exponential terms in (7) can be neglected, and we have Thus, with little or no error propagation, the global error effectively equals (a fraction of) the local error. To achieve tolerance proportionality, it is then sufficient to use the control objective h w ∞ ≤ TOL, corresponding to local error per step (EPS) control. The choice of EPUS or EPS determines the power q in the error model, but by suitably expressing the controller's parameters in terms of q, one can obtain similar overall dynamics for both control objectives as well as for methods of different orders.
Inserting the exact solution y(t) into (1), we obtain a residual, wherel n = c(ρ n ) h p+1 n y (p+1) (t n ). Compared to (6), the error constant c * has been replaced by a function c(ρ n ), which depends on the k − 1 most recent step ratios, collected in the vectorρ n = (ρ n−1 , . . . , ρ n−k+1 ). Letting u n = y n − y(t n ) denote the global error, we subtract (8) from (1) and linearize to obtain n . This is the multistep discretization of the variational equationu = J (t) · u + w. The function w is constant within each step, and where π(ρ n ) = c(ρ n )/b(ρ n ) introduces dynamics into the error model. In constant step size theory,ρ = (1, 1, . . . , 1) = 1. The method is often normalized so that b(1) = 1, implying that c(1) = c * , cf. [12, p. 373]. It then follows that π(1) = c * . Since we have used the normalization a k,n = 1 instead, the error constant is affected. Letπ(ρ) = π(ρ)/π(1). Then (9) becomes making the error constantĉ * = π(1) conform to the method normalization and error propagation model above. For the Adams methods, it always holds that b(ρ n ) ≡ 1, implying thatĉ * = c * even when our normalization is used. For general multistep methods, we have thus arrived at the error model with q = p for EPUS control and q = p + 1 for EPS. Letδ = {δ j } s 1 be given by the gradientδ Assume thatρ n = 1 +v n with v n ∞ 1. It follows that logρ n ≈v n , and Hence the approximationπ(ρ n ) ≈ s 1 ρ δj n−j is accurate to O( v 2 ∞ ), and establishes the multiplicative dynamic error model (5). This approximation enables the design of a controller that compensates the specific linear dynamics of the error. It is also a necessary simplification when the step number k is large, as the rational function π(ρ n ) consists of multinomials with a very large number of terms. For this reason, the determination of the powers δ j typically requires symbolic computation.
The error model is best illustrated by studying a simple example, and we choose the two-step Adams-Bashforth (AB2) of order p = 2. In the grid-independent formulation [1] the adaptive AB2 method reads where h n = t n+1 − t n . The local error r n = y(t n+1 ) − y n+1 at t n+1 is found by inserting exact data, i.e., y n = y(t n ), y n =ẏ(t n ) and y n−1 =ẏ(t n−1 ), and expanding in a Taylor series to obtain the asymptotic error expressioñ where higher order terms (here fourth order in the step size) have been neglected. For constant step size (ρ n−1 = 1) we recognize the classical asymptotic principal error function of the AB2 with its error constant, as ϕ n = 5 ... y /12. For variable steps, we now havê Consequently, the asymptotic local error model in multiplicative form (5) is While "dimensionless," the step ratio does not affect the order q, but modifies individual powers of the step sizes. This has little influence on error magnitude (as ρ n−1 ≈ 1) but alters step size dynamics and the subsequent control design. In one-step methods, the elementary controller is derived from (2) by assuming that ϕ is slowly varying along the solution to the differential equation. Thus, if the error r n deviates from TOL, the next step size can be chosen to make r n+1 ≈ TOL by solving for h n+1 , assuming that ϕ n+1 = ϕ n . Dividing (17) by (2) we obtain the control law (3), known as a deadbeat controller, [19]. This has finite impulse response (FIR), since an impulse in the sequence {ϕ j } ∞ 0 will be compensated by the controller in a finite number of steps, in this case a single step. In fact, it is easily shown that the control law yields ϕ n h q n+1 ≈ TOL, i.e., the proposed step size is matched to the previous value of the principal error function.
The same intuitive approach could be considered for multistep methods. Assuming that ϕ n+1 = ϕ n in the AB2 method and dividing the two equations TOL = ϕ n+1 h Apart from the unexpected power of TOL/r n in (18), the controller has an extra factor due to the error's dynamic character. Perhaps more surprisingly, (18) is not deadbeat, and has a less desirable behavior. There are also alternative ways to estimate the error, and a few will be described in detail in Section 5. Above, the error (16) was obtained by Taylor series expansion. It can be computed in practice, using a comparison method (at least) one order higher than AB2, say the AB3, operating in the asymptotic regime. Another common estimation technique is to compare the results from two methods of the same order. This is similar to the classical predictor-corrector technique. Thus, if the AB2 method is represented as a method constructing a polynomial of degree 2, the error can be estimated by comparing this polynomial at t n+1 to the result obtained by extrapolating the corresponding polynomial, constructed on the previous step. While crude, this difference can still be used to compute an asymptotically correct error. The dynamic asymptotic error model is then more complex, An attempt to derive control laws by the simple approach used above fails, and results in unstable controllers, both in EPUS and EPS modes. The elementary control law (3) fares a little better, and is stable in EPS mode in combination with (19), but unstable in EPUS mode. Thus, if methods, estimators and controllers are combined arbitrarily, performance is often underwhelming, see Fig. 1, where the elementary controller is compared to one of the dedicated dynamically compensated controllers we will develop in this paper. This motivates a thorough investigation of error estimation and step size control for adaptive multistep methods, focusing on the dynamics of the error model.

Tools from linear control theory
In the implementation of a discrete time system, it is important to synchronize variable indexation. When using the z transform, all variables, sequences and transforms are expressed in relation to the reference index n, corresponding to the point t n , and we use the following time domain numbering convention: 1. Given the numerical solution (t n , y n ), the method takes a step size h n to compute y n+1 at t n+1 = t n + h n 2. The method produces an error estimate r n at t n+1 , modeled in the asymptotic regime by an error model such as (5) 3. The controller computes the next step size ratio ρ n from r n 4. The controller then updates the step size according to h n+1 = ρ n h n 5. The computational process repeats from 1.
We shall analyze this process using linear control theory, which is based on the analysis of linear difference equations. The multiplicative forms of the error model and controller are necessary, since they can be converted to difference equations by taking logarithms, after which the z transform is used to analyze stability and frequency response. For a simple demonstration of this procedure, let us consider (5) for s = 2, i.e., Now, taking logarithms in (20), we obtain the difference equation where E is the forward shift operator, informally defined by Eu n = u n+1 for any discrete sequence. The rational function G is defined by Note that, due to the structure of (20), it follows that G(1) = q.
We let h = {h j } ∞ 0 denote the entire sequence of step sizes, and use the simplified notation log h = {log h j } ∞ 0 . A standard tool in discrete control theory is the z transform, defined for the sequence log h by Since the risk of misunderstanding is minimal, with a slight abuse of notation we follow [19] and let log h denote the z transform of the sequence log h. The z transform of the linear difference equation (21) is then The function G(z), referred to as the process model, is the z transform of the operator associated with the step size-error relation (21), while the additive term log ϕ is an external disturbance corresponding to the principal error function, which is to be "rejected" by the controller. This means that the controller adjusts the step size to counteract the influence of log ϕ.
In order to control the error magnitude by adjusting the step size, the error sequence log r is compared to the tolerance log TOL, and the sequence log TOL r j ∞ 0 = log TOL − log r is referred to as the control error. The controller aims to eliminate the control error, continually adjusting the step size according to the control law where C(z) is the z transform of the linear difference operator associated with the controller. The overall interaction between process and controller is usually represented by a block diagram, see Fig. 2.
Next we analyze the process-controller interaction. Combining the control law with the process model, we have the linear system , the computational process takes log h as input, producing an error estimate log r = G(z) log h + log ϕ. The principal error function log ϕ enters as an additive disturbance, to be compensated by the controller. The error estimate log r is fed back and compared to log TOL. The controller, represented by its transfer function C(z), constructs the next step size through log h = C(z)·(log TOL−log r). All sequences, as well as the controller and the process, are represented by their z transforms Solving for log h and log r in terms of the data log ϕ and log TOL, we obtain where we identify the four closed loop transfer functions, from log ϕ and log TOL respectively, to log h and log r. As the constant TOL is of little interest, we take TOL = 1 (log TOL := 0), meaning that the error is measured in units of TOL. Note that changing the constant TOL has no effect on the error, since R TOL (z) log TOL = 0 for any constant TOL. Likewise, it has little effect on the step size, other than adding a constant to the step size sequence log h.
The step size transfer function and error transfer function are, respectively, The asymptotic process model G(z) is found by analyzing the error estimator of a given linear multistep method, while the controller C(z) can be chosen to meet various design criteria for stability and smooth step size sequences. The closed loop transfer functions have the same poles, which correspond to the zeros of the characteristic equation and govern stability. They must remain well located inside the unit circle. It is sometimes possible to construct C(z) so as to move all poles to z = 0, corresponding to deadbeat control.
To model a bounded, quasi-periodic input one frequency at a time, we take log ϕ = {e iωn } for ω ∈ [0, π], with ω = π corresponding to the highest discrete frequency, log ϕ n = (−1) n . For such periodic inputs, it holds that where the complex value H ϕ (e iω ) reveals the phase and amplitude of the step size sequence compared to the input log ϕ = {e iωn }. Disregarding phase, the scaled step size frequency response and the error frequency response are given by, respectively, where the extra factor of q in |qH ϕ (e iω )| normalizes the step size attenuation so that A h (0) = 1. This factor is only a matter of convenience, reflecting that for the asymptotic process it always holds that G(1) = q.
Step size and error attenuations are dimensionless ratios, measured in the ISO unit of decibel (dB), defined by 20 · log 10 A h (ω). Thus an attenuation of +20 dB corresponds to an amplification by one order of magnitude. Likewise, an attenuation of −6 dB corresponds to a damping by a factor of 2.
The zeros (if any) of the step size transfer function |qH ϕ (z)| block signal transmission. Thus, selecting C(z) such that A h (π ) = 0, one can design low-pass digital filters that suppress noise and produce smooth step size sequences. This is covered in detail for the static process G(z) = q in [19,21].
In Fig. 3 we plot the step size frequency responses corresponding to the AB2 method operating in EPUS and EPS modes, using two different error estimators and the same two controllers as in Fig. 1, with a one-to-one correspondence between the plots. As before, the designation AB2(3) refers to a third order error estimator, With the third order estimator (top panels) the resonance is less pronounced. By contrast, a dynamically compensated controller (blue) eliminates the resonance and offers additional high frequency suppression, generating smooth step size sequences. Dashed black reference curve represents deadbeat control. Attenuation is measured in decibels, dB while AB2(2) refers to a second order estimator. For the latter, the frequency response reveals that the elementary controller has a damaging resonance near ω ≈ 1. This contaminates the transient performance of the elementary controller, even causing the unstable oscillatory behavior in EPUS mode observed in Fig. 1. Although EPUS mode is always worse than EPS, the resonance is much less prominent for the third order estimator. There are two conclusions to be drawn. First, as performance differs strongly, it is important to combine good error estimators with controllers so as to avoid resonance and instability. While the choice of error estimator is important, a dynamically compensated controller is always better than the elementary controller, which is unsuitable for multistep methods. Figure 3 demonstrates that a dynamically compensated controller cannot only manage the benign third order estimator, but also the second order estimator, with a nearly uniform behavior, irrespective of whether the objective is EPS or EPUS.

Control analysis of dynamic error models
Beginning with two-step methods, we assume that a computable error estimator can be modeled nearρ = 1 by The validity of this structure will be demonstrated in Section 5. The objective is to construct a controller that determines ρ n and updates the step size according to h n+1 = ρ n h n . The process has three parameters, (q, δ 1 , δ 2 ), and by combining this process with a controller of the following structure, we have three free parameters, (β, α 1 , α 2 ), to control the placement of the three poles of the closed loop. The control structure is general and follows the pattern of digital filters, [19], and we note that the elementary controller has β = 1 and α 1 = α 2 = 0. The general result is the following.

If this method is combined with the step size control law
taking α j = δ j ·β/q for all j , the closed loop has a single nontrivial pole at z = 1−β, with all other poles at z = 0. If in addition β = 1, the controller is deadbeat (finite impulse response).
Proof Taking logarithms in (30), we obtain log r n = (q + δ 1 ) log h n + (δ 2 − δ 1 ) log h n−1 + · · · − δ s log h n−s + log ϕ n implying that the transfer function is In a similar fashion, for the controller we obtain Factoring out the difference operator z − 1, this simplifies to Thus we find the control transfer function It follows that the closed loop error transfer function is As for the poles, starting at j = s, we recursively find that the coefficient of z s−j can be made zero by taking putting successively more poles at z = 0, until j = 1. If all α j have been chosen accordingly, z = 0 is a pole of multiplicity s and the denominator of R ϕ (z) is reduced to with the last pole located at z = 1 − β.

Remark 1
The factor z−1 in the numerator of R ϕ (z) represents a difference operator that removes any persistent disturbance in log ϕ. This first order adaptivity, [19], is a consequence of the controller's integral action (29), log h n+1 = log h n + log ρ n , which generates the step size by summing ("integrating") the step ratios. (3) in the case s = 2, for which α 1 = α 2 = 0 and β = 1. The poles are then determined by

Remark 2 Consider the elementary controller
In the static error model (like in all one-step methods) δ 1 = δ 2 = 0, implying that all roots are located at z = 0. This makes the elementary controller deadbeat for one-step methods. For multistep methods, however, the elementary controller is not deadbeat, as δ j = 0. In two-step methods, if the error estimator is one order higher than the method, there is a single δ coefficient. When δ 1 < 0 there are two complex conjugate poles, leading to the resonance peaks of the type already observed in Fig. 3. Although the elementary controller is stable in tandem with standard multistep methods such as AB2 and AM2, the matching is poor, inducing slowly damped oscillations in the step size sequence. As the modulus of the roots is |z| = √ |δ 1 /q|, the damping is weaker in EPUS mode than in EPS.

Remark 3 The theorem implies that there is a basic class of single-parameter controllers of the form
that can be applied in tandem with any linear multistep method, characterized by its asymptotic error parameters q and {δ j } s 1 . The closed loop is stable whenever β ∈ (0, 2), but in order to have a non-negative root, only β ∈ (0, 1] is of interest. It is necessary to include a step ratio compensator to achieve deadbeat control. While a deadbeat controller is not always preferable, its existence is of significance as pole locations depend continuously on the parameters. Thus there is a neighborhood of stable controllers near α j = δ j · β/q and β = 1, and the full set of parameters {α j } s 1 and β can be chosen to place the s + 1 poles at suitable locations.

Remark 4 Due to (36) the control problem is simplified by introducing the dynamically compensated error,r
which recovers the static model, asr n ≈ ϕ n h q n . In terms of the compensated error, the control system (36) above becomes This is merely a standard reduced gain integral controller, identical to the exponential forgetting used for one-step methods, [19]. By taking the gain parameter β = 2/3, which puts the nontrivial pole at z = 1/3, one obtains a good performance for all multistep methods operating in the asymptotic regime. This value is used in Figs. 1 and 3. Other types of controllers for one-step methods, including digital filters, [19], are also applicable to the dynamically compensated error model. As long as the step ratios are near 1, the dynamic compensator will not significantly affect the magnitude of the error, implying thatr n ≈ r n . However, the compensator will affect the dynamics and eliminate damaging resonances in the adaptive scheme.

Local asymptotic error estimators
Using the grid-independent representation of multistep methods in [1] the step from t n to t n+1 = t n + h n is based on constructing a polynomial P n+1 (t) satisfying a number of structural conditions and slack conditions, such that the solution to the differential equation y(t n+1 ) can be approximated by The error is defined by inserting the exact solution into the discretization, and satisfies y n+1 − y(t n+1 ) ≈ c y (ρ n ) · y (p y +1) (t n ) · h p y +1 n . However, this quantity is not directly computable. To obtain a reference value, we need a comparison method that produces another approximation, Although this quantity is not directly computable either, it follows that is computable. Here we need to distinguish two cases. Case 1. In the first case, the comparison method is one order higher, i.e., p z = p y +1. For example, if the method is AB2, one can use AB3 to compute the reference value. We shall refer to this method/error estimation combination as AB2 (3), where the number within paranthesis refers to the order p z of the error estimator. Now, if both methods operate in the asymptotic regime, the higher order result can be considered "exact," and the error estimate is effectivelỹ This implies that the error constant and the dynamicsπ l (ρ) are those from the lower order method. In effect, it emulates the situation where the actual error of the method has been computed.
Case 2. In the second case, the order of the comparison method is the same as that of the method, i.e., p z = p y . In the AB2 example, this method/error combination will be denoted AB2(2) to distinguish it from AB2(3) above. This computation can be arranged in different ways. If the solution at time t n+1 is represented by the value of a polynomial y n+1 = P n+1 (t n+1 ), one may e.g. use the polynomial from the previous step, P n (t), and compute the extrapolated reference value z n+1 = P n (t n+1 ). We then have obtaining a raw error estimate, whereπ e (ρ) accounts for the dynamics of the estimator, with which the controller interacts. This function is derived in the same way as before, by obtaining a variable step size formula for the method, as well as a formula for how the previous polynomial predicts z n+1 through extrapolation. When the difference between these formulas are expanded in a Taylor series, we obtain e n andπ e (ρ), which is approximated in the usual multiplicative way.
The raw estimatesl n and e n need further scaling to conform to the function w in (10). In Case 2 in EPUS mode, e n must be scaled back to the desired estimate r n of the method's actual error. According to (9), we want to control the magnitude of the quantity where the estimator's error constant is given by (1) .
In Case 1, the constant C e is replaced by C l = 1/b y (1) andπ e (ρ n ) is replaced bŷ π l (ρ n ). We can then describe the computational process as follows.
The complete system is illustrated in Fig. 4, where the classical use of the elementary controller corresponds to taking K(z) = 1 (no dynamic compensation) and F (z) = 1 (no digital filter). Note thatπ e (ρ n ) andπ l (ρ n ) are replaced by their respective multiplicative approximants to construct the dynamically compensated error estimate. Thenr n ≈ ϕ n h p y n in the asymptotic regime, enabling the use of standard controllers. Without the compensator, controller performance drops, sometimes exciting resonances or instability. Table 1 gives filter coefficients for the step size controllers available in the proof-of-concept software [2]. Fig. 4 A complete feedback control system with digital filter F (z) and dynamic compensator K(z), which transforms the error estimate log r using the step ratios log ρ. The controller is decomposed into three parts. The first is a simple gain of 1/q, converting the compensated control error into the scaled control error log c. The digital filter F (z) generates the step ratio log ρ. Finally, a simple integrator, 1/(z − 1), changes the step size, which enters the computational process G(z). Here, we have taken log TOL = 0 Table 1 Filter coefficients for various controllers, [2,21] For selected methods of step number k ≤ 4, we give the parameters C l , C e and {δ j } s 1 in Tables 2 and 3. The computation of these parameters is complex and was carried out in MATHEMATICA and MAPLE, while numerical emulations and experiments were done in MATLAB. Parameters for 5-step methods can also be provided on request. The method parameters are associated with the particular type of error estimation procedures described above, and may not work if other estimation procedures are used. Table 2 Error coefficients and parameters for selected 2-step, 3-step and 4-step methods In each method designation, the number within parenthesis refers to the order of the error estimator, which in this table is one order higher than the method order, i.e., p z = p y + 1 with s = k − 1. The parameters are independent of the actual choice of estimator as long as p z > p y . Method designations refer to those used in [1,2] Table 3 Error coefficients and parameters for selected 2-step, 3-step and 4-step methods In each method designation, the number within parenthesis refers to the order of the error estimator, which in this table is the same as the method order, i.e., p z = p y with s = k. The parameters are given for the specific estimation procedure obtained by using polynomial extrapolation. Method designations refer to those used in [1,2]

Experimental results
Four well-known standard problems, linear as well as nonlinear, were chosen to benchmark controller performance when different two-step methods were combined with two different types of error estimators. Three-step methods were also tried, and the controllers were either the elementary controller (the reference) and three dynamically compensated controllers: exponential forgetting with gain β = 2/3, and the digital filters H 211PI and H 211b. All combinations were run both in EPUS and EPS modes, and for three different tolerances, for the first three problems, to check step size sequences and step ratios, as well as tolerance proportionality. In all cases the step size was selected to control the absolute error, and care was taken to ensure that the methods operated in the asymptotic regime. Initial step sizes were equally spaced, and in implicit methods the nonlinear systems were solved by Newton iteration to full accuracy.
This provided a vast evaluation material, and only a few tests can be included to illustrate the new adaptivity in practical operation. The selected examples were chosen to illustrate various performance aspects. For example, a poor combination of method and error estimator may go unstable in tandem with the elementary controller, but with a dynamic compensator the problems can be rectified. But this relies crucially on knowing the δ j coefficients listed above for selected methods.
The first test problem is the two-compartment dilution process for t ∈ [0, 20], with exact solution The initial values were chosen as the exact solution values at the first steps.
The second test problem is the Lotka-Volterra equation, The initial values were taken as [1,1], using a second order Taylor expansion around this point. The third problem is the van der Pol equation, 20]. The initial values were [2,0], with a second order Taylor expansion providing additional initial values. We chose μ = 2 for nonstiff computation. Tolerance proportionality was checked by taking TOL = 10 −m for a suitable integer m in each test, as well as TOL = 10 −m±q . Thus, when the local error per (unit) step is controlled to keep ϕh q ∼ TOL, the step size should scale accordingly, like h ∼ TOL 1/q , along the entire solution. All computations confirm that the control system as well as the error estimators work in accordance with theory, producing three different step size sequences one order of magnitude apart, see Figs. 5 and 6.
We finally tried stiff computation, where the method operates in part outside standard asymptotic assumptions. Thus we solved the van der Pol equation for μ = 10 3 and t ∈ [0, 1000] with the BDF2 method. Here, too, the controllers work as expected, with the best performance obtained with the BDF2/AM2 combination and the H 211PI controller, see Figs. 7 and 8.
To also run a "large-scale" problem in moderately stiff computation, we ran a fourth problem, the 1D Brusselator reaction-diffusion equation,  The AB2(2) method, using polynomial extrapolation as error estimation, is combined with the elementary controller (upper panels) and the compensated exponential forgetting controller with gain 2/3 (lower panels), running in EPUS mode with TOL = 10 −5 , 10 −7 , 10 −9 . Linear problem (left), Lotka-Volterra (center), and van der Pol problem (right) all exhibit instabilities in the elementary controller in this setting, causing ringing, in particular in the Lotka-Volterra problem. The instability is overcome by the new controller, but is also less pronounced for sharp tolerances Fig. 6 The AM2(4) method, using the fourth order AB4 as error estimator, is combined with the elementary controller (upper panels) and the compensated exponential forgetting controller with gain 2/3 (lower panels), running in EPS mode with TOL = 10 −4 , 10 −8 , 10 −12 . Linear problem (left), Lotka-Volterra (center), and van der Pol problem (right) show good behavior with the elementary controller, but transient performance improves for the compensated controller, in particular in the linear and Lotka-Volterra problems Step sizes vary over seven orders of magnitude and scale correctly This problem was solved using the method of lines with a standard equidistant FDM space discretization of the diffusion term, where we chose Δx = 1/(N + 1) with N = 60. This setup deviates somewhat from the standard setup in the Bari test set, [15], the main differences being that we have chosen an initial condition for u that is richer in Fourier components, and we use N = 60 for somewhat higher resolution and increased stiffness, using a total of 120 equations. The problem was solved using the BDF2 method in EPS mode, with the error estimated by AB3, BDF2 extrapolation and AM2 methods, respectively. The absolute error was measured in the discrete L 2 norm, and the H 211P I controller was used at two tolerances, TOL = 10 −4 and TOL = 10 −7 , to verify order and tolerance proportionality. Equidistant starting data were generated by the explicit Cash-Karp 5th order Runge-Kutta method, [4]. In this moderately stiff setup, the compensated error control still works to exacting standards, with data for six runs presented in Table 4, with the corresponding step size sequences in Figs. 9 and 10. Step sizes are plotted on a linear scale for TOL = 10 −6 to reveal transient behavior and step size magnitude. Upper three curves use AM2 estimator together with the uncompensated elementary controller (blue); with compensated exponential forgetting (black); and with H 211PI (red). Lower three curves use polynomial error estimator, and the same color coding. The AM2 estimator is twice as efficient, while the compensated H 211PI offers the best dynamic response and smoothness

Conclusions
We have demonstrated that the asymptotic local error model for multistep methods has the form r n = ϕ n h q n ·π(ρ n ), whereρ n = {ρ n−j } s j =1 and the step ratios are defined by ρ n−j = h n−j +1 /h n−j , with t n+1 = t n + h n . The functionπ(ρ) accounts for step size history, defining the dynamically compensated error where the vectorδ = {δ j } s j =1 is characteristic of each linear multistep method and its error estimator. The multiplicative form of the compensator not only facilitates a linear control analysis, but is also far simpler to implement.  It is also shown that the conventional elementary controller may suffer damaging resonance phenomena, which are eliminated by instead controlling the compensated error. Another common issue is that many existing multistep implementations only allow the step size to be changed at fixed ratios, exciting transient behavior in the step size control. To overcome anomalies, one needs to continually adjust the step size and use a dynamic compensator.
Thus the implementation of multistep methods needs to be reconsidered. This includes specifying a method together with a dedicated error estimator and its δ j coefficients, following a practice similar to that of Runge-Kutta methods, where an "embedded form" is used. Due to the complexity of variable step size multistep methods, we believe that this is a substantial task with a potential to improve current software standards.
Moreover, with the dynamic compensator and a well selected controller, it is possible to achieve smooth step size sequences and tolerance proportionality in EPUS mode. The implementation becomes tolerance convergent for Lipschitz problems. The sharper the tolerance, the smoother is the step size sequence, cf. [22]. The adaptive numerical solution converges to the exact solution as TOL → 0. Asρ n → 1, the adaptive method gradually starts behaving like a constant step size method, reducing the need for the dynamic compensator.
Still, the choice of method order is not straightforward. At a given TOL, a higher order method can complete the integration in fewer steps. But fewer steps also mean a less regular step size sequence, possibly with a loss of tolerance proportionality. High Step sizes with H211PI controller BDF2AM2 BDF2AB3 BDF2extrapolation

Fig. 10
Brusselator problem is solved using BDF2 in EPS mode with dynamically compensated H 211PI control.
Step sizes are plotted on a logarithmic scale for t ∈ [0, 20]. Upper three curves use TOL = 10 −4 with AM2 estimator (blue); AB3 estimator (red) and BDF2 extrapolation (black). Lower three curves use the same color coding, but with TOL = 10 −7 . Notice that the first period differs from subsequent periods; compare Fig. 9. The AM2 estimator is again twice as efficient as the alternative error estimators as step sizes are twice as large. Tolerance proportionality is evident for all combinations order methods should therefore primarily be used when the accuracy requirement is high. For a smooth step size sequence and a predictable adaptive behavior, one typically needs several hundred steps to complete the integration, meaning that the method order should primarily be chosen with respect to TOL, and less so with an aim of minimizing total number of steps. It will still offer a good work/precision trade-off. Although the new techniques are also demonstrated to work as expected in stiff problems, further studies may be necessary for stiff computation, where the method at least in part may operate outside asymptotic conditions.