量子 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), \qquad \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)\rangle \propto \mathbf{u}(T)\)，或估计 \(\langle u(T)|M|u(T)\rangle\) 一类的全局泛函。

经典方法

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

\[N_{\mathrm{step}} = O\!\left(T\,\varepsilon^{-1/p}\right)\]

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

线性 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^\dagger\)）时，\(e^{At} = e^{iHt}\) 是酉演化，可直接用哈密顿量模拟（Trotter–Suzuki、截断 Taylor 级数、QSVT 等）实现，这与本系列 HHL 教程中受控酉 \(U = e^{iAt}\) 的约定一致。对一般的 \(A\)（例如耗散系统，\(A\) 的谱在左半平面），\(e^{At}\) 是收缩而非酉，需要后文的三类技术：Schrödingerization（情况二）、哈密顿量模拟的线性组合（情况三）或线性方程组方法（情况四）。

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

\[O\!\left((\|A\|\,t + \log(1/\varepsilon))\cdot C_U\right),\]

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

情况二：时变线性 ODE

当 \(A(t)\) 随时间变化时，\(\exp\!\big(\int_0^t A(s)\,ds\big)\) 一般不再是解，因为不同时刻的 \(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'\) 依次穿过各因子后重新求和得到）

\[\frac{d}{dt}e^{\Omega(t)} = \int_0^1 e^{s\Omega}\,\Omega'\,e^{(1-s)\Omega}\,ds\]

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

\[A(t) = \int_0^1 e^{s\Omega}\,\Omega'\,e^{-s\Omega}\,ds = \sum_{n=0}^{\infty}\frac{1}{(n+1)!}\,\mathrm{ad}_{\Omega}^{\,n}(\Omega'),\qquad \mathrm{ad}_{\Omega}(X):=[\Omega,X].\]

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

\[\Omega_2' = -\tfrac12[\Omega_1, A]\quad\Longrightarrow\quad \Omega_2(t) = \frac12\int_0^t ds\int_0^{s}[A(s),A(\tau)]\,d\tau.\]

综合起来，到二阶为止有

\[\mathbf{u}(t) = \exp\!\left(\int_0^t A(s)\,ds + \frac12\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^\dagger}{2},\qquad H_I := \frac{A-A^\dagger}{2i},\qquad A = H_R + i\,H_I,\]

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

\[\mathcal{K} \;=\; H_R\otimes P \;-\; H_I\otimes I \;=\; \frac{A+A^\dagger}{2}\otimes P \;+\; \frac{i}{2}\bigl(A-A^\dagger\bigr)\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\qquad\text{对支撑内一切 }\eta.\]

证明：把 \(\mathcal K\) 写到 \(\eta\) 表示（\(P\mapsto -i\partial_\eta\)），方程 \(i\partial_t\Phi = H_RP\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^{iH_I t}\,\Phi_0(\eta - H_R t)\)（代入验证：\(\partial_t\) 给出 \(iH_I\) 项加 \(-H_R\partial_\eta\) 项，与右端逐项吻合）。于是

\[e^{\eta}\Phi(t,\eta) = e^{iH_I t}\,e^{\eta}\,\mathbf u_0\,e^{-(\eta-H_Rt)}\,\Theta(\eta-H_Rt) = e^{(H_R+iH_I)t}\,\mathbf u_0 = e^{At}\mathbf u_0,\]

支撑条件 \(\Theta(\eta-H_Rt)\) 自动界定了恢复区域。\(\blacksquare\)

两个极端例子可以直接核对构造无误。其一，\(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_0e^{-\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_Rt)=\mathbf u_0e^{-(\eta+\gamma t)}\Theta(\eta+\gamma t)\)，恢复量 \(e^\eta\Phi = e^{-\gamma t}\mathbf u_0\) 正确给出衰减因子，且范数守恒 \(\|\Phi(t)\|^2 = \|\mathbf u_0\|^2\int_0^\infty e^{-2u}du = \tfrac12\|\mathbf u_0\|^2\) 与酉性一致。

对非交换情形，我们可以对生成元 \(L:=-H_R\partial_\eta + iH_I\) 的幂级数做归纳：由乘积法则，\(L\big[(e^{-\eta}\Theta)w\big]\) 在区域内部等于 \(e^{-\eta}\Theta\,(H_R+iH_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)\rangle\) 的成功概率正比于 \(\|\mathbf u(t)\|^2/\|\mathbf u_0\|^2\)：解指数衰减时需要振幅放大，代价约为衰减因子本身；解范数有界时后选择以常数概率成功。整体复杂度为

\[O\!\left(\|A\|_{\max}\,t\cdot \mathrm{polylog}(1/\varepsilon)\cdot C_U\right),\]

其中 \(\|A\|_{\max}\) 为矩阵元最大模。

情况三：非齐次线性 ODE

对 \(\dot{\mathbf{u}} = A\mathbf{u} + \mathbf{b}(t)\)（\(\mathbf b(t)\neq 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)\quad\Longrightarrow\quad \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'(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^{-As_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\widetilde H(\alpha)t}\)（\(\widetilde H(\alpha)\) 由 \(A\) 的 Hermitian 与反 Hermitian 部分组装）关于实参数 \(\alpha\) 的加权积分，积分用可积权重离散求积，从而使非齐次项与态制备的总代价从早期方法的 \(O(T^2/\varepsilon)\) 降至

\[O\!\left(T\cdot\mathrm{polylog}(T/\varepsilon)\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}}\bigl[A(s)\mathbf u(s) + \mathbf b(s)\bigr]\,ds .\]

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

\[\left(I - \tfrac{h}{2}A_{j+1}\right)\mathbf u_{j+1} - \left(I + \tfrac{h}{2}A_j\right)\mathbf u_j = \tfrac{h}{2}\big(\mathbf b_j + \mathbf b_{j+1}\big) + O(h^3),\]

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

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

\[\begin{split}L\mathbf U = \mathbf d,\qquad \mathbf d = \begin{pmatrix}\big(I+\tfrac h2 A_0\big)\mathbf u_0 + \tfrac h2(\mathbf b_0+\mathbf b_1)\\ \tfrac h2(\mathbf b_1+\mathbf b_2)\\ \vdots\\ \tfrac h2(\mathbf b_{M-1}+\mathbf b_M)\end{pmatrix},\end{split}\]

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

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

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

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

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

而对不稳定情形（\(a>0\)）有 \(|\rho|^M\approx e^{aT}\)，\(\kappa(L)\) 最高增长到 \(e^{\Theta(\|A\|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' \approx e^{(a-\lambda)h}\le 1\)，指数增长被压平；而稳定情形的 \(\Theta(M)\) 因子来自"小右端（\(O(h)\) 量级）映射到 \(O(1)\) 解"的固有范数比，是推动后续 LCHS/历史态方法取代朴素堆叠的动因之一。

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

\[|U\rangle = \sum_{j=1}^{M}\frac{\|\mathbf u_j\|}{\mathcal G}\,|j\rangle\otimes|u_j\rangle,\qquad |u_j\rangle := \frac{\mathbf u_j}{\|\mathbf u_j\|},\qquad \mathcal G := \Big(\sum_{j=1}^M\|\mathbf u_j\|^2\Big)^{1/2}.\]

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

\[P(M) = \frac{\|\mathbf u(T)\|^2}{\sum_j\|\mathbf u_j\|^2},\]

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

\[\|\mathbf u_j\| \;\le\; e^{\|A\|t_j}\,\|\mathbf u_0\| + \frac{\max\|\mathbf b\|}{\|A\|}\Big(e^{\|A\|t_j}-1\Big),\]

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

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

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

非线性 ODE 的量子求解

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

方法一：Carleman 线性化

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

我们先对标量二次 ODE 完整推导。设 \(\dot u = a u + 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+bu^2+c) = 2a\,u^2 + 2b\,u^3 + 2c\,u = 2a v_2 + 2b v_3 + 2cv_1,\]
\[\dot v_k = k\,v_{k-1}\,(a v_1 + b v_2 + c) = k\big(a v_k + b v_{k+1} + c\,v_{k-1}\big).\]

最后一步用了单项式乘法 \(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}\langle\mathbf u,A\mathbf u\rangle\le-\mu\|\mathbf u\|^2\)、\(\|B(\mathbf u,\mathbf u)\|\le\beta\|\mathbf u\|^2\)、\(\|\mathbf c\| = c\)，对 \(g(t):=\|\mathbf u(t)\|\) 求导并用三角不等式得

\[g' \;\le\; -\mu g + \beta g^2 + c .\]

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

\[g(t)\;\le\;\max\bigl(g(0),\,y_-^*\bigr)\qquad\forall t\ge0 .\]

于是 \(\|\mathbf w_{K+1}\| = \|\mathbf u\|^{K+1}\le R^{K+1}\) 随 \(K\) 几何衰减（\(R:=\sup_t g(t)\)），截断误差要降到 \(\varepsilon\) 只需

\[K \;=\; O\!\left(\frac{\log(1/\varepsilon)}{\log(1/R)}\right)\qquad(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\bigl(\mathrm{poly}(\kappa,\log(1/\varepsilon))\cdot T^{1+o(1)}\bigr),\]

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

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

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

\[f(\mathbf u) - f(\mathbf 0) = \int_0^1 \nabla f(s\mathbf u)\cdot\mathbf u\;ds\quad\Longrightarrow\quad f(\mathbf u) = A(\mathbf u)\,\mathbf u,\quad 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(\mathbf 0)=0\) 时非线性 ODE 严格等价于 \(\dot{\mathbf u} = A(\mathbf u)\mathbf u\)；若 \(f(\mathbf 0)\neq 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/\varepsilon^{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^\dagger)/2\)（由对角元 \(-k\) 与上下侧对角元 \(k/2\) 组成的对称三对角阵）、\(H_I = (M_K-M_K^\dagger)/(2i)\)（由反对称的侧对角部分 \(\mp ik/2\) 构成），并组装增广哈密顿量

\[\mathcal{K} = H_R\otimes P + \frac{i}{2}\bigl(M_K - M_K^\dagger\bigr)\otimes I,\]

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

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

第四步（后选择与测量）：在辅助维度的光滑区域上后选择，得到与 \(\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}\)，与右端同模式的系数相等，得

\[\bigl(A - i\omega_k I\bigr)\,\hat{\mathbf u}_k = \hat{\mathbf b}_k\qquad\forall k .\]

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

各方法复杂度对比

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

方法

ODE 类型

量子比特

门数

后选择

精度

哈密顿量模拟

线性（反 Hermitian \(A\)）

\(n\)

\(O((|A|T+\log(1/\varepsilon))\,C_U)\)

无

\(\varepsilon\)

Schrödingerization

线性（任意）

\(n+n_p\)

\(O(|A|_{\max}T\cdot\mathrm{polylog}(1/\varepsilon)\,C_U)\)

有

\(\varepsilon\)

线性方程组（HHL 类）

线性（任意）

\(\log M + n\)

\(O(\kappa(L)\,\mathrm{poly}(n,\log M))\)

无

\(\varepsilon\)

Carleman + 线性求解

多项式非线性

\(O(nK)\)

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

有

\(\varepsilon\)

线性化 + 黑盒求解

一般非线性

\(n\)

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

无

\(\varepsilon\)

Fourier ODE

周期线性 ODE

\(n\)

\(O(n\log n)\)

无

\(\varepsilon\)

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

具体例子

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

Lotka–Volterra 方程为

\[\dot{x} = \alpha x - \beta xy, \qquad \dot{y} = \delta xy - \gamma y ,\]

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

例子二：Van der Pol 振子

Van der Pol 振子为

\[\dot{x} = y, \qquad \dot{y} = \mu(1-x^2)y - x ,\]

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

例子三：Lorenz 系统（混沌）

Lorenz 系统为

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

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

局限性与开放问题

1. Carleman 爆炸

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

2. 混沌系统

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

3. 守恒量

量子模拟天然保持范数（酉演化），但 ODE 的解范数一般随时间变化。Schrödingerization 的后选择概率正比于 \(\|\mathbf u(t)\|^2/\|\mathbf u_0\|^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.