1 Introduction
Numerical methods for evolution equations are designed based on accuracy and stability. The theory of both is highly developed for constant time step and linear problems. Less is known for variable time steps and nonlinear problems.These cases are subtle. For example, for increasing time steps, the BDF2 method loses Astability and suffers nonphysical energy growth in the approximate solution [1]. Even the trapezoidal method is unstable when used with variable time steps, see e.g. [2], [3, pp. 181182]. Dahlquist, Liniger and Nevanlinna in [2] proposed a one parameter family of variablestep, oneleg, twostep methods (DLN), which are secondorder accurate, and variablestep, nonlinearly, longtime stable. Its detailed specification (given in Section 2), is sufficiently Gordian to deter its use in complex applications, in which a method with DLN’s excellent properties should be valued. Our preliminary work on adaptive timestepping for flow problems [4, 5] shows that (DLN) has promise, motivating the work herein.
Refactorization generally means restructuring of an existing algorithm without changing its behaviour. The goal of refactorization is to reduce complexity by creating a simple and clean logical structure, improving implementation, code readability, source maintainability and extensibility. Herein we show how (DLN) can be refactorized to be easily implemented in an intricate, possibly legacy/blackbox code, without modifying the ‘assemble and solve’ portion. While our refactorization can work for other base methods, to fix ideas for , we consider a method based on the fully implicit Euler method
(BE) 
Figure 1 illustrates the implementation of the (DLN) method in Algorithm 1, by adding a prefilter step to the data ahead of the nonlinear solver (BE), and a posterfilter step after the solver (BE). This algorithmic idea is our first contribution. In Section 3 we prove a new expression for the local truncation error, which simplifies the timeadaptive implementation of (DLN), and also recall its variable step stability property.
2 The DLN method and its refactorization
We consider a numerical approximation of the initial value problem
(2.1) 
at times , with the time step . To present the method of [2], let denote the step size variability and be an arbitrary parameter. The (DLN) method is a 2step method with coefficients:
(2.2) 
Note that are step size independent, while are step size dependent. Define the average step size as follows:
(2.3) 
The variable step DLN method of [2] as a oneleg ^{1}^{1}1 The ‘oneleg’ term was coined by Dahlquist in 1975 [6] to name the multistep methods which involve only one value of in each step. In particular, the leapfrog and BDF methods are oneleg multistep methods. method is
(DLN) 
Remark 2.1.
To reduce the complexity of implementing (DLN), we consider its implementation through
pre and postprocesses of an implicit (backward) Euler method,
described in Algorithm 1, and illustrated in
Figure 1.
// Define the refactorization coefficients
(2.4) 
Since the coefficients satisfy
Proof.
2.1 Related Work
The (DLN) method is variablestep stable outgrowth of a method of Liniger [9], which is nonautonomous stable (i.e. for ). The pre and postprocess steps in the Algorithm 1 are akin to time filters, highly developed as numerical methods in atmospheric science [10, 11, 12, 13, 14]. Recently it was noticed in [15] that this technique for adding stability can also increase accuracy. The idea of prefilter simple method postfilter was develloped in a different direction for constant time steps in [16].
The refactorization of an algorithm to reduce its cognitive complexity has been used in [7] to rearrange a family of oneleg onestep methods into a backward Euler code followed by postprocessing, and further applied for partitioning multiphysics problems [17, 18, 8, 19]. In [20], the authors describe the implementation of the (DLN) formulas in a Nordsieck formulation [21, 22] essentially identical to that of the backward differentiation formulas, facilitating to adapt Nordsieck formulation codes like DIFSUB [23, 24] to the (DLN) formulas.
3 Convergence analysis of (Dln)
While stability and consistency were already addressed in [2], we present complementary details on both which are useful for developing an adaptive (DLN) method.
3.1 Consistency error
In [25, 26], the variable timestep (DLN) method was implemented in an adaptive manner, using a local and global error estimator. Similar to [2], the authors of [20, 25, 26, 6] use the classical definition of the local truncation error
where . The above definition follows the approach taken in the analysis of linear multistep methods (see e.g. [27, page 27]), and involves both the differentiation defect , and the interpolation defect . Dahlquist raised in [28] the question of the appropriateness of this viewpoint: “We accept this definition, but we do not accept as the adequate local truncation error!”
Using the refactorized form (DLN2BE) and Theorem 2.1, we now prove that the local truncation error of the oneleg (DLN) method can be evaluated only by the differentiation defect (3.1), similarly to the midpoint rule [7] and the RungeKutta methods. The new expression (3.1) simplifies greatly the error estimation.
Proposition 3.1.
The local truncation error of (DLN) is the differentiation error and
Proof.
The consistency order and the coefficient of the leading term in (3.1) follow by Taylor expansions. On one hand, since (DLN) can be refactorized as the onestep method (DLN2BE), we further write the (DLN) method as follows
(3.1) 
On the other hand, when we integrate (2.1) on and on , multiply the results by and , respectively, and add, we obtain
Finally, approximating both integrals in the right handside with the chord quadrature rule, with as the point of evaluation on both intervals, gives
Remark 3.1.
In particular, for and , from (3.1) we have that
3.2 Gstability
Let and denote the inner product and norm in Euclidean space . For any pair of solutions of (2.1), a necessary and sufficient condition [29, page 384] for to be a nonincreasing function of is the contractivity (onesided Lipschitz) condition on :
(contractivity) 
The system (2.1) for which satisfies the (contractivity) condition is called dissipative, see e.g. the Definition in [30, page 268]. We recall that a RungeKutta method is Bstable, if the (contractivity) condition implies for any numerical solutions, see e.g. [31, page 359], or Definition 12.2 in [32]. Similarly, a 2step linear multistep method is called stable [33, 29, 34, 32] if there exists a real positive definite matrix such that its oneleg version is contractive, namely , where . In the case of the (DLN) method, there exists such a positive definite matrix (independent of the step size)
As pointed out by Dahlquist in [35], both stability and stability imply stability, and stability implies stability for constant time steps.
Proposition 3.2.
The ‘energy’ identity (3.2), implicit in [2], follows from algebraic manipulations, see e.g. [4]. The stability of (DLN), i.e. , follows from (3.2) and the (contractivity) assumption. The only (DLN) methods which yield the invariance of the solution are the symplectic (onestep midpoint) and (twostep midpoint) rules: the numerical dissipation vanishes if and only if .
References
 [1] R. D. Grigorieff, Stability of multistepmethods on variable grids, Numer. Math. 42 (3) (1983) 359–377. doi:10.1007/BF01389580.
 [2] G. G. Dahlquist, W. Liniger, O. Nevanlinna, Stability of twostep methods for variable integration steps, SIAM J. Numer. Anal. 20 (5) (1983) 1071–1085. doi:10.1137/0720076.

[3]
H. J. Stetter, Analysis of discretization methods for ordinary differential equations, SpringerVerlag, New YorkHeidelberg, 1973, Springer Tracts in Natural Philosophy, Vol. 23.

[4]
W. Layton, W. Pei, Y. Qin, C. Trenchea, Analysis of the variable step method of Dahlquist, Liniger and Nevanlinna for fluid flow, Numer. Methods Partial Differential Equations , accepted (2021).
 [5] Y. Qin, Y. Hou, W. Pei, J. Li, A variable timestepping algorithm for the unsteady Stokes/Darcy model, J. Comput. Appl. Math. 394 (2021) 113521, 14. doi:10.1016/j.cam.2021.113521.
 [6] G. Dahlquist, Error analysis for a class of methods for stiff nonlinear initial value problems, in: Numerical analysis (Proc. 6th Biennial Dundee Conf., Univ. Dundee, Dundee, 1975), 1976, pp. 60–72. Lecture Notes in Math., Vol. 506.
 [7] J. Burkardt, C. Trenchea, Refactorization of the midpoint rule, Appl. Math. Lett. 107 (2020) 106438. doi:10.1016/j.aml.2020.106438.
 [8] M. Bukač, A. Seboldt, C. Trenchea, Refactorization of Cauchy’s Method: A SecondOrder Partitioned Method for FluidThick Structure Interaction Problems, J. Math. Fluid Mech. 23 (3) (2021) 64. doi:10.1007/s0002102100593z.
 [9] W. Liniger, The contractive secondorder multistep formulas with variable steps, SIAM J. Numer. Anal. 20 (6) (1983) 1231–1238. doi:10.1137/0720093.
 [10] R. Asselin, Frequency filter for time integrations, Mon. Wea. Rev. 100 (6) (1972) 487–490. doi:10.1175/15200493(1972)100<0487:FFFTI>2.3.CO;2.
 [11] A. J. Robert, The integration of a spectral model of the atmosphere by the implicit method, Proc. WMOIUGG Symp. on NWP, Tokyo, Japan Meteorological Agency (VII) (1969) 19–24.
 [12] P. D. Williams, Achieving seventhorder amplitude accuracy in leapfrog integrations, Mon. Wea. Rev. 141 (2013) 3037–3051. doi:10.1175/MWRD1200303.1.
 [13] Y. Li, C. Trenchea, A higherorder Robert–Asselin type time filter, J. Comput. Phys. 259 (2014) 23–32. doi:10.1016/j.jcp.2013.11.022.
 [14] A. Guzel, C. Trenchea, The Williams step increases the stability and accuracy of the hoRA time filter, Appl. Numer. Math. 131 (2018) 158–173. doi:10.1016/j.apnum.2018.05.003.
 [15] A. Guzel, W. Layton, Time filters increase accuracy of the fully implicit method, BIT 58 (2) (2018) 301–315. doi:10.1007/s105430180695z.
 [16] V. DeCaria, S. Gottlieb, Z. J. Grant, W. J. Layton, A general linear method approach to the design and optimization of efficient, accurate, and easily implemented timestepping methods in cfd (2020). arXiv:2010.06360.
 [17] M. Bukač, C. Trenchea, Boundary update via resolvent for fluid–structure interaction, J. Numer. Math. 29 (1) (2021) 1–22. doi:10.1515/jnma20190081.
 [18] C. Trenchea, Partitioned conservative, variable step, secondorder method for magnetohydrodynamics in Elsässer variables, ROMAI J. 15 (2) (2019) 117–137.
 [19] M. Bukač, C. Trenchea, Adaptive, secondorder, unconditionally stable partitioned method for fluidstructure interaction, Tech. rep., University of Pittsburgh (2021).
 [20] D. S. Watanabe, Q. M. Sheikh, Oneleg formulas for stiff ordinary differential equations, SIAM J. Sci. Statist. Comput. 5 (2) (1984) 489–496. doi:10.1137/0905036.
 [21] A. Nordsieck, On numerical integration of ordinary differential equations, Math. Comp. 16 (1962) 22–49. doi:10.2307/2003809.
 [22] R. D. Skeel, Equivalent forms of multistep formulas, Math. Comp. 33 (148) (1979) 1229–1250. doi:10.2307/2006457.
 [23] C. W. Gear, Numerical initial value problems in ordinary differential equations, PrenticeHall, Inc., Englewood Cliffs, N.J., 1971.
 [24] C. W. Gear, Algorithm 407: Difsub for solution of ordinary differential equations [d2], Communications of the ACM 14 (3) (1971) 185–190. doi:10.1145/362566.362573.
 [25] G. Y. Kulikov, S. K. Shindin, Oneleg integration of ordinary differential equations with global error control, Comput. Methods Appl. Math. 5 (1) (2005) 86–96. doi:10.2478/cmam20050004.
 [26] G. Y. Kulikov, S. K. Shindin, Oneleg variablecoefficient formulas for ordinary differential equations and localglobal step size control, Numer. Algorithms 43 (1) (2006) 99–121. doi:10.1007/s1107500690435.
 [27] J. D. Lambert, Computational methods in ordinary differential equations, John Wiley & Sons, LondonNew YorkSydney, 1973, Introductory Mathematics for Scientists and Engineers.
 [28] G. G. Dahlquist, On oneleg multistep methods, SIAM J. Numer. Anal. 20 (6) (1983) 1130–1138. doi:10.1137/0720082.
 [29] G. G. Dahlquist, stability is equivalent to stability, BIT 18 (4) (1978) 384–401. doi:10.1007/BF01932018.
 [30] J. D. Lambert, Numerical methods for ordinary differential systems, John Wiley & Sons, Ltd., Chichester, 1991, The initial value problem.
 [31] J. C. Butcher, A stability property of implicit RungeKutta methods, BIT Numerical Mathematics 15 (4) (1975) 358–361. doi:10.1007/BF01931672.
 [32] E. Hairer, G. Wanner, Solving ordinary differential equations. II, Vol. 14 of Springer Series in Computational Mathematics, SpringerVerlag, Berlin, 2010, stiff and differentialalgebraic problems, Second revised edition. doi:10.1007/9783642052217.
 [33] G. G. Dahlquist, On the relation of Gstability to other stability concepts for linear multistep methods, Dept. of Comp. Sci. Roy. Inst. of Technology Report TRITANA7621 (1976).
 [34] G. G. Dahlquist, Positive functions and some applications to stability questions for numerical methods, in: Recent advances in numerical analysis (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1978), Vol. 41 of Publ. Math. Res. Center Univ. Wisconsin, Academic Press, New YorkLondon, 1978, pp. 1–29.
 [35] G. G. Dahlquist, On stability and error analysis for stiff nonlinear problems, Part I, Dept. of Comp. Sci. Roy. Inst. of Technology Report TRITANA7508 (1975).