Quantum Mathematics Learning Guide 63: Quantum Differential Equations, Linear ODEs, PDEs, Spectral Methods and Solution States

Quantum algorithms for differential equations do not usually print an entire numerical solution vector faster than a classical computer. Their central promise is different: under strong sparsity, conditioning and state-preparation assumptions, they can prepare a quantum state proportional to a discretised solution and estimate selected global observables with complexity that may scale polylogarithmically in the grid dimension.

This distinction is essential. A partial differential equation discretised on N grid points produces N unknown amplitudes. A quantum computer can encode N amplitudes in log₂N qubits, but extracting all N numbers requires Ω(N) information in the output. Quantum advantage is therefore plausible when the desired answer is itself compressed: an expectation value, overlap, norm ratio, probability, integral or another low-dimensional functional of the solution state.

Guide 19 owns the general quantum linear-systems problem. This guide owns the differential-equation pipeline that creates those systems: ODE time stepping and history states, finite-difference PDE matrices, Chebyshev/spectral collocation, conditioning, boundary conditions and physical readout.

Differential equation → discretise or embed time evolution → sparse/structured linear system → quantum linear-system or simulation primitive → solution state → estimate a compressed observable.

1. Linear ordinary differential equation

Consider

dx/dt=A(t)x(t)+b(t)

with x(t)∈ℂ^d.

If A and b are time independent, the exact solution is

x(t)=e^{At}x(0)+∫_0^t e^{A(t-s)}b ds.

The quantum problem is to encode information about this vector-valued function using oracles for A, b and the initial condition.

2. Homogeneous case

If b=0,

x(t)=e^{At}x_0.

If A=−iH with H Hermitian, this is ordinary unitary Schrödinger evolution and Guide 13/Guide 4 simulation methods apply.

General A is not anti-Hermitian. Norms can grow or decay, and e^At is not unitary. Quantum algorithms must therefore embed the nonunitary linear dynamics inside a larger unitary or linear-system construction.

3. Scalar decay example

Take

dx/dt=−λx, x(0)=1.

The exact solution is

x(t)=e^{-λt}.

For λ=2 and t=1.5,

x=e^{-3}≈0.049787.

This one-dimensional case is not a quantum-advantage problem; it is a calibration example for discretisation error.

4. Forward Euler as a linear recurrence

With step h=T/M, forward Euler gives

x_{k+1}=(I+hA)x_k+h b_k.

For the scalar decay equation,

x_{k+1}=(1−hλ)x_k.

The method is stable only in a restricted region of hλ. Quantum acceleration cannot rescue an unstable classical discretisation chosen before the quantum linear algebra begins.

5. History-state linear system

Stack all time slices into one vector

X=(x_0,x_1,…,x_M)^T.

The recurrence can be written as a block lower-bidiagonal linear system

L X = B.

For homogeneous Euler, L has identity blocks on the diagonal and −(I+hA) on the first subdiagonal.

A quantum linear-system algorithm can then prepare a state proportional to the history vector X, assuming L is sufficiently sparse/well conditioned and B can be prepared efficiently.

6. Why a history state is useful

The encoded state has schematic form

|X⟩∝Σ_{k=0}^{M}|k⟩|x_k⟩.

Measuring the time register can postselect a chosen time slice.

Alternatively, one can estimate time-averaged observables directly from the history state.

Padding repeated copies of the final-time solution can increase the probability of obtaining the terminal slice without changing its value.

7. Taylor-series ODE algorithms

Berry, Childs, Ostrander and Wang developed a quantum algorithm for sparse linear differential equations using high-order truncated Taylor expansions embedded in a linear system. [1]

High-order approximation reduces the number of time segments needed for a target precision compared with low-order Euler schemes.

The algorithmic complexity depends on sparsity, evolution time, norms/growth factors, condition numbers, precision and oracle costs. Quoting only logarithmic dependence on dimension hides these other resources.

8. Spectral methods

For sufficiently smooth solutions, polynomial spectral approximations can converge much faster than low-order finite differences.

Represent a scalar function on a finite interval using Chebyshev polynomials

x(t)≈Σ_{j=0}^{n} c_j T_j(τ(t))

after mapping time to τ∈[−1,1].

Differential equations enforced at collocation nodes become a structured linear system for coefficients or nodal values.

9. Chebyshev nodes

A common set is

τ_j=cos(jπ/n), j=0,…,n.

Nodes cluster near ±1, reducing Runge-type endpoint instability seen with equally spaced high-degree interpolation.

Childs and Liu developed quantum spectral methods with exponentially improved precision dependence for suitable smooth differential equations. [2]

10. Exponential convergence needs smoothness

Spectral methods can converge exponentially in polynomial degree when the solution is analytic in a sufficiently large complex neighbourhood.

Discontinuities, shocks, nonsmooth coefficients and boundary layers can destroy this advantage.

“Spectral means exponentially accurate” is therefore a regularity-dependent statement, not a universal property of all PDEs.

11. PDE discretisation

Consider Poisson’s equation

−∇²u(x)=f(x)

with boundary conditions.

Discretising space converts the PDE into

A u = f.

A is a sparse stiffness/Laplacian matrix for standard finite differences or finite elements.

At this stage the quantum task becomes a structured quantum linear-system problem.

12. One-dimensional finite-difference Laplacian

For −u”(x)=f(x) on [0,1] with u(0)=u(1)=0, use interior grid spacing h=1/(N+1).

The second derivative approximation gives

−u''(x_j)≈(−u_{j−1}+2u_j−u_{j+1})/h².

Thus

A=(1/h²) tridiag(−1,2,−1).

A is sparse and positive definite.

13. Worked three-point Poisson system

Take N=3 interior points, so h=1/4.

Then

A=16[[2,−1,0],[−1,2,−1],[0,−1,2]].

For f=(1,1,1)^T, solving gives

u=(3/32,1/8,3/32)^T=(0.09375,0.125,0.09375)^T.

A quantum linear-system algorithm would ideally prepare the normalised state proportional to this vector—not print the three entries automatically.

14. Condition number of the Laplacian

For the 1D Dirichlet finite-difference Laplacian, the smallest eigenvalue approaches a constant O(1) while the largest scales O(h⁻²)=O(N²) under physical scaling.

Thus the condition number grows as

κ(A)=O(N²).

This can erase apparent polylogarithmic dimension advantages unless preconditioning reduces κ.

15. Preconditioning

Classically, solve

M^{-1}Au=M^{-1}f

with preconditioner M chosen so M⁻¹A has a smaller condition number.

Quantumly, the challenge is stronger: M⁻¹, A and state preparation must remain efficiently implementable in the oracle/block-encoding model.

A mathematically excellent classical preconditioner may be unusable quantumly if applying it destroys sparsity or requires loading a dense inverse.

16. Boundary conditions

Dirichlet, Neumann, Robin and periodic boundary conditions change both the linear system and its spectrum.

For pure Neumann Poisson problems, the Laplacian has a zero eigenvalue corresponding to adding a constant. The equation is solvable only if f satisfies a compatibility condition, and one must fix a gauge/mean value to obtain uniqueness.

Quantum inversion cannot invert a singular operator without explicitly restricting to a compatible nonsingular subspace or regularising the problem.

17. Heat equation

For

∂_t u=α∇²u,

spatial discretisation gives

du/dt=−αA u

when A is the positive discrete −∇² operator.

This is exactly a large linear ODE. Eigenmodes with larger Laplacian eigenvalues decay faster.

The process is nonunitary, so algorithms use linear-system embeddings, linear-combination methods or dilation/block-encoding techniques rather than direct Hamiltonian simulation alone.

18. Wave equation

The second-order wave equation

∂_t²u=c²∇²u

can be converted to a first-order system by introducing velocity v=∂_tu:

d/dt [u;v]=[[0,I],[c²∇²,0]][u;v].

After spatial discretisation this becomes a structured linear ODE in twice the state dimension.

19. Schrödinger equation is the unitary special case

The time-dependent Schrödinger equation

i∂_tψ=Hψ

has A=−iH.

Because H is Hermitian, A is anti-Hermitian and the norm is preserved.

This case belongs primarily to Hamiltonian simulation rather than the general nonunitary differential-equation machinery here.

20. Nonlinear differential equations

For nonlinear dynamics

dx/dt=f(x,t),

the straightforward linear-system reduction disappears.

Research directions include Carleman linearisation, Koopman embeddings, perturbative expansions and variational quantum simulation.

These methods introduce truncation dimensions and convergence conditions that can grow badly with nonlinearity/time. Claims of generic exponential nonlinear-ODE speedup require particular caution.

21. Carleman linearisation

A polynomial nonlinear ODE can be lifted into an infinite-dimensional linear system over monomials of x.

For scalar

dx/dt=ax+bx²,

differentiate x², x³,… to obtain coupled linear equations among monomials, but each degree couples to higher degrees.

Truncating at finite degree creates an approximation whose error can grow with time or state magnitude.

22. Solution-state output

Suppose the discretised solution vector is u∈ℂ^N.

A quantum algorithm typically produces

|u⟩=u/||u||.

The norm ||u|| may need separate estimation if absolute scale matters.

The state makes quadratic functionals accessible, for example

⟨u|O|u⟩

for efficiently measurable O.

23. Useful compressed observables

Examples include:

• probability mass in a spatial region;
• overlap with a target mode;
• energy/norm estimates;
• integrated flux represented by a sparse observable;
• correlation with another prepared field;
• selected Fourier components.

If the scientific question requires all N grid values, quantum state compression offers little readout advantage.

24. State preparation cost

The right-hand side |f⟩, initial state |x₀⟩ and coefficient/operator oracles must be prepared efficiently.

If f is arbitrary classical data of length N and loading it costs O(N), that input cost can dominate an otherwise polylogarithmic linear-system solve.

Quantum differential-equation advantage is strongest when inputs have succinct structure: analytic functions, sparse expansions, generative circuits or data already produced quantumly.

25. Sparsity and block encoding

Finite-difference differential operators are often sparse.

Modern algorithms may assume oracle access to nonzero matrix entries or a block encoding of A/α inside a larger unitary.

The normalisation α enters complexity. A poor block encoding can erase an apparent theoretical gain.

26. Precision accounting

Total error includes several layers:

• physical-model error;
• spatial discretisation error;
• time discretisation or polynomial truncation error;
• linear-system approximation error;
• state-preparation error;
• observable-estimation shot error.

Reducing the quantum solver error far below the discretisation error wastes resources.

Guide-level error budgets should allocate tolerance across the whole pipeline.

27. Classical numerical baselines

Classical PDE solvers are extremely mature.

Strong baselines include:

• multigrid methods with near-linear complexity for elliptic PDEs;
• FFT/spectral solvers;
• preconditioned Krylov methods;
• adaptive finite elements;
• operator splitting;
• GPU sparse linear algebra.

Comparing a quantum algorithm only with dense Gaussian elimination can create a misleading speedup claim.

28. Where quantum advantage is most plausible

Promising mathematical regimes combine:

• very high-dimensional discretisations;
• sparse/structured well-conditioned operators;
• succinct state preparation;
• smooth solutions enabling high-order convergence;
• compressed output observables;
• many repeated solves where setup cost is amortised.

Breaking several of these assumptions can eliminate the advantage.

29. Common misconception: log N qubits means an N-point PDE solution is available in log N time

The quantum state can encode N amplitudes, but reading all N requires many measurements. Compression helps only when the desired output is also compressed.

30. Common misconception: HHL solves the PDE immediately

One must first choose a stable discretisation, implement boundary conditions, prepare the right-hand side, control condition number, construct the matrix oracle/block encoding and define a useful measurement of the solution state.

31. Common misconception: spectral methods always beat finite differences

Spectral convergence can be exceptional for smooth analytic solutions, but discontinuities, irregular domains and sharp boundary layers can favour local adaptive discretisations.

32. Worked synthesis problem

Discretise −u”=1 on [0,1] with three interior points and zero Dirichlet boundaries.

Step 1: Grid. h=1/4.

Step 2: Matrix. A=16[[2,−1,0],[−1,2,−1],[0,−1,2]].

Step 3: Solve. u=(3/32,1/8,3/32).

Step 4: Normalised quantum state. The squared norm is 2(3/32)²+(1/8)²=17/512. Therefore the amplitude ratios are proportional to 3:4:3 and the normalised state is (3|1⟩+4|2⟩+3|3⟩)/√34 after indexing the three grid points.

Step 5: Readout lesson. Measuring this state gives grid-index probabilities {9/34,16/34,9/34}, not the signed physical field values directly. Recovering absolute u requires separate norm/scale information and, for general signed/complex vectors, phase-sensitive observables.

33. Practice set

1. Write the general linear ODE used here.
2. What is special about A=−iH with H Hermitian?
3. How can time-stepping recurrences be turned into one linear system?
4. What is a history state?
5. Why can spectral methods converge rapidly for analytic solutions?
6. Write the 1D Dirichlet finite-difference Laplacian matrix pattern.
7. Why does Laplacian conditioning worsen with grid refinement?
8. What does a quantum linear-system algorithm usually output?
9. Why is reading every grid value expensive?
10. Name two compressed observables that may preserve quantum advantage.
11. Why can arbitrary classical state loading erase speedup?
12. Why should quantum PDE algorithms be compared with multigrid/Krylov/spectral classical baselines?

Answers

1. dx/dt=A(t)x+b(t).
2. The evolution is unitary Schrödinger-type evolution with norm preservation.
3. Stack all time slices into a block vector and encode each recurrence as a block row.
4. A superposition over a time register and corresponding solution vectors, schematically Σ_k|k⟩|x_k⟩.
5. Polynomial/Chebyshev coefficients decay rapidly when the function extends analytically beyond the real interval.
6. (1/h²)tridiag(−1,2,−1).
7. The largest discrete Laplacian eigenvalue scales as h⁻² while the smallest physical Dirichlet eigenvalue remains O(1).
8. A normalised state proportional to the solution vector, not a classical list of entries.
9. Extracting N independent classical values requires Ω(N) information/measurements in general.
10. Examples: overlap with a target mode, probability in a region, energy, selected Fourier coefficient or sparse quadratic observable.
11. Preparing an arbitrary N-dimensional input from classical memory can itself require O(N) work.
12. Those classical solvers exploit the same PDE structure and can be near-linear or better in practical regimes; dense elimination is not the relevant benchmark.

Sources and further study

[1] Dominic W. Berry, Andrew M. Childs, Aaron Ostrander and Guoming Wang, Quantum Algorithm for Linear Differential Equations with Exponentially Improved Dependence on Precision, Communications in Mathematical Physics 356, 1057–1081 (2017). A foundational high-order quantum linear-ODE algorithm.

[2] Andrew M. Childs and Jin-Peng Liu, Quantum spectral methods for differential equations, Communications in Mathematical Physics 375, 1427–1457 (2020). Quantum Chebyshev/spectral methods for smooth differential equations.

[3] Andrew M. Childs, Robin Kothari and Rolando D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. A modern quantum linear-system primitive relevant after PDE discretisation.

[4] Dong An, Di Fang and Lin Lin, Time-dependent unbounded Hamiltonian simulation with vector norm scaling. Related numerical-analysis ideas for structured differential evolution and discretisation-dependent quantum simulation.

Continue through Quantum Mathematics

Guide 61: Quantum Approximate Optimisation Algorithm, MaxCut, Cost Hamiltonians and Alternating Operators develops hybrid optimisation. Guide 62: Quantum Machine Learning, Feature Maps, Quantum Kernels and Variational Classifiers develops learning models. Guide 64: Quantum Simulation of Lattice Gauge Theories, Gauge Constraints, Wilson Loops and Kogut–Susskind Hamiltonians develops constrained field simulation.

Return to the BTT Mathematics Learning Hub.