Pith. sign in

REVIEW 2 major objections 7 minor 1 cited by

On the stability of IMEX BDF methods for DDEs and PDDEs

T0 review · 2 major / 7 minor · reviewed 2026-08-11 · deepseek-v4-flash

Pith's one-line read IMEX-BDF2 and BDF3 are stable for non-commuting delay systems when a scaled field of values lies in a disk, with explicit step bounds and unconditional stability inside radii $1/3$ and $1/7$.

desk verdict First field-of-values stability bounds for IMEX-BDF2/3 on delay systems with non-commuting matrices, but the key monotonicity claim is only numerically verified, so the headline step-size bound is not yet proven. read the letter →

arxiv 2412.12297 v1 pith:2ELESP3L submitted 2024-12-16 math.NA cs.NA

classification math.NAcs.NA MSC 65L2065M2065L06
keywords numericalstabilityIMEXBDFmethodsdelaydifferentialequationsunconditionalfieldofvaluespartialmethodlinesstep-sizebounds
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper studies when IMEX-BDF2 and IMEX-BDF3 remain stable for systems of delay differential equations of the form $y'(t)=-Ay(t)+By(t-\tau)$, where $A$ is positive definite Hermitian and $B$ is arbitrary. The classical scalar test equation covers only the case where $A$ and $B$ diagonalize simultaneously, so the paper extends the analysis to non-commuting matrices by replacing eigenvalues with the field of values of the scaled matrix $A^{p/2-1}BA^{-p/2}$. The main result is a sufficient condition: if that field of values fits inside a disk of radius $r$, then the methods are stable for $h\le|\chi(r)|/\lambda_d$ (BDF2) or $h\le|\tilde\chi(r)|/\lambda_d$ (BDF3), where $\lambda_d$ is the largest eigenvalue of $A$; inside disks of radius $1/3$ (BDF2) or $1/7$ (BDF3), stability holds for every step size. This gives a computable stability test for DDE systems and, after the method of lines, for parabolic partial delay differential equations.

What carries the argument

The central machinery is the field of values $F(X)=\{\langle x,Xx\rangle:\|x\|=1\}$, the set of all Rayleigh quotients of a matrix; for Hermitian $A$ it is a real interval, and the scaled field $F_p=F(A^{p/2-1}BA^{-p/2})$ plays the role of the eigenvalues of $A^{-1}B$ when the two matrices do not commute. The paper couples this with the scalar functions $\psi(z)$ and $\tilde\psi(z)$, which for each negative $z$ give the radius $\sigma_z$ of the largest disk around the origin contained in the scalar stability region of BDF2 or BDF3, respectively. Their inverses $\chi(r)$ and $\tilde\chi(r)$ convert a field-of-values radius $r$ into the maximal allowed $|z|$; dividing by $\lambda_d$ turns that into the step-size bound. The load-bearing monotonicity of $\psi$ and $\tilde\psi$ on $(-\infty,0)$ is what allows Lemma 18 to compare stability disks and Theorem 24 to use only the most negative $z$.

What would settle it

Directly compute the analytical derivative of $\psi_2(z)$ or $\tilde\psi_2(z)$ on the stated intervals and evaluate it at two points $z_2<z_1<0$; any negative value, or any pair with $\sigma_{z_2}>\sigma_{z_1}$ from a direct $\min_\theta|\mu(z,\theta)|$ calculation, would break Lemma 18 and invalidate the step-size bound in Theorem 24.

Watch

Extended reading notes

Core claim

On its own terms, the paper establishes Theorem 24: for the linear delay system $y'(t)=-Ay(t)+By(t-\tau)$ with $A$ Hermitian positive definite, let $\lambda_d$ be the largest eigenvalue of $A$, and suppose that for some real $p$ the set $F_p=F(A^{p/2-1}BA^{-p/2})$ is contained in the disk $D(0,r)$ with $0<r\le1$. Then IMEX-BDF2 is stable for $h\le|\chi(r)|/\lambda_d$ and IMEX-BDF3 for $h\le|\tilde\chi(r)|/\lambda_d$, with $\chi$ and $\tilde\chi$ the inverse functions built from the scalar stability radii in Section 3. Since the scalar radius functions tend to $1/3$ and $1/7$ as the scaled step $z$ tends to $-\infty$, the paper also obtains unconditional stability whenever $F_p\subseteq D(0,1/3)$ for BDF2 or $F_p\subseteq D(0,1/7)$ for BDF3. The proof passes through a scalar characteristic equation by multiplying with $A^{p-1}$ and forming Rayleigh quotients, then uses the monotonicity of the scalar radius functions to reduce the entire interval of possible $z$ values to the single most negative value $z_d=-\lambda_d h$. The simultaneous-diagonalization case is treated first and follows as a special case, and the numerical experiments apply the theorem to DDE systems and to parabolic PDDEs semidiscretized by the method of lines.

Load-bearing premise

The load-bearing premise is that the scalar functions $\psi(z)$ and $\tilde\psi(z)$ are non-decreasing for all negative $z$; this is checked in the paper by plotting numerical derivatives (Figure 4) rather than by an analytic proof, and the step-size bound would not follow if this monotonicity failed.

Editorial extensions

If this is right

  • Inside $D(0,1/3)$ for BDF2 or $D(0,1/7)$ for BDF3, no step-size restriction is needed for linear stability.
  • For non-commuting $A$ and $B$, the stability check becomes a numerical-radius computation for $A^{p/2-1}BA^{-p/2}$ instead of a simultaneous-diagonalization eigenvalue problem.
  • If $A^{-1}B$ is normal and the delay system is stable, the explicit bounds are $h\le1/(\lambda_d\sqrt2)$ for BDF2 and $h\le0.722965/\lambda_d$ for BDF3.
  • Semidiscretizing a parabolic PDDE by the method of lines produces a DDE system to which Theorem 24 applies directly, with $A$ the diffusion matrix and $B$ the delayed reaction matrix.
  • All the stability conditions are sufficient but not necessary, and the numerical examples show stability for some step sizes larger than the proved bounds.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Replacing the containing disk by the convex hull of $F_p$, or by a minimal enclosing ellipse, is a natural way to sharpen the step-size bound; the paper does not explore this.
  • The same construction of $\psi$ and $\chi$ functions could be repeated for higher-order IMEX-BDF methods, giving analogous sufficient conditions if the scalar stability regions are computed.
  • For the nonlinear delayed Burgers example, the stability analysis linearizes the delayed term at the initial history; a complete nonlinear stability proof remains open, and the observed stability for $h>h^*$ suggests the linear bound is conservative.
  • In the method-of-lines setting $\lambda_d$ grows like $(\Delta x)^{-2}$, so the proved bound should scale like $(\Delta x)^2$; testing that scaling would show how sharp the sufficient condition is for parabolic problems.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

2 major / 7 minor

Summary. The paper analyzes the stability of IMEX-BDF2 and IMEX-BDF3 methods applied to the linear delay test equation y'(t) = -A y(t) + B y(t-τ), where A is positive definite and B is arbitrary. It first studies the scalar characteristic equation, derives the functions ψ(z) and ψ̃(z) that give the minimal stability radius σ_z, and obtains sufficient stability conditions for the simultaneously diagonalizable case. The central contribution is Theorem 24 in Section 5, which uses the field of values to prove that if the scaled field of values F_p is contained in a disk of radius r, then both methods are stable for h ≤ |χ(r)|/λ_d (IMEX-BDF2) and h ≤ |χ̃(r)|/λ_d (IMEX-BDF3), with unconditional stability when F_p lies inside D(0,1/3) or D(0,1/7), respectively. Numerical experiments for DDEs and PDDEs illustrate the predicted step-size restrictions and convergence orders.

Significance. If Theorem 24 is fully established, the paper gives a genuinely useful and computable sufficient stability test for IMEX-BDF methods beyond the simultaneous-diagonalization setting; the bounds depend only on the largest eigenvalue of A and the radius of a field of values, and the paper provides reproducible symbolic and numerical codes for the constructions. The numerical experiments in Sections 5 and 6 support the qualitative predictions and show that the conditions are not vacuous. The main reservation is that a load-bearing monotonicity assertion in the proof of Theorem 24 is verified only numerically, which leaves the central step-size guarantee not completely proven as stated.

major comments (2)
  1. [Section 3.1, Theorems 15–16, Lemma 18, and Theorem 24] The monotonicity of ψ(z) and ψ̃(z) on (-∞,0) is asserted in Theorems 15 and 16, but the proofs check the sign of the derivatives only numerically ('Numerically, it is easy to observe...') with reference to Figure 4, and for IMEX-BDF3 the derivative of ψ̃2(z) is not even displayed. This monotonicity is exactly what Lemma 18 uses to conclude that D(0,σ_{z2}) ⊆ D(0,σ_{z1}) for z2 < z1 < 0, and Theorem 24 uses that conclusion to replace the entire interval F(hA) by its most negative endpoint z_d. If either function had a local minimum, the claimed bound h* would no longer be sufficient. Please provide an analytic proof of the derivative sign conditions, or an interval-arithmetic/computer-assisted certificate, or alternatively weaken the statement of Theorem 24 by explicitly assuming the numerically observed monotonicity.
  2. [Section 5.1, Eq. (29)] The displayed formula μ = ⟨v,Bv⟩/⟨v,A^p v⟩ omits the factor A^{p-1} before B. The derivation leading to Eq. (28) and the definition of F_p in Eq. (26) require μ = ⟨v,A^{p-1}Bv⟩/⟨v,A^p v⟩. This is not merely a typographical issue: the p=0 examples use F(A^{-1}B), which corresponds to ⟨v,A^{-1}Bv⟩/⟨v,v⟩, not ⟨v,Bv⟩/⟨v,v⟩. As written, the proof of Theorem 24 does not connect the field-of-values condition to the scalar characteristic equation. The formula and the surrounding derivation should be corrected and checked.
minor comments (7)
  1. [Section 2.1, Theorems 3 and 4] The estimates 'O(h)^3' and 'O(h)^4' should read 'O(h^3)' and 'O(h^4)', respectively.
  2. [Sections 4 and 5.1] There are inconsistent cross-references: the proof of Theorem 24 cites 'Proposition 11' where Corollary 11 is meant, and the proof of Corollary 25 cites 'Lemma 11' where Lemma 18 is meant.
  3. [Section 4, Theorem 19] The displayed bound contains a typo: '|χ(|µi)|' should be '|χ(|µi|)|'.
  4. [Section 3.1] In the formula for μ_{m,z,θ} of IMEX-BDF3, '3e2itheta' should be '3e^{2iθ}'.
  5. [Section 5.1] The remark after Definition 22 is numbered 'Remark 1', although earlier remarks are numbered 12 and 20; it should be renumbered.
  6. [Section 4.2, Table 1] In the h=0.05 row for IMEX-BDF2, the second error entry 4.9934·10^2 is inconsistent with the other entries and with the reported convergence rate; please check whether the exponent should be -2.
  7. [Section 3, Proposition 14] The proof of Proposition 14 appeals to Cauchy's integral theorem without presenting the argument in detail; please either expand the proof or cite the exact argument in [19], since the equality D_z = D̃_z is used implicitly in the figures and in the identification of σ_z.

Circularity Check

0 steps flagged · score 1.0 of 10

No circular derivation: the IMEX-BDF stability bounds are computed from explicit scalar characteristic equations and a field-of-values containment argument; overlapping citations are structural rather than definitional.

full rationale

The central claim (Theorem 24) is not equivalent to an input by construction. The scalar stability quantity σ_z is computed from the characteristic equations (15)-(16) by minimizing |µ_{m,z,θ}|; the functions ψ(z) and ψ̃(z) in (17)-(18) are explicit, and the inverse functions χ(r) and χ̃(r) in (19)-(20) are inverted from these scalar stability disks. The field-of-values argument maps each eigenvector v to z ∈ F(hA) = [−λ_d h, −λ_1 h] and to µ ∈ F_p, so h∗ = |χ(r)|/λ_d follows from the monotonicity lemma, not from any fit to the numerical experiments. The experiments in Tables 1-2 and Section 6 are genuine benchmarks and not used to determine h∗. There is reliance on prior work: Proposition 14 is justified by a sketch citing [19, p.8], and Proposition 23 is quoted from [19, Proposition 1] and [20, Proposition 2], one of which shares an author with the present paper. This is self-citation, but the cited result is a parameter-free matrix criterion that does not itself contain the target stability bound, so it does not make the derivation circular. The main caveat is correctness, not circularity: Theorems 15 and 16 assert monotonicity of ψ and ψ̃, but their proofs only state 'Numerically, it is easy to observe' and refer to Figure 4; Lemma 18 and hence the h∗ conversion depend on this monotonicity. There are also typographical slips, such as Eq. (29) writing µ = ⟨v, Bv⟩/⟨v, A^p v⟩ instead of the ⟨v, A^{p−1}Bv⟩ form used in the derivation. These gaps should be repaired, but they do not amount to a self-justifying equation or to a fitted parameter renamed as a prediction.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The stability functions ψ, ψ̃, χ, χ̃ are derived explicitly from the characteristic equations, not fit to data. The paper uses no free parameters; the experimental radii r are computed from fields of values, not chosen to match the stability boundary. The main unproved ingredients are the monotonicity claims and the D_z = innermost-region identification.

assumptions (4)
  • domain assumption Scalar P-stability characterization: for z in the A-stability region, the disk D(0, σ_z) is contained in the stability region D_z (Lemma 10, from [15, Theorem 1]).
    Imported from Koto [15]; underpins the disk-inclusion argument used throughout Sections 4 and 5.
  • ad hoc to paper Monotonicity of σ_z as a function of z for IMEX-BDF2 and IMEX-BDF3 (Theorems 15 and 16).
    The paper's proofs verify positivity of derivatives only numerically and reference plots; the result is treated as an assumption in Lemma 18 and Theorem 24.
  • standard math Standard field of values properties: convexity, spectral containment, and F(X) being a real segment for Hermitian X (Remark 1, from [9]).
    Used to identify F(hA) with the interval [-λ_d h, -λ_1 h] and to justify the field of values substitutions in Theorem 24.
  • ad hoc to paper Equality of the scalar stability region D_z with the innermost region enclosed by Γ_z (Proposition 14).
    The proof is sketched, relying on continuity and a reference to [19] for the Cauchy integral argument; the full argument is not reproduced.

how reviews work

0 comments
Cite this review

Pith. "Pith review of On the stability of IMEX BDF methods for DDEs and PDDEs." pith.science (2026). https://pith.science/paper/2ELESP3L

@misc{pith2026241212297,
  author       = {Pith},
  title        = {Pith review of: On the stability of IMEX BDF methods for DDEs and PDDEs},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/2ELESP3L}},
  note         = {Machine review of arXiv:2412.12297}
}
abstract

In this paper, the stability of IMEX-BDF methods for delay differential equations (DDEs) is studied based on the test equation $y'(t)=-A y(t) + B y(t-\tau)$, where $\tau$ is a constant delay, $A$ is a positive definite matrix, but $B$ might be any matrix. First, it is analyzed the case where both matrices diagonalize simultaneously, but the paper focus in the case where the matrices $A$ and $B$ are not simultaneosly diagonalizable. The concept of field of values is used to prove a sufficient condition for unconditional stability of these methods and another condition which also guarantees their stability, but according to the step size. Several numerical examples in which the theory discussed here is applied to DDEs, but also parabolic problems given by partial delay differential equations with a diffusion term and a delayed term are presented.

Figures

Figures reproduced from arXiv: 2412.12297 by the authors.

Figure 1
Figure 1. Γz with m = 0, u = 0, and z = −1, . . . , −50 (Γ−50 is the most interior, ..., Γ−1 is the most exterior region.) For m ≥ 1, let us set a concrete value z = −1. The regions, which we remember will correspond to the innermost region of the represented curves, are represented in [PITH_FULL_IMAGE:figures/full_fig_p005_1.png] view at source ↗
Figure 2
Figure 2. Γ−1 for different values of m. Proof. Let us first see that the mentioned region D˜ z exists, that is, the curve Γz is closed and bounded. Knowing that e iθ = cos(θ) + isin(θ), we will then have from the Definition 7 that any point on said curve is of the form µ = − (cos(θm) + isin(θm))Ps j=1(αj − zβj )(cos(θj) + isin(θj)) z Ps−1 j=1 β ∗ j (cos(θj) + isin(θj)) Multiplying by Ps−1 j=1(cos(θj) − isin(θj)) to the previ… view at source ↗
Figure 3
Figure 3. Representation of ψ(z) and ψ˜(z) functions. ψ1(z) ' = 0 ψ2(z) ' ψ3(z) ' = 4 3 z2 -14 -12 -10 -8 -6 -4 -2 0.05 0.10 0.15 0.20 0.25 (a) ψ ′ (z) ψ  1(z) ' = 0 ψ  2(z) ' ψ  3(z) ' = 20 21 z2 -14 -12 -10 -8 -6 -4 -2 0.2 0.4 0.6 0.8 1.0 (b) ψ˜′ (z) [PITH_FULL_IMAGE:figures/full_fig_p008_3.png] view at source ↗
Figures from the paper (7 more)
Figure 4
Figure 4. Figure 4: Representation of ψ ′ (z) and ψ˜′ (z) functions. IMEX BDF3 Reasoning in a similar way to the previous method, we arrive at the fact that the minimum function starts, in this case, from the equation of the curve µBDF3 = 2 − 9e iθ + 18e 2iθ − 11e 3iθ + 6e 3iθz 6(1 − 3e i…
Figure 5
Figure 5. Figure 5: Location of the numerical range F((−A) −1B). h2 = 0.005. We can observe numerically (see [PITH_FULL_IMAGE:figures/full_fig_p017_5.png]
Figure 6
Figure 6. Figure 6: Location of the numerical range F(A−1B) for n = 100 and different values of l. From this, it is possible to obtain the value of the step size constraint h ∗ for this problem. That is, taking into account Theorem 24 and the previous representation, we know that by impos…
Figure 7
Figure 7. Figure 7: Numerical solution of (31) given by the IMEX BDF2 method with h = 0.1 and te = 60((a) and (b)), and detail of the numerical solution with t ∈ [0, 2.5]. (a) u1 (b) u2 (c) Detail of u2, t ∈ [0, 2.5] [PITH_FULL_IMAGE:figures/full_fig_p020_7.png]
Figure 8
Figure 8. Figure 8: Numerical solution of (31) given by the IMEX BDF3 method with h = 0.1 and te = 60((a) and (b)), and detail of the numerical solution with t ∈ [0, 2.5]. It is remarkable that, even if we impose step sizes larger than the existing restriction for h in the IMEX BDF3, we f…
Figure 9
Figure 9. Figure 9: Location of the numerical range F(A−1B(z−m)) for n = 100. (a) IMEX BDF2 (b) IMEX BDF3 [PITH_FULL_IMAGE:figures/full_fig_p023_9.png]
Figure 10
Figure 10. Figure 10: Numerical solution of (36) given by each method with h = 0.1 and te = 20. 7 Conclusions In this paper, the stability of IMEX-BDF methods for DDEs is studied in detail for the test equation y ′ (t) = −Ay(t)+By(t−τ ), where τ is a constant delay, A is a positive definit…

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

  1. Sensitivity of ECG QRS Complexes to His-Purkinje Structure in Computational Heart Models

    q-bio.QM 2025-05 conditional novelty 6.0 of 10

    Variations in His-Purkinje structure have little individual effect on simulated QRS morphology, but parameter interactions can produce abnormal and premature QRS complexes.

Reference graph

Works this paper leans on

23 extracted references · 22 canonical work pages · cited by 1 Pith paper

  1. [1]

    Baker, G

    C. Baker, G. Bocharov, A. Filiz, N. Ford, C. Paul, F. Rihan, A. Tang, R. Thomas, H. Tian, and D. Wille. Numerical Modeling by Delay and Volterra Functional Differential Equations . 01 2006

  2. [2]

    Baker and E

    C. Baker and E. Buckwar. Numerical analysis of explicit one-step methods for stochastic delay differential equations. LMS Journal of Computation and Mathematics , 3, 01 2000

  3. [3]

    G. A. Bocharov and F. A. Rihan. Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics , 125(1):183–199, 2000

  4. [4]

    Calvo and T

    M. Calvo and T. Grande. On the asymptotic stability ofθ-methods for delay differential equations. Numerische Mathematik , 54(3):257–270, 1989

  5. [5]

    Calvo, J

    M. Calvo, J. I. Montijano, and L. R´ andez. A note on the stability of time–accurate and highly–stable explicit operators for stiff differential equations. J. Comput. Phys. , 436, 2021

  6. [6]

    Conte, J

    D. Conte, J. Martin-Vaquero, G. Pagano, and B. Paternoster. Stability theory of TASE-Runge- Kutta methods with inexact jacobian. SIAM Journal on Scientific Computing , pages A3628– A3657, 2024

  7. [7]

    T. A. Driscoll, N. Hale, and L. N. Trefethen. Chebfun Guide. Pafnuty Publications, Oxford, 2014

  8. [8]

    Gautschi

    W. Gautschi. Numerical analysis. Springer Science & Business Media, 2011

Show all 23 references
  1. [9]

    R. A. Horn and C. R. Johnson. Topics in matrix analysis . Cambridge university press, 1991

  2. [10]

    K. J. in’t Hour. The stability of θ-methods for systems of delay differential equations. Annals of Numerical Mathematics, 1:323–334, 1994

  3. [11]

    Jackiewicz, A

    Z. Jackiewicz, A. Bellen, and M. Zennaro. Stability analysis of one-step methods for neutral delay-differential equations. Numerische Mathematik , 52(6):605–620, 1987/88

  4. [12]

    C. R. Johnson. Numerical determination of the field of values of a general complex matrix. SIAM J. Numer. Anal. , 15(3):595–602, 1978

  5. [13]

    K¨ uchler and E

    U. K¨ uchler and E. Platen. Strong discrete time approximation of stochastic differential equations with time delay. Mathematics and Computers in Simulation (MATCOM) , 54(1):189–205, 2000

  6. [14]

    T. Koto. Stability of IMEX Runge–Kutta methods for delay differential equations. Journal of Computational and Applied Mathematics , 211(2):201–212, 2008

  7. [15]

    T. Koto. Stability of implicit-explicit linear multistep methods for ordinary and delay differential equations. Frontiers of Mathematics in China , 4:113–129, 2009

  8. [16]

    MacDonald, N

    N. MacDonald, N. MacDonald, C. Cannings, and F. Hoppensteadt. Biological Delay Systems: Linear Stability Theory . Cambridge Studies in Mathematical Biology. Cambridge University Press, 2008

  9. [17]

    M. C. Mackey and L. Glass. Oscillation and chaos in physiological control systems. Science, 197(4300):287–289, 1977

  10. [18]

    F. A. Rihan. Delay differential equations and applications to biology . Springer, 2021

  11. [19]

    Rodr ´ ıguez-Fern´ andez and J

    A. Rodr ´ ıguez-Fern´ andez and J. Mart ´ ın-Vaquero. On the stability ofθ-methods for DDEs and PDDEs. Applied Numerical Mathematics , 2024. 24

  12. [20]

    R. R. Rosales, B. Seibold, D. Shirokoff, and D. Zhou. Unconditional stability for multistep imex schemes: theory. SIAM Journal on Numerical Analysis , 55(5):2336–2360, 2017

  13. [21]

    Seibold, D

    B. Seibold, D. Shirokoff, and D. Zhou. Unconditional stability for multistep imex schemes: Practice. Journal of Computational Physics , 376:295–321, 2019

  14. [22]

    Takeuchi, W

    Y. Takeuchi, W. Ma, and E. Beretta. Global asymptotic properties of a delay sir epidemic model with finite incubation times. Nonlinear Anal., 42(6):931–947, nov 2000

  15. [23]

    Wanner, S

    G. Wanner, S. P. Norsett, and E. Hairer. Solving ordinary differential equations I. Nonstiff problems, volume 8. Springer, 2nd rev. ed. edition, 1993. 25

Pith tools

Reviewed August 11, 2026 · model on record in the stance chip above.