{"id":"0b710f43-948e-4555-9eba-81e3736440c2","arxiv_id":"2501.13515","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"The authors construct and test a new family of high-order implicit integrators for Hamiltonian systems based on the structural method, with strong energy preservation and applicability to non-separable problems.","lead":"This paper adapts the structural method for time integration to Hamiltonian systems, splitting the problem into physical equations (Hamilton's equations and their time derivatives) and stencil-only structural equations. The resulting implicit schemes show high order and very small long-time energy drift in benchmarks from pendulums to the outer solar system and charged particles.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Unconditional stability and fixed-point convergence are asserted rather than proved for the coupled Hamiltonian ZD/ZDS schemes, and the §3.2.2 stiffness caveat undercuts the large-time-step benefit claimed in §5.8 and the conclusions.","rationale":"The reader's weakest-assumption analysis identifies the same general region: scalar-method properties are assumed to carry over to the coupled Hamiltonian system. My stress-test narrows this to the most consequential unverified condition, namely the convergence of the fixed-point iteration and the accompanying unconditional-stability claim. The paper's numerical work is extensive and internally consistent, and the construction is transparent enough to be testable, but the advertised stability property is asserted rather than demonstrated in the Hamiltonian setting. This is a correctness gap rather than a demonstrated failure, so conditional acceptance remains the appropriate verdict. I mark agreement as partial because the reader also lists Az nonsingularity as a joint concern; I consider that less load-bearing, since it is shared with the scalar method and is numerically supported by the reported convergence orders.","tokens_in":38449,"tokens_out":4535,"duration_ms":44688,"concrete_test":"Run the Section 5.7.3 charged-particle benchmark with the ZDS R=4 scheme on a sequence of increasing time steps, e.g., N = 48T, 12T, 6T, 3T, and 1.5T, fixing tol = 1e-14, and record both the fixed-point residual after the published stopping criterion and the Hamiltonian error at T = 20000. If the iteration stalls or diverges at any of these Δt, or if the Hamiltonian error grows without bound where the published fine-grid results are flat, the 'unconditional stability' claim for the Hamiltonian scheme is not supported. If all coarse steps converge and keep the Hamiltonian bounded, the concern is resolved.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The load-bearing gap is that the properties advertised as the paper's main contribution—unconditional stability, high order, and practical efficiency via large time steps—are never proved for the Hamiltonian ZD/ZDS schemes. The stability and order statements appear only as remarks for the scalar ODE in §2.3.1 ('R = 2 provides a 4th-order unconditionally stable scheme') and §2.3.2, and are then carried over to the coupled x-p setting without new analysis. The Hamiltonian algorithm is not a direct application of the scalar analysis: it is a coupled fixed-point iteration in which Zx and Zp are updated from structural equations and then Dx and Dp are updated from the physical equations. Neither the convergence of this iteration nor the stability region of the resulting scheme is analyzed. The paper itself concedes the limitation in §3.2.2: 'for stiff problems, a more sophisticated fixed-point procedure may be implemented.' If this iteration fails to converge for large Δt, the 'unconditional stability' that drives the complexity comparison in Table 17 and the non-separable benchmark conclusions does not hold for the method as implemented. This is not a dispute about symplectic versus non-symplectic philosophy; it is a missing correctness condition for the central claim. The numerical tables are strong for small-to-moderate steps, but they contain no large-Δt stability test, which is exactly the regime the unconditional-stability claim concerns.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper adapts the structural method of Clain, Machado, and Malheiro to Hamiltonian systems, presenting ZD and ZDS variants for scalar, vector, and non-separable Hamiltonians. It derives physical and structural equations, proposes block fixed-point iterations, and benchmarks the resulting schemes against classical symplectic integrators on mass-spring, pendulum, Kepler, three-body, outer-solar-system, and charged-particle problems. The central advertised properties are high-order accuracy, unconditional stability, invariant preservation, and computational efficiency through large time steps.","tokens_in":38711,"tokens_out":5418,"duration_ms":51609,"significance":"If the claimed properties held, the method would be an attractive alternative to composition-based symplectic integrators, particularly for non-separable Hamiltonians and for raising order simply by increasing the block size R. The manuscript provides extensive numerical evidence: convergence orders are consistently reproduced, long-time energy and angular-momentum errors remain bounded in most benchmarks, and the ZDS variant achieves very high accuracy at moderate step sizes. The structural coefficients are derived from polynomial exactness rather than fitted to data, and the numerical studies use very high precision and a fair comparison with classical schemes. However, the central theoretical properties are asserted rather than proved for the Hamiltonian schemes, so the significance of the paper is currently contingent on analysis that is not supplied.","major_comments":[{"comment":"The order and unconditional-stability claims are stated only as remarks for the scalar ODE and then transferred to the Hamiltonian setting. The remark in §2.3.1 ('R = 2 provides a 4th-order unconditionally stable scheme while R = 4 reaches sixth-order accuracy') and the corresponding remark in §2.3.2 are not supported by a proof even for the scalar case in this manuscript, and no theorem establishes the order or stability of the coupled x-p system described in §3.2.1, §3.2.2, and §4.2. Since the Hamiltonian algorithm is a genuinely coupled fixed-point iteration, these properties cannot be taken as inherited from the scalar analysis.","section":"§2.3.1, §2.3.2, §3.2"},{"comment":"The algorithm requires the matrix Az to be non-singular, but the paper only assumes this ('Assuming the matrix Az ∈ RR×R is non-singular') and never proves it for the coefficient sets used for the reported values of R. If Az were singular for some R, the scheme would not be well-defined. The numerical experiments imply invertibility for the specific R tested, but the paper makes no general statement, and this assumption is load-bearing for the Hamiltonian schemes because they use the same structural equations (see §3.2.1, Eq. (14)-(15)).","section":"§2.3.1, Eq. (7)"},{"comment":"Convergence of the fixed-point iteration is not analyzed. The iteration alternates between the structural equations for Zx and Zp and the physical equations for Dx, Dp (and Sx, Sp in the ZDS case), and no contraction argument, step-size condition, or damping strategy is provided. The remark in §3.2.2 concedes that 'for stiff problems, a more sophisticated fixed-point procedure may be implemented,' which directly undercuts the unconditional-stability claim for the method as implemented. No large-Δt test is reported, so the regime in which unconditional stability is claimed is exactly the regime for which no numerical evidence is given.","section":"§3.2.2 and §4.2"},{"comment":"The complexity comparison and the conclusion that the method allows 'very large time steps' rest on unconditional stability and on the fixed-point iteration converging with few iterations. The iteration counts reported in Tables 18-21 are tied to the specific time steps used in the benchmarks; without a stability and convergence analysis, the complexity ratios in Table 17 do not establish the claimed advantage in the large-Δt regime. The statement in §5.8 'One of the key benefits of the structural method is its unconditional stability' is therefore not supported by the evidence in the paper.","section":"§5.8, Table 17 and bullet list"},{"comment":"The paper repeatedly labels the measured convergence rates as 'optimal' and attributes them to theoretical analysis (e.g., §5.1.2: 'the accuracy and order of convergence are the optimal ones, given by the theoretical analysis'). Since no Hamiltonian-specific theoretical analysis is supplied, these rates should be presented as numerical evidence for, not as consequences of, the claimed order. This distinction matters because the order properties are central to the paper's contribution.","section":"§5.1.2 and §5.2.2"}],"minor_comments":[{"comment":"The text 'the 3(R + 1)(R + 1) linear system' should read 'the 3(R + 1) linear system' (or 'the 3(R + 1) equations'), since there are 3(R + 1) unknown coefficients a_{r,d}.","section":"§2.2"},{"comment":"There are several typos: 'stuffer' should be 'suffer' in the Introduction; 'Hamtilonian' should be 'Hamiltonian' in §5.8.1; 'Not that this potential' should be 'Note that this potential' in §5.7.3; 'octa-precision' should be 'octuple precision' in the Conclusions; and 'Strömer-Verlet' is used instead of the standard 'Störmer-Verlet'.","section":"Throughout"},{"comment":"The remark 'very few iterations (3 to 12 in practice)' is inconsistent with the measured iteration counts in Tables 20 and 21, where the average number of fixed-point iterations reaches 28.4 for the ZDS scheme with R = 1 and 15.7 for the ZD scheme with R = 2 in the non-separable benchmark. Please reconcile the remark with the reported data.","section":"§3.2.2"},{"comment":"The classical scheme is labeled 'MA2' in some tables and 'MC2' in others (Tables 4-5 versus Tables 9 and 11), and the legend of Figure 2 also uses 'MC2'. Please use one name consistently.","section":"§5.4 and §5.5"},{"comment":"The manuscript does not state data or code availability. Given the very high precision and the reproducibility requirements of a numerical methods paper, a public implementation or detailed pseudocode for constructing the structural-equation coefficients would be helpful.","section":"§5"}],"recommendation":"major_revision","confidential_remarks":"The stress-test concern raised by the reader is valid and lands on a load-bearing point. The paper's own remarks concede the lack of analysis for stiff problems, and the absence of a large-Δt stability test is a real gap given the headline claim of unconditional stability. The numerical benchmarks are extensive and appear to provide strong evidence for the convergence orders and long-time boundedness at the tested steps, so the manuscript is not hopeless. I would recommend asking for a revised version that either provides a stability and convergence analysis for the Hamiltonian ZD/ZDS schemes, or clearly restates the claims as numerical observations and adds a large-Δt benchmark that tests the fixed-point iteration in the claimed regime."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here's my take. The paper is a real extension of the structural method to Hamiltonian systems, and the numerics are solid and extensive. What's new: the two-variable coupling (x and p), the body-space formulation for many-body problems, and the non-separable case. The structural coefficients come from polynomial exactness, not from fitting the test problems, and the convergence orders in the benchmarks match the claimed values. The comparison against standard symplectic integrators is honest, and the non-separable charged-particle results are genuinely useful: high-order ZDS variants track energy much better than Störmer-Verlet compositions.\n\nThe soft spot is exactly what the stress-test note says. The paper's headline properties—unconditional stability and high order—are stated as remarks for the scalar ODE, then carried over to the coupled Hamiltonian system without new analysis. The fixed-point iteration for the coupled x-p system has no convergence proof; §3.2.2 admits stiff problems may need a more sophisticated solver. The stability and complexity claims in §5.8 and the conclusions lean on this. The numerics all use small-to-moderate steps; there is no large-Δt stability test, which is precisely the regime the unconditional-stability claim concerns. So the empirical work supports high order and energy boundedness for the tested range, but not unconditional stability.\n\nAlso worth noting: the structural methods are not symplectic. They preserve energy to high accuracy on many problems, but in the Kepler and three-body tests the classical symplectic schemes preserve angular momentum to machine precision, while the structural schemes have errors that follow the method order. The Kepler LRL invariant drifts linearly unless you add a projection, and the projection degrades the other invariants. That's not fatal, but it tempers the 'better than symplectic' framing.\n\nFor a reader: people working on high-order implicit integrators, especially for non-separable Hamiltonian problems, will get a lot from the experiments. The paper deserves a serious referee—it's not a desk-reject—but the referee should push for either a proof of the stability/order transfer to the coupled system or a clear statement that these properties are observed numerically, not proven. The authors should also add a large-Δt stability experiment. With those changes it could be a solid journal paper.","headline":"A numerically strong extension of the structural method to Hamiltonian systems, but the headline stability claims are carried over from the scalar case without proof; worth a serious review if the gap is addressed.","tokens_in":39221,"tokens_out":2970,"would_cite":false,"duration_ms":27673,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65P10","37M15","65L06"],"pacs":[],"model":"deepseek-v4-flash","headline":"The structural method adapts to Hamiltonian systems, giving unconditionally stable integrators whose order is set by block size.","keywords":["structural method","Hamiltonian systems","unconditionally stable integrators","high-order accuracy","invariant preservation","compact schemes","fixed-point iteration","symplectic integrators"],"falsifier":"Compute $\\det(A_z)$ for a range of block sizes R: if any R gives a singular matrix, the scheme cannot be built for that R and the universal order/stability claim fails. Alternatively, take the non-separable charged-particle benchmark and increase the time step beyond the Stömer-Verlet limit; if the basic fixed-point iteration diverges for some R, the unconditional-stability claim for the Hamiltonian adaptation is refuted.","tokens_in":38264,"feed_emoji":"🪐","tokens_out":9162,"duration_ms":78166,"temperature":0.7,"pith_summary":"This paper claims that the structural method, a blockwise implicit integrator that separates the physical differential equation from a purely discretization-driven structural equation, can be adapted to Hamiltonian systems while keeping its two signature properties: unconditional stability and accuracy that increases with a block-size parameter R. The adaptation treats position x and momentum p as twin variables that obey the same structural equations, with coupling entering only through Hamilton's equations, and adds a second-derivative variant (ZDS) that is more compact for the same order. If the claim holds, the resulting schemes give a practical dial for accuracy: order R+2 for the ZD variant and 2(R+1) for ZDS, with numerical evidence showing Hamiltonian energy staying flat near 1e-25 over 100,000 time units and stable time steps roughly 20 times larger than Stömer-Verlet in a non-separable charged-particle problem. A sympathetic reader would care because it offers a route to very high-order, energy-preserving simulation without the exponential growth of composition sub-steps used by traditional symplectic methods.","feed_headline":"One block size sets the order of a stable Hamiltonian integrator","feed_subtitle":"Adapted structural schemes hit 4th to 10th order, hold energy flat over 100,000 time units, and use fewer calls than symplectic…","key_machinery":"The load-bearing object is the block of structural equations: for a block of R time steps, one needs R linearly independent relations among the unknown function values and derivatives at the R+1 grid points, with coefficients taken from the kernel of a polynomial-exactness matrix so that the relations are exact for polynomials up to degree R+1 (ZD) or 2(R+1) (ZDS). These coefficients depend only on the uniform time step and block size, not on the Hamiltonian, which is why the same equations can be applied unchanged to x and p and why they can be assembled in a preprocessing stage. The nonlinear Hamiltonian enters solely through the physical equations $D_x=\\partial_p H$, $D_p=-\\partial_x H$ (and their differentiated forms in ZDS), and the coupled block system is solved by a fixed-point iteration whose only matrix operation is the inverse of the block matrix $A_z$, assumed non-singular. Higher accuracy is obtained not by adding composition stages but by increasing R, which is the mechanism behind the stated order formulas and unconditional stability.","core_discovery":"The central discovery is that the structural decomposition, physical equations describing the Hamiltonian dynamics plus structural equations that only involve the grid, survives the transition from a scalar ODE to a coupled Hamiltonian system. For each block of R time steps the same R structural equations are written for position and for momentum, so the nonlinearity of the Hamiltonian lives entirely in the physical equations and the structural part can be precomputed once per grid. The ZDS formulation, which adds second-derivative physical equations, is the more efficient of the two: numerical benchmarks show order 4, 6, 8, and 10 accuracy for R=1, 2, 3, 4, with position errors three orders of magnitude smaller than the classical symplectic schemes at matched order, and Hamiltonian deviation flat in time rather than growing. The paper further claims unconditional stability, demonstrated in practice by running the charged-particle benchmark at time steps roughly 20 times larger than the Stömer-Verlet stability limit, and argues from a complexity analysis that high-order structural schemes need fewer function calls than composition-based symplectic methods.","pith_inferences":["Because the structural coefficients depend only on the grid, the same precomputed matrices could be reused for a Hamiltonian whose parameters change during the simulation, enabling cheap parameter sweeps or time-dependent Hamiltonians; the paper does not explore this.","If the proposed extension to third- and higher-order derivatives is implemented, the method's accuracy should scale far beyond order 10; a direct test would be to compare a hundred-digit benchmark against the analytical two-spring solution, which would verify whether the fixed-point iteration converges at extreme accuracy.","The benchmarks suggest a hybrid use: classical symplectic methods when angular momentum is the invariant of interest (they reach machine precision), structural ZDS when total energy and high accuracy over very long times matter, and projection only when a specific non-energy invariant must be bounded.","A practical implementation detail not fully investigated in the paper is acceleration of the fixed-point iteration by warm-starting each block from the previous block's converged values; the paper notes that a better initialization would cut iterations, and this is directly testable in the non-separable benchmarks."],"forward_implications":["For separable Hamiltonian systems the order of the ZDS scheme is raised simply by increasing the block size, with the numerical tests reaching order 10 at R=4, so no new implementation is needed to change accuracy.","Hamiltonian error in long runs (up to 100,000 time units) stays essentially constant rather than drifting, at levels near 1e-25 in quadruple precision, matching or exceeding what the compared symplectic schemes deliver in the separable benchmarks.","In the non-separable charged-particle test, the structural schemes remain stable with time steps about 20 times larger than the Stömer-Verlet limit, cutting total function calls substantially despite the extra fixed-point iterations.","Projection onto an additional invariant manifold (the Laplace-Runge-Lenz vector in the Kepler problem) bounds that invariant's error but converts the flat error curves of the other two invariants into linear growth, a trade-off the paper documents.","The fixed-point iteration typically needs 1 to 3 iterations per block on the outer solar system at tolerance 1e-15, so the per-step cost of the high-order structural schemes is competitive with composition-based symplectic methods."],"supporting_citations":[{"why":"Supplies the original structural method: the physical/structural equation split, the kernel construction of structural coefficients, and the scalar stability and order results this paper adapts.","marker":"Clain et al. [2023a]"},{"why":"Establishes geometric numerical integration and composition methods that serve as both background and the baseline framework for the classical schemes.","marker":"Hairer et al. [2006]"},{"why":"Provides the Stömer-Verlet integrator, the second-order baseline and the starting point for the composition methods compared in the benchmarks.","marker":"Verlet [1967]"},{"why":"Supplies the second-order symplectic scheme (MC2) used as a comparison method.","marker":"McLachlan and Atela [1992]"},{"why":"Supplies the fourth-order symplectic scheme (CS4) used as a comparison method.","marker":"Sanz-Serna and Calvo [1993]"},{"why":"Supplies the composition constants for the sixth- and eighth-order symplectic schemes (KL6, KL8) used as comparisons.","marker":"Kahan and Li [1997]"},{"why":"Provides the quad-double arithmetic library that lets the paper resolve the very small errors of high-order schemes.","marker":"Hida et al. [2000]"},{"why":"Supplies the DifferentialEquations.jl environment through which the classical symplectic baselines are run.","marker":"Rackauckas and Nie [2017]"}],"fun_headline_variants":["Block size tunes Hamiltonian solver from 4th to 10th order","Precomputed structural equations give stable high-order integrators","Energy holds flat as structural scheme hits order 10","Hamiltonian structural method: fewer calls, energy conserved","One grid stencil sets accuracy for Hamiltonian dynamics"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The paper carries the scalar method's stability and accuracy guarantees over to the coupled position-momentum system without a fresh proof, and the construction needs the block matrix $A_z$ (the matrix that determines the block of new positions) to be invertible and the fixed-point iteration to converge; if either condition fails for some block size or Hamiltonian, the claimed unconditional stability and order do not hold.","fun_headline_variants_meta":{"raw":{"variants":["Block size tunes Hamiltonian solver from 4th to 10th order","Precomputed structural equations give stable high-order integrators","Energy holds flat as structural scheme hits order 10","Hamiltonian structural method: fewer calls, energy conserved","One grid stencil sets accuracy for Hamiltonian dynamics"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000189,"raw_usage":{"total_tokens":1333,"prompt_tokens":939,"completion_tokens":394,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":555,"completion_tokens_details":{"reasoning_tokens":314}},"tokens_in":555,"tokens_out":394,"duration_ms":4130,"temperature":1.0,"reasoning_tokens":314,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T15:51:58.105129+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute $\\det(A_z)$ for a range of block sizes R: if any R gives a singular matrix, the scheme cannot be built for that R and the universal order/stability claim fails. Alternatively, take the non-separable charged-particle benchmark and increase the time step beyond the Stömer-Verlet limit; if the basic fixed-point iteration diverges for some R, the unconditional-stability claim for the Hamiltonian adaptation is refuted.","supporting_citations":[],"review_version":1}