# 量子 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}$，已知量移到右端，得到块下双对角线性系统 $$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},$$ 矩阵 $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$，否则右端恒正、解无界），向量场在 $y4\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. --- > 返回目录：[量子计算算法教程系列](https://chenzhaoyun.com/index.php/archives/54/)