量子 PDE 求解详解：复杂度演进与最新突破

量子偏微分方程（partial differential equation, PDE）求解是量子计算最具变革性潜力的应用之一。PDE 描述了物理世界的基本规律——流体力学、电磁学、量子场论、金融定价——而经典求解的计算代价随空间维度指数增长。本文系统梳理量子 PDE 求解的复杂度演进史，从 2014 年前后的初步方案到 2024 年的最优结果。记号与本系列 HHL 教程保持一致：条件数（condition number）记 \(\kappa\)，酉演化采用 \(e^{iHt}\) 约定，\(n=\lceil\log_2 N\rceil\) 为编码 \(N\) 维离散系统所需的量子比特数；为避免与条件数混淆，热方程的扩散系数记 \(\alpha\)（而非部分文献的 \(\kappa\)）。

问题定义

PDE 的一般形式

我们考虑发展型方程

\[\frac{\partial u}{\partial t} = \mathcal{L}[u] + f(\mathbf{x}, t), \qquad u(\mathbf{x}, 0) = u_0(\mathbf{x}),\]

其中 \(\mathcal{L}\) 是空间微分算符（如 Laplace 算符 \(\Delta\)、散度 \(\nabla\cdot\)、旋度 \(\nabla\times\)），\(f\) 为源项；椭圆型方程（如 Poisson 方程）没有时间变量，直接化为边值问题。主要类型如下表所示。

类型

代表方程

物理意义

抛物型

热方程 \(u_t = \alpha\Delta u\)

扩散、热传导

双曲型

波方程 \(u_{tt} = c^2\Delta u\)

声波、电磁波

椭圆型

Poisson 方程 \(\Delta u = f\)

静电场、引力

混合型

Navier-Stokes 方程

流体力学

经典复杂度

对 \(d\) 维空间、每维 \(N\) 个格点（共 \(N^d\) 个自由度）、时间精度 \(\varepsilon\) 的问题，经典方法的代价为：\(p\) 阶时间格式需要 \(O(T\,\varepsilon^{-1/p})\) 步（局部误差 \(O(h^{p+1})\) 累积为全局误差 \(O(h^p)\)），每步代价由空间离散算子的应用或求解决定，如下表。

方法

每步复杂度

步数

总复杂度

有限差分（显式）

\(O(N^d)\)

\(O(T\varepsilon^{-1/p})\)

\(O(N^d\,T\,\varepsilon^{-1/p})\)

有限元（隐式稀疏求解）

\(O(N^d)\)

\(O(T\varepsilon^{-1/p})\)

\(O(N^d\,T\,\varepsilon^{-1/p})\)

谱方法（FFT）

\(O(N^d\log N)\)

\(O(T\varepsilon^{-1/p})\)

\(O(N^d\log N\,T\,\varepsilon^{-1/p})\)

维度灾难（curse of dimensionality）体现在自由度数 \(N^d\) 上：\(d=3\) 时 \(N^3\) 已十分巨大，\(d>10\) 时经典方法完全不可行。量子编码把空间自由度压缩到 \(n = d\log_2 N\) 个量子比特，这是量子 PDE 求解的根本动机；但门数是否也随之改善，取决于下文讨论的条件数与输出方式。

量子 PDE 求解的复杂度演进

第一代：量子模拟直接方法（2014-2019）

Berry 的高阶 ODE 算法 [1] 与 Lloyd 等人的早期方案确立了基本范式：对线性 PDE，先做空间离散化得到线性 ODE 系统 \(\dot{\mathbf{u}} = A\mathbf{u}\)，其中 \(A\) 由空间差分矩阵构成，再用 Trotter 分解逐步推进 \(e^{A\Delta t}\) 类的演化。

对 Hermitian 部分（可写成 \(i\dot{\mathbf u} = H\mathbf u\) 的系统，如波方程的适当变形，见下文），问题直接归结为稀疏哈密顿量模拟。\(A\) 的稀疏度为 \(L=O(d)\)（每个格点只与 \(2d\) 个邻居耦合），编码需要 \(n = \lceil d\log_2 N\rceil\) 个量子比特。用一阶 Trotter 分解实现时间 \(T\) 的演化，门数为

\[O\!\left(L\cdot n\cdot\frac{T}{\Delta t}\cdot\varepsilon^{-1/(2p)}\right),\]

其中 \(p\) 阶 Suzuki 公式把每步误差压到 \(O(\Delta t^{p+1})\)。这一代方法的关键限制是：空间复杂度对 \(d\) 指数改善，但时间复杂度中 \(\Delta t\) 受稳定性与精度双重约束，整体对 \(1/\varepsilon\) 只有多项式依赖（早期约为 \(O(T^2/\varepsilon)\) 量级），且只适用于能直接写成酉演化的系统。

第二代：高阶方法（2020-2022）

Costa–An–Jordan–Liu 的波方程量子算法与同期工作把时间积分系统地高阶化：用高阶 Magnus 展开处理时变部分、用高阶 Suzuki 公式或截断 Taylor 级数降低分解误差，将复杂度改进为

\[O\!\left(L\cdot n\cdot\big(T^{1+o(1)} + \log(1/\varepsilon)\big)\right),\]

即对模拟时间近线性、对精度对数。这一步把"对 \(1/\varepsilon\) 多项式依赖"的瓶颈彻底移除，代价是对 \(T\) 的 \(T^{1+o(1)}\) 超线性因子。

第三代：Schrödingerization（2022-2024）

Jin–Liu–Yu 提出的 Schrödingerization 框架 [2] 是突破性的：它不要求算符 Hermitian、不依赖 Trotter 分解，适用于任意线性 PDE。我们给出其构造的完整推导，这与本系列 ODE 教程中的推导一致。

空间离散化后，线性 PDE 化为 \(\dot{\mathbf w} = A\mathbf w\)（\(A\) 为 \(N^d\times N^d\) 矩阵，可以是非 Hermitian 的）。把 \(A\) 分解为 Hermitian 与反 Hermitian 部分：

\[H_R := \frac{A+A^\dagger}{2},\qquad H_I := \frac{A-A^\dagger}{2i},\qquad A = H_R + iH_I,\]

引入辅助变量 \(\eta\in\mathbb{R}\)（动量算符 \(P=-i\partial_\eta\)，自伴）与初值 \(\Phi(0,\eta) = \mathbf w_0\,e^{-\eta}\Theta(\eta)\)，则增广哈密顿量取为

\[\mathcal{K} \;=\; \frac{A+A^\dagger}{2}\otimes P \;+\; \frac{i}{2}\bigl(A - A^\dagger\bigr)\otimes I .\]

\(\mathcal K\) 是 Hermitian 的（两个直积项的每个因子都 Hermitian），因此 \(e^{-i\mathcal K t}\) 是酉演化，可用任意标准模拟方法实现。在 \(\eta\) 表示下薛定谔方程化为输运方程 \(\partial_t\Phi = -H_R\partial_\eta\Phi + iH_I\Phi\)（速度 \(-H_R\) 叠加相位旋转）；当 \([H_R,H_I]=0\) 时沿特征线求解给出 \(\Phi(t,\eta) = e^{iH_It}\,\Phi_0(\eta-H_Rt)\)，于是恢复量 \(e^{\eta}\Phi(t,\eta) = e^{iH_It}e^{H_Rt}\mathbf w_0 = e^{At}\mathbf w_0\) 精确成立；非交换情形由生成元幂级数的逐项归纳同样给出精确恢复（严格表述见 [2]）。

我们用热方程逐模式验证构造。\(A = \alpha\Delta\) 是 Hermitian 的（\(H_R = \alpha\Delta\)、\(H_I=0\)），故 \(\mathcal K = \alpha\Delta\otimes P\)，\(\eta\) 表示下方程为 \(\partial_t\Phi = -\alpha\Delta\,\partial_\eta\Phi\)。对特征模式 \(e^{ikx}\)（\(\Delta e^{ikx} = -k^2 e^{ikx}\)，见下文谱方法一节的推导），输运速度为 \(\alpha k^2\)，解为

\[\Phi_k(t,\eta) = \Phi_k(0,\,\eta+\alpha k^2 t)\quad\Longrightarrow\quad e^{\eta}\Phi_k(t,\eta) = e^{-\alpha k^2 t}\,\hat u_k(0),\]

恰好给出热方程的衰减因子 \(e^{-\alpha k^2 t}\)。整体复杂度为

\[O\!\left(\mathrm{poly}(n)\cdot\big(\|A\|\,t + \log(1/\varepsilon)\big)\right),\]

突破点在于：对局域差分算符 \(\|A\| = O(d/h^2)\) 与格点总数 \(N^d\) 无关（只与尺度相关），因此总代价是 \(\mathrm{poly}(n)\cdot T\) 而非 \(\mathrm{poly}(N^d)\cdot T\)。后选择概率正比于解的范数平方与初值范数平方之比；解指数衰减（长期热演化）时需要振幅放大，这一代价与问题的物理收缩率相关，无法由算法侧消除。

第四代：黑盒 PDE 求解（2023-2024）

An–Liu–Lin 的 LCHS（linear combination of Hamiltonian simulation）框架 [3] 针对非齐次线性系统 \(\partial_t u = \mathcal Lu + f\) 给出了近最优的态制备方案。出发点仍是 Duhamel 公式（推导见本系列 ODE 教程情况三，由常数变易法得到）：

\[u(T) = e^{\mathcal LT}u_0 + \int_0^T e^{\mathcal L(T-s)}\,f(s)\,ds .\]

量子实现分为三步：第一步，把积分离散化为 \(M\) 个时间片上的求积 \(\sum_m w_m\,e^{\mathcal L(T-s_m)}f(s_m)\)；第二步，对每个 \(m\)，用 LCHS 把非酉块 \(e^{\mathcal L(T-s_m)}\) 分解为一族可模拟酉演化 \(e^{-i\widetilde H(\alpha)t}\) 的加权积分（\(\widetilde H(\alpha)\) 由 \(\mathcal L\) 的 Hermitian 与反 Hermitian 部分组装，积分权重可解析给出）；第三步，用受控操作把所有时间片叠加成历史态，末时刻分量即解。相比早期方法，总代价从 \(O(T^2/\varepsilon)\) 降至

\[O\!\left(T\cdot\mathrm{polylog}(T/\varepsilon)\right),\]

这是已知的近最优时间依赖。

主要 PDE 类型的具体处理

热方程（抛物型）

热方程为 \(\partial_t u = \alpha\Delta u\)。我们先在傅里叶基下对角化（完整推导见下文谱方法一节）：周期域上 \(u(x,t) = \sum_k \hat u_k(t)e^{ikx}\)，逐项代入得 \(\dot{\hat u}_k = -\alpha k^2\hat u_k\)，解为

\[\hat u_k(t) = e^{-\alpha k^2 t}\,\hat u_k(0),\]

即每个模式的模长严格衰减——传播子 \(e^{\alpha\Delta t}\) 是收缩半群而非酉算符。这正是量子模拟的挑战所在：\(e^{i(-i\alpha\Delta)t} = e^{\alpha\Delta t}\) 需要"虚时间演化"，不是酉操作。可行的量子路线有三条。其一，Schrödingerization：按第三节的构造 \(\mathcal K = \alpha\Delta_h\otimes P\)（\(\Delta_h\) 为离散 Laplacian），通过辅助维度与后选择实现收缩；解衰减越多，后选择概率越低，需配合重缩放（见 ODE 教程情况四）。其二，QSVT：对块编码的 \(\alpha\Delta_h\) 施加多项式变换 \(P(x)\approx e^{\alpha x t}\)，直接以多项式代价（随 \(\alpha\|\Delta_h\|t\) 与 \(\log(1/\varepsilon)\) 多项式增长）实现收缩半群。其三，隐式时间步进（如 Crank–Nicolson）把每步化为线性方程组 \((I-\tfrac{\alpha\Delta t}{2}\Delta_h)u^{(m+1)} = (I+\tfrac{\alpha\Delta t}{2}\Delta_h)u^{(m)}\)，再调用线性求解器，代价由该系统的条件数（\(O(N^2)\)，见 Poisson 一节）决定。整体复杂度 \(O(\mathrm{poly}(n)\cdot T)\)。

波方程（双曲型）

波方程为 \(\partial_{tt}u = c^2\Delta u\)。降为一阶系统的朴素做法是令 \(v = \partial_t u\)，得到

\[\begin{split}\partial_t\begin{pmatrix}u\\ v\end{pmatrix} = \begin{pmatrix}0 & I\\ c^2\Delta & 0\end{pmatrix}\begin{pmatrix}u\\ v\end{pmatrix}.\end{split}\]

需要指出：这个生成矩阵并不是 Hermitian 的——按分块伴随的计算，\(\begin{pmatrix}0&I\\ c^2\Delta&0\end{pmatrix}^\dagger = \begin{pmatrix}0 & c^2\Delta^\dagger\\ I & 0\end{pmatrix}\)，它与原矩阵相等当且仅当 \(c^2\Delta = I\)，一般并不成立。波方程守恒的是能量范数而非标准范数，正确的 Hermitian 化必须引入平方根算符。记 \(B := c\,(-\Delta)^{1/2}\)（\(\Delta\) 的谱非正，故 \(B\) 是 Hermitian 正算符），做变量替换 \(\psi := (Bu,\ v)^{\mathsf T}\)，我们逐行求导：

\[\partial_t(Bu) = Bv,\qquad \partial_t v = c^2\Delta u = -c^2(-\Delta)u = -B^2u = -B(Bu),\]

即 \(\partial_t\psi = M\psi\)，\(M = \begin{pmatrix}0 & B\\ -B & 0\end{pmatrix}\)。由 \(B = B^\dagger\) 直接验算

\[\begin{split}M^\dagger = \begin{pmatrix}0 & (-B)^\dagger\\ B^\dagger & 0\end{pmatrix} = \begin{pmatrix}0 & -B\\ B & 0\end{pmatrix} = -M,\end{split}\]

即 \(M\) 反 Hermitian，因此 \(e^{Mt}\) 是酉演化，\(\|\psi(t)\|^2 = \|Bu(t)\|^2 + \|v(t)\|^2\) 守恒（这正是能量）。写成薛定谔形式 \(i\partial_t\psi = H\psi\)，哈密顿量为 \(H = iM = \begin{pmatrix}0 & iB\\ -iB & 0\end{pmatrix}\)，它是 Hermitian 的（验证：\(H^\dagger = -iM^\dagger = -i(-M) = iM = H\)）。在傅里叶基下 \(B\) 对角（\(Be^{ikx} = c|k|e^{ikx}\)），每个模式以频率 \(ck\) 谐振；单模式核对：\(u = \cos(ckt)\cos(kx)\) 给出 \(Bu = ck\cos(ckt)\cos(kx)\)、\(v = -ck\sin(ckt)\cos(kx)\)，故 \(\|\psi\|^2 = c^2k^2(\cos^2+\sin^2) = c^2k^2\) 为常数。离散化后 \(B = c(-\Delta_h)^{1/2}\) 一般是稠密的，但恰好可用谱方法（下节）经 QFT 对角化实现，或者保留非 Hermitian 的一阶形式改用 Schrödingerization/LCHS。复杂度 \(O(\mathrm{poly}(n)\cdot T)\)。

Poisson 方程（椭圆型）

Poisson 方程 \(\Delta u = f\)（等价地 \(-\Delta u = -f\)，我们解正定的 \(-\Delta u = g\) 形式）没有时间变量，有限差分离散化直接得到线性方程组。我们把整条推导链逐步展开。

第一步（中心差分格式及其误差）。在 \([0,1]\) 上取网格 \(x_i = ih\)（\(i=0,\dots,N+1\)，\(h = 1/(N+1)\)），内部节点为未知量。由 Taylor 展开

\[u(x\pm h) = u(x) \pm h u'(x) + \frac{h^2}{2}u''(x) \pm \frac{h^3}{6}u'''(x) + \frac{h^4}{24}u^{(4)}(x) \pm O(h^5),\]

两式相加时奇次项相消，得 \(u(x+h)+u(x-h) = 2u + h^2u'' + \frac{h^4}{12}u^{(4)} + O(h^6)\)，移项即

\[u''(x_i) = \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} - \frac{h^2}{12}\,u^{(4)}(x_i) + O(h^4),\]

因此离散算符 \(L_hu_i := \frac{2u_i - u_{i-1}-u_{i+1}}{h^2}\) 以 \(O(h^2)\) 精度逼近 \(-u''\)。

第二步（边界条件进入矩阵）。齐次 Dirichlet 条件 \(u(0)=u(1)=0\) 意味着 \(u_0 = u_{N+1} = 0\)，未知量是 \(u_1,\dots,u_N\)，系统为 \(\mathbf A\mathbf u = \mathbf g\)，其中

\[\begin{split}\mathbf A = \frac{1}{h^2}\begin{pmatrix} 2 & -1 & & \\ -1 & 2 & -1 & \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 \end{pmatrix}.\end{split}\]

非齐次 Dirichlet 条件 \(u(0)=g_0\)、\(u(1)=g_1\) 时，边界值是已知量，从矩阵行移到右端：第 \(1\) 行的方程 \(\frac{2u_1-u_0-u_2}{h^2} = g(x_1)\) 化为 \(\frac{2u_1-u_2}{h^2} = g(x_1) + \frac{g_0}{h^2}\)（第 \(N\) 行同理加 \(g_1/h^2\)）。也就是说，Dirichlet 条件不改变矩阵、只改变右端。Neumann 条件 \(u'(0)=0\) 则改变矩阵：引入镜像点 \(u_{-1}\) 并用中心差分 \(\frac{u_1-u_{-1}}{2h}=0\) 得 \(u_{-1}=u_1\)，把它代入第 \(0\) 行的离散方程 \(\frac{2u_0-u_{-1}-u_1}{h^2} = g(x_0)\)，得 \(\frac{2u_0-2u_1}{h^2} = g(x_0)\)——首行变为 \((2,-2)/h^2\)（此朴素形式不对称，对称化需改用半格点格式 \((1,-1)/h^2\)）。原则上：边界条件通过修改边界行（Neumann）或右端（Dirichlet）进入离散系统。

第三步（谱与条件数）。我们推导 \(\mathbf A\) 的全部特征值。对候选特征向量 \(v_i = \sin(i\theta)\)（\(i=1,\dots,N\)），用和差恒等式 \(\sin((i\pm1)\theta) = \sin(i\theta)\cos\theta \pm \cos(i\theta)\sin\theta\) 计算

\[\frac{-v_{i-1}+2v_i-v_{i+1}}{h^2} = \frac{2-2\cos\theta}{h^2}\,\sin(i\theta),\]

故 \(v\) 是特征向量、特征值为 \(\lambda(\theta) = \frac{2-2\cos\theta}{h^2} = \frac{4\sin^2(\theta/2)}{h^2}\)。边界条件要求 \(v_0 = \sin 0 = 0\)（自动满足）且 \(v_{N+1} = \sin((N+1)\theta) = 0\)，后者给出 \(\theta = \frac{k\pi}{N+1}\)，\(k=1,\dots,N\)。于是

\[\lambda_k = \frac{4}{h^2}\sin^2\!\frac{k\pi}{2(N+1)},\qquad k = 1,\dots,N .\]

两端估计：最小特征值用 \(\sin x\approx x\)（\(x\to0\)），并注意 \(h = 1/(N+1)\)：

\[\lambda_{\min} = \lambda_1 \approx \frac{4}{h^2}\cdot\frac{\pi^2}{4(N+1)^2} = \frac{\pi^2}{\big(h(N+1)\big)^2} = \pi^2,\]

与连续算符 \(-\frac{d^2}{dx^2}\) 在 \((0,1)\) 上的第一特征值 \(\pi^2\) 一致，这是对推导的自洽核对。最大特征值用恒等式 \(\frac{N\pi}{2(N+1)} = \frac\pi2 - \frac{\pi}{2(N+1)}\) 与 \(\sin(\frac\pi2-x) = \cos x\)：

\[\lambda_{\max} = \lambda_N = \frac{4}{h^2}\cos^2\frac{\pi}{2(N+1)} \;\approx\; \frac{4}{h^2} = 4(N+1)^2 .\]

因此条件数

\[\kappa(\mathbf A) = \frac{\lambda_{\max}}{\lambda_{\min}} \approx \frac{4(N+1)^2}{\pi^2} = O(N^2).\]

第四步（高维与量子求解）。\(d\) 维情形的离散算符是 Kronecker 和 \(\mathbf A_d = \sum_{j=1}^d\, I^{\otimes(j-1)}\otimes\mathbf A\otimes I^{\otimes(d-j)}\)，其特征向量为一维正弦函数的张量积、特征值为 \(d\) 个一维特征值之和，故 \(\kappa\) 仍为 \(O(N^2)\)；矩阵稀疏度为 \(2d+1\)（中心加每维两个邻居），未知量 \(N^d\) 个，编码需 \(n = \lceil d\log_2N\rceil\) 个量子比特。用量子线性求解器（HHL/QSVT）求 \(\mathbf A\mathbf u = \mathbf g\) 的门数为

\[O\big(\kappa\cdot\mathrm{poly}(n)\cdot\log(1/\varepsilon)\big) = O\big(N^2\,\mathrm{polylog}(N^d)\,\log(1/\varepsilon)\big).\]

与经典比较：稀疏 Cholesky/共轭梯度为 \(O(N^d\sqrt\kappa\log(1/\varepsilon))\) 量级，因此量子优势要求 \(N^d \gg N^2\mathrm{polylog}\)，即 \(d\ge3\) 才有明确的渐近优势（\(d=1,2\) 时条件数抵消对数编码收益）；此外输出必须是泛函而非全场。缓解 \(\kappa\) 的途径包括预条件化（preconditioning，Clader–Jacobs–Sprouse 的早期分析）与下节的谱方法。

谱方法：傅里叶基下的对角化

对周期边界问题，谱方法（spectral method）用傅里叶基代替差分，二阶导数算符在该基下严格对角。我们逐项推导。取规范正交基 \(\phi_k(x) = \frac{1}{\sqrt{2\pi}}e^{ikx}\)（\(k\in\mathbb Z\)，域 \([0,2\pi)\)）。正交性由直接积分验证：

\[\begin{split}\langle\phi_j,\phi_k\rangle = \frac{1}{2\pi}\int_0^{2\pi}e^{i(k-j)x}\,dx = \begin{cases}1, & k=j,\\[2pt] \dfrac{e^{2\pi i(k-j)}-1}{2\pi i(k-j)} = 0, & k\ne j,\end{cases}\end{split}\]

其中用了 \(e^{2\pi i m} = 1\)（\(m\) 为整数）。对 \(\phi_k\) 求两次导数，每次求导产生因子 \(ik\)：

\[\phi_k'' = \frac{(ik)^2}{\sqrt{2\pi}}\,e^{ikx} = -k^2\,\phi_k .\]

因此微分算符 \(\frac{d^2}{dx^2}\) 在傅里叶基下是对角的，本征值为 \(-k^2\)（同样地 \(\Delta\) 在多维情形的本征值为 \(-|\mathbf k|^2\)）。把 \(u(x) = \sum_k\hat u_k\phi_k\) 代入线性 PDE，由正交性可逐项比对系数：对热方程得 \(\dot{\hat u}_k = -\alpha k^2\hat u_k\)（上节已用）；对 Poisson 方程 \(-\Delta u = g\) 得

\[k^2\,\hat u_k = \hat g_k\qquad\Longrightarrow\qquad \hat u_k = \frac{\hat g_k}{k^2}\quad(k\ne0).\]

\(k=0\) 模式必须单独讨论：它对应右端的积分 \(\hat g_0 = \frac{1}{\sqrt{2\pi}}\int_0^{2\pi}g\,dx\)，方程 \(0\cdot\hat u_0 = \hat g_0\) 有解当且仅当 \(\int g = 0\)（相容性条件，物理上即总电荷为零），此时 \(\hat u_0\) 任意（解可差常数）。离散版本：格点 \(x_j = 2\pi j/N\) 上的差分算符是循环矩阵，被离散傅里叶变换 \(F\) 对角化为 \(D^2 = F^\dagger\,\mathrm{diag}(-\tilde k^2)\,F\)（\(\tilde k\in\{-N/2,\dots,N/2-1\}\) 为折叠频率）。量子实现恰好匹配：QFT 以 \(O(n^2)\) 门实现 \(F\)，对角因子 \(\mathrm{diag}(1/k^2)\) 或 \(\mathrm{diag}(e^{-\alpha k^2t})\) 是对角预言机（每行一个非零元），逆 QFT 返回物理空间。整条链路不含条件数因子 \(\kappa\)，代价由态制备与 QFT 决定；其代价是仅适用于周期边界，且谱精度要求解光滑。这与有限差分形成互补：差分方法通用但受 \(\kappa = O(N^2)\) 制约，谱方法对角化彻底但受周期性与光滑性制约。

对流方程（一阶双曲型）

对流方程为 \(\partial_tu + c\,\partial_xu = 0\)。周期边界下用中心差分离散一阶导数：\(D_x := \frac{1}{2h}(S_+-S_-)\)，其中 \(S_\pm\) 为平移算符。由 \(S_+^\dagger = S_-\) 直接验证

\[D_x^\dagger = \frac{1}{2h}(S_+^\dagger - S_-^\dagger) = \frac{1}{2h}(S_- - S_+) = -D_x,\]

即 \(D_x\) 反 Hermitian，其特征值为纯虚数 \(\pm\frac{i}{h}\sin(kh)\)（对波数 \(k\)，由 \(S_\pm e^{ikx} = e^{\pm ikh}e^{ikx}\) 直接算出）。因此半离散系统 \(\dot{\mathbf u} = A\mathbf u\)（\(A = -cD_x\)）的生成元是反 Hermitian 的，传播子 \(e^{At} = e^{-iHt}\)（\(H := iA = -icD_x\) Hermitian）是酉演化，可用标准哈密顿量模拟直接实现，复杂度 \(O(\mathrm{poly}(n)\cdot T)\)，且中心差分的纯虚谱意味着无数值耗散（中性稳定）。需要说明的是，稳定化的迎风格式（upwind）会破坏反对称性，届时应改用 Schrödingerization 或 LCHS 处理非酉传播子。

复杂度演进总结

下表汇总各代方法的标度（\(n = d\log_2N\)，\(T\) 为模拟时间，\(K\) 为 Carleman 截断阶数）。

年份

方法

PDE 类型

总复杂度

关键突破

2014

Berry（高阶方法）[1]

线性

\(O(\mathrm{poly}(n)\,T^2/\varepsilon)\)

对数空间编码

2019

Lloyd et al.

线性（Hermitian）

\(O(\mathrm{poly}(n)\,T)\)

直接量子模拟

2020

Berry et al.（高阶积分）

线性

\(O(\mathrm{poly}(n)\,T^{1+o(1)})\)

高阶时间积分

2022

Schrödingerization [2]

任意线性

\(O(\mathrm{poly}(n)\,T)\)

打破 Hermitian 限制

2023

An–Liu–Lin（LCHS）[3]

非齐次线性

\(O(\mathrm{poly}(n)\,T\,\mathrm{polylog}(T/\varepsilon))\)

近最优态制备

2023

Costa et al.（Carleman）[4]

非线性（二次）

\(O(\mathrm{poly}(n,K)\,T)\)

非线性系统

2024

Jin et al.（改进 Schrödingerization）[2]

一般线性 PDE

\(O(\mathrm{poly}(n)\,T)\)

统一框架

与经典方法的对比

量子优势明确的场景有三类。其一，高维 PDE（\(d\ge3\)）：量子编码需 \(n = d\log_2N\) 个量子比特对比经典的 \(N^d\) 存储，且门数为 \(\mathrm{poly}(n)\cdot T\) 对比经典的 \(N^d\,T\,\varepsilon^{-1/p}\)（对适定的演化问题）；对椭圆型问题还需 \(\kappa = O(N^2)\) 与 \(N^d\) 的权衡（见 Poisson 一节）。其二，长时间模拟：量子方法对 \(T\) 近线性且对精度对数。其三，多查询场景：同一 PDE 在多个初值或多个右端下的批量求解，可共享离散化与模拟电路。

量子优势不明确的场景同样有三类。其一，低维 PDE（\(d=1,2\)）：经典方法已经高效，条件数与读出开销抵消量子收益。其二，强非线性 PDE：Carleman 截断阶数 \(K\) 与条件数可能爆炸。其三，输出问题：完整读出 \(u(\mathbf x,T)\) 需要 \(O(N^d)\) 次测量，量子优势只保留在泛函估计（如 \(\langle u|M|u\rangle\)）与采样场景。

开放问题

  1. 非线性 PDE 的最优量子算法：Carleman 截断阶数与条件数的普适控制仍是核心难题，耗散占优之外的 regime 缺乏多项式保证。

  2. 守恒律 PDE：激波与间断解缺乏光滑性，谱方法失效、差分方法的稳定格式（迎风、通量限制）又是非酉的，量子处理方式尚不明确。

  3. 多尺度 PDE：含快慢尺度的问题（如大气模型）需要自适应网格，量子编码的灵活性不足。

  4. 输出效率：如何以少于 \(O(N^d)\) 的代价提取宏观物理量（平均速度、最大涡量、能谱）仍是瓶颈，现有方案限于特定泛函形式。

  5. 实际验证：目前量子 PDE 求解的实验验证限于 \(O(10)\) 量子比特的小规模演示，与实用规模之间仍有数量级差距。

总结

量子 PDE 求解的复杂度从 2014 年的 \(O(\mathrm{poly}(n)\,T^2/\varepsilon)\) 演进到 2024 年的 \(O(\mathrm{poly}(n)\,T)\)。Schrödingerization 框架的引入是最重要的突破——它以增广哈密顿量 \(\mathcal K = \frac{A+A^\dagger}{2}\otimes P + \frac{i}{2}(A-A^\dagger)\otimes I\) 把任意线性演化嵌入酉动力学，打破了"只能模拟 Hermitian 系统"的限制；LCHS 进一步把非齐次项的代价做到近最优。离散化层面，有限差分给出稀疏但条件数 \(O(N^2)\) 的系统，谱方法给出严格对角化但限于周期光滑问题，二者构成量子 PDE 工具箱的两端。非线性问题的 Carleman 线性化在 2023 年取得实质改进，但其普适性与输出效率仍是开放难题；量子 PDE 求解的最终实用化取决于这两点的进展。

参考文献：

  1. Berry, D. W. (2014). High-order quantum algorithm for solving linear differential equations. Journal of Physics A, 47, 105301.

  2. Jin, S., Liu, N., & Yu, Y. (2022). Quantum simulation of partial differential equations via Schrödingerization. arXiv:2212.13969.

  3. An, D., Liu, J., & Lin, L. (2023). Linear combination of Hamiltonian simulation for nonunitary dynamics with optimal state preparation cost. Physical Review Letters, 131, 150601.

  4. Costa, P. C. S., Jordan, S. P., & Ostrander, A. (2023). Improved quantum algorithms for linear and nonlinear differential equations. Quantum, 7, 913.