Central differences, space only

What it fills

• Replace a spatial domain with a finite grid

  Approximate the spatial derivatives of a PDE on finitely many points, leaving time continuous, so that what remains is a system of ordinary differential equations in the grid values. The method of lines: the continuum is gone, the clock is not.

When it applies

Linden, Montanaro and Shao derive the stencil from Taylor's theorem with remainder and require the solution to be four times differentiable in space for the bound to hold. Their analysis fixes the fourth spatial derivatives by assumption rather than reading them off the problem, and they note the constraint applies to the solution of the heat equation rather than to the initial condition — though a bound on the initial condition implies one at later times, because the discrete time-evolution operator cannot increase the infinity norm and commutes with the discretized partial-derivative operators.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  heat equation  du/dt = alpha * (d^2u/dx_1^2 + ... + d^2u/dx_d^2)  on the
       hypercubic region  x_i in [0, L],  i = 1..d, with periodic boundary
       conditions in each x_i -- but NOT in t                       (Sec. I A)
       initial condition u0, diffusivity alpha > 0, n grid points per
       dimension, so Delta x = L / n;  e_i = unit vector along x_i

requires  u four times differentiable in each x_i, on the region swept out
    by t in [0, T] -- the stencil error bound below only holds to that
    order                                                           (Sec. I A)
    # the paper states the 4th-derivative bound as a constraint on the
    #   SOLUTION u, not on the initial condition u0             (Sec. I A, p.7)
    # it adds that a bound on u0 implies one on u at later t. CAREFUL: the
    #   justification it gives is imported, not internal to this method --
    #   it appeals to the discretisation argument of Theorem 1, where a
    #   discrete time-evolution operator L cannot increase the infinity
    #   norm and commutes with the discretised partial-derivative operators.
    #   Theorem 1 is the FULLY discretized (forward-time centred-space)
    #   scheme and needs Delta t <= Delta x^2 / (2*d*alpha) for L to be
    #   stochastic at all                                    (Thm. 1, Eq. 14)
    # this construction takes no time step, so it builds no such L; the
    #   transfer of the bound is borrowed from that separate analysis

# --- stencil each second spatial derivative, one x_i at a time ------------
d^2u/dx_i^2(x)  ~=  ( u(x + Delta x * e_i) + u(x - Delta x * e_i)
                       - 2*u(x) ) / Delta x^2                          (Eq. 7)
    # error, from Taylor's theorem with remainder:
    #   |d^2u/dx_i^2(x) - stencil(x)| <= (Delta x^2 / 12) * sup_x |d^4u/dx_i^4(x)|
    #                                                                  (Eq. 9)
    # the time derivative du/dt is left EXACT here -- discretizing it too
    #   (Eq. 6, Eq. 8) is a separate step (forward-time stepping) this
    #   method does not take

# --- close the stencil into a matrix, using periodicity in space ----------
H = n x n matrix: -2 on the diagonal, 1 on the two neighbouring
    off-diagonals, and 1 in the two corners -- the wraparound entries the
    periodic boundary condition in x_i licenses                       (Eq. A3)
    # H carries only the NUMERATOR of the one-variable stencil above: its
    #   entries are -2, 1, 1 with no 1/Delta x^2. That factor is supplied
    #   separately, outside A, in the ODE below. So H*u~ by itself is not
    #   an approximation to d^2u/dx_i^2                        (Eq. A3 vs Eq. 7)

# --- tensor one copy of H into each of the d spatial dimensions -----------
A = sum_{j=1}^{d}  kron( I_n^{(j-1)}, H, I_n^{(d-j)} )                 (Eq. A2)
    # sparsity Theta(d); A is diagonalized by the d-fold tensor product of
    #   the quantum Fourier transform                                (Lemma 23)
    eigenvalues:  lambda_{j_1,...,j_d} = sum_{i=1}^{d} -4*sin^2(j_i*pi/n)
                                                                        (Eq. A4)
    ||(alpha/Delta x^2) * A|| = 4*alpha*d / Delta x^2
        # this is the norm the runtime of a downstream ODE solver scales
        #   with, not a claim proved about the stencil itself           (Eq. A8)
        # the paper writes this norm as ||A||, reusing the symbol of Eq. A2
        #   for the SCALED coefficient matrix; written out here to keep the
        #   two objects apart
        # Eq. A8 then rewrites it as Theta(alpha^2*d^2*zeta*T/epsilon) --
        #   that further step holds only under the Delta x chosen in
        #   Corollary 2, not for arbitrary grids                    (Cor. 2)

return  the semi-discretized system
    d u~/dt = (alpha / Delta x^2) * A * u~,   u~(0) = u0 restricted to the grid
                                                                       (Eq. A1)
# NOT part of this method: choosing a time-stepping scheme or ODE solver
#   for the returned system. Appendix A calls the result "a system of
#   ODEs" precisely because time is handed on untouched, to be discretized
#   or solved by whatever comes next                                  (App. A)

Cost, as the source states it

Linden, Montanaro and Shao's Eq. (9) bounds the second-derivative stencil itself: the difference from the exact derivative is at most $\frac{h^2}{12}\sup_x |d^4u/dx^4(x)|$ at spacing $h$. Applied to the heat equation and stated in their Appendix A, discretizing only the spatial variables gives the ODE system $d\tilde{u}/dt = (\alpha/\Delta x^2) A\tilde{u}$, whose generator is a sum of $d$ tensor factors of one circulant $n \times n$ matrix, has sparsity $\Theta(d)$, and carries eigenvalues $\lambda_j = -4\sin^2(j\pi/n)$ so that $\lVert A\rVert = 4\alpha d/\Delta x^2$.

Implementations

• `MethodOfLines.jl`'s centered-difference scheme, with the time variable left undiscretized

  MethodOfLines.jl is the SciML organisation's Julia package for turning a symbolically stated PDE into something a time-stepping solver can take, described by its own README as "a package for automated finite difference discretization of symbolically-defined PDEs in N dimensions". A `PDESystem` from ModelingToolkit.jl goes in — equations, boundary conditions, and `Interval` domains from DomainSets.jl — and a problem ready for `solve` comes out. Not an `ODEProblem`, though, and the distinction is the package's own: since the v1 change the time-dependent path returns a `DAEProblem`, and the `discretize` docstring warns that "Explicit Runge–Kutta methods such as `Tsit5()` solve `ODEProblem`s, not the `DAEProblem` returned by this method". What makes this the semi-discrete method rather than a full discretization is one argument: `MOLFiniteDifference` takes the grid spacings first and, second, a variable its docstring's Fields section calls "The independent variable left undiscretized, or `nothing` for a fully discretized system". Passing `t` grids space alone and hands time on; passing nothing discretizes everything and the docstring says the result is then a `NonlinearProblem`. The package is the successor to the same organisation's DiffEqOperators.jl — a comment in `src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl` calls that package "the previous home of this package" — and carries a `DerivativeOperator` struct, defined at line 1 of `src/discretization/derivative_operator.jl`, as its stencil container.

  `MOLFiniteDifference(dxs, time = nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align = CenterAlignedGrid(), kwargs...)` is the entry point, and the default `approx_order = 2` is what makes the second-derivative stencil the three-point one: `CompleteCenteredDifference(derivative_order, approximation_order, dx)` sets `stencil_length = derivative_order + approximation_order - 1 + (derivative_order + approximation_order) % 2`, which at `derivative_order = 2` and `approximation_order = 2` evaluates to $3$, and fills it with `(1 / dx^derivative_order) * calculate_weights(derivative_order, zero(T), dummy_x)` over the symmetric index range `dummy_x = (-div(stencil_length, 2)):div(stencil_length, 2)` — Fornberg's recurrence in `fornberg_calculate_weights.jl`, not a hard-coded $(1, -2, 1)$. Note that unlike the corpus record's $H$, which carries only the numerator, these coefficients already include the $1/\Delta x^{2}$. Two qualifiers decide whether a given term reaches this scheme at all. First, the centered ruleset is generated only for even orders — `generate_cartesian_rules` iterates `orders[iseven.(orders)]` and its docstring says "Any even ordered derivative may be adequately approximated by these" — so a first-order spatial derivative goes to `UpwindScheme()` by default, and there are separate schemes for the nonlinear and spherical Laplacians. Second, the operator is only the wraparound matrix of the record's periodic case when the boundary supplies one: the two near-boundary branches of `central_difference_weights_and_stencil` are guarded by `& !haslower` and `& !hasupper`, so a point near a boundary that has a lower or upper entry falls through to the interior branch, whose taps are `[bwrap(II + i * I1, bs, s, jx) for i in half_range(D.stencil_length)]`; without one, those rows instead take `low_boundary_coefs`/`high_boundary_coefs`, one-sided stencils of length `derivative_order + approximation_order`. `bwrap` is the wraparound itself, defined at `src/discretization/interface_boundary.jl` line 26 as a fold over `boundary_index`, whose periodic arm reaches `_wrapperiodic` at line 33 under the comment "shift l-1: u[1] ~ u[end]". `SciMLBase.discretize` then assembles the result, and since the v1 change its docstring reads "For a time-dependent system this builds a `DAEProblem`", because "MethodOfLines emits residuals of the form `D(u) - f ~ 0`, which are already implicit-DAE form". The `ODEProblem` is the fallback, and the docstring is specific about what does what: the systems that cannot be posed as a first-order DAE "fall back to `mtkcompile` plus an `ODEProblem`, which scalarizes the array equations", and separately "Supplying `analytic` selects the compiled `ODEProblem` path", which the code confirms at line 252 with `ode_path = analytic !== nothing`.

  No dataset; everything is symbolic. The inputs are a ModelingToolkit `PDESystem` and the grid spacings, and the package's own examples supply closed-form initial conditions rather than files — the heat tutorial poses `eq = Dt(u(t, x)) ~ Dxx(u(t, x))` with `bcs = [u(0, x) ~ cos(x), u(t, 0) ~ exp(-t), u(t, 1) ~ exp(-t) * cos(1)]`, both `t` and `x` given the domain `Interval(0.0, 1.0)` at `dx = 0.1`, and the 1-D diffusion tests use the method of manufactured solutions against $e^{-t}\cos x$. Grid spacings may be given as a number per variable or as an explicit non-uniform vector, and the non-uniform path has its own `CompleteCenteredDifference` method taking `x::AbstractVector` at line 77 of `centered_diff_weights.jl`. That method asserts nothing itself; the restriction sits one file over, in the `DX <: AbstractVector` method of `central_difference_weights_and_stencil` at `centered_difference.jl` line 37, whose `@assert length(bs) == 0` carries the message "Interface boundary conditions are not yet supported for nonuniform dx dimensions". That assertion, together with the matching one at `half_offset_centred_difference.jl` line 44 saying the same of periodic conditions, is what forbids the wraparound on a non-uniform grid.

  github.com/SciML/MethodOfLines.jl, `master` at `806d7a85` (2026-08-26), read 2026-08-27; Julia; MIT, whose `LICENSE` (21 lines) opens "MIT License" and "Copyright (c) 2022 SciML Open Source Scientific Machine Learning Organization". The files a reader opens are `src/interface/MOLFiniteDifference.jl` (97 lines; docstring lines 1-52, `struct MOLFiniteDifference{G}` at line 53, the constructor with `approx_order = 2` at lines 66-71, last touched by `82874237`, 2026-08-21, whose subject is "feat!: make array-form DAE discretization the v1 default" followed by the number of the pull request that landed it, pull request 650); `src/discretization/schemes/centered_difference/centered_diff_weights.jl` (153 lines; `CompleteCenteredDifference` for scalar `dx` at line 4, the vector-`x` method at line 77); `src/discretization/schemes/centered_difference/centered_difference.jl` (227 lines; `central_difference_weights_and_stencil` at lines 5 and 32, `central_difference` at lines 59 and 177, `generate_cartesian_rules` at lines 68 and 99; both this file and `centered_diff_weights.jl` were last touched by `613428f5`, 2026-01-04, a formatting-only commit); `src/discretization/schemes/fornberg_calculate_weights.jl` (86 lines; `calculate_weights` at line 20); `src/discretization/interface_boundary.jl` (`bwrap` at line 26, `_wrapperiodic` at line 33); and `src/MOL_discretization.jl` (290 lines; the `discretize` docstring at lines 212-243 and `function SciMLBase.discretize(` at line 244). The newest GitHub release and tag is `v0.11.13` (2026-06-25) while `Project.toml` line 3 on `master` reads `version = "1.1.1"`, so the v1 line described above is unreleased on GitHub at the time of reading; the Julia General registry, which is where a Julia release actually lands, was not checked.

  The package does publish accuracy studies, but none of them touches this scheme. `docs/src/tutorials/weno_showcase.md` makes the weighted relative $L^{2}$ error its stated figure of merit and reports a clustered grid at 129 points holding the layer to a relative error of $\sim 10^{-4}$, which it calls two orders of magnitude better than the uniform grid with twice the points; `docs/src/tutorials/nonuniform_weno.md` adds a manufactured-solution study with an estimated order of convergence, and reports the maximum error near the front dropping from $1.5 \times 10^{-2}$ on the uniform grid to $1.1 \times 10^{-3}$ on the clustered one; `benchmark/README.md` opens "MethodOfLines.jl WENO Benchmarks" and times rather than measures error. All of it is WENO on advection — the convergence study's equation is `Dt(u(t, x)) ~ -Dx(u(t, x))` under `advection_scheme = WENOScheme()` — and that tutorial's own limitations section says `WENOScheme` discretizes "first-order spatial derivatives only", while "Even-order terms such as diffusion are handled by the standard centered schemes". So nothing published measures the centered second-derivative scheme, and its only numbers are in the repository's tests, all of them classical CPU runs of Julia with no quantum device or simulator anywhere. `test/Diffusion/MOL_1D_Linear_Diffusion.jl` (909 lines) is the relevant file. Its "Test 00: Dt(u(t,x)) ~ Dxx(u(t,x))" at line 25 grids $[0, \pi]$ at 30 points, builds four discretizations — the default, `grid_align = edge_align`, an explicit `approx_order = order` with `order = 2`, and `approx_order = 4` — and asserts each against $e^{-t}\cos x$ at every saved time with `isapprox` at `atol = 0.01`, dropping the two endpoint columns first. Test 03 at line 172 runs homogeneous Neumann conditions on 300 points and adds a conservation check asserting that the composite trapezoidal integral stays within `atol = 1.0e-9` of its value at the first saved time. Its testset name says "order 8", and line 203 does set `order = 8`, but that variable never reaches either `MOLFiniteDifference` call, so both arms in fact run at the constructor default `approx_order = 2` — the name is as stale as the comment a few lines below it, which still attributes the test's $\sim 4 \times 10^{-15}$ conservation figure to `Tsit5`, the explicit solver the v1 `DAEProblem` default no longer selects. Both tests call `solve(prob; ...)` without naming an integrator, so which algorithm actually ran is whatever that path defaults to and not a stated choice, and no error is attributed to the spatial stencil separately from the time integration in either test.

• The same three-point Laplacian assembled twelve times over in `HelloQuantum`'s heat-equation benchmark

  • Quantum PDE Solvers in Practice: Application-Driven Benchmarking of the Heat Equation

  Brightskies Inc.'s benchmark of eleven quantum kernels plus one classical baseline on the 1-D Dirichlet heat equation, written up as arXiv:2607.12688 and released as a repository that its own description calls "A simple introduction to quantum computing and programming through an easy to follow tutorial for the heat equation simulation." The paper's opening complaint is precisely about the layer this method occupies: "Quantum PDE solvers are difficult to evaluate in practice because published studies use different discretizations, output models, reconstruction rules, and hardware assumptions", and its answer is to hold the discretization fixed and vary only the solver — "eleven kernels under the same problem instances and readout contract". What is held fixed is the object recorded here. Every one of the twelve kernels starts from the same tridiagonal $-2/1/1$ stencil divided by $\Delta x^{2}$, and the reference all eleven quantum kernels are scored against is the exact continuous-time solution of that semi-discrete system, which the README's kernel table names a "Semi-discrete spectral reference computed in the Laplacian eigenbasis; the baseline every quantum kernel is scored against." Starting from is not the same as consuming, and the difference is what keeps this entry here rather than one capability over: the four linear-solver kernels — HHL, QSVT, QLS-Fourier and VQLS — go on to form $A = I - \alpha\,\Delta t\,L$ on top of the stencil, an implicit Euler step that belongs to the full-discretization slot, and the time discretization is per-kernel throughout: Trotter steps for the real-time kernels, none at all for the classical baseline. The stencil assembly is the only part common to all twelve. The README's own tutorial derivation, by contrast, is a full space-and-time scheme it names "Forward-Time, Central-Space (FTCS)", which is not what the executed classical kernel does.

  There is no shared discretization module. Each of the twelve kernel files builds the same stencil for itself, under eight different names, seven of them private: `_build_spectral_cache` in `classical_kernel.py`, `_build_discrete_laplacian` in `hhl_kernel.py` and `vqls_kernel.py`, `_build_laplacian` in `qsvt_kernel.py`, `qls_fourier_kernel.py` and `schrodingerization_kernel.py`, `_construct_laplacian` in `qsm_kernel.py` and `schade_hamiltonian_kernel.py`, `_heat_hamiltonian_dense` in `avqds_kernel.py`, `_laplacian_pauli_sum` in `var_qite_kernel.py`, `_generate_hamiltonian_terms` in `qite_kernel.py`, and the one public one, the module-level `generate_laplacian_pauli` in `hamiltonian_simulation_kernel.py`. The twelve reach the same matrix by three routes. Seven write an explicit index loop — `lap[i, i] = -2.0`, `lap[i, i - 1] = 1.0`, `lap[i, i + 1] = 1.0` — with no wraparound entries, because the boundary condition is Dirichlet rather than periodic; the array is called `laplacian` in the classical kernel, `L` in HHL, VQLS and Schrödingerisation, `lap` in QSM, Schade-Hamiltonian and AVQDS. Two, QSVT and QLS-Fourier, skip the loop entirely and assemble it vectorised, `return (np.diag(diag) + np.diag(off, 1) + np.diag(off, -1)) / dx ** 2`. Three never form the matrix at all, recursing on the qubit count in the Pauli basis from the single-qubit base case $-2I + X$. The scaling differs too: most builders return the bare stencil over $\Delta x^{2}$ and apply $\alpha$ later — the classical kernel at `eigh(-alpha * laplacian)`, Hamiltonian simulation at `scale = -alpha / (dx ** 2)` when it caches the Pauli sum — while Schrödingerisation's returns `alpha * L / dx ** 2` and AVQDS's writes `-alpha * lap` into its padded block. Three things about the copies are the interesting part. `ClassicalKernel` never steps in time at all: `_build_spectral_cache` calls `eigh(-alpha * laplacian)` once, and `compute_step`, typed to return `None`, accumulates `self._elapsed_time`, re-evaluates the exact solution from the cached initial condition with `decay = np.exp(-self._eig_vals * self._elapsed_time)` and writes the result into `grid.current_temperature[start:end]` rather than returning it; its class docstring calls this "an exact (up to floating-point precision) solution to the semi-discrete heat equation for a single time step", a phrase that reads oddly beside a kernel that recomputes from $t = 0$ on every call. `SchrodingerizationKernel` skips the numerical diagonalisation entirely and writes the eigenvalues in closed form, `eig_vals_N = alpha / dx ** 2 * (2.0 * np.cos(np.pi * k / (N + 1)) - 2.0)` with DST-I eigenvectors — the Dirichlet counterpart of the record's periodic $\lambda_j = -4\sin^{2}(j\pi/n)$, since $2\cos\theta - 2 = -4\sin^{2}(\theta/2)$. And `HamiltonianSimulationKernel._ensure_cache` refuses to run unless $r = \alpha\Delta t/\Delta x^{2} < 0.5$, raising "CFL violation" — a stability condition inherited from the README's FTCS framing rather than from the semi-discretization, and one that constrains the time step this method does not take.

  No external dataset. `ModelHandler.initialize_grid` reads every parameter from a JSON configuration and allocates `np.zeros(grid.domain_size, dtype=np.float32)`; the initial condition is one of three closed forms the paper names pulse, Gaussian and bimodal, applied by `_apply_pulse`, `_apply_gaussian` and `_apply_bimodal` in `src/component/source_injector.py`. The repository ships eight preset directories, `config/n3` through `config/n10`, each with its own `model.json` and `system.json`, plus a top-level `config/model.json` the README calls the legacy onboarding default and which reads `"size": 23`, `"alpha": 0.2`, `"dt": 0.2`, `"total_time": 120.0`. The paper's benchmark uses four of the eight, which the README names as `config/n4/` through `config/n7/`. All eight obey the same rule — `size` is $2^{n} + 2$ so the interior is exactly $N = 2^{n}$ points with one Dirichlet point each side, and `dt` is chosen to hold $r = \alpha\Delta t/\Delta x^{2}$ fixed. `config/n4/model.json` reads `"size": 18`, `"alpha": 0.01`, `"dt": 0.1384083044982699`, `"total_time": 1.0`, `"sampling": 0.058823529411764705`, `"boundary_length": 1`, `"half_length": 0`, `"num_shots": 100000` — that is $N = 16$ interior points, $\Delta x = 1/17$, and $r = 0.4$ exactly. Computed from their own values, $r$ is $0.4$ at all eight, from `n3` at $N = 8$ to `n10` at $N = 1024$. The boundary value is a separate file, `config/n4/system.json`, holding `{"boundary": 0.0}`.

  github.com/brightskiesinc/HelloQuantum, `main` at `6cdc6da3` (2026-07-24), read 2026-08-27; Python (`requirements.txt` pins `qiskit>=2.3.0`, `qiskit-aer>=0.17.2`, `numpy>=1.24.0`, `scipy>=1.10.0`); Apache-2.0, its 201-line `LICENSE` matched by a per-file header reading "Copyright 2026 Brightskies Inc." and "Licensed under the Apache License, Version 2.0 (the License); you may not use this file except in compliance with the License.", where the header itself spells License in escaped double quotes, reproduced here without them. The stencil-bearing files are `src/kernels/concrete/classical/classical_kernel.py` (105 lines; `class ClassicalKernel(BaseKernel)` at line 25, `_build_spectral_cache` at line 46, the loop at lines 55-62, `eigh` at line 64, `compute_step` at line 68); `src/kernels/concrete/quantum/hamiltonian_simulation_kernel.py` (589 lines; `generate_laplacian_pauli` at line 176, the CFL guard at lines 348-353, `scale = -alpha / (dx ** 2)` at line 357); `src/kernels/concrete/quantum/schrodingerization_kernel.py` (774 lines; `_build_laplacian` at line 200, the closed-form eigenvalues at line 268); `src/kernels/concrete/quantum/hhl_kernel.py` (622 lines; `_build_discrete_laplacian` at line 129, `_build_heat_equation_matrix` at line 141, `A = np.eye(N) - alpha_diff * dt * L` at line 151); `src/kernels/concrete/quantum/qsvt_kernel.py` (569 lines; the `np.diag` builder at lines 120-124, `A = np.eye(N) - alpha * dt * L` at line 131); and `src/kernels/concrete/quantum/qls_fourier_kernel.py` (479 lines; the same builder at lines 119-122, the same matrix at line 128). The marching driver is `src/engine/simulator.py` (118 lines), which computes `num_steps = int(total_time / dt) if dt > 0 else 0` at line 90 and calls each kernel's `compute_step(grid)` in a loop. `6cdc6da3` is both the repository's `main` HEAD and the most recent commit touching each of those files; the repository carries no tags or releases, and `CITATION.cff` declares `version: "1.0.0"` and `date-released: "2026-04-26"`, capitalising the same four author names differently from the arXiv page ("ElKarargy", "AbdelAziz", "AbdelRahman") and ordering Hatem before ElSayed.

  The numbers are published as tables, not only as figures. Under "Terminal Relative L2 Error" the README carries two per-kernel tables, eleven kernels by three initial conditions by three backends, one at $n = 4$ ($N = 16$) and one at $n = 7$ ($N = 128$), every cell a terminal relative $\ell_2$ error measured against this method's own output — the figure index names the reference "the semi-discrete classical reference". The three backends are all simulated — statevector ("Exact double-precision amplitude evolution"), an ideal shot-based `AerSimulator` at $10^5$ shots per step whose sampling floor the README puts at $\sim 1/\sqrt{N_{\text{shots}}} \approx 3 \times 10^{-3}$, and a noisy backend under "a portable depolarising/readout model" — with no run on quantum hardware. On the pulse initial condition at $n = 4$, statevector: QSM and Schade-Hamiltonian both read `<1e-15`, Schrödingerisation `4.23e-6`, QITE `6.13e-3`, VQLS `2.77e-2`, QSVT `8.90e-2`, HHL `1.36e-1`, var-QITE `1.64e-1`, Hamiltonian simulation `4.60e-1`, AVQDS `7.08e-1`, QLS-Fourier `8.04e-1`. At $n = 7$ the two transform kernels hold at `<1e-15` and Schrödingerisation rises to `5.64e-4`, while the rest degrade — HHL to `2.13e-1`, QLS-Fourier and AVQDS to `1.00e+0`. The noisy column is where it collapses: Schade-Hamiltonian goes from `<1e-15` on statevector to `9.02e-1` noisy at $n = 4$ and `1.00e+0` at $n = 7$, and HHL is `7.15e-1` noisy at $n = 4$ and a dash at $n = 7$ — a dash the README says marks a run "not performed for an implementation reason, not zero error", since "VQLS, Schrödingerisation, and QLS-Fourier are statevector-only by construction" and HHL on the noisy backend is restricted to $n = 4$ because "larger circuits fail Qiskit two-qubit synthesis". A third table, the norm-mismatch ablation at $n = 7$, attributes the smooth-IC plateau shared by Hamiltonian simulation, AVQDS and QLS-Fourier to reconstruction normalisation rather than to any algorithm: the norm-only error is 94%, 100% and 100% of the residual on the pulse initial condition and, for all three kernels alike, 23% on Gaussian and 29% on bimodal. Nothing in the benchmark measures the truncation error of the spatial stencil itself — the semi-discrete solution is the reference, so the $O(\Delta x^{2})$ term this method pays is common to every arm and invisible in every cell of all three tables.

What it needs

Nothing below this — it bottoms out here.

Other ways to fill the same slot

Different approaches

• Graph-Laplacian finite differences

  Discretize the domain onto a lattice and read the discrete Laplacian off the resulting graph: off-diagonal entries minus one between neighbours, each diagonal entry the degree of its vertex. Higher-order stencils are obtained by factorizing the operator through hypergraph incidence matrices, which is what lets the error fall faster than the second power of the spacing while keeping a form a simulator can consume.

In the Atlas

No record in the Atlas covers this yet. The catalogue is circuits and primitives; this part of the literature is not in it.

Sources

• Quantum vs. classical algorithms for solving the heat equation↗Noah Linden, Ashley Montanaro, Changpeng Shao · 2020theory yes · simulation ? · hardware no