Refactorization of a variable step, unconditionally stable method of Dahlquist, Liniger and Nevanlinna

by   William Layton, et al.
University of Pittsburgh

The one-leg, two-step time-stepping scheme proposed by Dahlquist, Liniger and Nevanlinna has clear advantages in complex, stiff numerical simulations: unconditional G-stability for variable time-steps and second-order accuracy. Yet it has been underutilized due, partially, to its complexity of direct implementation. We prove herein that this method is equivalent to the backward Euler method with pre- and post arithmetic steps added. This refactorization eases implementation in complex, possibly legacy codes. The realization we develop reduces complexity, including cognitive complexity and increases accuracy over both first order methods and constant time steps second order methods.


page 1

page 2

page 3

page 4


Numerical analysis of an efficient second order time filtered backward Euler method for MHD equations

The present work is devoted to introduce the backward Euler based modula...

Walking into the complex plane to 'order' better time integrators

Most numerical methods for time integration use real time steps. Complex...

Analysis of the variable step method of Dahlquist, Liniger and Nevanlinna for fluid flow

The two-step time discretization proposed by Dahlquist, Liniger and Neva...

On the Analysis of the Second Order Time Filtered Backward Euler Method for the EMAC formulation of Navier-Stokes Equations

This paper considers the backward Euler based linear time filtering meth...

Proper Selection of Obreshkov-Like Numerical Integrators Used as Numerical Differentiators

Criteria for Obreshkov-like numerical integrators to be used as numerica...

Doubly-Adaptive Artificial Compression Methods for Incompressible Flow

This report presents adaptive artificial compression methods in which th...

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 A-stability and suffers non-physical 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. 181-182]. Dahlquist, Liniger and Nevanlinna in [2] proposed a one parameter -family of variable-step, one-leg, two-step methods (DLN), which are second-order accurate, and variable-step, nonlinearly, long-time 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 time-stepping 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/black-box 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


Figure 1 illustrates the implementation of the (DLN) method in Algorithm 1, by adding a pre-filter step to the data ahead of the nonlinear solver (BE), and a poster-filter 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 time-adaptive 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


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 2-step method with coefficients:


Note that are step size independent, while are step size dependent. Define the average step size as follows:


The variable step DLN method of [2] as a one-leg 111 The ‘one-leg’ 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 one-leg multistep methods. method is

Remark 2.1.

The (DLN) methods are indexed by the free parameter . When , the (DLN) method becomes the (implicit) midpoint rule [7, 8]

(one-step midpoint)

while for , the (DLN) method is the (implicit) midpoint rule with double time step

(two-step midpoint)

To reduce the complexity of implementing (DLN), we consider its implementation through pre- and post-processes of an implicit (backward) Euler method, described in Algorithm 1, and illustrated in Figure 1.

Input: , and ;
// Pre-process : interpolation
// Evaluate quantities in (2.2) and (2.3)

// Define the refactorization coefficients
// Evaluate the time-step for BE
// Set the BE time interval: , and
;   ;
// backward Euler
Solve for :    // Post-process : extrapolation
  // the DLN solution

,  If desired: Estimate Error and adapt

Algorithm 1 Refactorization of the (DLN) method

Since the coefficients satisfy







Backward Euler


Figure 1: Refactorization of the (DLN) method as a pre- and post-processed (BE) method
Theorem 2.1.

Algorithm (1) is equivalent to the (DLN) method.


First using the notations (2.4), the post-processing step writes

Using also the pre-processing relations, the backward Euler step in Algorithm 1


translates to

Finally, by (2.4), this shows that the backward-Euler based Algorithm 1 yields the solution of the (DLN) method. ∎

2.1 Related Work

The (DLN) method is variable-step -stable outgrowth of a method of Liniger [9], which is non-autonomous -stable (i.e. for ). The pre- and post-process 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 one-leg one-step methods into a backward Euler code followed by post-processing, and further applied for partitioning multi-physics 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 time-step (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 one-leg (DLN) method can be evaluated only by the differentiation defect (3.1), similarly to the midpoint rule [7] and the Runge-Kutta 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


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 one-step method (DLN2BE), we further write the (DLN) method as follows


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 hand-side with the chord quadrature rule, with as the point of evaluation on both intervals, gives

which by (2.2)-(2.3) yields the (3.1) method. ∎

Remark 3.1.

In particular, for and , from (3.1) we have that

3.2 G-stability

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 non-increasing function of is the contractivity (one-sided Lipschitz) condition on :


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 Runge-Kutta method is B-stable, if the (contractivity) condition implies for any numerical solutions, see e.g. [31, page 359], or Definition 12.2 in [32]. Similarly, a 2-step linear multistep method is called -stable [33, 29, 34, 32] if there exists a real positive definite matrix such that its one-leg 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 (DLN) method is unconditionally -stable, and


where the -coefficients are

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 (one-step midpoint) and (two-step midpoint) rules: the numerical dissipation vanishes if and only if .


  • [1] R. D. Grigorieff, Stability of multistep-methods on variable grids, Numer. Math. 42 (3) (1983) 359–377. doi:10.1007/BF01389580.
  • [2] G. G. Dahlquist, W. Liniger, O. Nevanlinna, Stability of two-step 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, Springer-Verlag, New York-Heidelberg, 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 time-stepping algorithm for the unsteady Stokes/Darcy model, J. Comput. Appl. Math. 394 (2021) 113521, 14. doi:10.1016/
  • [6] G. Dahlquist, Error analysis for a class of methods for stiff non-linear 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 Second-Order Partitioned Method for Fluid-Thick Structure Interaction Problems, J. Math. Fluid Mech. 23 (3) (2021) 64. doi:10.1007/s00021-021-00593-z.
  • [9] W. Liniger, The -contractive second-order 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/1520-0493(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. WMO-IUGG Symp. on NWP, Tokyo, Japan Meteorological Agency  (VII) (1969) 19–24.
  • [12] P. D. Williams, Achieving seventh-order amplitude accuracy in leapfrog integrations, Mon. Wea. Rev. 141 (2013) 3037–3051. doi:10.1175/MWR-D-12-00303.1.
  • [13] Y. Li, C. Trenchea, A higher-order Robert–Asselin type time filter, J. Comput. Phys. 259 (2014) 23–32. doi:10.1016/
  • [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/s10543-018-0695-z.
  • [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 time-stepping 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/jnma-2019-0081.
  • [18] C. Trenchea, Partitioned conservative, variable step, second-order method for magneto-hydrodynamics in Elsässer variables, ROMAI J. 15 (2) (2019) 117–137.
  • [19] M. Bukač, C. Trenchea, Adaptive, second-order, unconditionally stable partitioned method for fluid-structure interaction, Tech. rep., University of Pittsburgh (2021).
  • [20] D. S. Watanabe, Q. M. Sheikh, One-leg 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, Prentice-Hall, 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, One-leg integration of ordinary differential equations with global error control, Comput. Methods Appl. Math. 5 (1) (2005) 86–96. doi:10.2478/cmam-2005-0004.
  • [26] G. Y. Kulikov, S. K. Shindin, One-leg variable-coefficient formulas for ordinary differential equations and local-global step size control, Numer. Algorithms 43 (1) (2006) 99–121. doi:10.1007/s11075-006-9043-5.
  • [27] J. D. Lambert, Computational methods in ordinary differential equations, John Wiley & Sons, London-New York-Sydney, 1973, Introductory Mathematics for Scientists and Engineers.
  • [28] G. G. Dahlquist, On one-leg 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 Runge-Kutta 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, Springer-Verlag, Berlin, 2010, stiff and differential-algebraic problems, Second revised edition. doi:10.1007/978-3-642-05221-7.
  • [33] G. G. Dahlquist, On the relation of G-stability to other stability concepts for linear multistep methods, Dept. of Comp. Sci. Roy. Inst. of Technology Report TRITA-NA-7621 (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 York-London, 1978, pp. 1–29.
  • [35] G. G. Dahlquist, On stability and error analysis for stiff non-linear problems, Part I, Dept. of Comp. Sci. Roy. Inst. of Technology Report TRITA-NA-7508 (1975).