Pith. sign in

REVIEW 3 major objections 6 minor 18 references

New stable, explicit, first order method to solve the heat conduction equation

T0 review · 3 major / 6 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read A new explicit 'constant-neighbour' update solves the spatially discretized heat equation with unconditional stability and first-order accuracy in time, avoiding matrix solves.

desk verdict The 'constant-neighbour' method is exactly exponential Euler/ETD1 with a diagonal splitting, so the central novelty claim doesn't hold, though the stability and error analysis are correct for the stated assumptions. read the letter →

arxiv 1908.09500 v1 pith:TI46NSI5 submitted 2019-08-26 physics.comp-ph physics.flu-dyn

Endre Kovács This is my paper · ORCID
classification physics.comp-phphysics.flu-dyn MSC 65M0665M1265L20 PACS 44.10.+i
keywords heatconductionexplicitmethodsunconditionallystablestiffequationsconstant-neighbourmethodfirstorderconvergenceequationirregulargrids
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

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

The reading

The paper claims that a new explicit update rule, called the 'constant-neighbour' method, solves the spatially discretized heat equation with unconditional stability and first-order accuracy in time. During a timestep the method holds every cell's neighbours fixed, solves the resulting one-variable linear differential equation exactly, and advances each temperature with an exponential weight. Because the update needs only local sums of neighbour temperatures and no matrix factorization, it can take timesteps far beyond the explicit Euler limit on stiff heat-conduction meshes. The paper supports the claim with a two-cell analytic comparison, a proof in the appendix, and timing tests on a severely stiff 5000-cell random system. A sympathetic reader should care because this combination would remove the usual tradeoff between small timesteps and matrix solves.

What carries the argument

The load-bearing object is the constant-neighbour update: freeze each cell's neighbours at the start of the step and solve the resulting scalar linear ODE exactly. Its work is to replace the coupled stiff system by independent exponential steps. For cell $i$, the step is governed by the characteristic time $\tau_i = C_i / \sum_j U_{ij}$; the old temperature is damped by $e^{-h/\tau_i}$, while the neighbour-and-source average enters with weight $1 - e^{-h/\tau_i}$. Because these two weights are nonnegative and sum to one, the update obeys a discrete maximum principle for the homogeneous problem, keeping the new temperature inside the old range of cell temperatures when sources are absent. The same expansion that gives the two-cell exact comparison supplies the first-order convergence proof in the appendix: the constant-neighbour step matches the exact evolution operator through terms of order $h$, with a local remainder of order $h^2$.

What would settle it

Take a one-dimensional uniform heat-conduction problem with a smooth initial profile and an exact Fourier solution, and run the constant-neighbour update with a timestep far above the explicit Euler stability limit, for example $h = 100 \cdot 2 C / \sum U$; if the numerical temperatures grow, oscillate, or leave the predicted range, the unconditional-stability claim is false. The same experiment with a sequence of step sizes settles the order claim: the global error should fall linearly with $h$ if the method is truly first order.

Watch

Extended reading notes

Core claim

The central discovery is the update formula (3), derived from the spatially discretized heat equation \[ C_i \frac{dT_i}{dt} = \sum_{j \in \mathrm{neigh}(i)} U_{ij}(T_j - T_i) + Q_i. \] If each neighbour temperature $T_j$ is held constant during the step, the remaining equation for $T_i$ is linear and can be solved exactly, giving \[ T_i(t+h) = T_i(t) $e^{{-h/\tau_i}}$ + \left(1 - $e^{{-h/\tau_i}}$\right) \frac{\sum_j U_{ij} T_j(t) + Q_i}{\sum_j U_{ij}}, \quad \tau_i = \frac{C_i}{\sum_j U_{ij}}. \] The coefficients multiplying the temperatures are nonnegative and sum to one, so the homogeneous update is a convex average of old cell values; the paper argues this makes the method unconditionally stable for heat conduction. Comparing the Taylor expansion of this frozen-coefficient step with the exact matrix-exponential solution of the system shows the local error is second order in $h$, hence first-order convergence. The paper presents exact two-cell and stiff 5000-cell tests as numerical support.

Load-bearing premise

The load-bearing premise is that all couplings between cells are nonnegative, all heat capacities are positive, and the conductances and source terms are held constant during each timestep; if any of these fails, the convex-weight stability argument and the first-order error proof no longer apply.

Editorial extensions

If this is right

  • Timesteps can be chosen from accuracy needs alone, because the explicit update has no stability restriction for heat conduction problems.
  • Each step is a single local pass over cells and their neighbour lists, so the method is easy to parallelize and cheap on large meshes.
  • The scheme applies to any space dimension and to irregular or inhomogeneous grids, since only per-cell capacities and per-edge conductances are needed.
  • For quick approximate solutions of very large stiff problems, the reported 5000-cell test shows the constant-neighbour method can be far faster than standard implicit solvers, at the price of first-order accuracy.
  • The first-order error scale means the method suits exploratory or warm-start computations; the paper notes a second-order version is being developed.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Although the paper does not develop the connection, the constant-neighbour step is a form of local linearization for exponential integration; viewing it that way suggests the stability analysis could be extended to stiff reaction-diffusion systems whose frozen linearization still has nonnegative off-diagonal couplings.
  • If thermal conductivity depends on temperature, freezing the conductances at the start of each step would preserve the convex-weight stability argument, but the error analysis would need an extra term for the changing coefficients; this is a testable extension the paper does not make.
  • For problems with negative effective couplings, such as a linearized reaction term that amplifies temperature deviations, the unconditional-stability proof does not apply; one could derive a step-size restriction from the requirement that the update weights stay nonnegative, turning the method into a conditionally stable explicit scheme for those cases.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 6 minor

Summary. The manuscript proposes an explicit one-step scheme, Eq. (3), for the spatially discretized heat equation (1). The update is derived by freezing the neighbor temperatures and source terms within each timestep, solving the resulting uncoupled scalar linear ODEs exactly, and then updating all cells in parallel. The authors claim the method is new, unconditionally stable, first-order in time, and faster than standard explicit or implicit solvers on stiff heat-conduction problems. They verify first-order convergence on a two-cell problem with an exact solution and compare errors and runtimes on a random 100x50 stiff lattice against Matlab's ode15s/ode23s and a Crank-Nicolson-type implementation. The Appendix gives a local truncation error calculation showing O(h^2) local error, which the authors use to infer first-order global convergence.

Significance. If the novelty claim were supported, an explicit unconditionally stable first-order method for heat conduction would be practically useful for large stiff systems. The derivation is parameter-free, and the stability intuition via convex combinations is clear and accessible. However, the central claimed novelty is not supported: Eq. (3) is algebraically identical to the exponential Euler/ETD1 method applied to a diagonal splitting of the semi-discrete operator, a standard method in the exponential-integrator literature that the paper itself cites as [3]. The numerical experiments provide some evidence of practical performance, but they do not compare against ETD1. After correcting the novelty claim, the contribution would be a compact pedagogical derivation and a performance study, not a new algorithm.

major comments (3)
  1. [Section 2, Eq. (3)] The update formula Eq. (3) is algebraically identical to the exponential Euler/ETD1 method applied to the diagonal splitting of Eq. (1). Let D=diag(-1/tau_i) be the diagonal part of M, and R=M-D the off-diagonal part plus the source term. The standard ETD1 update T^{n+1}=e^{hD}T^n + h*phi_1(hD)(R T^n + Q), with phi_1(z)=(e^z-1)/z, reduces componentwise exactly to Eq. (3) after substituting tau_i=C_i/sum_j U_ij. Since ETD1 is a well-known method and is cited in the manuscript as [3], the abstract's claim that this is a 'novel explicit and stable numerical algorithm' and the introduction's claim of 'fundamentally new explicit methods' are not supported. The stability and first-order analyses in Section 2 and the Appendix are properties of this known scheme, not of a new method. The authors must either identify a genuinely new algorithmic component or reframe the paper as an application/exposition of ETD1 to heat conduction.
  2. [Section 2 and Appendix] The stability and order claims are proved only under restrictive hypotheses that are not stated in the abstract or Section 5. The convex-combination stability argument in Section 2 requires nonnegative conductances U_ij and positive capacities C_i, and the Appendix's local error calculation assumes that M and Q are constant over the timestep. For temperature-dependent conductivities or negative coupling terms (e.g., linearized convection or sources with negative slope), the convex-weight argument fails and the frozen-coefficient local error analysis does not apply. The authors should explicitly restrict the claims to the linear, constant-coefficient, nonnegative-coupling case, or supply proofs for the broader claims made in the abstract. Additionally, the Appendix proves only local truncation error O(h^2); the global first-order claim also requires a zero-stability/consistency argument, which is not stated.
  3. [Section 4, Table 1] The numerical comparison does not include an ETD1/exponential-Euler baseline, even though Eq. (3) is exactly that method. The abstract claims the method can give approximate results for extremely large systems faster than commonly used explicit or implicit methods, but the experiments compare only against Matlab's ode15s/ode23s and a Crank-Nicolson-type implementation. Since ETD1 is the relevant reference method, the authors should compare against it directly and state whether their implementation differs in any way. Otherwise the performance claim is not demonstrated for the comparison class that would establish the method's advantage. There are also inconsistencies in the reported runtimes: the text states ode15s needed 677s, while Table 1 lists ode15s runtimes of 181, 248, and 400s.
minor comments (6)
  1. [Section 3 and Section 4] Figure 2 is used twice: once in Section 3 for error versus stepsize and once in Section 4 for a spatial temperature snapshot. The figures should be renumbered.
  2. [Table 1] The notation 'CN, 42 10−⋅' is unclear; the authors should explicitly write the timestep sizes, e.g., h=4e-2, in the table or caption.
  3. [Section 4] The random parameters are generated with MATLAB's rand without a fixed seed, so the numerical experiment is not exactly reproducible. Providing a seed or code would strengthen the reproducibility of the reported results.
  4. [Section 4, Fig. 3] The statement that 'the infamous phenomenon of order reduction has never been observed' is too broad, as it is based on only two test problems; this claim should be qualified or supported by additional experiments.
  5. [Section 3] The exact solution formulas are badly garbled by typesetting; they should be re-set cleanly so the reader can verify the comparison.
  6. [Section 2] The sentence 'the result is always bounded' is not strictly true in the presence of a nonzero constant source, since temperatures grow linearly in time; the intended statement is about the homogeneous part being non-increasing in norm.

Circularity Check

0 steps flagged · score 0.0 of 10

Derivation is self-contained: Eq. (3) is obtained by solving a stated frozen-coefficient ODE, and the stability/order analyses compare against the exact ODE solution rather than importing the conclusion.

full rationale

The paper derives its update (3) directly from the spatially discretized heat equation (1) by freezing all neighbour temperatures and source/conductance coefficients during the step, solving the resulting scalar linear ODE (2) exactly, and then writing the solution as an explicit one-step formula. No fitted parameters enter this construction: tau_i is defined from C_i and the U_ij conductances of the discretized problem, and the update coefficients are explicit functions of h/tau_i and the U_ij/C_i ratios. The stability argument is an internal convex-combination property of the coefficients in (3), and the first-order claim is checked in the Appendix by expanding both the exact matrix-exponential solution of (1) and the proposed formula in powers of h and showing that the difference is O(h^2); this is a direct local-truncation-error computation, not an appeal to a prior conclusion. The single self-citation [5] only states that an earlier no-source version was already published; it is not used to justify stability, order, or novelty. Even if Eq. (3) were algebraically identical to exponential Euler / ETD1 with a diagonal splitting, that would be a novelty-overlap or correctness concern about the word 'new', not a circularity: the derivation does not assume ETD1 as an input, and nothing in the paper's analysis is equivalent by construction to its conclusions. No step reduces the claimed result to its own inputs.

Assumptions & free parameters 0 free parameters · 5 assumptions · 0 invented entities

No fitted constants are used: the timestep h is a user choice and the random test distributions are benchmark inputs. The assumptions are standard linear heat-diffusion modeling, power-series linear algebra, the frozen-coefficient step, and trust in the MATLAB reference solution. No new physical entities are introduced.

assumptions (5)
  • domain assumption The spatially discretized heat equation is a linear ODE dT/dt = M T + Q whose conductance matrix has C_i > 0 and U_ij >= 0, making M an M-matrix with negative diagonal and nonnegative off-diagonals.
    Used in Section 2 to define the ODE system and to guarantee the convex-combination weights in Eq. (3).
  • standard math The exact solution of the ODE system (1) can be written with a matrix exponential and expanded in a power series.
    Appendix Section 8 expands e^{Mt} T0 plus the source term to compare local errors.
  • domain assumption Heat sources and conductance coefficients are constant over each timestep h.
    The frozen-neighbour derivation in Section 2 and the first-order proof in the Appendix assume coefficients do not change during the step.
  • domain assumption MATLAB ode15s with strict tolerances (RelTol=AbsTol=1e-8) provides the reference solution used for all error metrics.
    Used in Section 4 and Table 1 as the ground truth for MaxD, SumD, and SumEnD.
  • domain assumption The method does not suffer from order reduction on the stiff benchmark, so the global error remains O(h).
    The paper states in Section 4 that order reduction 'has never been observed' but supplies no proof, so the uniform convergence is an empirical assumption.

how reviews work

0 comments
Cite this review

Pith. "Pith review of New stable, explicit, first order method to solve the heat conduction equation." pith.science (2026). https://pith.science/paper/TI46NSI5

@misc{pith2026190809500,
  author       = {Pith},
  title        = {Pith review of: New stable, explicit, first order method to solve the heat conduction equation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/TI46NSI5}},
  note         = {Machine review of arXiv:1908.09500}
}
read the original abstract

We introduce a novel explicit and stable numerical algorithm to solve the spatially discretized heat or diffusion equation. We compare the performance of the new method with analytical and numerical solutions. We show that the method is first order in time and can give approximate results for extremely large systems faster than the commonly used explicit or implicit methods.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

18 extracted references · 18 canonical work pages

  1. [3]

    constant- neighbour method

    We make a simplification: when we calculate the new value of a variable Ti, we neglect that other variables are also changing during the timestep. It means that we consider Tj a constant if ji≠ , thus we can call it “constant- neighbour method”. So we have to solve uncoupled, linear ODEs: i i ii i'T a mT= + (2) where 0,i ij j ijia mT Q≠= +∑ and ii i 1/m τ...

    work page

  2. [1]

    The PDEs for real -life problems can rarely be solved analytically

    Introduction and the studied problem It is well known that the simplest Fourier-type heat conduction phenomena are described by the heat equation, which is a second-order parabolic partial differential equation (PDE), with the following form: T Tqt α∂ =∆+∂ , where /( )kcαρ= is the thermal diffusivity, q, k, c, and ρ is the volumetric intensity of heat sou...

    work page

  3. [2]

    The proposed method Now we introduce the core method to solve the ODE system dT MT Qdt = +   (1) through the following two steps:

    work page

  4. [4]

    The analytical solution of eq

    We solve the obtained equations analytically. The analytical solution of eq. (2) at the end of the timestep is the following: ii ii i i ii ij i 0,j i i i i ( ) () 1 () 1 ji ijj neigh hh hh Tt h Tt e a e U CTt e T Q eCU ττ ττ τ τ −− −− ≠ =  += ⋅ + ⋅− =     = ⋅+ + −   ∑ ∑ Thus we suggest the following simple formula to obtain the val...

    work page

  5. [5]

    It also implies that the process is easily parallelizable and even vectorizable

    It is obviously explicit, one can calculate the new values without solving a system of equations or even without using matrices. It also implies that the process is easily parallelizable and even vectorizable

    work page

  6. [6]

    It is unconditionally stable for the heat conduction equation

    work page

  7. [7]

    We performed numerical tests on several sys tems, but here we present only two different examples

    It can be easily applied regardless of the number of space dimension, grid irregularity and inhomogeneity of the heat conduction medium. We performed numerical tests on several sys tems, but here we present only two different examples

    work page

  8. [8]

    Verification: comparison with an exact result If we have only two variables and one heat source, eq. (e1) has the following analytical solution: 2 1 0,1 a 1 2 0,2 a ( ) exp 1 exp 1 exp ( ) exp 1 exp 1 exp Ctt tT t T St S C tt tT t T St S T T τττ τ τττ τ −− −  = − ++ −      −− −  = − +− −       + ...

    work page

Show all 18 references
  1. [9]

    constant-neighbour

    Comparison with numerical results for a stiff system The second system is a rectangle -shaped lattice, N x=100, Ny=50. A value (3 7 )10 rand−⋅ were given to the capacities iC , and to the inverse conductances xi i,1/ 1/ yUU (the resistances), where rand is a random number gene...

  2. [10]

    This method is explicit, unconditionally stable and first order in time

    Summary We presented a new numerical algorithm to solve the heat conduction or diffusion equation with external sources. This method is explicit, unconditionally stable and first order in time. We illustrated the performance of the method for a simple, analytically soluble cas...

  3. [11]

    Acknowledgement I would like to thank András Gilicz for introducing me to this way of thinking and for the long years of collaboration

  4. [12]

    Skvortsov, L.: Explicit stabilized Runge-Kutta Methods, Comput. Math. & Math. Phys., Vol. 51, No. 7, pp. 1153, (2011)

    2011

  5. [13]

    Özişik, M. N. et al.: Finite Difference Methods in Heat Transfer, ISBN 69781482243468, eBook, CRC Press, pp224. (2017)

    2017

  6. [14]

    Hochbruck, M

    M. Hochbruck, M. and Ostermann, A.: Exponential integrator s, Acta Numerica, 19, pp.209, (2010)

    2010

  7. [15]

    R., McKee, S

    Gourlay, A. R., McKee, S. The construction of hopscotch methods for parabolic and elliptic equations in two space dimensions with a mixed derivative, J. Comput. Appl. Math., 3, No. 3. pp. 201, (1977)

    1977

  8. [16]

    Kovács, E.; Gilicz, A.: New Stable Method t o Solve Heat Conduction Problems in Ex tremely Large Systems , Design of Machines and Structures, Vol. 8, No. 2., pp. 30. (2018)

    2018

  9. [17]

    et al.: Numerical Solution Of Ordinary Differential Equations, Wiley, Chapter 9.3

    Atkinson, K. et al.: Numerical Solution Of Ordinary Differential Equations, Wiley, Chapter 9.3. (2009)

    2009

  10. [18]

    Constant Neighbour

    Appendix: The proof that the method is first order in time Using the power series form of the exponential function, the exact solution of (1) is the following: ( ) 1 0 23 23 23 2 0 () 1 1 ... ...2 3! 2 3! Mt MtTt e T e M Q tt ttM tM M T tM M Q −= +− =   =+ + + + ++ + + ...

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.