PaperlayerView exact revision on arXiv

Source-linked AI summary

Optimized linear multistep methods for stiff problems

Tyler Jensen, David I Ketcheson, Lajos Lóczi

arXiv:2610.22685v1math.NAcs.NA

TL;DR

A-stable linear multistep methods cannot be more than second-order accurate, motivating the search for methods that also damp stiff modes while keeping error constants small. The paper optimizes L-stable and prescribed-damping methods; its five-step BDFL52 error constant lies within 10^-6 of the theoretical global optimum.

  • Problem

    A-stable linear multistep methods cannot exceed second-order accuracy, while stiff problems require methods that damp stiff modes.

    Source: §1 · p. 2, §4 · p. 14

  • Method

    The paper uses symbolic techniques and stability-verification algorithms to optimize A-stable second-order methods with full or prescribed stiff-mode damping.

    Source: §1 · p. 2, §4 · p. 14, §4.2.1 · p. 18

  • Results

    The BDFL52 error constant differs by at most 10^-6 from the theoretical global optimum.

    Source: §3.3 · p. 13

  • Takeaways & Limitations

    The optimized methods include L-stable three- and four-step schemes, an essentially optimal five-step scheme, and methods with prescribed nonzero damping.

    Source: §3.3 · p. 13, §4 · p. 14, §4.2.1 · p. 18

    In the Van der Pol example, overcoming the stated issue would require an extremely small step size, making the proposed methods very costly.

    Source: §5.2 · p. 25

Abstract

from arXiv · show

We design and test several new optimized linear multistep methods for stiff problems. Each method is A-stable, implicit, and second-order accurate. We focus on achieving the smallest possible error constant while also ensuring that stiff modes are damped. We recover some known optimal L-stable three- and four-step methods, and provide new and improved five-step methods. We also develop new methods that are not L-stable but provide a prescribed rate of damping in the infinitely stiff limit. Variable-step-size formulae are also provided for some of the new methods. Numerical tests demonstrate that the improved error constants lead to improved accuracy.

1 Background and motivation

The paper targets second-order A-stable linear multistep methods that balance small error constants against strong damping of stiff modes. It develops optimized L-stable and partially damped methods, alongside symbolic verification tools and variable-step implementations.

1 Background and motivation

  • Because no A-stable linear multistep method can exceed second-order accuracy, the paper focuses on implicit methods of order two.The restriction follows from Dahlquist’s second barrier [5, Section V.1].Source: §1 · p. 2
  • The authors seek A-stable second-order methods with small error constants and strong stiff-mode damping, recognizing that these are conflicting optimization criteria.Relaxing L-stability can reduce the error constant, while stronger damping can increase it.Source: §1 · p. 2
  • For three to five steps, symbolic techniques identify methods with the smallest error-constant modulus among A-stable L(0)-stable methods; the five-step result is within 10^-6 of the theoretical optimum.The five-step method’s optimality is not proved, and the error constant approaches the trapezoidal lower bound as the step count increases.Source: §1 · p. 3
  • The work also recovers Vatsa et al.’s optimal BDF2OPT(4) and BDF2OPT(5) methods as BDFL32 and BDFL42 within the corresponding special classes [13, Section III].The paper takes a more general approach by considering broader classes of linear multistep methods.Source: §1 · p. 3
  • For prescribed damping factors 0 < R∞ < 1, the paper optimizes two-step methods and three restricted families of three-step methods, and supplies variable-step implementations for optimal BDFL methods.Exploring all three-step methods would be computationally demanding.Source: §1 · pp. 3–4

1.1 Notation and preliminaries

  • The preliminaries establish polynomial-root terminology and recursive Schur–Cohn tests used to reason about roots inside disks and certify stability.The paper also presents algorithms to test A-stability and L(κ)-stability.Source: §1.1 · p. 4

2 A-stable LMMs with small error constant and damping factor

This section defines the accuracy and stability criteria used to design implicit linear multistep methods, then formulates optimization of the error constant under those criteria. It also develops tests for A-stability and prescribed asymptotic damping.

2 A-stable LMMs with small error constant and damping factor

  • The section sets up the initial-value ODE and defines a k-step linear multistep method through its α- and β-coefficients; βk≠0 makes it implicit.The methods considered are irreducible, since reducible k-step methods are essentially equivalent to methods with fewer steps.Source: §2 · pp. 5–6

2.1 Local order of consistency and the error constant

  • The error constant is defined for an order-p method using its consistency residual and the second characteristic polynomial, and is well-defined under consistency and zero-stability.Second-order accuracy requires the local truncation error to be O(h^2).Source: §2.1 · p. 6

2.2 Some stability concepts

  • For implicit methods, the damping factor as z approaches infinity equals the modulus of the dominant root of the second characteristic polynomial.The asymptotic damping factor defines L(κ)-stability by requiring R∞≤κ; κ=0 gives L-stability when the method is also A-stable.Source: §2.2 · pp. 7–8

2.3 Characterization of A-stable methods

  • A-stability is non-local, and a root-locus curve contained in the closed right half-plane is not sufficient to establish it.The section gives a second-order counterexample and notes that the root-locus curve generally differs from the boundary of the absolute-stability region.Source: §2.3 · p. 8

2.4 Characterizations of L(κ)-stability

  • For κ>0, testing L(κ)-stability uses the reduced polynomial of Q(ζ)=σ(κζ) and a recursive decision process.The recursive test reduces polynomial degree until a linear polynomial remains; parameterized β-coefficients can make successive polynomials increasingly complicated.Source: §2.4 · p. 10

2.5 Optimization

  • The optimization seeks the smallest error constant for second-order methods subject to A-stability, implicitness, and structural coefficient constraints.For larger step counts, the authors may restrict the search to lower-dimensional subsets, which can yield useful methods without guaranteeing global optimality.Source: §2.5 · p. 11

3 Optimal 2nd-order A-stable and L(0)-stable BDFL methods with k ∈{3, 4, 5} steps

The section optimizes second-order A-stable, L(0)-stable BDFL methods for three to five steps, minimizing their error constants while checking stability and irreducibility. The five-step method is essentially optimal, while the three-step method halves BDF2’s error constant.

3.1 Three-step methods

  • For three steps, β3=3/5 minimizes |C| at 1/6 among A-stable candidates, half BDF2’s error constant, and recovers a method proposed in.The method is irreducible when β3 ≠ 2/3, and the A-stable parameter range is [3/5, 1].Source: §3.1 · p. 11

3.3 Summary

  • Table 1 collects the optimal three- and four-step methods and essentially optimal five-step method with their error constants; Figures 1–4 show error comparisons and stability regions.Figure 1 compares error constants and step counts; Figures 2–4 display absolute-stability regions.Source: Table 1 · p. 12, §3.3 · p. 12, Figure 1 · p. 13

3.2 Four- and five-step methods

  • For four- and five-step searches, the authors minimize over a necessary-condition superset and verify the candidate with full A-stability checks; the five-step search uses a 1/500 grid.The grid produces 1001 linear inequalities defining a superset of the A-stable methods.Source: §3.2 · p. 12, §3.3 · p. 13

3.3 Summary

  • The BDFL52 error constant is within 10^-6 of the theoretical global optimum, but the method is established only as essentially optimal.The five-step parameter space has more than 200 components, and some were discarded when minimization did not finish in reasonable time.Source: §3.3 · p. 13

4 Optimal A-stable and L(κ)-stable methods

The section characterizes optimal A-stable methods with prescribed stiff-mode damping, first for two-step methods and then for restricted three-step families. Across methods with a fixed step count, reducing the damping factor and reducing the error constant trade off.

4.1 2-step 2nd-order methods

  • For irreducible two-step methods, Corollary 4.2 identifies the optimal A-stable, L(κ)-stable family for κ∈[0,1); the strict-class search yields only an infimum.The section presents the resulting optimal methods and notes that the earlier strict-class minimum is not attained.Source: §4.1 · pp. 16–17

4.2 Some optimal 3-step 2nd-order A-stable and L(κ)-stable methods

  • Because the full three-step parameter search is more demanding, the analysis restricts attention to selected structural families.The paper studies three restricted cases, rather than the full three-step space.Source: §4.2 · p. 17

4.2.1 Case 1: β0 = β1 = 0

  • With β0=β1=0, the authors identify optimal irreducible three-step methods for prescribed κ; the κ=0 case recovers BDFL32, while κ=1 gives a reducible method equivalent to trapezoidal.The methods are among those selected for presentation in Table 2.Source: §4.2.1 · p. 18

4.2.2 Case 2: β0 = β2 = 0

  • With β0=β2=0, necessary conditions characterize the A-stable region exactly, and the error constant is minimized on its lower curved boundary.The resulting optimal error-constant graph is presented in Figure 9.Source: §4.2.2 · p. 20

4.2.3 Case 3: β1 = β2 = 0

  • With β1=β2=0, the optimal irreducible three-step methods share the error constant |C|=1/6 for every κ∈[0,1].Source: §4.2.3 · p. 20

4.3 Some optimized methods

  • For each fixed step count, reducing κ and reducing the error constant |C| trade off.The paper summarizes this tradeoff for the optimized two- and three-step methods.Source: §4.3 · p. 20

5 Numerical tests

The tests show that optimized methods can outperform BDF2 after stiff transients are damped, while fixed-step methods lose accuracy on the van der Pol problem. The results also motivate the variable-step-size extension developed next.

5 Numerical tests

  • The tests use TR-BDF2 for the k −1 starting values and reference solutions where analytical solutions are unavailable.Table 2 identifies selected optimal A-stable, L(κ)-stable methods and their error constants; the optimal step count depends on κ.Source: §5 · pp. 20, 22, Table 2 · p. 22

5.1 Prothero–Robinson problem

  • After the initial transient, the L-stable BDFL methods outperform BDF2, and their errors are roughly proportional to their error constants.Errors spike at steps 1 and k, when the methods use initial data far from the slow manifold.Source: §5.1 · p. 22, Figure 10 · p. 23
  • Methods with smaller κ damp the initial oscillations more rapidly; after sufficient damping, even intermediate-κ methods yield smaller errors than BDF2.For the trapezoidal rule, κ = 1, damping is insignificant on the displayed time scale.Source: §5.1 · pp. 22, 24, Figure 11 · p. 24, Figure 12 · p. 25

5.2 Van der Pol oscillator

  • With h = 10^-3, all optimal BDFL methods lose accuracy in under one van der Pol period.The authors state that making the step extremely small would be costly, motivating adaptive step sizes.Source: §5.2 · pp. 24–25

6 Variable-step-size formulation

The paper gives proof-of-concept variable-step-size extensions for the optimal BDFL methods, with stability checks and an error estimator. Tests show improved accuracy for van der Pol, while the authors identify ad hoc coefficient choices as a limitation.

6 Variable-step-size formulation

  • The VSS extension is a proof of concept; more sophisticated coefficients and step-size control are left for future work.The stated goal is to show that a reasonably efficient extension is possible.Source: §6 · p. 25

6.1 Setup for VSS

  • The variable-step formulation indexes each step by its ending interval and uses step-size ratios to express its second-order consistency conditions.With constant step size, the local truncation-error expression reduces to the fixed-step formula.Source: §6.1 · pp. 26–27

6.2 Adaptive BDFL methods

  • For adaptive BDFL methods, the authors fix some fixed-step coefficients and solve the remaining consistency conditions, rather than optimizing error over step ratios.They give a BDFL32 construction and provide BDFL42 and BDFL52 coefficients in Appendix A.Source: §6.2 · p. 27

6.3 Stability

  • Frozen-coefficient root checks give a practical indication of admissible step-ratio variations, but do not establish zero-stability for arbitrary ratio sequences.The exact stable set is determined for BDFL32; for BDFL42 and BDFL52, the stable sets include a neighborhood of equal step ratios.Source: §6.3 · pp. 28–29, Figure 14 · p. 29

6.4 Error estimation and step-size control

  • A cubic Hermite extrapolation estimates local error from the two most recent accepted solution values and derivatives, without additional right-hand-side evaluations.Its O(h^4) error makes the difference an asymptotically correct estimate of the leading O(h^3) local solution error.Source: §6.4 · pp. 29–30

6.5 Variable-step-size methods: Robertson’s problem

  • On Robertson’s problem, errors for the VSS methods are roughly proportional to their error constants under the chosen logarithmically spaced step sequence.The sequence uses 2000 points from 10^-6 to 10^8 seconds, with successive-step ratio about 1.0162569.Source: §6.5 · p. 30, Figure 15 · p. 30, Figure 16 · p. 31

6.6 Van der Pol oscillator revisited

  • Adaptive BDFL methods remain accurate over multiple van der Pol periods and take about 3.38×10^5 steps, versus 2×10^6 fixed-size steps.In the work-precision study, proposed methods improve accuracy over BDF2 at the same cost across most tested tolerances; BDFL4 and BDFL5 reject up to 30% of steps at tight tolerances.Source: §6.6 · pp. 31–32, Figure 17 · p. 32, Figure 18 · p. 32

A Variable-step-size coefficients for the optimal BDFL42 and BDFL52 methods

Statements and declarations

The authors report no relevant competing interests. They make the code used for numerical results and figures available and disclose limited AI assistance in code generation.

Statements and declarations

Data availability statement

Use of artificial intelligence

Loading 2610.22685v1…