REVIEW 4 major objections 3 minor 1 cited by
An implicit two-stage fourth-order time discretization is derived for Lax-Wendroff-type solvers and shown to be A-stable (hence L-stable) for a specific parameter interval.
Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →
T0 review · deepseek-v4-flash
2026-08-03 19:08 UTC pith:IBUWS3UX
load-bearing objection Genuine two-stage fourth-order implicit scheme with correct order conditions, but the A-stability interval is inferred from a numerical scan rather than proven, so the 'L-stable' headline claim is not yet established. the 4 major comments →
An L-Stable Implicit Two-Stage Fourth-Order Temporal Discretization Scheme for Lax-Wendroff-Type Solvers Applied to Stiff Problems
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
Core claim
The paper constructs an implicit two-stage fourth-order (TSFO) time discretization for Lax-Wendroff-type solvers and claims it is A-stable (hence L-stable) when the free parameter D is set to -C and C lies in [0.018824, 0.045589], with C=0.018924 recommended. The order conditions are derived by Taylor expansion, giving unique stage-1 coefficients and a two-parameter family for stage 2; Lemma 3.1 forces G(∞)=0 via D=-C, and Lemma 3.2 locates the poles of the amplification factor in the right half-plane. A systematic numerical scan of |G(iy)| over C∈[0,0.10], y∈[10^-8,10^4] supplies the A-stability interval. Newton iteration is provided to solve the nonlinear implicit stages.
What carries the argument
The amplification factor G(z) for the scalar model L(u)=λu, with half-step factor R(z)=(1+z/4+z²/48)/(1-z/4+z²/48) that is unimodular on the imaginary axis. After enforcing D=-C, G(z) takes the form (3.11) with denominator factoring into the half-step poles and a quadratic; pole locations follow from (3.22)-(3.23). The maximum modulus principle reduces A-stability to three conditions: analyticity in the closed left half-plane, |G(iy)|≤1 on the imaginary axis (checked numerically), and |G(∞)|≤1 (guaranteed by D=-C). Newton iteration with the Jacobians in (4.6) and (4.12) is the practical solver.
Load-bearing premise
The claim that |G(iy)|≤1 for all real y in the stated C-interval rests on a finite numerical scan, not on an analytic proof.
What would settle it
For C=0.018824, compute the maximum of |G(iy)| over y∈R using high-precision arithmetic or a much finer scan; if any value exceeds 1, the claimed A-stability interval is invalid.
If this is right
- For stiff balance laws, the scheme can be combined with GRP-based spatial reconstruction to form an implicit Lax-Wendroff-type solver that embeds source terms in time derivatives, avoiding operator splitting.
- Two-stage fourth-order accuracy implies roughly half the stages of a four-stage Runge-Kutta method, so to reach a given error the implicit TSFO should require fewer stages per step.
- The parameter C provides a tuning knob; values near 0.018924 give strong damping, and the interval endpoints indicate the effective stability range.
- On the tested benchmarks (Robertson, ozone, van der Pol with ν=100), the scheme remains stable at larger time steps than classical implicit RK4 and yields errors about an order of magnitude smaller at fixed step size.
- The Newton iteration accelerates convergence to the implicit solution, making the method practical for multidimensional problems.
Where Pith is reading between the lines
- Because the A-stability interval comes from a scan over a finite grid, the true stability region may extend further; in particular Lemma 3.2 also admits C<-1/39, so a negative-C stable branch is plausible and worth checking.
- The scheme's reliance on ∂t L suggests it is naturally suited for implicit-explicit (IMEX) splitting where fast stiff reactions are solved implicitly while slow convection remains explicit; the paper's framework could be extended to that hybrid.
- The sharp endpoints 0.018824 and 0.045589 likely mark a transition where |G(iy)| just touches 1 at some finite y; proving that boundary analytically would turn the numerical scan into a theorem.
- If the temporal derivatives are supplied by a GRP solver, the full implicit TSFO-GRP method would unify flux and source coupling, potentially improving the accuracy of reactive compressible-flow simulations near discontinuities.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes an implicit two-stage fourth-order (TSFO) temporal discretization for ODEs and Lax-Wendroff-type solvers. The scheme is derived by Taylor expansion and method of undetermined coefficients, yielding Stage 1 with fixed coefficients and Stage 2 with free parameters C and D. Setting D = -C makes the amplification factor vanish at infinity. The paper invokes the maximum modulus principle to reduce A-stability to analyticity in the closed left half-plane and boundedness of |G(iy)| on the imaginary axis. Lemma 3.2 establishes the analyticity condition, but the imaginary-axis check is replaced by a finite numerical scan over C in [0,0.10] and y in [10^-8,10^4]. This leads to Conclusion 3.3 asserting A-stability (hence L-stability) for C in [0.018824,0.045589] with optimal C=0.018924. A Newton iteration is then formulated, and numerical experiments on stiff benchmarks are reported, claiming fourth-order accuracy in two stages and smaller errors than implicit RK4.
Significance. If the stability and accuracy claims were rigorously established, the proposed scheme would be a useful contribution to implicit Lax-Wendroff/GRP-type time discretizations, potentially offering a two-stage fourth-order implicit alternative with L-stability and no operator splitting. The derivation is transparent, the Newton implementation is concrete, and the source code is made publicly available. However, the central L-stability claim is not proven in the manuscript: the A-stability interval rests on a finite numerical scan, and the numerical experiments contain errors below machine precision without any statement of arithmetic precision. The significance of the contribution is therefore not yet established.
major comments (4)
- [Section 3, Conclusion 3.3] The A-stability interval (3.33) is inferred from a finite numerical scan over C in [0,0.10] and y in [10^-8,10^4], as explicitly stated in the paragraph preceding Figure 3.1: 'Analytically proving |G(iy)|≤1 for all y∈R is complex. Instead, we employ a systematic numerical scanning algorithm.' This is not a proof. A finite grid cannot certify sup_y |G(iy)|≤1; there is no discretization-error bound for the scan, no asymptotic bound for y outside the sampled range, and no coverage of the region C<-1/39 that Lemma 3.2 permits. Since A-stability is the basis for the title's 'L-Stable' claim and for all subsequent stability statements, the central load-bearing result is unverified.
- [Section 5, Tables 1-12] The reported errors are far below double-precision roundoff for O(1) solutions. For example, Table 1 reports an L2 error of 9.618514177224e-19 at Δt=1.953125e-3. Such values are not meaningful unless the computations use extended or arbitrary precision, which is not stated. Moreover, the observed convergence orders fluctuate at intermediate step sizes, e.g., Table 3 shows orders of 1.755 and 5.636 at consecutive refined steps, indicating roundoff or non-asymptotic contamination. This undermines the empirical verification of fourth-order accuracy and the claimed comparison with implicit RK4.
- [Section 5.1, Eq. (5.4)] The first numerical test is the affine system du1/dt=-1000u1+1, du2/dt=-u2+1, which contains a constant source. The stability analysis in Section 3 is performed for the homogeneous linear test equation L(u)=λu, i.e., without source terms. The paper does not justify transferring the A-stability conclusion to this affine problem, for example by analyzing the error equation or incorporating the inhomogeneity into the amplification factor. The observed stability at large Δt is therefore not direct evidence for the theoretical stability claim.
- [Section 3, Eq. (3.34)] The 'optimal' parameter C=0.018924 is selected from the same numerical scan that defines the stability interval, and no quantitative criterion for 'superior numerical damping characteristics' is given. The final method is therefore not parameter-free in the sense claimed; the parameter C is tuned to the scan. This is a secondary but real issue, as the paper presents (3.34) as a definite recommendation without a derived optimality condition.
minor comments (3)
- [Section 5, Eqs. (5.2)-(5.3)] The definition of Error_{L2}(Δt) is repeated twice; the second occurrence should presumably be Error_{L∞}(Δt). The text also refers to 'reference solution' in two different ways; please clarify which definition is used for each table.
- [Section 3, Figure 3.1] The manuscript references 'the figure 3.1' but no actual figure content appears in the text. Either include the figure or remove the reference.
- [Throughout] There are numerous typos and grammatical errors, e.g., 'there is a clear and pressing need to developing', 'the follows', and inconsistent use of 'Chapter' versus 'Section'. A careful editorial pass is needed.
Circularity Check
No circularity found: the scheme is derived by Taylor matching, fourth-order accuracy is checked against external/reference solutions, and the stability scan is an incompleteness concern, not a circular one.
full rationale
The paper's derivation chain is self-contained and does not reduce to its own inputs. Stage coefficients (2.18) and (2.21) are obtained by matching Taylor expansions of the exact evolution (2.5)-(2.6) against the proposed ansatz (2.7)-(2.8), which is the standard constructive derivation of a temporal discretization, not a prediction fitted to itself. The claimed fourth-order accuracy is then verified in Section 5 against exact solutions and a highly refined RK4 reference solution (Tables 1-12), so the accuracy result has independent external content. The L-stability claim in Conclusion 3.3 is built from the scheme's own amplification function: condition G(∞)=0 is derived in Lemma 3.1, pole location in the closed left half-plane is analyzed in Lemma 3.2, and the remaining requirement sup_y |G(iy)|≤1 is checked by direct evaluation of |G(iy)|. That verification is finite and numerical rather than analytic — the paper itself states 'Analytically proving |G(iy)|≤1 for all y∈R is complex. Instead, we employ a systematic numerical scanning algorithm' — but this is a completeness/rigor gap, not circularity: the scanned quantity is exactly the stability condition, and the conclusion is not assumed or implied by construction. Likewise, the 'optimal' C=0.018924 is selected from the scan, but the fourth-order accuracy is independent of C and is confirmed numerically, so the tuning does not force the headline accuracy result. Citations to TSFO methods [1-13] are contextual background and do not carry the load of the derivation; no uniqueness theorem or central premise is imported solely through self-citation. Overall, no step in the derivation reduces to an assumed conclusion.
Axiom & Free-Parameter Ledger
free parameters (1)
- C (and D=-C) =
0.018924 (optimal); interval [0.018824,0.045589]
axioms (3)
- domain assumption L(u) is sufficiently smooth that all required derivatives exist.
- domain assumption Scalar test equation L=lambda u with Re(lambda)<=0 is sufficient to assess stability for systems.
- ad hoc to paper The numerical scan grid C in [0,0.10], y in [10^-8,10^4] is sufficient to certify sup_y |G(iy)|<=1.
read the original abstract
The explicit two-stage fourth-order (TSFO) temporal-spatial coupling method is efficient and compact but suffers severe time-step restrictions for stiff problems with multiple scales. To address Professor Jiequan Li's call for an implicit extension, this paper first constructs an implicit TSFO time discretization scheme using the method of undetermined coefficients and Taylor expansion. Second, using a model equation and the maximum modulus principle, sufficient conditions for L-stability are derived. Third, a Newton iteration accelerates convergence. Numerical experiments on classical stiff benchmarks show that the proposed implicit scheme achieves fourth-order temporal accuracy in two stages. Compared to the classical fourth-order implicit Runge-Kutta method, it allows larger stable time steps and reduces convergence errors by an order of magnitude. More importantly, this implicit scheme can be extended to construct an implicit TSFO temporal-spatial coupling method that captures flow-field correlations and handles strong discontinuities, fundamentally contrasting with method-of-lines approaches. Additionally, it unlocks Lax-Wendroff-type solvers to naturally and synchronously embed both stiff source terms and flow transport into time derivatives, thereby avoiding operator-splitting errors.
Figures
Forward citations
Cited by 1 Pith paper
-
Arbitrary-Order Pad\'e-Closed Anchored Two-Derivative Time Discretizations: $s$ Active Stages, Order $2s$, and $L$-Stability
Anchored two-derivative one-step methods with s active stages attain global order 2s and exact second-subdiagonal Padé L-stability for every s, with no hidden stage poles.
Reference graph
Works this paper leans on
-
[1]
Kolgan, Application of the principle of minimum derivatives to the construction of diference schemes for computing discontinuous solutions of gas dynamics (in Russian), Uch
V. Kolgan, Application of the principle of minimum derivatives to the construction of diference schemes for computing discontinuous solutions of gas dynamics (in Russian), Uch. Zap. TsaGI, Russia 3.6 (1972) 68-77
1972
-
[2]
van Leer, Towards the ultimate conservative difference scheme I
B. van Leer, Towards the ultimate conservative difference scheme I. the quest for monotonicity, Lect. Notes Phys. 18 (1973) 163-168
1973
-
[3]
van Leer, Towards the ultimate conservative difference scheme II
B. van Leer, Towards the ultimate conservative difference scheme II. mononicity and conservation combined in a second order scheme, J. Comput. Phys. 14 (1974) 361-370
1974
-
[4]
Toro and V
E. Toro and V. Titarev, Solution of the Generalised Riemann Problem for Advection-Reaction Equations, Proc. Roy. Soc. London A 458 (2002) 271-281
2002
-
[5]
Titarev and E
V. Titarev and E. Toro, ADER: Arbitrary high order Godunov approach, J. Scientific Computing , 17 (2002) 609-618
2002
-
[6]
Toro and V
E. Toro and V. Titarev, ADER Schemes for Scalar Hyperbolic Conservation Laws with Source Terms in Three Space Dimensions, J. Comput. Phys. 202.1 (2005) 196-215
2005
-
[7]
Balsara, T
D. Balsara, T. Rumpf, M. Dumbser, et al, Efficient, high accuracy ADER-WENO schemes for hydrohynamics and divergencefree magnetohydrodynamics, J. Comput. Phys. 228.7 (2009) 2480-2516
2009
-
[8]
E. Toro, A. Hidalgo, ADER finite volume schemes for nonlinear reaction-diffusion equations, Appl. Numer. Math. 59 (2009) 73-100
2009
-
[9]
Hidalgo, M
A. Hidalgo, M. Dumbser, ADER schemes for nonlinear systems of stiff advection-diffusion-reaction equations, J. Sci. Comput., 48 (2011) 173-189
2011
-
[10]
Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes with time-accurate local time stepping for hyperbolic conservation laws
M.Dumbser. Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes with time-accurate local time stepping for hyperbolic conservation laws. Computer methods in applied mechanics and engineering, 2014, 280: 57-83
2014
-
[11]
Hidlgo, O
M.Dumbser, A. Hidlgo, O. Zanotti. High order space-time adaptive ADER-WENO finite volume schemes for non-conservative hyperbolic systems. Computer methods in applied mechanics and engineering, 2014, 268: 359-387
2014
-
[12]
Boscheri, M
W. Boscheri, M. Dumbser. An efficient quadrature-free formulation for high order Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes on unstructured meshes. Journal of scientific computing, 2014, 48: 173-189
2014
-
[13]
Boscheri, M
W. Boscheri, M. Dumbser. A direct Arbitrary-Lagrangian Eulerian ADER-WENO finite volume scheme on unstructured tetrahedral meshes for conservative and non-conservative hyperbolic systems in 3D. Journal of computational physics, 2014, 275: 484-523
2014
-
[14]
Boscheri, R
W. Boscheri, R. Loubere, M. Dumbser. Direct ArbitraryLagrangian-Eulerian ADER-MOOD finite volume schemes for multidimensional hyperbolic conservation laws. Journal of computational physics, 2015, 292: 56-87
2015
-
[15]
Dumbser, R
W.Boscheri, M. Dumbser, R. Loubere. Cell centered direct Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes for nonlinear hyperelasticity. Computers and fluids, 2016, 134/135: 111-129
2016
-
[16]
Boscheri, M
W. Boscheri, M. Dumbser. High order accurate direct Arbitrary Lagrangian-Eulerian ADER-WENO finite volume schemes on moving curvilinear unstructured meshes. Computers and fluids, 2016, 136: 48-66
2016
-
[17]
Zanotti, F
O. Zanotti, F. Fambri, M. Dumbser, et al. Space-time adaptive ADER discontinuous Galerkin finite element schemes with a posteriori sub-cell finite volume limiting. Computers and fluids, 2015, 118: 204-224
2015
-
[18]
Fambri, M
F. Fambri, M. Dumbser, O. Zanotti. Space-time adaptive ADER-DG schemes for dissipative flows-compressible Navier-Stokes and resistive MHD equations. Computer physics communications, 2017, 220: 297-318
2017
-
[19]
Ben-Artzi and J
M. Ben-Artzi and J. Falcovitz, A Second Order Godunov-Type Scheme for Compressible Fluid Dynamics, J. Comput. Phys. 55 (1984) 1-32
1984
-
[20]
Ben-Artzi and J
M. Ben-Artzi and J. Falcovitz, A High Resolution Upwind Scheme for Quasi 1-D Flows, In Numerical Methods for the Euler Equations of Fluid Dynamics, INRIA, SIAM (1985) 66-83
1985
-
[21]
Ben-Artzi and J
M. Ben-Artzi and J. Falcovitz, Generalized Riemann Problems in Computational Gas Dynamics, Cambridge Univ. Press (2003)
2003
-
[22]
X. Lei, J. Li. Transversal effects of high order numerical schemes for compressible fluid flows. Appl. Math. Mech., 40(3): 343-354, 2019
2019
-
[23]
Ben-Artzi and J
M. Ben-Artzi and J. Li. Hyperbolic Balance Laws: Riemann Invariants and Hyperbolic Balance Laws. Numer. Math., 106: 69-425, 2007
2007
-
[24]
J. Li, T. Liu, Z. Sun. Implementation of the GRP scheme for computing radially symmetric compressible fluid flows. J. Comput. Phys., 228: 5867-5887, 2009
2009
-
[25]
J. Li, Y. Wang, Thermodynamical effects and high resolution methods for compressible fluid flows, J. Comput. Phys. 343 (2017) 340-354
2017
-
[26]
Ben-Artzi and J
M. Ben-Artzi and J. Falcovitz, An upwind second-order scheme for compressible duct flows, SIAM J. Sci. Comput. 7 (1986) 744-768
1986
-
[27]
Ben-Artzi, The generalized Riemann problem for reactive flows, J
M. Ben-Artzi, The generalized Riemann problem for reactive flows, J. Comput. Phys. 81.1 (1989) 70-101
1989
-
[28]
J. Li, G. Chen, The generalized Riemann problem method for the shallow water equations with bottom topography, Int. J. Numer. Methods Eng. 65.6 (2006) 834-862
2006
-
[29]
F. Zhou, G. Chen, S. Noelle, et al, A well-balanced stable generalized Riemann problem scheme for shallow water equations using adaptive moving unstructured triangular meshes, Numer. Meth. Fluids, 73.3 (2013) 266-283
2013
-
[30]
J. Qi, J. Li. A fully discrete ALE method over untwisted time-space control volumes. Numerical methods in fluids, 2017, 83(8): 625-641
2017
-
[31]
X. Lei, J. Li, A staggered-projection Godunov-type method for the Baer-Nunziato two-phase model, J. Comput. Phys. 437 (2021) 110312
2021
-
[32]
TW. Chen, Z. Du, Generalized Riemann problem method for the Kapila Model of compressible multiphase flows, Phys. Fluids 37.7 (2025) 076149
2025
-
[33]
Z. Huo, Z. Jia, A GRP-based tangential effects preserving, high resolution and efficient ghost fluid method for the simulation of two-dimensional multi-medium compressible flows, Comput. Fluids 276 (2024) 106261
2024
-
[34]
Zhang, W
Q. Zhang, W. Sheng, The generalized Riemann problem scheme for a laminar two-phase flow model with two-velocities, J. Comput. Phys. 506 (2024) 112929
2024
-
[35]
Hazra, D
A. Hazra, D. Balsara, P. Chandrashekar, et al, Multidimensional generalized Riemann problem solver for Maxwell's equations, J. Sci. Comput. 96 (2023) 26
2023
-
[36]
K. Wu, H. Tang, A direct Eulerian GRP scheme for spherically symmetric general relativistic hydrodynamics, SIAM J. Sci. Comput. 38 (2016) B458-B489
2016
-
[37]
Z. Yang, P. He, H. Tang, A direct Eulerian GRP scheme for relativistic hydrodynamics: one-dimensional case, J. Comput. Phys. 230.22 (2011) 7964-7987
2011
-
[38]
Z. Yang, H. Tang, A direct Eulerian GRP scheme for relativistic hydrodynamics: Two-dimensional case, J. Comput. Phys. 231.4 (2012) 2116-2139
2012
-
[39]
Sheng, Q
W. Sheng, Q. Zhang, Y. Zheng, A direct Eulerian GRP scheme for a blood flow model in arteries, SIAM J. Sci. Comput. 43 (2021) A1975-A1996
2021
-
[40]
Kuang, H
Y. Kuang, H. Tang, Second-order direct Eulerian GRP schemes for radiation hydrodynamical equations, Comput. Fluids 179 (2019) 163-177
2019
-
[41]
Ben-Artzi and J.Q
M. Ben-Artzi and J.Q. Li, Consistency of finite volume approximations to nonlinear hyperbolic balance laws, Math. Comp. 90.327 (2020) 141-169
2020
-
[42]
Ben-Artzi , J.Q
M. Ben-Artzi , J.Q. Li, Regularity of fluxes in nonlinear hyperbolic balance laws, Commun. Appl. Math. Comput. 5.3 (2023) 1289-1298
2023
-
[43]
Runge, U ber die numerisehe Aufl o sing yon Differentialgleichungen, Math
C. Runge, U ber die numerisehe Aufl o sing yon Differentialgleichungen, Math. Ann. 46 (1895) 167-178
-
[44]
Kutta, Beitrag zur n a herungsweisen Integration totaler Differentialgleichungen, Z
W. Kutta, Beitrag zur n a herungsweisen Integration totaler Differentialgleichungen, Z. Math. Phys. 46 (1901) 435-453
1901
-
[45]
Kuntzmann, Neure Entwicklungen der Methoden yon Runge und Kutta, Z
J. Kuntzmann, Neure Entwicklungen der Methoden yon Runge und Kutta, Z. Angew. Math. Mech. 41 (1961) T28-T31
1961
-
[46]
Butcher, Implicit Runge-Kutta processes, Math
J.C. Butcher, Implicit Runge-Kutta processes, Math. Comp. 18 (1964) 50-64
1964
-
[47]
Li and Z
J. Li and Z. Du, A two-stage fourth order time-accurate discretization for Lax Wendroff type flow solvers I. Hyperbolic conservation laws, SIAM J. Sci. Comput. 38 (2016) A3046-A3069
2016
-
[48]
Li, Two-stage fourth order: Temporal-spatial coupling in computational fluid dynamics (CFD), Aero
J. Li, Two-stage fourth order: Temporal-spatial coupling in computational fluid dynamics (CFD), Aero. Adv. 1.3 (2019)
2019
-
[49]
Cheng, Z
J. Cheng, Z. Du, X. Lei, Y. Wang and J. Li, A two-stage fourth-order discontinuous Galerkin method based on the GRP solver for the compressible Euler equations, Comput. Fluids 181 (2019) 248-258
2019
-
[50]
Du and J.Q
Z. Du and J.Q. Li, A Hermite WENO reconstruction for fourth order temporal accurate schemes based on the GRP solver for hyperbolic conservation laws, J. Comput. Phys. 355 (2018) 385-396
2018
-
[51]
F. Zhao, X. Ji, W. Shyy and K. Xu, Compact higher-order gas-kinetic schemes with spectral-like resolution for compressible flow simulations, Adv. Aerodyn. 1.1 (2019) 34
2019
-
[52]
Yuan and H
Y. Yuan and H. Tang, On the explicit two-stage fourth-order accurate time discretizations, J. Comput. Math. 41.2 (2023) 305-324
2023
-
[53]
L. Pan, K. Xu, Q. Li and J. Li, An efficient and accurate two-stage fourth-order gas-kinetic scheme for the Euler and Navier-Stokes equations, J. Comput. Phys. 326 (2016) 197-221
2016
-
[54]
Y. Chen, C. Chen, F. Xiao, X. Li and X. Shen, A two-stage fourth-order multi-moment global shallow water model on cubed sphere, Mon. Wea. Rev. 148.10 (2020) 4267-4279
2020
-
[55]
Zhang, Q
C. Zhang, Q. Li, Z. Wang, J. Li and S. Fu, A two-stage fourth-order gas kinetic CPR method for Navier-Stokes equations on triangular meshes, J. Comput. Phys. 451 (2022) 110830
2022
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.