Graph-Laplacian finite differences

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

The truncation bound assumes the exact solution has the derivatives the stencil's Taylor expansion uses: Costa, Jordan and Ostrander state that a discretization with error $O(a^k)$ of an $m$-th derivative 'is only justified if the exact solution is $(k+m)$-times differentiable'. A scatterer is modelled as removed lattice points, and the presence of one breaks translational invariance, so the Laplacian can no longer simply be diagonalized by a Fourier transform. The hypergraph incidence-matrix factorizations the higher orders rest on are, in the authors' own words, not known to appear elsewhere in the literature.

Requires

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

Example

given  a linear PDE reducible to the wave equation d^2 phi/dt^2 = grad^2 phi
       on a domain Omega in R^D, initial data phi(x,0) and dphi/dt(x,0),
       a lattice spacing a, a boundary condition -- Dirichlet (phi = 0) or
       Neumann (grad phi . n = 0) -- and a target order k >= 2 (even)
       -- k > 2 is supported only with Dirichlet boundary conditions      (Sec. 7.2)

requires  the exact solution to be (k+2)-times differentiable, since the
    O(a^k) discretization of the Laplacian (a second derivative, so m=2
    here) below is only justified to that order of smoothness            (Sec. 9)
requires  if a scatterer is present, that it be modelled as removed
    lattice points -- this breaks translational invariance, so the
    resulting Laplacian can no longer be diagonalized by a Fourier
    transform the way the periodic-lattice construction below is         (Sec. 1, fn. 1)

if Neumann and k > 2: stop -- unsupported, see the Neumann branch below   (Sec. 7.2)

# --- lay out the lattice graph -----------------------------------------------
Ga = cubic grid of spacing a covering Omega, edges between nearest neighbours
    # a scatterer is modelled by deleting its interior vertices and their
    #   incident edges before anything below runs

# --- k = 2: read the Laplacian directly off the graph -------------------------
# both boundary conditions at k=2 are built straight off Ga; no periodic or
#   circulant matrix is constructed anywhere in this branch                (Sec. 3)
if k = 2:
    L_ij = -1            if vertices i, j are adjacent
    L_ii = deg(i)         # number of vertices i is connected to
    L_ij = 0              otherwise                                       (Eq. 2)
        # -(1/a^2) L approximates the continuum Laplacian; at finite a the
        #   truncation error is O(a^2)                                    (Eq. 2)
    if Neumann:
        L <- exactly the above, on Ga with the scatterer already removed;
             add no self-loops                                            (Sec. 3)
    else if Dirichlet:
        for each boundary vertex v: add a weighted self-loop at v, weight
            = (edges v is missing relative to an interior vertex)
            # forces every diagonal entry of L to the interior degree, so
            #   phi is pinned to 0 outside the boundary                   (Sec. 3, Fig. 2)

# --- k > 2 (Dirichlet only): raise the order via Lagrange interpolation -------
else:
    N = k/2
    take the second derivative of the Lagrange interpolation formula and
        expand to radius N; this step assumes only a UNIFORM lattice,
        x_j = j*a -- not a periodic one                                   (Eq. 40, 41)
    assume in addition periodic boundary conditions, so that the same
        formula holds at every point and not just at x0; only then can L
        be written as a circulant matrix in the cyclic-shift S,
        e.g. for k=4:  L = -(1/a^2)( (5/2)*1 -(4/3)(S+S^T) +(1/12)(S^2+(S^T)^2) )
                                                                          (Eq. 43)
        # the paper introduces periodicity here and only here, to turn the
        #   pointwise formula of Eq. 41 into a single matrix              (Sec. 6)
    # coefficients for orders up to k=10 are tabulated, not re-derived     (App. C)
    L <- principal submatrix of that periodic L, restricted to the
         interior vertices
        # this restriction IS the Dirichlet condition at higher order, and
        #   it is the only branch here that consumes the periodic L       (Eq. 44)
    # Neumann is excluded above because this higher-order L is not
    #   symmetric once restricted to the interior vertices, and the
    #   decoupled second-order dynamics below need a symmetric (Hermitian)
    #   generator                                                         (Sec. 7.2)

# --- factor L = B B^T through an incidence matrix -----------------------------
if k = 2:
    B_ij = sqrt(W_j)   if j is a self-loop of i, or an edge with i as source
    B_ij = -sqrt(W_j)  if j is an edge with i as sink
    B_ij = 0           otherwise                                          (Eq. 8)
    # holds for any graph, weighted or unweighted, with or without
    #   self-loops; the arbitrary edge orientation changes B but not
    #   B B^T, which always equals L of the undirected graph              (Eq. 8)
else:
    ansatz B = sum_j a_j S^j over the hyperedges of radius N, each
        touching up to N+1 vertices; solve the resulting system of
        degree-2 polynomial equations in the a_j for the entries
            # e.g. at k=4: B = c*S - (c+b)*1 + b*S^T solves to
            #   c ~= 1.07735, b ~= 0.07735                                (Eq. 49-51)
    # for D > 1 dimensions: build one such (L_d, B_d) pair per axis d
    #   from the disconnected path/cycle subgraph in that direction, then
    #   vertically concatenate the B_d into one incidence matrix          (Sec. 7.4)

# --- assemble the block Hamiltonian and read off the linear IVP --------------
H = (1/a) * [[0, B], [B^T, 0]]
    # Hermitian by construction for ANY B, regardless of which
    #   factorization of L produced it                                    (Eq. 4)
psi(0) = (phi_V(0), phi_E(0))  with  phi_V(0) = phi(x,0) sampled on Ga
    # phi_E(0) is left unconstructed here -- exact preparation requires
    #   inverting B, which is a separate state-preparation step, not part
    #   of this discretization                                            (Sec. 4.3)

return  generator H and initial vector psi(0) for the linear IVP
        d(psi)/dt = -i H psi, together with the truncation error O(a^k)
        this replacement carries
# differentiating the IVP once more gives d^2(phi_V)/dt^2 = -(1/a^2) B B^T
#   phi_V = -(1/a^2) L phi_V, i.e. exactly the discretized wave equation on
#   the vertex block alone; the edge block obeys the analogous equation for
#   B^T B and is not itself the quantity of interest                      (Eq. 6, 7)
# the order is not free of what the resulting Hamiltonian costs to simulate:
#   a D-dimensional order-k Laplacian has a D(k/2+1)-sparse incidence
#   matrix, so an s-sparse H forces k = 2(s/D) - 2                        (Sec. 8)

Cost, as the source states it

Costa, Jordan and Ostrander give the truncation error of the second-order stencil as $O(a^2)$ at finite lattice spacing $a$, and generalise it: a $k$-th order Laplacian gives truncation errors of order $a^k$. The order is not free of what the stencil feeds — a $D$-dimensional Laplacian of order $k$ has a $D(k/2+1)$-sparse incidence matrix, so an $s$-sparse Hamiltonian corresponds to $k = 2(s/D) - 2$, which is the order-against-sparsity trade in the form the paper states it.

Implementations

• CERFACS qaths — the lattice graph, its Laplacian and its incidence matrix compiled into oracles

  • Practical Quantum Computing: solving the wave equation using a quantum approach

  Built because the algorithm had never been run. Suau, Staffelbach and Calandra open by saying that despite the large number of algorithms available, it is hard to find an actual implementation of a quantum differential equation solver, Hamiltonian simulation being the unique exception by solving the time-dependant Schrödinger equation, and that to the best of their knowledge their work is the first to analyse experimentally the characteristics of a quantum PDE solver. Their framing is cost rather than correctness: an experimental study of the costs, in gate number and in execution time on an idealised hardware created from realistic gate data, of one of the direct quantum algorithms, namely the wave equation solver of Costa, Jordan and Ostrander. The instance is what section II calls a simplified version of the wave equation — the one-dimensional line $[0,1]$ with constant propagation speed $c = 1$, their Eq. (1), with no assumption made on the initial data; the paper claims no minimality for it, and the construction it borrows is stated in any number of dimensions. The boundary condition their Eq. (2) writes as $\frac{\partial}{\partial x}\phi(0,t) = \frac{\partial}{\partial x}\phi(1,t) = 0$ is called Dirichlet throughout the paper; what they then build is the self-loop graph Costa, Jordan and Ostrander assign to Dirichlet, and Appendix B 2 drops the two boundary nodes on the ground that their value is always equal to 0, so the graph rather than Eq. (2) is what the implementation follows.

  The discretization is Appendix B and it is the graph-Laplacian recipe read straight off the source paper. Space becomes a graph $G_{\delta x}$ whose vertices are the discretisation points and whose edges join nearest neighbours; Eq. (B1) defines $L(G_{\delta x})_{i,j}$ as $\deg(v_i)$ when $i = j$, as $-1$ when $i \neq j$ and $v_i$ is adjacent to $v_j$, and as $0$ otherwise; Eq. (B3) sets $A = -\frac{1}{\delta x^2}L(G_{\delta x})$, which Eq. (B4) checks against the second-order stencil $(\phi_{i-1,t} - 2\phi_{i,t} + \phi_{i+1,t})/\delta x^2$ of Eq. (B2). Eq. (B6) is the block Hamiltonian carrying $B$ above the diagonal and $B^\dagger$ below it, and Eq. (B7) is the factorization condition $BB^\dagger = L(G_{\delta x})$. The incidence matrix is assembled in three steps that the paper states as an algorithm: index the vertices arbitrarily in $[0, N_d-1]$, orient and index each edge arbitrarily in $[0, N_d-2]$, then read Eq. (B8) — $1$ for a self-loop of $i$ or an edge with $i$ as source, $-1$ for an edge with $i$ as sink, $0$ otherwise. The paper makes the arbitrariness load-bearing rather than incidental: changing the orientation or either ordering changes $B$ but not $BB^\dagger$, and that freedom is spent picking the ordering and orientation that produce an easy-to-implement $B$, which for Dirichlet is the bidiagonal $B_d$ of Eq. (B9) whose $B_dB_d^\dagger$ is the tridiagonal matrix with $2$ on the diagonal and $-1$ beside it, Eq. (B10). From there the work is compilation, not mathematics: $H$ is decomposed into two $1$-sparse Hermitian matrices, both of spectral norm $1$, and each is handed to three oracles — $M$ returning the column index of the first non-zero element of a row, $V$ its absolute value, $S$ its sign — implemented in oracles.py as get_oracle_dirichlet1_1d_wave_equation and get_oracle_dirichlet2_1d_wave_equation, bundled with the arithmetic routines by linking_sets/basic.py's get_linking_set, and stepped by evolve_1D_dirichlet.py's evolve_1d_dirichlet. One substitution against the source paper is deliberate and stated: Costa, Jordan and Ostrander chose a quantum-walk Hamiltonian simulation for its nearly optimal asymptotics, and this implementation uses a Trotter-Suzuki product formula instead, at order $k = 1$, for its good experimental results and simpler implementation.

  No dataset. Everything is synthetic and generated in the repository: the validation initial condition is a sine sampled on the lattice, exposed as the console script qaths.simulation.1D_dirichlet.sin over _cli/solving/dirichlet_1D_sin.py, and the reference against which the circuits are checked is a classically assembled pair of sparse matrices built by construct_Dirichlet_Hamiltonians_1D in tests/utils/matrices.py, whose docstring says it returns 2 Hamiltonian matrices that sum to the Hamiltonian described in the paper and builds them as scipy.sparse.coo_matrix incidence blocks. The only external data in the whole study are hardware numbers, not inputs: IBM Q Melbourne gate times, with the GF pulse time approximated via arithmetic mean to 347ns, the GD pulse time 100ns and the buffer time 20ns. Parameter settings are stated per figure rather than swept blind, and the paper's two benchmarks are held apart: the validation run is $N_d = 32$, $t = 0.4$, $\varepsilon = 10^{-3}$; the solver sweeps of FIG. 2 and FIG. 4, in section IV B, hold $t = 1$, $\varepsilon = 10^{-5}$ and $k = 1$ while $N_d$ runs; the Hamiltonian-simulation sweeps of FIG. 1, in section IV A, instead fix two of 32 discretisation points, $t = 1$ and $\varepsilon = 10^{-5}$ and sweep the third, so $\varepsilon$ is the swept axis in FIG. 1(c) rather than a held constant there; and the automated tests parametrise discretisation_points_number over the two values $4$ and $6$ only.

  The paper's Supplementary Material names one artefact — the implementation of the quantum wave equation solver is available at https://gitlab.com/cerfacs/qaths — and that GitLab project (id 10961444, default branch master, public, created 2019-02-21, last activity 2023-06-20) is the entry's subject. It is Python 3.6 against qat, the Atos Python library shipped with the Quantum Learning Machine, whose setup.py declares install_requires as numpy, scipy, matplotlib, myqlm, qprof and qprof-myqlm. Licence is CeCILL-B: the LICENSE file leads with the note that all files with the exception of those in third_party/ are under the CeCILL-B license agreement, though no third_party/ directory exists on master. The header block sits beneath the attribution rather than above it — the copyright and contributor lines come first and the sentence This software is governed by the CeCILL-B license under French law and abiding by the rules of distribution of free software. follows — and it is not universal: of the 95 Python files on master read 2026-08-27, 93 carry it, in 68 of them under a Copyright CERFACS line with a month, in one under Copyright CERFACS / LIRMM (12/2019) and in 24 under Copyright TOTAL / CERFACS / LIRMM with a 2020 date, that last set including both linking_sets modules named below, every one of the 93 naming Adrien Suau as contributor; src/qaths/qram.py and docs/conf.py carry no header at all. Reading the tree the same day, the discretization lives in four files under src/qaths/applications/wave_equation/: oracles.py at 1533 lines, which carries the whole graph-to-matrix-to-oracle derivation in its module docstring and draws it with the four figures that docstring embeds by an image directive — docs/images/non_oriented_dirichlet_graph.png, non_oriented_neumann_graph.png, oriented_dirichlet_graph_numbered.png and oriented_neumann_graph_numbered.png, four of the seven files checked in under docs/images/; evolve_1D_dirichlet.py at 399 lines, whose unconditional entry points are evolve_1d_dirichlet, evolve_1d_dirichlet_no_repetition, evolve_1d_dirichlet_no_time_adjustment and solve_1d_dirichlet_stationary, with a fifth, initialise_1d_dirichlet_stationary, defined inside a try importing StatePreparation from qat.linalg.oracles, so that under the myqlm setup.py actually declares the except ImportError branch fires and logs that it is disabling the initialise_1d_dirichlet_stationary function; utils.py at 78 lines, whose compute_qubit_number_from_considered_points_1d fixes the arity at $n = \lceil\log_2(2N_d-3)\rceil$, its docstring defining $N_d$ as the number of discretisation points and the function's own argument as $N_d-2$ because the two boundary points are not considered; and linking_sets/basic.py at 79 lines beside linking_sets/arithmetic_adder.py at 80 lines. The last commit touching oracles.py is 65feb48a of 2020-04-07, Adapted oracles comparator AbstractGate to a normal comparator; the newest commit on master is 3a510af4 of 2021-11-02; the only tag is 1.0.0a, dated 2019-05-29, so the repository has never been tagged at the version 1.0.0 its own setup.py declares. The README is explicit that this repository is not maintained on a daily basis anymore and lists Implementing higher-order Laplacian discretisation. among features not built — which is the honest reading of what the artefact covers, since Appendix D derives the fourth-order factorization on paper and the paper states that the results shown in it have all been generated using the second-order discretisation formula.

  Nothing ran on quantum hardware and the paper says why: even if IBM Q Melbourne has 14 qubits, the quantum circuits constructed in this paper are not runnable because they require more qubits, so Melbourne enters only as a gate set and a table of pulse times. The runs are on qat simulators and the Atos QLM, access to which the acknowledgements credit to Reims University, the ROMEO HPC center, Total, the CCRT and Atos. The correctness result is the one number worth carrying: at $N_d = 32$, $t = 0.4$ and a requested $\varepsilon = 10^{-3}$, FIG. 3's caption reports the two solutions as too close to be able to notice a difference, they overlap on the graph, and the error between them is of the order of $10^{-7}$, which the paper reads as 4 orders of magnitudes smaller than the error we asked for. The cost results are fits rather than proofs, and the paper labels them as such: for the solver, section IV B's FIG. 2 gives the qubit count as $11 + 3\log_2(N_d)$, FIG. 4(a) the gate count against discretisation size as $\lambda N_d^{3/2}\log_2(N_d)^2$ at $\lambda = 300\,000$, and FIG. 4(b) the IBM Q Melbourne execution-time estimate at that same shape with $\delta = 0.06$, every constant chosen arbitrarily to fit the experimental data. The precision fit belongs to the other benchmark and should not be read as the solver's: FIG. 1(c), in section IV A on Hamiltonian simulation, plots the gates needed to simulate the Hamiltonian of Appendix B using the oracles of Appendix C, at 32 discretisation points and $t = 1$, against $\varepsilon$, and fits $\alpha\varepsilon^{-1/2}$ at $\alpha = 130\,000$; the paper prints no precision panel for the solver itself, though the repository ships the generator for one, _cli/data_generation/ibmq_gate_number_vs_precision_wave_equation_solver.py, beside its Hamiltonian-simulation twin. The derived complexity for the solver at $k = 1$ is Eq. (23), $O(N_d^{3/2}\log_2(N_d)^2 t^{3/2}/\sqrt{\varepsilon})$, and it carries no numerical constant inside the $O$: Eq. (22) states the general case as $O(5^{2k} t N_d \log_2(N_d)^2 (tN_d/\varepsilon)^{1/(2k)})$, whose coefficient at $k = 1$ would be $25$, and the paper drops it. The conclusion the authors draw against their own algorithm is the part a reader should not lose: even though the asymptotic scaling is better than classical algorithms, they found out that the constants hidden in the big-O notation were huge enough to make the solver less efficient than classical solvers for reasonable discretisation sizes.

• QuDiffEq.jl — the incidence matrix written out by hand, checked against OrdinaryDiffEq

  A Julia package from the QuantumBFS organisation, the same organisation as the Yao.jl simulator it is built on, written during Julia's Season of Contribution 2019. Its README describes it as quantum algorithms for solving differential equations and lists four references, of which the fourth is Costa, Jordan and Ostrander's wave-equation paper; the other three are the linear- and nonlinear-ODE algorithms the package's src/ implements. The wave equation is therefore not the package's subject but its worked example, and the example exists to exercise the package's Taylor-truncation Hamiltonian simulation on a generator that means something physical. What makes it an artefact for this method rather than for the simulation underneath it is that the discretization is not imported from anywhere: the incidence matrix is typed out in the example itself, and the Laplacian is obtained from it rather than written down directly.

  The whole thing is 81 lines. Vertices and edges are fixed at vertx = 7 and ege = 8 — seven interior points of a segment, six oriented edges between them and the two self-loops that make it Dirichlet — and the incidence matrix is built by the loop B[i,i] = -1 followed by B[i,i+1] = 1 over the seven rows, then by the single override B[1,1] = 1. That override changes no sparsity pattern: column 1 already held exactly one entry, because no row 0 exists to carry its partner, so it was a self-loop before the line ran and all the line does is flip its sign from the $-\sqrt{W_j}$ that Eq. (8) of the source paper assigns a sink to the $+\sqrt{W_j}$ it assigns a self-loop — the entry is squared on the way to $BB^{T}$, which is identical either way. The second self-loop needs no override at all: column 8 falls out of the loop's own last row, B[7,8] = 1. Columns 2 to 7 are the six oriented path edges, and the eight columns together give a $BB^{T}$ that is tridiagonal with $2$ on the diagonal and $-1$ beside it, which is the Dirichlet graph Laplacian of Eq. (8) reduced to its one-dimensional case. make_hamiltonian, taking $B$, its transpose and $a$, pads to nextpow(2, vertx+ege), so a $15 \times 15$ operator becomes $16 \times 16$ and the register is four qubits wide, writes $B$ into the upper-right block and $B^{T}$ into the lower-left, and returns -im/a*H, which is $-iH$ for the source paper's Eq. (4) $H$, namely $\frac{1}{a}$ times the block matrix carrying $B$ above the diagonal and $B^{\dagger}$ below it; the $1/a$ in the code is the one Eq. (4) already carries, not a second division. do_pde then calls taylorsolve, from src/TaylorTrunc.jl, ten times at Taylor order $k = 2$ and step $t = 10^{-2}$, rescaling each output by the inverse-probability factor the solver returns and keeping the first vertx entries as the vertex block. The check is the interesting half: the Laplacian is recovered from that same $B$ by the constant D1, assigned $-\frac{1}{a^2}$ times $B$ times its transpose and commented Laplacian Operator, and integrated classically by OrdinaryDiffEq's Tsit5() at fixed step dt=0.01 as the first-order system $du_1/dt = D u_2$, $du_2/dt = u_1$, whose $u_1$ obeys $d^2u_1/dt^2 = -\frac{1}{a^2}BB^{T}u_1$. So what the example certifies is the quantum stepper against the same discrete operator, not the discretization against the continuum — the truncation error $O(a^2)$ is common to both arms and cancels out of the comparison.

  No dataset. The initial field is a sine sampled by sn = sin.((0.0:de:(vertx-1)*de)*2*pi/((vertx-1)*de)) at de = 0.5; de cancels out of that expression, which reduces to $\sin(2\pi j/6)$ for $j = 0,\dots,6$, one full period sampled at the seven vertices and vanishing at the first and the last, which is what the Dirichlet graph wants. The initial velocity is zero, marked by the comment Intial conditions (stationary to begin with), the typo the file's. The spacing inside the operator is a separate constant, a = 1e-1, commented spactial discretization, and it appears identically in the Hamiltonian and in D1, so both arms carry the same $a$ and therefore the same discrete operator — which is what makes them comparable at this $a$. The example fixes $a$ once and never varies it, so nothing here measures how the agreement would hold at another value; $a$ sets the norm of the generator, and both the fixed Taylor order and the fixed classical step are sensitive to that.

  GitHub, QuantumBFS/QuDiffEq.jl, file Examples/wave_equation_1D.jl, Julia, 81 lines read 2026-08-27. Licence MIT, the LICENSE file opening Copyright (c) 2019 Divyanshu Gupta followed by the standard permission grant; the README repeats that the project is licensed under the MIT License. The file has exactly two commits and has not been edited since the later of them, a3243855 of 2019-09-12, wave equations sim, by Divyanshu Gupta; master has taken 45 commits since that date, the newest being 83dd4b4a of 2024-08-24, update API. Tags are v0.1.0 and v0.1.1 only. The functions a reader would open are the example's own make_hamiltonian and do_pde, and taylorsolve in src/TaylorTrunc.jl, whose docstring reads Simulates a Hamiltonian using the Taylor truncation method. Returns the state register and inverse probability of finding it. One caveat belongs in this paragraph rather than in a footnote: test/runtests.jl includes six test files and Examples/ is not among them, so the assertion at the end of the example is never executed by the repository's CI workflow.

  The example's own assertion is the only figure it produces: isapprox.(res1,s1, atol = 1e-2) over all entries, that is, the eleven stored vertex vectors — the initial one plus ten steps — agree with the eleven classical ones to an absolute tolerance of $10^{-2}$. The author's Nextjournal write-up runs exactly this cell and records the outcome as Test Passed at a wall time of 3.1s. The simulator is Yao.jl's state-vector simulator running in-process; no hardware, no noise model and no gate count are reported anywhere in the example or the post. The scale is small enough to state plainly — seven vertices, eight edges, a four-qubit vertex-plus-edge register before the Taylor ancillas, ten steps of $10^{-2}$ over $t \in [0, 0.1]$ and eleven stored states — so this artefact is evidence that the factorization is implementable and self-consistent, not evidence about cost.

What it needs

Nothing below this — it bottoms out here.

Other ways to fill the same slot

Different approaches

• Central differences, space only

  Replace each second spatial derivative by the three-point central difference on a uniform grid and leave the time derivative alone. What comes out is the plainest form the method of lines takes: one generator, assembled from a single stencil repeated at every interior point, acting on the vector of grid values.

In the Atlas

• Quantum algorithm for simulating the wave equation

  Simulate the wave equation under Dirichlet and Neumann boundary conditions on a quantum computer, using Hamiltonian simulation and quantum linear system algorithms as subroutines.

Sources

• Quantum Algorithm for Simulating the Wave Equation↗Pedro C.S. Costa, Stephen Jordan, Aaron Ostrander · 2017theory yes · simulation ? · hardware no