# 量子 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$，得到

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

需要指出：这个生成矩阵并不是 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$ 直接验算

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

即 $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$，其中

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

非齐次 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)$）。正交性由直接积分验证：

$$\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}$$

其中用了 $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 处理非酉传播子。

### Navier-Stokes 方程

Navier-Stokes 方程为

$$\partial_t\mathbf u + (\mathbf u\cdot\nabla)\mathbf u = -\nabla p + \nu\Delta\mathbf u + \mathbf f,\qquad \nabla\cdot\mathbf u = 0 .$$

困难集中在两点：非线性对流项 $(\mathbf u\cdot\nabla)\mathbf u$，以及压力–速度耦合（压力由不可压缩约束隐式决定）。现有量子路线有三类。其一，Carleman 线性化加 Schrödingerization：把二次非线性展开为无穷维线性系统后截断，再用第三节的框架模拟（复杂度 $O(\mathrm{poly}(n,K)\cdot T)$，$K$ 为截断阶数，其可控制性取决于耗散 $\nu$ 是否压制对流）。其二，迭代线性化：每个时间片内冻结非线性项求解线性 PDE，经典外环更新，类似 ODE 教程的方法二。其三，格子 Boltzmann 方法（lattice Boltzmann）的量子化：把流体动力学改写为格点上的松弛–迁移线性动力学加非线性碰撞项，后者在低马赫数下可截断为低次多项式，再行线性化。

## 复杂度演进总结

下表汇总各代方法的标度（$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.
---

> 返回目录：[量子计算算法教程系列](https://chenzhaoyun.com/index.php/archives/54/)