{"id":"dac8df99-ac55-46a0-8b89-b123a20a2a19","arxiv_id":"2506.04188","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A kernel-compression plus linear-chain-trick reformulation lets standard stiff ODE/DAE codes such as Radau5 solve fractional differential equations with memoryless, variable-step, fast linear algebra.","lead":"This paper shows how to solve fractional differential equations by first replacing the memory-heavy fractional kernel with a sum of exponentials, turning the problem into a large but structured system of stiff ordinary differential equations that existing solvers such as Radau5 can handle.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Algorithm 1 does not deliver condition (5.9) as stated: (5.13) and (5.14) give at best Cε(1+t^α) with C>1, so Theorem 1's εu(t) bound is not established.","rationale":"Reader's verdict is already CONDITIONAL, and I agree with that. My analysis sharpens the weakest assumption. The theorem's condition (5.9) is not met by Algorithm 1 because of two concrete quantitative mismatches: the 3/Γ(α+1) factor on δ≤t≤T and the high-tail behaviour near 0. This is an internal inconsistency, not a difference from consensus, and it is directly load-bearing because Theorem 1 is the only theoretical justification for the kernel-replacement step. The rest of the paper—linear chain trick, Jacobian structure, fast linear algebra—is standard or clearly explained. The numerical experiments in Tables 7.1–8.2 are consistent with an error proportional to max(ε, Tol), so I do not doubt the practical soundness. The missing comparison with [25] and absent code hash are real but secondary; they do not affect the correctness of the theoretical claim. Therefore no verdict change is needed; the condition for acceptance should include either proving (5.9) with the constants from Algorithm 1 or restating the theorem with the appropriate constant.","tokens_in":20441,"tokens_out":13551,"duration_ms":125337,"concrete_test":"Implement Algorithm 1 for α=1/2, ε=10^{-5}, T=1 (or a Table 7.1 row) and form k̃(s)=T_N^M(s,h) from (5.6). Use adaptive quadrature on a fine t-grid to compute R(t)=[∫_0^t |s^{α−1}/Γ(α)−k̃(s)| ds]/[ε(1+t^α)] and report max R(t) over [0,T]. If max R≤1, (5.9) holds and the concern is resolved. If max R>1, rerun the Table 7.1 experiment with ε_alg=ε/5 in Algorithm 1; if the observed error scales with ε_alg rather than with the original ε, the theorem needs only the constant correction and the conditional verdict stands.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Conditional on (5.9), Theorem 1 is correct and the reduction in §3 is sound. The gap is that Algorithm 1 is claimed to produce a kernel approximation satisfying (5.9), but the bounds displayed in §5.3 do not show this. For δ≤t≤T, (5.13) gives ∫_δ^t |k−k̃| ds ≤ 3ε t^α/Γ(α+1); since Γ(α+1)<1 for 0<α<1, this exceeds ε(1+t^α) for large t. On (0,δ), the high-tail term in (5.5) is not bounded by εk(s): the factor Γ(1−α,se^{Nh})/Γ(1−α) tends to 1 as s→0, so only its integral is O(ε), via (5.14). Combining the two regions yields a bound Cε(1+t^α) with C>1, not ε(1+t^α). Consequently, running Algorithm 1 with the user's ε does not, as stated, put the hypotheses of Theorem 1 in force; the rigorous error statement needs either a smaller ε_alg=ε/C in Algorithm 1 or a constant C in Theorem 1. This is a repairable correctness gap, but it is the load-bearing step between the construction and the accuracy claim.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper proposes to solve fractional ODEs (Caputo or Riemann-Liouville) by rewriting them as Volterra integral equations, approximating the kernel t^{alpha-1}/Gamma(alpha) by a sum of exponentials via the Beylkin-Monzon quadrature, and converting the convolution integral into an augmented system of ODEs. The augmented system is stiff and is integrated with Radau5, using a purpose-built structured linear solver that reduces the linear-algebra cost from O((d+D)^3) to O(d^3+D). The paper gives explicit formulas and an algorithm for choosing the quadrature parameters, a theorem showing that an O(epsilon) weighted L1 kernel error implies an O(epsilon) solution error, and numerical experiments with known exact solutions, including scalar tests, a Brusselator system, multi-term problems, and 1D fractional PDEs.","tokens_in":20617,"tokens_out":13794,"duration_ms":129869,"significance":"If the claimed accuracy and efficiency hold, the paper offers a practically useful, memory-efficient alternative to quadrature-based fractional ODE solvers, with the notable advantage of variable step sizes inherited from Radau5. The kernel-compression construction is explicit and depends only on the accuracy parameter epsilon and a quadrature parameter; Theorem 1 provides a clean propagation of the kernel approximation error into the solution; and the structured linear solver is a genuine algorithmic contribution, demonstrated by large speedups in the experiments. The driver examples are publicly available, which aids reproducibility. The main theoretical caveat, namely the constant mismatch between Algorithm 1 and the hypothesis of Theorem 1, does not affect the substance of the method and is fixable.","major_comments":[{"comment":"The construction as stated does not put the hypothesis (5.9) of Theorem 1 in force for the user-supplied epsilon. From (5.6), the bound on [delta,T] is int_delta^t |k(s)-k_tilde(s)| ds <= 3 epsilon t^alpha / Gamma(alpha+1), and since Gamma(alpha+1) < 1 for 0<alpha<1, the coefficient 3/Gamma(alpha+1) is larger than 3, so for all sufficiently large t this exceeds epsilon(1+t^alpha). On (0,delta), equation (5.14) only controls the integral of k(s), not the pointwise high-tail contribution in (5.5), which adds a further O(epsilon) term. Combining the two regions gives at best C epsilon (1+t^alpha) with C>1, not epsilon(1+t^alpha). Consequently Algorithm 1 with the user's epsilon does not, as stated, satisfy the assumption of Theorem 1. This is the load-bearing link between the construction and the error bound, but it is repairable either by scaling epsilon in Algorithm 1 by 1/C or by stating Theorem 1 with an explicit constant on the right-hand side.","section":"5.3, Eqs. (5.13)-(5.14), Algorithm 1"}],"minor_comments":[{"comment":"The factor e^{alpha i h} in the displayed sums appears inconsistent with the integrand e^{(1-alpha)s} in (5.1)-(5.2); the incomplete-gamma substitution suggests that it should be e^{(1-alpha) i h}. Please correct and re-verify (5.6).","section":"Eq. (5.5)"},{"comment":"The sentence stating that Algorithm 1 chooses parameters so that the kernel approximation satisfies (5.9) is stronger than what (5.13)-(5.14) show; please rephrase it to say that the weighted L1 error is O(epsilon), or that a version of (5.9) with an explicit constant is satisfied.","section":"5.3, before Algorithm 1"},{"comment":"No error estimate is given for the alpha>1 constructions (6.3)-(6.4); a sentence explaining that the same weighted-L1 argument extends to the monomial-times-exponential kernel, or explicitly noting that the proof is not carried out, would clarify the scope of the theoretical guarantee.","section":"Section 6"},{"comment":"The reported CPU times for the present method at T=220 differ (0.051 s in Table 7.4 and 0.8 s in the text) and correspond to different tolerances and accuracies; please state the settings used for the 0.8 s figure.","section":"7.2, Table 7.4 and text"},{"comment":"There are minor typos, including 'frcational' in Section 2.2 and 'developped' in the Introduction, that should be corrected.","section":"2.2 and Introduction"},{"comment":"The role of the two dummy variables added per integral block in dc_sumexp.f is described operationally but not explained; a sentence on why the Radau5 interface requires them would improve reproducibility.","section":"4.3"}],"recommendation":"major_revision","confidential_remarks":null},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear colleague,\n\nHere is my read on Guglielmi–Hairer. The paper delivers something real: a complete, memoryless, variable-step recipe for fractional ODEs by kernel compression with Beylkin–Monzón exponential sums, the linear chain trick, and a purpose-built linear solver that lets Radau5 run on the augmented stiff system. The extension to alpha > 1 and the banded treatment for 1D fractional PDEs are useful too. The benchmarks against Garrappa's codes show dramatic speedups, and the test problems have exact solutions, so the demonstrations are credible.\n\nThe central error analysis (Theorem 1) is clean and parameter-free, and the kernel approximation bounds in Section 5.1 are standard. But there is a genuine soft spot in the bridge between the algorithm and the theorem. Algorithm 1 claims to choose parameters so that the kernel error satisfies condition (5.9). The displayed bounds do not deliver that. From (5.13), the integral over [delta, T] is at most 3ε t^α / Γ(α+1), and since Γ(α+1) < 1 for 0 < α < 1, this exceeds ε t^α. Adding the (0, delta) contribution only gives C ε (1 + t^α) with C > 1. So running Algorithm 1 with the user's epsilon does not put the hypotheses of Theorem 1 in force as stated. The fix is straightforward: either shrink epsilon in Algorithm 1 by a known constant or state Theorem 1 with a slightly larger constant. This is a repairable correctness gap, but it is load-bearing between the construction and the claimed accuracy bound. The authors even say (5.13) is \"up to a constant compatible with (5.9),\" so they know; they just didn't propagate the constant.\n\nOther issues are smaller. The drivers are said to be available from the authors' homepages, but no repository or hash is given, which makes reproduction harder than it should be. The paper also claims to know no variable-step code for fractional problems, yet cites [25] (López-Fernández–Lubich–Schädle), an adaptive fast oblivious convolution method; a sentence justifying why that does not count would be enough. The coupling between the Radau5 tolerance and the kernel error is empirical, but the experiment in Table 7.1 addresses it directly, so I do not see that as a flaw.\n\nOverall this is a solid numerical-methods contribution. The core reduction is sound and the numerical evidence is convincing. I would send it to a serious referee, with the expectation that the constant factor in the error bound get fixed and the code location be made explicit. The paper is worth reading for anyone working on fractional ODE solvers, and I would cite it after the fix.\n\nBest,\n[Your name]","headline":"A genuinely useful memoryless fractional-ODE solver built on known pieces, but the paper's rigorous error claim has a constant-factor gap between Algorithm 1 and Theorem 1.","tokens_in":21246,"tokens_out":1801,"would_cite":true,"duration_ms":17480,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["26A33","34A08","65L06","45D05","65F05"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that fractional differential equations can be solved as stiff augmented ODE systems after replacing the fractional kernel $t^{\\alpha -1}/\\Gamma(\\alpha)$ by an explicit sum of exponentials, with an…","keywords":["fractional differential equations","Caputo derivative","Riemann-Liouville derivative","Volterra integral equations","exponential sum approximation","stiff ODE solvers","memoryless convolution","structured linear systems"],"falsifier":"Take any $\\alpha\\in(0,1)$ and the parameters from Algorithm 1, compute the left side of (5.9) numerically on a mesh that extends to $t=0$, and check whether it stays below $\\varepsilon(1+t^\\alpha)$; if it exceeds that bound, the proof of Theorem 1 no longer covers the method. A direct experiment would run the scalar test problem (7.1) with $\\varepsilon$ set well below the integrator tolerance and compare the observed error against $\\varepsilon\\,u(T)$ from (5.11); an observed error much larger than predicted would falsify the claimed $O(\\varepsilon)$ accuracy.","tokens_in":20129,"feed_emoji":"🧮","tokens_out":10653,"duration_ms":90257,"temperature":0.7,"pith_summary":"Fractional derivatives are nonlocal, so standard stiff ODE codes cannot be applied directly to fractional differential equations; this paper turns that obstacle into an augmented but standard ODE system. The key reduction replaces the fractional kernel $t^{\\alpha -1}/\\Gamma(\\alpha)$ by an explicit finite sum of exponentials, so each convolution integral becomes a few scalar ODEs with no history to store. The paper proves (Theorem 1) that if the kernel approximation satisfies a weighted $\\ell^1$-type error bound, the resulting solution error is $O(\\varepsilon)$, and it gives Algorithm 1 to choose the exponential-sum parameters so that bound holds. It then shows how the stiff code Radau5 can solve the augmented system, with a linear solver that exploits the arrow-shaped Jacobian to run in $O(d^3 + D)$ time instead of $O((d+D)^3)$. If the method works as claimed, fractional initial-value problems become solvable with variable steps, constant memory, and the reliability of mature stiff integrators.","feed_headline":"Exponential sums erase the memory in fractional differential equations","feed_subtitle":"Replace the fractional kernel by exponentials, and a standard stiff integrator handles the rest — no history to store.","key_machinery":"The linear chain trick is the object that carries the argument: replacing the convolution kernel by exponentials converts the fractional integral into the output of ODEs $\\dot z_i=-\\gamma_i z_i+G(t,y)$, so the memory lives in the current state rather than in stored history. The trapezoidal-rule exponential-sum approximation, applied after the change of variables $z=e^s$ to the integral representation (5.1), produces explicit $c_i$ and $\\gamma_i$ with the exponentially convergent bound (5.3); Algorithm 1 picks $\\delta,h,M,N$ to make the weighted $\\ell^1$ kernel error in (5.9) small, and Theorem 1 is the bridge from kernel error to solution error. The arrow-shaped Jacobian (3.6), with diagonal blocks $\\mathrm{diag}(-\\gamma_i)$ and rank-one off-diagonal blocks, makes the linear algebra fast: eliminating the $z$-variables costs $O(D)$ per step and leaves only a $d\\times d$ system (4.1), giving total $O(d^3+D)$ complexity.","core_discovery":"The central claim is that the nonlocal memory of a fractional differential equation is an artifact of the kernel representation: once $k(t)=t^{\\alpha -1}/\\Gamma(\\alpha)$ is replaced by $\\sum_{i=1}^n c_i e^{-\\gamma_i t}$, the fractional integral $J^\\alpha f(t)$ becomes a weighted sum of variables $z_i(t)$ that satisfy $\\dot z_i = -\\gamma_i z_i + G(t,y)$, and the original problem (2.8) becomes the stiff ODE system (3.4)-(3.5) with the arrow-shaped Jacobian (3.6). For $0<\\alpha<1$ the coefficients $c_i,\\gamma_i$ come from a trapezoidal-rule discretization of the real integral representation (5.1), with parameters $\\delta,h,M,N$ chosen by Algorithm 1 so that the kernel error satisfies (5.9); Theorem 1 converts that condition into the solution bound $\\|y-\\tilde y\\| \\le \\varepsilon\\, u(t)$, where $u$ solves the scalar Abel-type equation (5.11). For $\\alpha>1$ the kernel is split as $t^{m-1}/((\\alpha-1)\\cdots(\\alpha-m+1))$ times $t^{\\alpha_0-1}/\\Gamma(\\alpha_0)$, and the same chain trick yields a bi-diagonal block structure, so the whole range $\\alpha>0$ is covered. The paper supports the claim with experiments on a scalar problem with known exact solution, a fractional Brusselator, a multi-term fractional equation, and one-dimensional fractional PDEs, including comparisons with existing constant-step and memory-intensive codes.","pith_inferences":["The same kernel-replacement machinery should transfer to other weakly singular convolution kernels, such as logarithmic kernels or Abel kernels, provided a weighted $\\ell^1$ error bound analogous to (5.9) can be verified; Theorem 1 itself is not specific to the fractional power kernel.","A practical risk not settled by the paper is the behavior of the approximation on the singular strip $0<t<\\delta$; a natural test would be to replace the analytic $\\delta$ bound by an explicit exponential fit on $(0,\\delta)$ and see whether the observed global error still tracks $\\varepsilon$.","Because the augmented system is stiff regardless of the original problem, the method's value hinges on the structured linear solver; for very large $D$ the $O(D)$ elimination is the enabling step, so improving that block elimination would extend the method to large systems with many integral terms.","The demonstrated CPU gains over step-size-fixed fractional codes suggest that variable-step fractional solvers are the main practical payoff; a systematic benchmark across tolerances and $\\alpha$ values would show where the memoryless approach overtakes quadrature-based splitting methods."],"forward_implications":["Any stiff or differential-algebraic ODE solver can be pointed at a fractional problem by feeding it the augmented system (3.4), so fractional ODEs inherit variable step-size control and tolerance-based adaptivity.","Memory use is bounded and independent of the integration horizon, since no past solution values are stored; the paper demonstrates runs to $T=5000$ for a multi-term equation and $T=1000$ for PDE examples at constant memory.","For semi-discretized one-dimensional fractional PDEs, the banded structure of the spatial Jacobian is preserved in the Schur complement (Lemma 8.1), so the per-step cost scales linearly with the number of spatial points.","The method detects stability changes in fractional problems: the multi-term example reproduces the instability onset near $\\alpha\\approx 0.654298$, so the exponential-sum approximation does not mask the underlying dynamics."],"supporting_citations":[{"why":"Supplies the trapezoidal-rule exponential-sum coefficients $c_i,\\gamma_i$ used to approximate the fractional kernel.","marker":"[6]"},{"why":"Establishes the linear chain trick and the $O(d^3+D)$ structured linear algebra that this paper extends to fractional kernels.","marker":"[15]"},{"why":"Provides the Radau5 stiff/differential-algebraic integrator that the proposed method is built on and adapted to.","marker":"[17]"},{"why":"Gives the exponentially convergent trapezoidal-rule error estimate used in (5.3) to bound the kernel approximation.","marker":"[33]"},{"why":"Provides the scalar fractional test problem with known exact solution used in the accuracy and timing experiments.","marker":"[12]"},{"why":"Supplies the fractional Brusselator test problem and comparison software used for the CPU and memory experiments.","marker":"[13]"},{"why":"Provides the multi-term benchmark problem whose exact solution and stability threshold are used to validate the method.","marker":"[34]"}],"fun_headline_variants":["Exponential kernel sums make fractional ODEs memoryless","No history: exponentials replace fractional kernel","Fractional equations solved by local exponential sums","Kernel to exponentials: memory becomes ODEs"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The accuracy guarantee assumes the exponential-sum kernel error stays small all the way down to time zero, but the proof only checks it from a small cutoff $\\delta$ onward and trusts that the singular interval before $\\delta$ contributes a negligible error.","fun_headline_variants_meta":{"raw":{"variants":["Exponential kernel sums make fractional ODEs memoryless","No history: exponentials replace fractional kernel","Fractional equations solved by local exponential sums","Kernel to exponentials: memory becomes ODEs"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000304,"raw_usage":{"total_tokens":1870,"prompt_tokens":1192,"completion_tokens":678,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":808,"completion_tokens_details":{"reasoning_tokens":617}},"tokens_in":808,"tokens_out":678,"duration_ms":6059,"temperature":1.0,"reasoning_tokens":617,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T10:45:23.213220+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take any $\\alpha\\in(0,1)$ and the parameters from Algorithm 1, compute the left side of (5.9) numerically on a mesh that extends to $t=0$, and check whether it stays below $\\varepsilon(1+t^\\alpha)$; if it exceeds that bound, the proof of Theorem 1 no longer covers the method. A direct experiment would run the scalar test problem (7.1) with $\\varepsilon$ set well below the integrator tolerance and compare the observed error against $\\varepsilon\\,u(T)$ from (5.11); an observed error much larger than predicted would falsify the claimed $O(\\varepsilon)$ accuracy.","supporting_citations":[{"cited_title":"Beylkin and L","cited_arxiv_id":null,"evidence_quote":"Supplies the trapezoidal-rule exponential-sum coefficients $c_i,\\gamma_i$ used to approximate the fractional kernel."},{"cited_title":"Guglielmi and E","cited_arxiv_id":null,"evidence_quote":"Establishes the linear chain trick and the $O(d^3+D)$ structured linear algebra that this paper extends to fractional kernels."},{"cited_title":"Hairer and G","cited_arxiv_id":null,"evidence_quote":"Provides the Radau5 stiff/differential-algebraic integrator that the proposed method is built on and adapted to."},{"cited_title":"Trefethen and J.A.C","cited_arxiv_id":null,"evidence_quote":"Gives the exponentially convergent trapezoidal-rule error estimate used in (5.3) to bound the kernel approximation."},{"cited_title":"Diethelm, N","cited_arxiv_id":null,"evidence_quote":"Provides the scalar fractional test problem with known exact solution used in the accuracy and timing experiments."},{"cited_title":"Garrappa","cited_arxiv_id":null,"evidence_quote":"Supplies the fractional Brusselator test problem and comparison software used for the CPU and memory experiments."},{"cited_title":"Xue and L","cited_arxiv_id":null,"evidence_quote":"Provides the multi-term benchmark problem whose exact solution and stability threshold are used to validate the method."}],"review_version":1}