This paper presents an efficient quantum spectral-Galerkin method for constant-coefficient second-order linear elliptic equations. The method employs a boundary-adapted Jacobi tensor-product basis that incorporates the boundary condition into the approximation space and produces Galerkin systems with polynomially bounded condition numbers and sparsity. For analytic solutions satisfying a suitable uniform regularity condition, spectral convergence allows the normalized Galerkin coefficient state to be prepared with complexity $\mathrm{poly}(d,\log(1/ε))$, where $d$ is the spatial dimension and $ε$ is the error tolerance, under standard sparse-access and state-preparation assumptions. At the state-preparation level, this improves the precision dependence over fixed-order finite-difference and finite-element discretizations. We further show how to estimate prescribed linear functionals of the solution, with an additional $O(1/ε)$ factor arising from norm and inner-product estimation.
This paper presents an efficient quantum spectral-Galerkin method for constant-coefficient second-order linear elliptic equations. The method employs a boundary-adapted Jacobi tensor-product basis that incorporates the boundary condition into the approximation space and produces Galerkin systems with polynomially bounded condition numbers and sparsity. For analytic solutions satisfying a suitable uniform regularity condition, spectral convergence allows the normalized Galerkin coefficient state to be prepared with complexity poly(d,log(1/ε)), where d is the spatial dimension and ε is the error tolerance, under standard sparse-access and state-preparation assumptions. At the state-preparation level, this improves the precision dependence over fixed-order finite-difference and finite-element discretizations. We further show how to estimate prescribed linear functionals of the solution, with an additional O(1/ε) factor arising from norm and inner-product estimation.
Quantum computers can produce a quantum encoding of the solution of a system of differential equations exponentially faster than a classical algorithm can produce an explicit description. However, while high-precision quantum algorithms for linear ordinary differential equations are well established, the best previous quantum algorithms for linear partial differential equations (PDEs) have complexity p o l y ( 1 / ϵ ) , where ϵ is the error tolerance. By developing quantum algorithms based on adaptive-order finite difference methods and spectral methods, we improve the complexity of quantum algorithms for linear PDEs to be p o l y ( d , log ( 1 / ϵ ) ) , where d is the spatial dimension. Our algorithms apply high-precision quantum linear system algorithms to systems whose condition numbers and approximation errors we bound. We develop a finite difference algorithm for the Poisson equation and a spectral algorithm for more general second-order elliptic equations.
Spatial derivatives pose a basic challenge for quantum algorithms for partial differential equations: their standard discretizations have operator norms that grow as the mesh is refined, leading to large block-encoding normalization factors. Even when preconditioning improves access to the solution, estimating derivative-dependent physical observables can reintroduce polynomial dependence on the inverse mesh width. We develop a multilevel quantum algorithm for estimating linear and quadratic observables of elliptic equations that addresses this readout cost. The algorithm decomposes the output into corrections across nested discretizations and realizes fine--coarse cancellation coherently, before measurement. A Ritz-complement factorization connects the O(h2) size of the corrected response with a block encoding at the same scale, allowing this decay to offset growing readout normalizations. Under efficient quantum access to the corrected coordinates, input data, and readouts, observables with readout scale O(h−χ), 0≤χ≤2, can be estimated to additive accuracy ε at cost O(ε−1) using amplitude estimation or O(ε−2) using direct sampling, with only polylogarithmic dependence on the inverse finest mesh width hL−1. We give explicit constructions of the required operator oracles for one-dimensional piecewise linear finite elements with diffusion piecewise constant on a fixed coarsest partition, and for Fourier spectral hierarchies. Output-specific bias estimates connect the discrete results to continuum accuracy.
Quantum linear-system algorithms (QLSAs) have rigorous worst-case complexity guarantees, but their runtimes are often chosen from spectral information assumed in advance. What is largely lacking is an a posteriori progress flag: most QLSA workflows, unlike the classical counterparts, do not provide a built-in mechanism to signal whether a particular instance has already converged.
For discretizations of elliptic PDEs −∇⋅(a(x)∇u(x))=f(x), with divergence--gradient structure −∇⋅(a(x)∇)≈Ah=Gh†Gh, we formulate a stable first-order ODE whose limiting solution block is the desired Galerkin solution. The PDE-dependent scale is then \(\norm{G_h}=\bigO(h^{-1})\), comparable to factorized QLSA constructions with square-root condition-number scaling. We design an augmented dynamics with residual variables, in which measuring a residual register gives an on-the-fly convergence indicator without reconstructing the solution vector. For smooth right-hand sides, dynamic stopping can reduce the evolution time and gate count relative to a fixed worst-case schedule, and may also reduce exposure to accumulated hardware errors. Numerical experiments for a two-dimensional finite element Poisson problem show that the residual-register probability follows the actual error decay and, for some right-hand sides, can stop the quantum circuit well before a conservative worst-case runtime estimate is reached.
We develop a Schrödingerization algorithm for quantum linear systems problems in two higher dimensions. The earlier LC-Schrödingerization approach represents the solution as a time integral of a homogeneous convection response and implements this integral by a linear combination of unitaries (LCU) over evolution times. We instead use Duhamel's principle to incorporate the integral into an inhomogeneous convection equation with zero initial data. Fourier projection in the convection variable and Schrödingerization in a second auxiliary variable then realize the solution without the separate LCU step. A joint choice of the kernel and recovery procedure gives an evolution time independent of the target accuracy, a uniformly bounded L2 kernel normalization, and accurate recovery on a fixed interval. We establish the periodization and discretization bounds and analyze the interval recovery probability. With block preconditioning, we retain a single logarithmic precision factor and obtain linear condition-number dependence without variable-time amplitude amplification. Under exact oracle access and given a valid constant-factor solution-norm estimate, the algorithm uses O(κAlog(1/ε)) queries to each original input oracle, with constant success probability and ℓ2 state error at most ε. The matrix-query bound matches the standard worst-case scaling when the supplied norm bounds are tight. UnitaryLab simulations on positive-definite and indefinite systems demonstrate the solution accuracy and probability gain from interval recovery.
This paper presents a quantum algorithm for the solution of prototypical second-order linear elliptic partial differential equations discretized by d d -linear finite elements on Cartesian grids of a bounded d d -dimensional domain. An essential step in the construction is a BPX preconditioner, which transforms the linear system into a sufficiently well-conditioned one, making it amenable to quantum computation. We provide a constructive proof demonstrating that, for any fixed dimension, our quantum algorithm can compute suitable functionals of the solution to a given tolerance t o l \mathtt {tol} with an optimal complexity of order t o l − 1 \mathtt {tol}^{-1} up to logarithmic terms, significantly improving over existing approaches. Notably, this approach does not rely on the regularity of the solution and achieves a quantum advantage over classical solvers already in two dimensions, whereas prior quantum methods required at least four dimensions for asymptotic benefits. We further detail the design and implementation of a quantum circuit capable of executing our algorithm, present simulator results, and report numerical experiments on current quantum hardware, confirming the feasibility of preconditioned finite element methods for near-term quantum computing.
Can a quantum PDE algorithm avoid the polynomial cost of resolving a fine spatial mesh? For standard fixed-order discretizations, direct classical methods require work polynomial in h−1, or equivalently in the number of spatial degrees of freedom Nh=Θ(h−d). Direct quantum implementations of a parabolic semigroup still have coherent complexity O(T/h), and gradient-dependent observables such as heat flux and dissipation introduce additional mesh dependence. Decay of the solution norm will further suppress the postselection probability for preparing a normalized final state.
We develop a multilevel quantum algorithm that estimates linear and quadratic observables directly and places the fine--coarse cancellation inside the circuit before measurement. A contour-based LCU reconstructs each target-time correction from a coherent family of shifted resolvent differences. Rather than block encoding the fine and coarse inverses separately, we encode their difference through a shifted Ritz--Schur factorization, exposing its O(hℓ2) two-grid normalization. For Fourier hierarchies, the corresponding SELECT oracle consists of a quantum Fourier or sine transform, a spectral-band selector, and reversible diagonal arithmetic. We also give a non-Fourier realization based on energy-orthogonal dyadic midpoint details in one dimension, together with structured tensor-product extensions under fixed-rank coefficient and access assumptions.
For readouts with derivative order 0≤χ≤2, optimized amplitude estimation removes all polynomial dependence on the finest mesh size. Under the stated access assumptions, both linear and quadratic observables can be estimated with complexity O(1+(Tε)−1), with only polylogarithmic dependence on h−1.