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 →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
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.
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
- 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.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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)
- [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.
- [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)
- [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.
- [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.
- [Section 4, Theorem 19] The displayed bound contains a typo: '|χ(|µi)|' should be '|χ(|µi|)|'.
- [Section 3.1] In the formula for μ_{m,z,θ} of IMEX-BDF3, '3e2itheta' should be '3e^{2iθ}'.
- [Section 5.1] The remark after Definition 22 is numbered 'Remark 1', although earlier remarks are numbered 12 and 20; it should be renumbered.
- [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.
- [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
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
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]).
- ad hoc to paper Monotonicity of σ_z as a function of z for IMEX-BDF2 and IMEX-BDF3 (Theorems 15 and 16).
- standard math Standard field of values properties: convexity, spectral containment, and F(X) being a real segment for Hermitian X (Remark 1, from [9]).
- ad hoc to paper Equality of the scalar stability region D_z with the innermost region enclosed by Γ_z (Proposition 14).
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 from the paper (7 more)
Forward citations
Cited by 1 Pith paper
-
Sensitivity of ECG QRS Complexes to His-Purkinje Structure in Computational Heart Models
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
- [1]
-
[2]
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
work page 2000
-
[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
work page 2000
-
[4]
M. Calvo and T. Grande. On the asymptotic stability ofθ-methods for delay differential equations. Numerische Mathematik , 54(3):257–270, 1989
work page 1989
- [5]
- [6]
-
[7]
T. A. Driscoll, N. Hale, and L. N. Trefethen. Chebfun Guide. Pafnuty Publications, Oxford, 2014
work page 2014
- [8]
Show all 23 references
-
[9]
R. A. Horn and C. R. Johnson. Topics in matrix analysis . Cambridge university press, 1991
1991
-
[10]
K. J. in’t Hour. The stability of θ-methods for systems of delay differential equations. Annals of Numerical Mathematics, 1:323–334, 1994
1994
-
[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
1987
-
[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
1978
-
[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
2000
-
[14]
T. Koto. Stability of IMEX Runge–Kutta methods for delay differential equations. Journal of Computational and Applied Mathematics , 211(2):201–212, 2008
2008
-
[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
2009
-
[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
2008
-
[17]
M. C. Mackey and L. Glass. Oscillation and chaos in physiological control systems. Science, 197(4300):287–289, 1977
1977
-
[18]
F. A. Rihan. Delay differential equations and applications to biology . Springer, 2021
2021
-
[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
2024
-
[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
2017
-
[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
2019
-
[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
2000
-
[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
1993
Reviewed August 11, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.