量子 ODE 求解

量子常微分方程（ordinary differential equation, ODE）求解是量子科学计算的核心方向之一。它的基本格局可以概括为一句话：线性 ODE 已经有一整套成熟的量子算法——哈密顿量模拟、Schrödingerization、以及把时间离散化成线性方程组后调用 HHL/QSVT——而非线性 ODE 必须先线性化（最常用的是 Carleman 线性化），再回到线性工具箱。本词条按”问题定义—线性方法—非线性方法—复杂度对比—具体例子”的顺序展开，对每一条推导链都逐步给出所用恒等式、代入过程与化简结果；记号与本百科 HHL 教程保持一致（条件数记 $\kappa$，酉演化采用 $e^{iHt}$ 约定，$n=\lceil {\log }_{2}N\rceil$ 为量子比特数）。

问题定义

一般 ODE

我们考虑初值问题（initial value problem）

$\frac{d\mathbf{u}}{dt}=f(\mathbf{u},t),\mathbf{u}(0)={\mathbf{u}}_0\in {\mathbb{R}}^d,$

其中 $f:{\mathbb{R}}^d\times [0,T]\to {\mathbb{R}}^d$ 关于两个变量均 Lipschitz 连续，此时解局部存在且唯一。按照右端对 $\mathbf{u}$ 的依赖方式，我们区分两类：

• 线性 ODE：$\dot{\mathbf{u}}=A(t)\mathbf{u}+\mathbf{b}(t)$，其中 $A(t)\in {\mathbb{C}}^{d\times d}$，$\mathbf{b}(t)\in {\mathbb{C}}^d$；
• 非线性 ODE：$f(\mathbf{u},t)$ 含 $\mathbf{u}$ 的二次或更高次项，例如 Lotka–Volterra 方程。

量子算法的输出约定与 HHL 算法相同：我们不要求读出全部 $d$ 个分量（完整读出需要 $O(d)$ 次测量，会抵消量子加速），而是制备与解成比例的量子态 $|u(T)⟩\propto \mathbf{u}(T)$，或估计 $⟨u(T)|M|u(T)⟩$ 一类的全局泛函。

经典方法

我们先回顾经典格式的精度规律，因为它们决定了与量子算法对比的基线。以 $p$ 阶单步法为例，局部截断误差为 $O(h^{p+1})$，把 $[0,T]$ 分成 $N_{\mathrm{step}}=T/h$ 步后误差累积为全局误差 $O(h^p)$，因此达到精度 $\epsilon$ 需要

$N_{\mathrm{step}}=O\left(T{\epsilon }^{-1/p}\right)$

步。对显式与隐式 Euler（$p=1$）这给出 $O({\epsilon }^{-1})$ 步；对四阶 Runge–Kutta（$p=4$）给出 $O({\epsilon }^{-1/4})$ 步；Adams 类多步法同理得到 $O({\epsilon }^{-1/p})$。每步的主要代价是对一般稠密右端 $f$ 的 $O(d)$ 次分量求值，因此总代价为 $O(dT{\epsilon }^{-1/p})$。信息基复杂度（information-based complexity）的经典结果表明，对仅假设 Lipschitz 连续的一般右端，达到全局精度 $\epsilon$ 需要 $\Theta (1/\epsilon )$ 次右端求值，因此一阶格式的步数已是阶数最优；更高的阶数需要相应的高阶光滑性假设。

线性 ODE 的量子求解

情况一：常系数线性 ODE

我们先看最简单的 $\dot{\mathbf{u}}=A\mathbf{u}$，其中 $A$ 为常数矩阵。直接对矩阵指数求导（用级数定义逐项验证，$\frac{d}{dt}{A}^k{t}^k/k!=A\cdot A^{k-1}{t}^{k-1}/(k-1)!$，求和后两项相消）可知解为

$\mathbf{u}(t)=e^{At}{\mathbf{u}}_0.$

量子实现的关键在于 $e^{At}$ 是否为酉算符。当且仅当 $A$ 反 Hermitian（即 $A=iH$、$H=H^{†}$）时，$e^{At}=e^{iHt}$ 是酉演化，可直接用哈密顿量模拟（Trotter–Suzuki、截断 Taylor 级数、QSVT 等）实现，这与本百科 HHL 教程中受控酉 $U=e^{iAt}$ 的约定一致。对一般的 $A$（例如耗散系统，$A$ 的谱在左半平面），$e^{At}$ 是收缩而非酉，需要后文的三类技术：Schrödingerization（情况二）、哈密顿量模拟的线性组合（情况三）或线性方程组方法（情况四）。

在块编码（block-encoding）模型下，设 $A$ 的块编码代价为 $C_U$，则以精度 $\epsilon$ 在归一化态上实现非酉变换 $e^{At}$ 的调用次数为

$O\left((\Vert A\Vert t+\log (1/\epsilon ))\cdot C_U\right),$

另需付出与范数放大 $e^{\Vert A\Vert t}$ 相关的后选择或振幅放大代价。这一 $\Vert A\Vert t$ 的线性依赖是最优哈密顿量模拟复杂度直接移植的结果。

情况二：时变线性 ODE

当 $A(t)$ 随时间变化时，$\exp ({\int }_0^{t}A(s)ds)$ 一般不再是解，因为不同时刻的 $A$ 不可交换。记方程的基本解（fundamental solution）为矩阵 $\Phi (t)$，它满足 $\dot{\Phi }=A(t)\Phi$、$\Phi (0)=I$，则 $\mathbf{u}(t)=\Phi (t){\mathbf{u}}_0$。处理时变性的标准方式有两种。

第一种方式是 Magnus 展开（Magnus expansion）。我们做拟设 $\Phi (t)=e^{\Omega (t)}$ 并推导 $\Omega$ 满足的方程。对恒等式（可由 $e^{\Omega }={\sum }_n{\Omega }^n/n!$ 逐项求导并将 ${\Omega }^{\prime }$ 依次穿过各因子后重新求和得到）

$\frac{d}{dt}{e}^{\Omega (t)}={\int }_0^1{e}^{s\Omega }{\Omega }^{\prime }{e}^{(1-s)\Omega }ds$

两边右乘 $e^{-\Omega }$，得

$A(t)={\int }_0^1{e}^{s\Omega }{\Omega }^{\prime }{e}^{-s\Omega }ds={\sum }_{n=0}^{\infty }\frac{1}{(n+1)!}{\mathrm{ad}}_{\Omega }^n({\Omega }^{\prime }),{\mathrm{ad}}_{\Omega }(X):=[\Omega ,X].$

这是一个关于 $\Omega$ 的隐式方程，我们按 $t$ 的幂次迭代求解：零阶给出 ${\Omega }^{\prime }=A$，即 ${\Omega }_1(t)={\int }_0^{t}A(\tau )d\tau$（量级 $O(t)$）；把 $\Omega ={\Omega }_1+{\Omega }_2+\cdots$ 代回上式，交叉项 $\frac{1}{2}[{\Omega }_1,{\Omega }_1^{\prime }]=O(t^2)$ 必须被 ${\Omega }_2$ 抵消，于是

${\Omega }_2^{\prime }=-\frac{1}{2}[{\Omega }_1,A]⟹{\Omega }_2(t)=\frac{1}{2}{\int }_0^{t}ds{\int }_0^s[A(s),A(\tau )]d\tau .$

综合起来，到二阶为止有

$\mathbf{u}(t)=\exp \left({\int }_0^{t}A(s)ds+\frac{1}{2}{\int }_0^{t}ds{\int }_0^s[A(s),A(\tau )]d\tau +O(t^3)\right){\mathbf{u}}_0.$

$[A(s),A(\tau )]=0$ 时对易子项消失，退化回常系数情形；该项就是”时变性”的全部低阶修正。截断到含 $O(t^p)$ 项的 Magnus 展开给出 $p$ 阶方法（局部误差 $O(t^{p+1})$），每一项指数因子可再由 Trotter–Suzuki 分解实现。

第二种方式是 Schrödingerization，它把 $\dot{\mathbf{u}}=A\mathbf{u}$（$A$ 任意）嵌入一个增广的薛定谔方程，我们完整推导其构造。先把 $A$ 分解为 Hermitian 部分与反 Hermitian 部分：

$H_R:=\frac{A+A^{†}}{2},H_I:=\frac{A-A^{†}}{2i},A=H_R+i{H}_I,$

其中 $H_R=H_R^{†}$、$H_I=H_I^{†}$。引入辅助变量 $\eta \in \mathbb{R}$，其上的动量算符 $P=-i{\partial }_{\eta }$ 是自伴的，在动量基 $\{|p⟩\}$ 下满足 $P|p⟩=p|p⟩$。定义增广 Hilbert 空间 ${\mathbb{C}}^d\otimes L^2(\mathbb{R})$ 上的增广哈密顿量

$\mathcal{K}=H_R\otimes P-H_I\otimes I=\frac{A+A^{†}}{2}\otimes P+\frac{i}{2}(A-A^{†})\otimes I.$

由于 $H_R\otimes P$ 与 $H_I\otimes I$ 的每个直积因子都是 Hermitian 的，$\mathcal{K}$ 是 Hermitian 的，因此演化 $e^{-i\mathcal{K}t}$ 是酉的，可用标准哈密顿量模拟实现。初值取在半直线上的指数剖面

$\Phi (0,\eta )={\mathbf{u}}_0{e}^{-\eta }\Theta (\eta ),$

其中 $\Theta$ 为阶跃函数。下面的命题给出解的恢复方式。

命题（Schrödingerization 的精确性，交换情形）：若 $[H_R,H_I]=0$，则增广薛定谔方程 $i{\partial }_t\Phi =\mathcal{K}\Phi$ 的解满足

$e^{\eta }\Phi (t,\eta )=e^{At}{\mathbf{u}}_0\text{对支撑内一切 }\eta .$

证明：把 $\mathcal{K}$ 写到 $\eta$ 表示（$P\mapsto -i{\partial }_{\eta }$），方程 $i{\partial }_t\Phi =H_{R}P\Phi -H_I\Phi$ 两边乘 $1/i=-i$ 并展开 $P$，化为输运方程

${\partial }_t\Phi =-H_R{\partial }_{\eta }\Phi +i{H}_I\Phi .$

这是速度为 $-H_R$ 的输运叠加一个相位旋转，交换情形下其解为 $\Phi (t,\eta )=e^{i{H}_{I}t}{\Phi }_0(\eta -H_{R}t)$（代入验证：${\partial }_t$ 给出 $i{H}_I$ 项加 $-H_R{\partial }_{\eta }$ 项，与右端逐项吻合）。于是

$e^{\eta }\Phi (t,\eta )=e^{i{H}_{I}t}{e}^{\eta }{\mathbf{u}}_0{e}^{-(\eta -H_{R}t)}\Theta (\eta -H_{R}t)=e^{(H_R+i{H}_I)t}{\mathbf{u}}_0=e^{At}{\mathbf{u}}_0,$

支撑条件 $\Theta (\eta -H_{R}t)$ 自动界定了恢复区域。$■$

两个极端例子可以直接核对构造无误。其一，$A=i\omega$ 为纯虚标量（$H_R=0$，$H_I=\omega$）：增广哈密顿量只剩 $\mathcal{K}=-\omega I\otimes I$，方程化为 ${\partial }_t\Phi =i\omega \Phi$，解为 $\Phi =e^{i\omega t}{\mathbf{u}}_0{e}^{-\eta }\Theta$，恢复量 $e^{\eta }\Phi =e^{i\omega t}{\mathbf{u}}_0$ 正是酉演化。其二，$A=-\gamma$ 为耗散标量（$H_R=-\gamma$，$H_I=0$）：输运速度为 $-H_R=+\gamma$，按命题的解 $\Phi (t,\eta )={\Phi }_0(\eta -H_{R}t)={\mathbf{u}}_0{e}^{-(\eta +\gamma t)}\Theta (\eta +\gamma t)$，恢复量 $e^{\eta }\Phi =e^{-\gamma t}{\mathbf{u}}_0$ 正确给出衰减因子，且范数守恒 $\Vert \Phi (t){\Vert }^2=\Vert {\mathbf{u}}_0{\Vert }^2{\int }_0^{\infty }{e}^{-2u}du=\frac{1}{2}\Vert {\mathbf{u}}_0{\Vert }^2$ 与酉性一致。

对非交换情形，我们可以对生成元 $L:=-H_R{\partial }_{\eta }+i{H}_I$ 的幂级数做归纳：由乘积法则，$L[(e^{-\eta }\Theta )w]$ 在区域内部等于 $e^{-\eta }\Theta (H_R+i{H}_I)w=e^{-\eta }\Theta Aw$，而 ${\partial }_{\eta }\Theta$ 产生的 $\delta (\eta )$ 项沿特征线传播、只影响支撑边界；逐项归纳得 $L^n{\Phi }_0$ 在内部恒等于 $e^{-\eta }{A}^n{\mathbf{u}}_0$，因此 $e^{\eta }\Phi (t,\eta )=e^{At}{\mathbf{u}}_0$ 在未被边界扰动波及的区域上精确成立（严格表述见 Jin–Liu–Yu [3]）。注意不同文献中 $P$ 与其符号的约定会交换两个直积项的位置，本词条采用上述可逐项验算的约定。

从增广态中提取归一化解 $|u(t)⟩$ 的成功概率正比于 $\Vert \mathbf{u}(t){\Vert }^2/\Vert {\mathbf{u}}_0{\Vert }^2$：解指数衰减时需要振幅放大，代价约为衰减因子本身；解范数有界时后选择以常数概率成功。整体复杂度为

$O\left(\Vert A{\Vert }_{\max }t\cdot \mathrm{polylog}(1/\epsilon )\cdot C_U\right),$

其中 $\Vert A{\Vert }_{\max }$ 为矩阵元最大模。

情况三：非齐次线性 ODE

对 $\dot{\mathbf{u}}=A\mathbf{u}+\mathbf{b}(t)$（$\mathbf{b}(t)\ne 0$，$A$ 常数），我们用常数变易法（variation of constants）推导 Duhamel 公式。做拟设 $\mathbf{u}(t)=e^{At}\mathbf{v}(t)$，代入方程：

$A{e}^{At}\mathbf{v}+e^{At}\dot{\mathbf{v}}=A{e}^{At}\mathbf{v}+\mathbf{b}(t)⟹\dot{\mathbf{v}}=e^{-At}\mathbf{b}(t).$

从 $0$ 到 $t$ 积分并注意 $\mathbf{v}(0)={\mathbf{u}}_0$，得 $\mathbf{v}(t)={\mathbf{u}}_0+{\int }_0^t{e}^{-As}\mathbf{b}(s)ds$，代回拟设并用 $e^{At}{e}^{-As}=e^{A(t-s)}$（常数矩阵可交换），即得 Duhamel 公式：

$\mathbf{u}(t)=e^{At}{\mathbf{u}}_0+{\int }_0^t{e}^{A(t-s)}\mathbf{b}(s)ds.$

量子实现的关键是把卷积积分离散化。我们把 $[0,T]$ 分成步长 $h$ 的时间片，并在每片上推导精确的单步公式。定义 $\mathbf{w}(s):=e^{-A(t_j+s)}\mathbf{u}(t_j+s)$，求导得 ${\mathbf{w}}^{\prime }(s)=e^{-A(t_j+s)}(-A\mathbf{u}+\dot{\mathbf{u}})=e^{-A(t_j+s)}\mathbf{b}(t_j+s)$；从 $0$ 到 $h$ 积分、再左乘 $e^{A(t_j+h)}$，并用常数 $A$ 下的因子分解 $e^{A(h-s)}=e^{Ah}{e}^{-As}$，得到

$\mathbf{u}(t_{j+1})=e^{Ah}\left(\mathbf{u}(t_j)+{\int }_0^h{e}^{-As}\mathbf{b}(t_j+s)ds\right).$

被积的积分可用 $p$ 阶求积（即量子求积，quantum quadrature）近似为 ${\sum }_l{w}_l{e}^{-A{s}_l}\mathbf{b}(t_j+s_l)$，误差 $O(h^{p+1})$；于是每步的作用是”一个酉演化块 $e^{iHh}$ 的线性组合”。An–Liu–Lin 的 LCHS（linear combination of Hamiltonian simulation）框架 [4] 恰好把 $e^{At}$ 严格分解为一族可模拟的酉演化 $e^{-i\tilde{H}(\alpha )t}$（$\tilde{H}(\alpha )$ 由 $A$ 的 Hermitian 与反 Hermitian 部分组装）关于实参数 $\alpha$ 的加权积分，积分用可积权重离散求积，从而使非齐次项与态制备的总代价从早期方法的 $O(T^2/\epsilon )$ 降至

$O\left(T\cdot \mathrm{polylog}(T/\epsilon )\right).$

情况四：离散化为线性方程组

早期量子 ODE 算法（Lloyd 等人的方案及其后 Clader 等人的发展）走的是另一条路：把 $\dot{\mathbf{u}}=A(t)\mathbf{u}+\mathbf{b}(t)$ 整体离散化为一个线性方程组，再调用 HHL/QSVT。我们把这条推导链逐步展开。

第一步（时间网格与积分形式）。取网格 $t_j=jh$（$j=0,\dots ,M$，$t_M=T$），对方程在 $[t_j,t_{j+1}]$ 上积分：

$\mathbf{u}(t_{j+1})-\mathbf{u}(t_j)={\int }_{t_j}^{t_{j+1}}[A(s)\mathbf{u}(s)+\mathbf{b}(s)]ds.$

第二步（梯形离散化）。对光滑被积函数 $g$，梯形法则的误差为 ${\int }_{t_j}^{t_{j+1}}g=\frac{h}{2}(g(t_j)+g(t_{j+1}))+O(h^3)$，因此

$\left(I-\frac{h}{2}{A}_{j+1}\right){\mathbf{u}}_{j+1}-\left(I+\frac{h}{2}{A}_j\right){\mathbf{u}}_j=\frac{h}{2}({\mathbf{b}}_j+{\mathbf{b}}_{j+1})+O(h^3),$

其中 $A_j:=A(t_j)$。这是隐式格式（梯形格式是 A 稳定的），全局误差 $O(h^2)$，故取 $M=O(T/\sqrt{\epsilon })$ 即可达精度 $\epsilon$。

第三步（堆叠成线性方程组）。把未知量堆叠为 $\mathbf{U}=({\mathbf{u}}_1;\dots ;{\mathbf{u}}_M)\in {\mathbb{C}}^{Md}$，已知量移到右端，得到块下双对角线性系统

$L\mathbf{U}=\mathbf{d},\mathbf{d}=\left(\begin{array}{c}(I+\frac{h}{2}{A}_0){\mathbf{u}}_0+\frac{h}{2}({\mathbf{b}}_0+{\mathbf{b}}_1)\\ \frac{h}{2}({\mathbf{b}}_1+{\mathbf{b}}_2)\\ ⋮\\ \frac{h}{2}({\mathbf{b}}_{M-1}+{\mathbf{b}}_M)\end{array}\right),$

矩阵 $L$ 的对角块为 $I-\frac{h}{2}{A}_{j+1}$、次对角块为 $-(I+\frac{h}{2}{A}_j)$。$L$ 是稀疏的（每行至多 $2s$ 个非零元，$s$ 为 $A(t)$ 的稀疏度），维数 $Md$ 只需 ${\log }_2(Md)$ 个量子比特编码——这是对维度 $d$ 的指数内存优势的来源。

第四步（条件数分析）。$L$ 并不接近单位阵（其次对角块是 $-I$ 量级的），其条件数需要单独估计。我们以标量常系数 $A=a<0$（稳定情形）、$\mathbf{b}=0$ 为例做显式计算：此时 $L$ 是下双对角矩阵，$L_{jj}=1-\frac{ha}{2}$、$L_{j,j-1}=-(1+\frac{ha}{2})$，其逆的元素为 $(L^{-1}{)}_{ji}={\rho }^{j-i}/(1-\frac{ha}{2})$，其中

$\rho =\frac{1+ha/2}{1-ha/2}\approx e^{ha}<1.$

逆矩阵的列范数为几何级数之和 ${\sum }_{j\ge i}|\rho {|}^{j-i}/(1-\frac{ha}{2})$；利用 $1-\rho =\frac{-ha}{1-ha/2}$，该和恰为 $\frac{1}{|a|h}=\frac{M}{|a|T}$。又 $\Vert L\Vert \le 2$，故

$\kappa (L)=\Theta (M)(\text{稳定情形}),$

而对不稳定情形（$a>0$）有 $|\rho {|}^M\approx e^{aT}$，$\kappa (L)$ 最高增长到 $e^{\Theta (\Vert A\Vert T)}$。这告诉我们：朴素堆叠系统的条件数随步数线性增长，不稳定时甚至指数增长。针对后者的标准补救是重缩放（rescaling）：令 ${\mathbf{u}}_j={\gamma }_j{\mathbf{w}}_j$（${\gamma }_j=e^{\lambda t_j}$、$\lambda \ge$ 谱增长上界），代入第二步的方程并除以 ${\gamma }_{j+1}$，次对角块获得因子 ${\gamma }_j/{\gamma }_{j+1}=e^{-\lambda h}<1$，等效传播子从 $\rho \approx e^{ah}$ 变为 ${\rho }^{\prime }\approx e^{(a-\lambda )h}\le 1$，指数增长被压平；而稳定情形的 $\Theta (M)$ 因子来自”小右端（$O(h)$ 量级）映射到 $O(1)$ 解”的固有范数比，是推动后续 LCHS/历史态方法取代朴素堆叠的动因之一。

第五步（量子求解与编码、归一化因子）。对 $L\mathbf{U}=\mathbf{d}$ 调用 QSVT 线性求解器，得到历史态（history state）

$|U⟩={\sum }_{j=1}^M\frac{\Vert {\mathbf{u}}_j\Vert }{\mathcal{G}}|j⟩\otimes |u_j⟩,|u_j⟩:=\frac{{\mathbf{u}}_j}{\Vert {\mathbf{u}}_j\Vert },\mathcal{G}:=({\sum }_{j=1}^M\Vert {\mathbf{u}}_j{\Vert }^2{)}^{1/2}.$

这里出现了量子 ODE 求解特有的编码问题：线性求解器输出的是归一化态，各时间片按自身范数加权。测量时间寄存器得到 $|M⟩$（即末时刻）的概率为

$P(M)=\frac{\Vert \mathbf{u}(T){\Vert }^2}{{\sum }_j\Vert {\mathbf{u}}_j{\Vert }^2},$

成功后系统寄存器坍缩为 $|\mathbf{u}(T)⟩/\Vert \mathbf{u}(T)\Vert$——相对分量正确，但绝对范数 $\Vert \mathbf{u}(T)\Vert$ 尚未获得。范数信息可以用两种方式补齐。其一是经典上界：由递推 $\Vert {\mathbf{u}}_{j+1}\Vert \le (1+h\Vert A\Vert )\Vert {\mathbf{u}}_j\Vert +h\max \Vert \mathbf{b}\Vert$ 与归纳法，并利用 $1+x\le e^x$ 和几何级数求和 ${\sum }_{i=0}^{j-1}(1+h\Vert A\Vert {)}^{j-1-i}hb\le \frac{b}{\Vert A\Vert }(e^{\Vert A\Vert t_j}-1)$，得到离散 Grönwall 不等式

$\Vert {\mathbf{u}}_j\Vert \le e^{\Vert A\Vert t_j}\Vert {\mathbf{u}}_0\Vert +\frac{\max \Vert \mathbf{b}\Vert }{\Vert A\Vert }(e^{\Vert A\Vert t_j}-1),$

它给出 $\mathcal{G}$ 与 $P(M)$ 的先验多项式界，可用于判定后选择代价。其二是量子估计：$\Vert \mathbf{u}(T){\Vert }^2={\mathbf{d}}^{†}(L^{-1}{)}^{†}{\Pi }_M{L}^{-1}\mathbf{d}$（${\Pi }_M$ 为末时间片投影）可用重叠估计以 $\mathrm{polylog}$ 代价近似。进一步地，对解范数随时间剧烈变化的系统，可做代换 $\mathbf{u}(t)=e^{\lambda t}\mathbf{w}(t)$：代入原方程得

$\dot{\mathbf{w}}=(A(t)-\lambda I)\mathbf{w}+e^{-\lambda t}\mathbf{b}(t),$

选择 $\lambda$ 落在 $A$ 的谱实部区间内即可把 $\Vert \mathbf{w}(t)\Vert$ 的动态范围压平，使各时间片权重接近、后选择概率保持在多项式水平——这正是后文 Carleman 算法中重缩放技术的同一思想。

非线性 ODE 的量子求解

非线性 ODE 是量子科学计算的核心挑战：量子力学中的演化是线性的，不能直接模拟非线性动力学。现有方法可以分为四大类。

方法一：Carleman 线性化

Carleman 线性化（Carleman linearization）由 Liu–Kolda–Sun 引入量子算法领域 [1]，并由 Costa–Jordan–Ostrander 加以改进 [2]。核心思想是：把非线性 ODE 通过引入高阶单项式新变量化为无穷维线性 ODE，再截断到有限维。

我们先对标量二次 ODE 完整推导。设 $\dot{u}=au+b{u}^2+c$，引入 $v_k:=u^k$（$k=1,2,\dots$）。由链式法则 ${\dot{v}}_k=k{u}^{k-1}\dot{u}$，逐阶代入：

${\dot{v}}_1=a{v}_1+b{v}_2+c,$ ${\dot{v}}_2=2u(au+b{u}^2+c)=2a{u}^2+2b{u}^3+2cu=2a{v}_2+2b{v}_3+2c{v}_1,$ ${\dot{v}}_k=k{v}_{k-1}(a{v}_1+b{v}_2+c)=k(a{v}_k+b{v}_{k+1}+c{v}_{k-1}).$

最后一步用了单项式乘法 $v_{k-1}{v}_1=v_k$、$v_{k-1}{v}_2=v_{k+1}$。这是一个有限带宽的无穷维线性 ODE：$v_k$ 只与 $v_{k-1},v_k,v_{k+1}$ 耦合。向量情形的推导遵循同样的 Leibniz 规则：对 $\dot{\mathbf{u}}=A\mathbf{u}+B(\mathbf{u},\mathbf{u})+\mathbf{c}$（$B$ 双线性）与 ${\mathbf{w}}_k:={\mathbf{u}}^{\otimes k}$，对 $k$ 个因子中每一个槽位求导并替换该槽位的 $\dot{\mathbf{u}}$，线性项保持次数 $k$（给出 ${\sum }_j{I}^{\otimes (j-1)}\otimes A\otimes I^{\otimes (k-j)}$ 作用在 ${\mathbf{w}}_k$ 上）、二次项升到 $k+1$ 次、常数项降到 $k-1$ 次，故

${\dot{\mathbf{w}}}_k=A_k{\mathbf{w}}_k+B_k{\mathbf{w}}_{k+1}+{\mathbf{c}}_k{\mathbf{w}}_{k-1},$

即截断前系统是关于阶数指标 $k$ 的块三对角（标量情形即三对角）线性系统。截断到 $K$ 阶（丢弃含 ${\mathbf{w}}_{K+1}$ 的项）得到有限维线性 ODE $\dot{\mathbf{V}}=M_K\mathbf{V}+\mathbf{c}$，随后即可用情况二至四的任何线性方法处理。

截断何时精确？关键在于解的一致有界性。设线性部分耗散、非线性有界：$\mathrm{Re}⟨\mathbf{u},A\mathbf{u}⟩\le -\mu \Vert \mathbf{u}{\Vert }^2$、$\Vert B(\mathbf{u},\mathbf{u})\Vert \le \beta \Vert \mathbf{u}{\Vert }^2$、$\Vert \mathbf{c}\Vert =c$，对 $g(t):=\Vert \mathbf{u}(t)\Vert$ 求导并用三角不等式得

$g^{\prime }\le -\mu g+\beta g^2+c.$

与标量方程 $y^{\prime }=-\mu y+\beta y^2+c$ 比较：其平衡点为 $y_{\pm }^{\ast }=\frac{\mu \pm \sqrt{{\mu }^2-4\beta c}}{2\beta }$（要求 ${\mu }^2\ge 4\beta c$，否则右端恒正、解无界），向量场在 $y<y_{-}^{\ast }$ 时为正、在 $y\in (y_{-}^{\ast },y_{+}^{\ast })$ 时为负，故由比较原理

$g(t)\le \max (g(0),y_{-}^{\ast })\forall t\ge 0.$

于是 $\Vert {\mathbf{w}}_{K+1}\Vert =\Vert \mathbf{u}{\Vert }^{K+1}\le R^{K+1}$ 随 $K$ 几何衰减（$R:={\sup }_{t}g(t)$），截断误差要降到 $\epsilon$ 只需

$K=O\left(\frac{\log (1/\epsilon )}{\log (1/R)}\right)(R<1).$

耗散占优时（${\mu }^2>4\beta c$ 且初值小）$R$ 可以通过整体缩放压小：令 $\mathbf{u}=\gamma \tilde{\mathbf{u}}$，方程变为 $\dot{\tilde{\mathbf{u}}}=A\tilde{\mathbf{u}}+\gamma \beta B(\tilde{\mathbf{u}},\tilde{\mathbf{u}})+\mathbf{c}/\gamma$，有效非线性强度从 $\beta$ 降为 $\gamma \beta$，这就是 Costa–Jordan–Ostrander 重缩放技术的机制 [2]。在该机制下他们的算法复杂度为

$O(\mathrm{poly}(\kappa ,\log (1/\epsilon ))\cdot T^{1+o(1)}),$

其中 $T$ 为模拟时间、$\kappa$ 为（截断后）Carleman 矩阵的条件数。方法的局限也由同一分析暴露：非线性占优（$R\ge 1$）时所需截断阶数 $K$ 随模拟时间增长，可能出现指数爆炸；$\kappa$ 本身也可能指数增长；且方法原则上只适用于多项式非线性（一般非线性需先做多项式逼近，把逼近误差计入 $\epsilon$）。2023 年的改进 [2] 正是针对这三点：更高阶的 Carleman 截断把 $K$ 从指数级压到多项式级、重缩放控制条件数、多项式逼近扩展到非多项式右端。

方法二：线性化 + 黑盒 ODE 求解

An–Liu–Lin [4] 提出把非线性右端逐点改写为参数化线性形式。由微积分基本定理（沿从 $0$ 到 $\mathbf{u}$ 的线段积分梯度），

$f(\mathbf{u})-f(0)={\int }_0^1\nabla f(s\mathbf{u})\cdot \mathbf{u}ds⟹f(\mathbf{u})=A(\mathbf{u})\mathbf{u},A(\mathbf{u}):={\int }_0^1\nabla f(s\mathbf{u})ds,$

其中用到 $\frac{d}{ds}f(s\mathbf{u})=\nabla f(s\mathbf{u})\cdot \mathbf{u}$。当 $f(0)=0$ 时非线性 ODE 严格等价于 $\dot{\mathbf{u}}=A(\mathbf{u})\mathbf{u}$；若 $f(0)\ne 0$，可增广一个恒等于 $1$ 的分量 $\tilde{\mathbf{u}}=(\mathbf{u};1)$ 把常数项并入线性部分。在每个时间片上，我们用当前估计的 $\mathbf{u}$ 冻结出 $A(\mathbf{u})$，用量子线性求解器（QSVT）局部求解线性系统，再用经典外环更新参数。由于 $A(\mathbf{u})$ 依赖 $\mathbf{u}$，这是一个量子–经典混合的迭代格式，复杂度为 $O(\mathrm{poly}(n)\cdot T/{\epsilon }^{1/p})$（$p$ 为积分阶数），其优势是原则上不限于多项式非线性。

方法三：Schrödingerization + 线性化

Jin–Liu–Yu 的框架 [3] 把前两种思想串联起来：先用 Carleman（或 Taylor）线性化把非线性 ODE 化为高维线性 ODE，再对线性系统做 Schrödingerization，最后用标准量子电路模拟增广哈密顿量。以 $\dot{u}=-u+u^2$ 为例，我们逐步写出完整流程。

第一步（Carleman 化）：这是方法一中标量情形取 $a=-1$、$b=1$、$c=0$ 的特例，递推为 ${\dot{v}}_k=-k{v}_k+k{v}_{k+1}$，堆叠成 $\dot{\mathbf{V}}=M_K\mathbf{V}$，其中 $M_K$ 是对角元 $-1,-2,\dots ,-K$、上次对角元 $1,2,\dots ,K-1$ 的双对角矩阵。

第二步（Schrödingerization）：按情况二中推导的构造，计算 $H_R=(M_K+M_K^{†})/2$（由对角元 $-k$ 与上下侧对角元 $k/2$ 组成的对称三对角阵）、$H_I=(M_K-M_K^{†})/(2i)$（由反对称的侧对角部分 $\mp ik/2$ 构成），并组装增广哈密顿量

$\mathcal{K}=H_R\otimes P+\frac{i}{2}(M_K-M_K^{†})\otimes I,$

按情况二的命题，它在初值 $\Phi (0,\eta )=\mathbf{V}(0)e^{-\eta }\Theta (\eta )$ 下以 $e^{\eta }\Phi (t,\eta )$ 恢复 $e^{M_{K}t}\mathbf{V}(0)$。

第三步（量子模拟）：$\mathcal{K}$ Hermitian 且稀疏（$M_K$ 双对角 ⟹ $\mathcal{K}$ 的每行非零元数为常数），用 Trotter 或 QSVT 实现 $e^{-i\mathcal{K}t}$。

第四步（后选择与测量）：在辅助维度的光滑区域上后选择，得到与 $\mathbf{V}(T)$ 成比例的态，其前 $d$ 个分量即原解；高阶分量携带的是截断信息，读出时丢弃。

方法四：量子 Fourier ODE 求解器

最近的工作（arXiv:2504.10218）利用量子傅里叶变换在频率空间中求解 ODE，其原理可用谐波平衡（harmonic balance）严格表述。设解是周期为 $T_p$ 的周期函数（极限环、周期轨道均属此类），展开 $\mathbf{u}(t)={\sum }_k{\hat{\mathbf{u}}}_k{e}^{i{\omega }_{k}t}$（${\omega }_k=2\pi k/T_p$），右端同样展开 $\mathbf{b}(t)={\sum }_k{\hat{\mathbf{b}}}_k{e}^{i{\omega }_{k}t}$。对常系数线性 ODE $\dot{\mathbf{u}}=A\mathbf{u}+\mathbf{b}(t)$，逐项代入：左边 $\frac{d}{dt}$ 对第 $k$ 个模式给出 $i{\omega }_k{\hat{\mathbf{u}}}_k{e}^{i{\omega }_{k}t}$，与右端同模式的系数相等，得

$(A-i{\omega }_{k}I){\hat{\mathbf{u}}}_k={\hat{\mathbf{b}}}_k\forall k.$

于是每个频率模式是一个独立的 $d$ 维线性方程组，可堆叠为一个块对角系统交给 HHL/QSVT 求解，再用逆 QFT 回到时域。该方法对周期解特别高效（无需时间网格，复杂度约 $O(\mathrm{polylog})$ 每模式），其代价是只适用于解具有已知周期或可逼近为周期的问题。

各方法复杂度对比

下表汇总各方法的关键指标，其中 $n$ 为编码维数所需的量子比特数、$n_p$ 为辅助维度所需比特数、$K$ 为 Carleman 截断阶数、$C_U$ 为相应矩阵块编码的单次代价、$T$ 为模拟时间。

方法	ODE 类型	量子比特	门数	后选择	精度
哈密顿量模拟	线性（反 Hermitian $A$）	$n$	$O((\Vert A\Vert T+\log (1/\epsilon ))C_U)$	无	$\epsilon$
Schrödingerization	线性（任意）	$n+n_p$	$O(\Vert A{\Vert }_{\max }T\cdot \mathrm{polylog}(1/\epsilon )C_U)$	有	$\epsilon$
线性方程组（HHL 类）	线性（任意）	$\log M+n$	$O(\kappa (L)\mathrm{poly}(n,\log M))$	无	$\epsilon$
Carleman + 线性求解	多项式非线性	$O(nK)$	$O(\mathrm{poly}(nK,\kappa )T^{1+o(1)})$	有	$\epsilon$
线性化 + 黑盒求解	一般非线性	$n$	$O(\mathrm{poly}(n)T/{\epsilon }^{1/p})$	无	$\epsilon$
Fourier ODE	周期线性 ODE	$n$	$O(n\log n)$	无	$\epsilon$

表中需要注意两点。其一，线性方程组方法的 $\kappa (L)$ 按情况四的分析为 $\Theta (M)$（稳定）至 $e^{\Theta (\Vert A\Vert T)}$（不稳定，重缩放后可压平），这是该方法被后继框架取代的主要原因。其二，所有方法的门数都还要乘以态制备代价；对 Carleman 方法，解范数随时间的变化由重缩放控制，否则后选择概率指数衰减。

具体例子

例子一：Lotka-Volterra 方程（捕食者-被捕食者）

Lotka–Volterra 方程为

$\dot{x}=\alpha x-\beta xy,\dot{y}=\delta xy-\gamma y,$

右端是二次的，属于 Carleman 方法的直接适用范围。我们先核对不动点：令右端为零，由 $\alpha x=\beta xy$ 得 $y^{\ast }=\alpha /\beta$，由 $\delta xy=\gamma y$ 得 $x^{\ast }=\gamma /\delta$；对参数 $\alpha =1$、$\beta =0.1$、$\delta =0.075$、$\gamma =1.5$，正平衡点为 $(x^{\ast },y^{\ast })=(20,10)$，轨道围绕它闭合并保持有界，因此截断所需的 $R=\sup \Vert \mathbf{u}\Vert$ 一致有界。数值实验表明，对模拟时间 $T=10$，截断阶数 $K=10$ 已足以达到典型精度要求。需要说明的是，Carleman 变量是全部单项式 $x^a{y}^b$（$a+b\le K$），状态空间维数为 $\left(\genfrac{}{}{0ex}{}{d+K}{K}\right)$ 量级（此处 $d=2$）。

例子二：Van der Pol 振子

Van der Pol 振子为

$\dot{x}=y,\dot{y}=\mu (1-x^2)y-x,$

其非线性项 $x^{2}y$ 是三次的。Carleman 递推中三次项使相邻阶数的耦合从 $k\leftrightarrow k\pm 1$ 变为 $k\leftrightarrow k\pm 2$，且振幅以 $\mu$ 增长，因此需要比二次系统更高的截断阶数才能达到同样精度。

例子三：Lorenz 系统（混沌）

Lorenz 系统为

$\dot{x}=\sigma (y-x),\dot{y}=x(\rho -z)-y,\dot{z}=xy-\beta z.$

右端同样是二次的，但解对初值极度敏感（混沌、蝴蝶效应）：初始误差 $\delta$ 经时间 $t$ 放大为 $\delta e^{\lambda t}$（$\lambda$ 为最大 Lyapunov 指数）。量子算法在有限时间内可以精确预测（其复杂度含 $\mathrm{poly}(e^{\lambda T})$ 型因子），但长期行为所需的精度随 $T$ 指数增长——这是信息论层面的障碍，与计算模型无关。

局限性与开放问题

1. Carleman 爆炸

对强非线性（大系数或高次非线性），由方法一的分析，当 $R\ge 1$ 时截断阶数 $K$ 必须随 $T$ 增长，最坏情形下指数爆炸；同时条件数 $\kappa$ 可能随之指数增长。重缩放与高阶截断 [2] 缓解但未彻底消除该问题。

2. 混沌系统

对混沌 ODE（如 Lorenz），量子算法不能超越经典算法的长期预测能力。蝴蝶效应是信息论的而非计算复杂度的障碍：任何计算模型都必须指数精确地知道初值。

3. 守恒量

量子模拟天然保持范数（酉演化），但 ODE 的解范数一般随时间变化。Schrödingerization 的后选择概率正比于 $\Vert \mathbf{u}(t){\Vert }^2/\Vert {\mathbf{u}}_0{\Vert }^2$，解衰减时需要振幅放大；线性方程组方法则需要按情况四的方式追踪归一化因子。

4. 边界条件

对初值问题，初始条件通过态制备进入；但边值问题（条件分布在区间两端）需要额外的打靶或伴随型编码，其量子实现仍是活跃的研究方向。

总结

量子 ODE 求解已从线性 ODE 的哈密顿量模拟，发展出时变系统的 Magnus 展开与 Schrödingerization、非齐次系统的 LCHS、以及离散化为线性方程组的 HHL 类方法；非线性系统则经 Carleman 线性化回到同一工具箱。2023 年前后的改进（重缩放、高阶 Carleman、最优态制备）把多个复杂度指标从指数级压到多项式级。核心的未解难题依然是：Carleman 截断的普适控制、病态系统的条件数、以及解的读出效率。

---

参考文献：

1. Liu, J., Kolda, T. G., & Sun, R. (2021). Quantum-inspired classical algorithm for nonlinear differential equations. arXiv:2105.02124.
2. Costa, P. C. S., Jordan, S. P., & Ostrander, A. (2023). Improved quantum algorithms for linear and nonlinear differential equations. Quantum, 7, 913.
3. Jin, S., Liu, N., & Yu, Y. (2022). Quantum simulation of partial differential equations via Schrödingerization. arXiv:2212.13969.
4. 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.

---

---

参考资料

• 本词条整理自《量子计算算法教程》原文：ch06-scientific-computing/quantum-ode-solver-tutorial
• 内容遵循 CC BY-NC-SA 4.0 许可协议