pith. sign in

arxiv: 2605.19942 · v1 · pith:ET6NSYYTnew · submitted 2026-05-19 · 🧮 math.NA · cs.NA

A second-order product-type implicit-explicit Runge-Kutta method preserving unit length and energy dissipation structures for gradient flows of vector fields

Jianan Li , Shuang Liu , Tao Tang , Jiang Yang This is my paper

Pith reviewed 2026-05-20 04:14 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords gradient flows of vector fieldsunit length constraintenergy dissipationIMEX Runge-Kutta methodsharmonic map heat flowprojection methodsstructure-preserving schemessecond-order accuracy
0
0 comments X

The pith

A linear second-order scheme preserves unit length and Dirichlet energy dissipation for gradient flows of unit vector fields

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

Gradient flows of unit vector fields model systems like harmonic map heat flows where vectors must remain exactly unit length while energy dissipates over time. Standard high-order schemes often lose one or both properties after projection steps are added to enforce the constraint. The paper builds a product-type IMEX Runge-Kutta construction that stays linear at second order and pairs it with projection to keep both the length constraint and energy decrease intact at every step. This matters because it supports reliable long-time simulations without artificial growth in energy or drift away from unit length in physical applications.

Core claim

The authors develop a general methodology for constructing product-type IMEX-RK schemes and apply it to obtain a linear second-order method that, together with projection, simultaneously preserves the unit length constraint and dissipates the original Dirichlet energy for gradient flows of unit vector fields. They verify the properties through numerical experiments on the harmonic map heat flow.

What carries the argument

Product-type IMEX-RK methods, which arrange implicit and explicit stages multiplicatively to support structure preservation, combined with a projection step onto the unit sphere.

If this is right

  • The scheme applies to related models such as nematic liquid crystals and magnetization dynamics while retaining the same preservation properties.
  • Long-time integrations remain stable without artificial energy increase or constraint violation.
  • The method achieves second-order accuracy while staying linear and avoiding nonlinear solves at each step.
  • Numerical stability holds under the verified conditions for the target gradient-flow equations.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The product-type construction could extend to other geometric constraints or energies in nonlinear evolution equations.
  • Testing the scheme on three-dimensional domains or with adaptive time stepping would check practical scalability.
  • Comparison against existing nonlinear or fully implicit structure-preserving methods could quantify efficiency advantages.

Load-bearing premise

The product-type IMEX-RK construction combined with projection enforces both the unit-length constraint and strict energy dissipation at second order without order reduction or instability.

What would settle it

A numerical test on the harmonic map heat flow in which the computed vector length deviates from one by more than machine precision or the discrete energy increases between steps would disprove the preservation claims.

Figures

Figures reproduced from arXiv: 2605.19942 by Jianan Li, Jiang Yang, Shuang Liu, Tao Tang.

Figure 1
Figure 1. Figure 1: Absolute stability region S π 2 of the Butcher tableau (3.14) for PRK sheme (3.4) and the Butcher tableau (4.3) for PRK scheme (4.2). 4.2. Convergence Rate Test. Example 4.1. To investigate the convergence rates, we consider the gradient flow (1.1) in the computational domain Ω = (− 1 2 , 1 2 ) 2 with parameters α = β = 1. The initial condition is prescribed as follows: m1(x, y, 0) = 0.3 sin(πx) sin(πy), m… view at source ↗
Figure 2
Figure 2. Figure 2: Example 4.2: Numerical solutions m at t = 0, 0.02, 0.04, 0.048, 0.049, 0.1 using the PRK2 scheme with τ = 10−4 . (a) Initial condition. (b) t = 0.049. (c) t = 0.1 [PITH_FULL_IMAGE:figures/full_fig_p022_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Example 4.2: Numerical solutions m for nodals close to the origin at t = 0, 0.049, 0.1 using the PRK2 scheme with τ = 10−4 . higher accuracy at larger time steps while preserving the correct dissipation behavior. Although a rigorous proof of energy dissipation for scheme (4.2) is left for future work, its numerical performance is virtually identical to that of the proven PRK2 scheme, underscoring the intri… view at source ↗
Figure 4
Figure 4. Figure 4: Example 4.2: Comparison of energy evolution and computational efficiency for SIP1, LM2, and the proposed PRK scheme. The work-precision diagrams are generated using a fixed spatial mesh size h = 1/48 and a sequence of time steps τ = 2 × 10−3/2 j for j = 1, 2, . . . , 7 [PITH_FULL_IMAGE:figures/full_fig_p023_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Example 4.3: Numerical solutions m at t = 0, 0.015, 0.05, 0.6 computed using the PRK2 scheme with τ = 10−3 . Neumann boundary conditions are imposed on the remaining boundaries. Specifically: on the lower slab {(x, y, z) ∈ Ω : z = 0}, the nematic rods are aligned parallel to the x-axis with m = (1, 0, 0)⊤. For the upper slab {(x, y, z) ∈ Ω : z = 0.5}, the rods are uniformly aligned along the y-axis with m … view at source ↗
Figure 6
Figure 6. Figure 6: Example 4.3: Numerical solutions m at {(x, y, z) ∈ Ω : z = 0.5} for t = 0, 0.015, 0.05, 0.6 computed using the PRK2 scheme with τ = 10−3 . REFERENCES [1] J. H. Adler, T. J. Atherton, D. B. Emerson, and S. P. MacLachlan, An energy￾minimization finite-element approach for the Frank–Oseen model of nematic liquid crys￾tals, SIAM J. Numer. Anal., 53 (2015), pp. 2226–2254. [2] G. Akrivis, B. Li, and D. Li, Energ… view at source ↗
Figure 7
Figure 7. Figure 7: Example 4.4: Numerical solutions m at t = 0, 0.05, 0.15, 0.5 computed using the PRK2 scheme with τ = 5 × 10−3 . (a) t = 0. (b) t = 0.05. (c) t = 0.15. (d) t = 0.5 [PITH_FULL_IMAGE:figures/full_fig_p027_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Example 4.4: Numerical solutions m at {(x, y, z) ∈ R 3 : x = y = 0.5} for t = 0, 0.05, 0.15, 0.5 computed using the PRK2 scheme with τ = 5 × 10−3 [PITH_FULL_IMAGE:figures/full_fig_p027_8.png] view at source ↗
read the original abstract

Gradient flows of unit vector fields arise in a wide range of physical models such as harmonic map heat flows, nematic liquid crystals, and magnetization dynamics. Designing numerical schemes that simultaneously preserve the unit length constraint and dissipate energy is essential for reliable simulations of such systems. Although projection methods can effectively enforce the unit length constraint, ensuring energy dissipation under projection, especially in high-order schemes, remains challenging. Unlike traditional implicit-explicit Runge-Kutta (IMEX-RK) methods, in this work we propose a general methodology for constructing product-type IMEX-RK schemes that offers greater adaptability to various models with the goal of designing structure-preserving numerical schemes. For gradient flows of unit vector fields with Dirichlet energy, we design a linear and second-order numerical scheme that simultaneously preserves energy dissipation and the unit length constraint by using product-type IMEX-RK methods and projection techniques. Numerical experiments verify the accuracy, stability, and structure-preserving properties of the scheme. According to our best knowledge, this is the first second-order linear scheme that can preserve both the unit length and the original Dirichlet energy for harmonic map heat flows.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

1 major / 1 minor

Summary. The paper proposes a second-order product-type implicit-explicit Runge-Kutta (IMEX-RK) scheme, combined with a projection step, for gradient flows of unit vector fields such as the harmonic map heat flow. The method is designed to be linear while exactly enforcing the unit-length constraint ||u^{n+1}|| = 1 and preserving the Dirichlet energy dissipation structure. Numerical experiments are presented to confirm second-order accuracy, stability, and the structure-preserving properties, with the authors claiming this is the first such linear second-order scheme for these models.

Significance. If the claimed exact preservation of both the unit length and the original energy dissipation law holds at full second order without order reduction or instability, the result would advance structure-preserving methods for constrained gradient flows arising in physics. The general product-type IMEX-RK construction could offer a flexible framework for other models, enabling reliable long-time simulations that maintain physical invariants.

major comments (1)
  1. [Scheme construction (product-type IMEX-RK with projection)] The central claim requires that the product-type IMEX-RK stages plus final projection simultaneously enforce ||u^{n+1}||=1 exactly and E(u^{n+1}) ≤ E(u^n) - dissipation term at full second-order accuracy. The product structure is intended to cancel cross terms so that the energy law survives the explicit treatment of the nonlinear map and the subsequent projection, but the projection is a nonlinear operation whose first-order perturbation can re-introduce positive contributions to the discrete energy balance unless the stage weights and the specific form of the product-type Butcher tableau are shown to annihilate those contributions up to O(Δt^3). No independent verification (e.g., machine-checked energy identity or parameter-free derivation) is cited, so the claim rests on an algebraic identity whose validity after projection is not obviously robust.
minor comments (1)
  1. [Abstract] The abstract states that numerical experiments verify the properties but does not specify the test problems, time-step ranges, or quantitative measures (e.g., discrete energy decay rates or observed convergence orders) used to confirm energy dissipation and absence of order reduction.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for the positive evaluation of the significance of our work and for the detailed feedback on the scheme construction. We address the major comment point by point below, providing clarifications on the energy identity and projection step.

read point-by-point responses

  1. Referee: [Scheme construction (product-type IMEX-RK with projection)] The central claim requires that the product-type IMEX-RK stages plus final projection simultaneously enforce ||u^{n+1}||=1 exactly and E(u^{n+1}) ≤ E(u^n) - dissipation term at full second-order accuracy. The product structure is intended to cancel cross terms so that the energy law survives the explicit treatment of the nonlinear map and the subsequent projection, but the projection is a nonlinear operation whose first-order perturbation can re-introduce positive contributions to the discrete energy balance unless the stage weights and the specific form of the product-type Butcher tableau are shown to annihilate those contributions up to O(Δt^3). No independent verification (e.g., machine-checked energy identity or parameter-free derivation) is cited, so the claim rests on an algebraic identity whose validity after projection

    Authors: We appreciate the referee's emphasis on rigorously establishing the post-projection energy balance at second order. In the manuscript (Section 3 and Theorem 3.1), the product-type IMEX-RK construction is derived so that the Butcher tableau weights satisfy both the second-order accuracy conditions and the algebraic cancellation of cross terms arising from the explicit nonlinear map and implicit linear operator; this yields an exact discrete energy dissipation E(v^{n+1}) ≤ E(u^n) − Δt · (positive term) for the pre-projection stage value v^{n+1}. The final normalization projection u^{n+1} = v^{n+1}/||v^{n+1}|| enforces ||u^{n+1}|| = 1 exactly. To control the nonlinear perturbation, we note that consistency implies ||v^{n+1}|| − 1 = O(Δt^2). A direct Taylor expansion of the Dirichlet energy around v^{n+1} shows that E(u^{n+1}) − E(v^{n+1}) = O(Δt^3) because the first-order variation vanishes on the sphere and the quadratic remainder is multiplied by the O(Δt^2) deviation; the specific weights annihilate any O(Δt) or O(Δt^2) positive contributions that could violate the inequality. The derivation is parameter-free and follows directly from the order conditions and product structure; we have added an explicit expansion and remark in the revised proof to make this step self-contained. While we do not supply a machine-checked formal proof, the algebraic steps are elementary and can be verified by direct substitution. revision: yes

Circularity Check

0 steps flagged

No significant circularity; scheme construction and preservation claims are self-contained

full rationale

The paper introduces a product-type IMEX-RK construction combined with projection to enforce unit-length and energy dissipation for the target gradient flows. The derivation relies on algebraic cancellation in the Butcher tableau and projection step, presented as a direct consequence of the chosen weights and stages rather than a fitted parameter or self-citation chain. No load-bearing step reduces to a prior result by the same authors that itself assumes the target property. Numerical experiments are cited only for verification, not as the source of the preservation claim. The central assertion of being the first such second-order linear scheme follows from the explicit construction and does not loop back to its own inputs by definition.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The central claim rests on the assumption that a product-type IMEX-RK discretization plus projection can be arranged to inherit both the length constraint and energy dissipation from the continuous gradient flow.

axioms (1)
  • domain assumption The continuous gradient flow of unit vector fields with Dirichlet energy admits a discretization via product-type IMEX-RK and projection that inherits both structures at second order.
    This modeling assumption underpins the scheme construction and is invoked to justify the preservation properties.

pith-pipeline@v0.9.0 · 5740 in / 1166 out tokens · 36234 ms · 2026-05-20T04:14:09.366505+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

What do these tags mean?
matches
The paper's claim is directly supported by a theorem in the formal canon.
supports
The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends
The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses
The paper appears to rely on the theorem as machinery.
contradicts
The paper's claim conflicts with a theorem or certificate in the canon.
unclear
Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

57 extracted references · 57 canonical work pages

  1. [1]

    Akrivis and B

    G. Akrivis and B. Li and D. Li , title =. SIAM J. Sci. Comput. , volume =

    work page

  2. [2]

    François Alouges , title =. Disc. Cont. Dyn. Syst. Ser. S , volume =

    work page

  3. [3]

    Alouges , title =

    F. Alouges , title =. SIAM J. Numer. Anal. , volume =

    work page

  4. [4]

    Tang and X

    T. Tang and X. Wu and J. Yang , title =. J. Sci. Comput. , volume =

    work page

  5. [5]

    An and H

    R. An and H. Gao and W. Sun , title =. SIAM J. Numer. Anal. , volume =

    work page

  6. [6]

    Cohen and R

    R. Cohen and R. Hardt and D. Kinderlehrer and S.-Y. Lin and M. Luskin , title =. Theory and Applications of Liquid Crystals , publisher =

    work page

  7. [7]

    Ramage and E

    A. Ramage and E. C. Gartland , title =. SIAM J. Sci. Comput. , volume =

    work page

  8. [8]

    L. A. Landau and E. Lifshitz , title =. Perspectives in Theoretical Physics , publisher =

    work page

  9. [9]

    Bartels and A

    S. Bartels and A. Prohl , title =. SIAM J. Numer. Anal. , volume =

    work page

  10. [10]

    Bartels and J

    S. Bartels and J. Ko and A. Prohl , title =. Math. Comp. , volume =

    work page

  11. [11]

    Alouges and E

    F. Alouges and E. Kritsikis and J. Steiner and J.-C. Toussaint , title =. Numer. Math. , volume =

    work page

  12. [12]

    Cheng and J

    Q. Cheng and J. Shen , title =. SIAM J. Sci. Comput. , volume =

    work page

  13. [13]

    J. H. Adler and T. J. Atherton and D. B. Emerson and S. P. MacLachlan , title =. SIAM J. Numer. Anal. , volume =

    work page

  14. [14]

    Cheng and J

    Q. Cheng and J. Shen , title =. SIAM J. Numer. Anal. , volume =

    work page

  15. [15]

    Du and S

    Q. Du and S. Liu and J. Yang , title =. J. Sci. Comput. , volume =

    work page

  16. [16]

    Brezis and J.-M

    H. Brezis and J.-M. Coron and E. H. Lieb , title =. Commun. Math. Phys. , volume =

    work page

  17. [17]

    Hardt and D

    R. Hardt and D. Kinderlehrer and F.-H. Lin , title =. Commun. Math. Phys. , volume =

    work page

  18. [18]

    T. Rivi. Everywhere discontinuous harmonic maps into spheres , journal =

    work page

  19. [19]

    C. W. Oseen , title =. Trans. Faraday Soc. , volume =

    work page

  20. [20]

    F. C. Frank , title =. Discuss. Faraday Soc. , volume =

    work page

  21. [21]

    E. G. Virga , title =

    work page

  22. [22]

    L. D. Landau and E. Lifshitz , title =. Phys. Z. Sowjetunion , volume =

    work page

  23. [23]

    Guo and S

    B. Guo and S. Ding , title =

    work page

  24. [24]

    Wang , title =

    C. Wang , title =. Arch. Rational Mech. Anal. , volume =

    work page

  25. [25]

    Bartels , title =

    S. Bartels , title =. SIAM J. Numer. Anal. , volume =

    work page

  26. [26]

    Gui and B

    X. Gui and B. Li and J. Wang , title =. SIAM J. Numer. Anal. , volume =

    work page

  27. [27]

    Pistella and V

    F. Pistella and V. Valente , title =. Numer. Methods Partial Differ. Equ. , volume =

    work page

  28. [28]

    Wang and C

    X.-P. Wang and C. J. Garc. A. J. Comput. Phys. , volume =

    work page

  29. [29]

    2020 , issn =

    Second-order semi-implicit projection methods for micromagnetics simulations , journal =. 2020 , issn =

    work page 2020

  30. [30]

    Chen and C

    J. Chen and C. Wang and C. Xie , title =. Appl. Numer. Math. , volume =

    work page

  31. [31]

    Li and J

    B. Li and J. Yang and Z. Zhou , title =. SIAM J. Sci. Comput. , volume =

    work page

  32. [32]

    E and X.-P

    W. E and X.-P. Wang , title =. SIAM J. Numer. Anal. , volume =

    work page

  33. [33]

    Fu and J

    Z. Fu and J. Yang , title =. J. Comput. Phys. , volume =

    work page

  34. [34]

    Fu and J

    Z. Fu and J. Shen and J. Yang , title =. Sci. China Math. , year =

    work page

  35. [35]

    Fu and T

    Z. Fu and T. Tang and J. Yang , title =. Math. Comp. , year =

    work page

  36. [36]

    Li and Z

    J. Li and Z. Xu and J. Yang and X. Yang , title =. 2026 , journal =

    work page 2026

  37. [37]

    U. M. Ascher and S. J. Ruuth and R. J. Spiteri , title =. Appl. Numer. Math. , volume =

    work page

  38. [38]

    The scalar auxiliary variable (

    Shen, Jie and Xu, Jie and Yang, Jiang , journal=. The scalar auxiliary variable (. 2018 , publisher=

    work page 2018

  39. [39]

    SIAM Rev

    A new class of efficient and robust energy stable schemes for gradient flows , author=. SIAM Rev. , volume=. 2019 , publisher=

    work page 2019

  40. [40]

    and Jin, Shi and Russo, Giovanni , title =

    Caflisch, Russel E. and Jin, Shi and Russo, Giovanni , title =. SIAM J. Numer. Anal. , volume =

    work page

  41. [41]

    Additive. Appl. Numer. Math. , volume =. 2003 , issn =

    work page 2003

  42. [42]

    2016 , publisher=

    Numerical methods for ordinary differential equations , author=. 2016 , publisher=

    work page 2016

  43. [43]

    Boscarino, Sebastiano and Filbet, Francis and Russo, Giovanni , title =. J. Sci. Comput. , volume =

    work page

  44. [44]

    Boscarino, Sebastiano , title =. SIAM J. Numer. Anal. , volume =

    work page

  45. [45]

    Yan, Gui and Rui, Du and Wang, Cheng , title =. Numer. Math. Theory Methods Appl. , year =

    work page

  46. [46]

    High order asymptotic preserving. J. Comput. Phys. , volume =. 2015 , issn =

    work page 2015

  47. [47]

    Cao, Weichen and Yang, Hengli and Chen, Wenbin , title =. J. Sci. Comput. , year =

    work page

  48. [48]

    Du, Qiang and Ju, Lili and Li, Xiao and Qiao, Zhonghua , title =. SIAM J. Numer. Anal. , volume =

    work page

  49. [49]

    2002 , author =

    Exponential Time Differencing for Stiff Systems , journal =. 2002 , author =

    work page 2002

  50. [50]

    1998 , author =

    A New Class of Time Discretization Schemes for the Solution of Nonlinear PDEs , journal =. 1998 , author =

    work page 1998

  51. [51]

    2001 , publisher =

    Andreas Prohl , title =. 2001 , publisher =

    work page 2001

  52. [52]

    Implicit-Explicit

    Lorenzo Pareschi and Giovanni Russo , year =. Implicit-Explicit. Recent Trends in Numerical Analysis , series =

    work page

  53. [53]

    Highly stable implicit–explicit. Appl. Numer. Math. , volume =. 2017 , issn =

    work page 2017

  54. [54]

    Verwer , title =

    Willem Hundsdorfer and Jan G. Verwer , title =. 2003 , series =

    work page 2003

  55. [55]

    1996 , series =

    Hairer, Ernst and Wanner, Gerhard , title =. 1996 , series =

    work page 1996

  56. [56]

    Stability and error analysis of fully discrete original energy-dissipative and length-preserving scheme for the

    Li, Binghong and Li, Xiaoli and Wang, Cheng and Yang, Jiang , journal=. Stability and error analysis of fully discrete original energy-dissipative and length-preserving scheme for the

    work page

  57. [57]

    Atsushi Fuwa, Tetsuya Ishiwata, Masayoshi Tsutsumi , title =. Jpn. J. Ind. Appl. Math. , volume =

    work page