数值微积分#

Simone Montangero

本章将回顾利用有限精度数学执行微分和积分计算的数值标准方法。事实上，调用数值方法强大帮助所需付出的代价是引入一个有限精度截断，必须仔细考虑该截断，因为它会作为数值误差出现。然而，只要这种误差被控制住并正确界定，例如达到机器精度（见附录A），我们就说最终结果是数值精确的。接下来，我们将回顾求解函数积分和偏微分方程最成功的数值方法，并将它们具体应用于求解少体系统的薛定谔方程，即当系统的波函数及其上作用的算子可以显式写出时的情形。这些方法构成了基于张量网络方法的数值方法的基础，我们将在本书后续部分引入这些方法。

3.1 经典求积#

计算一维实函数 $f ( x )$ 积分的最直观方法直接遵循黎曼积分，即，在函数 $f$ 的支撑空间上定义网格 $\{x_{i} \}_{i = 1}^{N}$，并用 $N - 1$ 个梯形近似函数 $f$：

如果我们假设格点间距是等距的，$h = x_{i} - x_{i + 1}$，则得到梯形法则

其中 $f ( x_{i} ) \equiv f_{i}$，$w_{i} = 1 , 1 / 2$。因此，根据梯形法则，积分表示为在点 $x_{i}$ 处求值的函数值 $f_{i}$ 用系数 $w_{i}$ 加权后的和，再乘以 h。根据定义，结果在格点间距 h 趋于零的极限下是精确的。然而，了解每个有限 $h$ 引入的误差至关重要。误差的阶次可以通过泰勒展开来评估：误差源于用梯形近似函数 $f$ 之下每个底为 $h$ 的面积，而每个梯形的求和并不引入额外误差。因此，只需计算单个梯形的误差，不失一般性，我们假设该梯形中心在零点。然后，我们利用在 $- h / 2$ 和 $h / 2$ 处的两个函数值计算梯形的面积，并将结果与围绕零点泰勒展开得到的解析积分进行比较，即

$f ( 0 )$ 的值可以用两个端点 $f ( - h / 2 )$ 和 $f ( h / 2 )$ 来表示，再次对它们进行泰勒展开并将两个表达式相加：$\begin{array} {r} {f ( 0 ) = \frac{1} {2} ( f ( - h / 2 ) + f ( h / 2 ) ) - h^{2} / 8 f^{\prime \prime} ( 0 ) + \dots} \end{array}$。最后，我们得到积分值作为端点函数的表达式：

与式 (3.2) 比较表明，梯形法则在 $h$ 的二阶以内是精确的，因为它忽略了项 $\frac{5 h^{2}} {12} f^{\prime \prime} ( 0 )$ 及更高阶项。这一结果毫不奇怪，因为将积分近似为梯形之和等价于在每个区间内对函数 $f ( x )$ 进行线性近似。

显然，只需在 $f ( x )$ 的泰勒展开中保留更高阶项，就可以改进积分算法的精度。简单地重复前面的步骤，并考虑三个函数求值点（这是对函数进行二阶近似所必需的），我们得到

这就是辛普森法则（Simpson's rule），它规定$w_{1} = w_{N} = 1 / 3$，对于偶数项$w_{i} = 4 / 3$，对于奇数项$w_{i} = 2 / 3$。这个规定的精度可以像之前一样再次验证：高阶展开式为

代入关系式$f^{\prime \prime} ( 0 ) = ( f ( h ) - 2 f ( 0 ) + f ( - h ) ) / h^{2} + \mathcal{O} \bigl ( h^{2} \bigr )$，并忽略$\mathcal{O} \big ( h^{4} \big )$项，即可得到辛普森法则。注意，近似原则上应该是二阶的，然而，问题的对称性保证了额外一阶的精度。

这种处理方式可以推广到更高阶。但这样很快就会变得非常冗长和复杂。更方便的做法是将问题重新表述为一个线性方程组，用标准工具即可求解。事实上，我们可以直接强加式(3.2)中引入的条件，即，函数的积分与$N$个加权函数值$w_{i} f_{i}$之和在高达某阶$M + 1 \leq N$时相等。总之，对于高达$f ( x ) = x^{M}$阶的函数，我们可以强加以下条件：

可以很容易地验证，前两个方程给出了一个关于系数$w_{i}$的线性方程组，其解（在假设$a = - h / 2$且$b ~ = ~ h / 2$下）就是梯形法则。类似地，前三个方程得到辛普森法则。这种方法也强调了这种数值近似的全部潜力并未被充分挖掘。确实，先验地设定格点间距h并不方便：为了提高数值积分精度，我们可以将位置$x_{i}$也作为变量包含进来。因此，式(3.9)不仅可以求解权重$w_{i}$，还可以求解评估点$x_{i}$：变量数量加倍后，就可以用相同数量的函数求值次数来求解更多的方程。例如，强加$x_{1} = - x_{2} $（但$x_{1} \neq - h / 2$）和$| a | = | b |$，由四个方程组成的方程组

可以求解得到$\begin{array} {l l l l} {w_{1}} & {=} & {w_{2}} & {=} & {( b \mathrm{~ - ~} a ) / 2} \end{array}$和$\begin{array} {r l r l} {x_{1}} & {{} =} & {- x_{2}} & {{} =} \end{array}$ $\sqrt{( b^{3} - a^{3} ) / ( 3 ( b - a ) )}$。因此，只需改变函数采样点的位置，近似精度就能从$\mathcal{O} \big ( h^{2} \big )$提高到$\mathcal{O} \big ( h^{4} \big )$。这种适用于任意阶M的推广过程被称为高斯积分（Gaussian integration），我们将在下一节中回顾。

3.2 高斯积分#

在上一节中，我们看到可以针对$N$个函数采样点$f_{i}$的权重$w_{i}$建立并求解一个$N$元线性方程组。在$N = 2$的情况下，我们还表明，通过恰当选择采样点——即放弃标准的等间距格点——可以在不增加计算成本的情况下使积分的精度翻倍。接下来，我们将证明由于多项式基函数（如勒让德多项式）的性质，这总是可能的。确实，通过采样$N$个点上的函数值并恰当选择权重$w_{i}$，我们已经证明了可以求解前$N$个方程(3.7)–(3.10)。一般来说，以下方程定义了变量$x_{i}$的一组多项式$Z_{j}$：

对于这些方程，一般而言，在恰当选择权重 $w_{i}$ 的情况下，有 $Z_{j} ( x_{1} , x_{2} , \ldots x_{N} ) \neq 0 \forall j > N - 1$ 且 $Z_{j} ( x_{1} , x_{2} , \dots x_{N} ) =$ $0 \forall j \leq N - 1$。我们的任务是找到一组位置 $\{x_{1} , x_{2} , \ldots x_{N} \}$，使得尽可能多的多项式 $Z_{j}$ 为零，从而有效提高数值积分的近似阶数。从现在开始，我们将积分限重新缩放为 $a = - 1$ 和 $b = 1$，并将利用勒让德多项式（Legendre polynomial）$\begin{array} {r} {L_{n} ( x ) \equiv \sum_{m = 0}^{n} p_{n , m} x^{m}} \end{array}$ 及其性质，特别是：

1. $L_{n}$ 是一个 $n$ 阶多项式，具有 $n$ 个零点 $x_{i}$，满足 $L_{n} ( x_{i} ) = 0$，且 $L_{0} \equiv 1$

2. $L_{n}$ 构成 $[ - 1 ~ : ~ 1 ]$ 上函数空间的正交基：任何其他多项式都可以用这组基展开，特别地，$x^{m} =$ $\scriptstyle \sum_{n = 0}^{m} q_{m , n} L_{n} ( x )$，其中 $q_{m , n}$ 是 $p_{n , m}$ 的逆矩阵

3. $L_{n}$ 与 $x^{m}$（$m < n$）正交，因为 $x^{m} = \sum_{j=0}^{m} q_{m,j} L_{j} (x)$，且 $\langle L_{j}(x) | L_{n}(x) \rangle = 0$ 对所有 $j < n$ 成立。

现在，我们准备寻找本节开头所提问题的解：找到一组 $\{\bar{x}_{i} \}$，使得例如 $Z_{N} ( {\bar{x}}_{i} ) = 0$。为此，我们将前 $N + 1$ 个方程（3.12）—（3.15）分别乘以系数 $p_{n , j}$ 并求和，得到以下关系式

(3.19)

方程（3.19）的第二项恒为零，因为 $\begin{array} {r} {\int_{- 1}^{1} L_{n} ( x ) d x = \langle L_{0} ( x ) | L_{n} ( x ) \rangle =} \end{array}$ 0。最后，多项式 $L_{n} ( x )$ 恰好有 $N$ 个零点，因此如果选择这些零点作为 $\{\bar{x}_{i} \}$，则得到 $Z_{N} ( x_{1} , x_{2} , . . . x_{N} ) = 0$。类似地，可以证明后续的 $Z_{N + j}$ 多项式对于每个 $j \le N - 1$ 也恒为零。实际上，为了证明 $Z_{N + j}$ 为零，只需将 $Z_{j}$ 到 $Z_{N + j}$ 的方程乘以与之前相同的系数，即可得到一般关系式

和之前一样，右侧第一项对于相同的 $x_{i}$ 选择恒为零，而第二项等于 $\left. L_{j} ( x ) \middle | L_{n} ( x ) \right.$ ，对于任意 $j < n$ 均为零。总之，利用 2N 个自由参数（权重和采样位置），我们求解了 2N 个方程，从而保证了数值积分的精度为 $\mathcal{O} \big ( h^{2 N} \big )$。

3.3 不含时薛定谔方程#

在本章余下部分，我们将回顾求解偏微分方程数值解的最常见且公认的方法。特别地，我们专注于不含时和含时薛定谔方程的精确数值解。

量子力学的基本问题是寻找描述系统的哈密顿量的本征值和本征矢。一旦哈密顿量被对角化，系统的所有信息都可以被计算出来。实现这一目标的一般方法是写出哈密顿量矩阵表示，并利用第2章介绍的数值方法进行求解。下面，我们介绍将连续空间中定义的量子力学问题转化为离散化版本的标准方法，这对于将哈密顿量写成离散且有限的矩阵表示至关重要。主要方法分为两类：如第3.3.1节所述，通过引入格点对空间进行离散化；以及如第3.3.2节所示，引入合适的基函数来展开哈密顿量。

3.3.1 有限差分法#

在存在势能 $V ( \vec{x} )$ 的三维空间中，求解单粒子量子力学问题的不含时薛定谔方程可以形式化地写成

其中 $\hbar = 1$，$\vec{x}$ 和 $\vec{p}$ 分别是位置和动量的三维矢量。为了数值求解该方程，可以人为地在空间 $\{\Delta \vec{x} \}$ 上引入一个格点，从而定义一个离散的本征值问题，然后在 $\Delta \vec{x} \to 0$ 的极限下，事后检查解向连续解的收敛性。对于上一章介绍的积分，作为 $\Delta \vec{x}$ 函数的误差估计对于问题的准确离散化至关重要（见附录 B）。下面，我们将简要回顾最简单的选择，即在每个空间方向引入具有恒定间距 h 的格点。离散化方程 (3.21) 的第一个必要步骤是利用有限差分以 h 的某阶精度来近似二阶导数。与章节 B.3 中对一阶导数的处理类似，函数 $f ( x )$ 的二阶导数的离散版本可以通过函数在给定点附近的泰勒展开来计算

将上述两个方程相加并求解函数的二阶导数，我们得到

总之，在一维区间 $[ a : b ]$ 上（其中 $h = ( b - a ) / N$），四阶精度的不含时薛定谔方程的离散版本由下式给出

其中 $V_{i} = V ( x_{i} )$，我们设 $2 m = 1$。将式 (3.24) 推广到更多维度和更高阶是直接的。正如读者无疑已经注意到的，上述矩阵是三对角的（且通常是稀疏的）；因此，它可以被高效地对角化，见第2.3节。然而，特别是在维度大于一的情况下，本征值问题的规模通常会迅速变得不可处理，因为它随单一维度的格点数的维度次幂增长。例如，对每个维度取一千个点，在三维空间将导致一个十亿规模的矩阵，即使使用第2章回顾的强大方法也难以处理。因此，在这种情况下强烈推荐采用更高效的离散化方法，例如有限元方法（finite elements methods）[129]。

3.3.2 变分法#

处理高维问题的第二种方法是避免引入格点，而使用一组有限的基函数进行离散化。确实，总是可以选择一组合适的函数来展开系统的波函数：

其中 $a_{i} ~ = ~ \langle \phi_{i} | \psi \rangle$，而 $| \phi_{i} \rangle$ 可以通过物理直觉来选择，或者基于这样的假设：增加它们的数量 K，任何一组独立的波函数最终都将覆盖整个希尔伯特空间。稍后我们将看到，函数 $| \phi_{i} \rangle$ 最好是正交归一化的，但这并非提供问题精确解的必要条件。引入基函数 $| \phi_{i} \rangle$ 后，可以直接计算哈密顿量在该基下的矩阵表示，例如，对于一个处于一维外势中的粒子：

以及重叠矩阵：

选择一组正交归一基显然会得到 $S ~ = ~ \mathbb{1}$，但一般情况下并非如此。注意，这正是 Hartree-Fock 方法（结合其他近似）所采用的相同方法，它目前是用于计算分子构型最强大的方法之一 [17, 134]；也是引入诸如 Hubbard 模型和 Bose-Hubbard 模型等方法的基础，这些模型分别用于描述例如固体中的电子和光学晶格中的冷原子 [96, 98]。

现在可以计算任意波函数 $| \psi \rangle$ 的能量，同时需确保归一化，即计算：

并最小化该表达式以求得基态能量。如果波函数 $| \phi_{i} \rangle$ 是哈密顿量本征函数集，则容易证明任意波函数的能量满足 $E \ge E_{0}$，其中 $E_{0}$ 是系统的基态能量 [135]。相反，给定一组由 K 个波函数 $| \phi_{i} \rangle$ 构成的一般有限基，可以在所张的子空间内寻找系统基态的从上方近似（任何波函数的能量期望值大于或等于 $E_{0}$）。一旦完成最小化，总是可以增加基函数的数量来检查收敛性。

式 (3.28) 中的表达式被任何满足以下等式（以所选基函数 $| \phi_{i} \rangle$ 表示）且能量 E 取最小值的（归一化）波函数所最小化：

若 $S ~ = ~ \mathbb{1}$，则这又是一个维度为 $K$ 的本征值问题，可通过第2章所述的标准方法求解。反之，它表示一个广义本征值问题，可通过将其转化为标准本征值问题来求解，转化方式如下：重叠矩阵 $S$ 由定义为厄米矩阵，因此可通过酉矩阵 $W$ 对角化，使得 $W^{\dagger} \hat{S} W = D$，其中 $D$ 是某个实对角矩阵。对每个本征矢进行适当缩放后，可以找到正交变换使得 $V^{\dagger} S V = \mathbb{1}$，并定义波函数 $\left| \psi \right. = V \left| \xi \right.$，从而将方程(3.29)的本征问题重写为

因此，可以对算子 $V^{\dagger} \hat{H} V$ 和波函数 $\xi$ 定义一个标准本征值问题：求解后，可通过关系 $\left| \psi \right. = V \left| \xi \right.$ 找到原始问题的解。

我们以一个应用此方法的例子来结束本节，即如何计算非谐振子的解：一个受三阶势 $V ( x ) = x^{2} + x^{3}$ 作用的一维量子粒子。在这个简单例子中，自然选择谐振子解作为基函数集 $| n \rangle$，并在该基中展开势能。显然，结果是一个矩阵，对角线上是谐振子本征能量，非对角项由下式给出：

其中 $a^{\dagger}$ 和 $a$ 分别是产生算符和湮灭算符。因此哈密顿量是一个带宽为七的带状矩阵，需对其进行对角化。

3.4 含时薛定谔方程#

在本节中，我们回顾求解含时薛定谔方程的最常用方法，这些方法基于我们目前介绍的线性代数方法：我们的目标是求解问题，即根据下式找到量子系统的演化态 $| \psi ( t = T ) \rangle$：

给定边界条件 $| \psi ( t = 0 ) \rangle = | \psi_{0} \rangle$。此类数值方法的可能应用非常广泛，例如，从量子比特或自旋1/2系统的演化，到量子线的输运性质，或化学反应或光捕获过程的描述。

3.4.1 谱方法（Spectral Method）#

对于不含时哈密顿量，最直接的方法遵循标准量子力学步骤：将初始态分解为哈密顿量的本征函数 $\left| E_{i} \right.$，并精确计算时间演化态：

其中 $H \left| E_{i} \right. = E_{i} \left| E_{i} \right.$ 且 $c_{i} = \langle E_{i} | \psi_{0} \rangle$。这确实是最精确的方法，并且允许系统演化几乎无限长的时间 $T$，因为误差主要来自本征态 $| E_{i} \rangle$ 的计算。然而，这一优良特性也埋下了局限性的种子：计算系统的整个谱是扩展性方面的限制因素。事实上，尽管第2章介绍了强大的方法，但这种方法主要限于希尔伯特空间尺寸较小的系统，即少体离散系统或低维连续系统。

3.4.2 劈裂算子方法（Split Operator Method）#

另一个在实空间中求解不含时薛定谔方程的有效方法是劈裂算子方法。与大多数偏微分方程的数值解法一样，第一步是离散化时间变量，使得 $T = n \Delta t$（另见附录B）：第3.4节中薛定谔方程在单时间步 $\Delta t$ 后波函数的正式解为：

其中我们设 $\hbar = 1$。现在考虑一个定义在连续空间上、由动能和势能部分组成的通用一维哈密顿量：

劈裂算子方法源于以下观察：哈密顿量的这两部分都是单个变量（$\hat{x}$ 或 $\hat{p}$）的函数，因此它们分别在坐标表象和动量表象中对角化。于是，尽管这两部分不对易，利用 Baker-Campell-Hausdorff 关系 [135]，我们可以将系统从 $t$ 传播到 $t + \Delta t$ 的时间演化算子拆分为：

现在可以写出一个简单的算法来计算时间演化：式 (3.36) 右边的每一项要么在坐标表象中对角化，要么在动量表象中对角化。位置表象和动量表象之间的基变换就是傅里叶变换 $\mathcal{F}$，对此存在一种高效的算法——快速傅里叶变换（Fast Fourier Transform, FFT）。FFT 的复杂度为 $\mathcal{O}(N \log N)$，明显优于一般的基变换（矩阵向量乘法，$\mathcal{O}(N^2)$）。将基变换插入式 (3.36) 中，我们得到：

总之，只要在应用 $U_T(p)$ 或 $U_V(\boldsymbol{x})$ 之前进行适当的基变换，所有算子都以其对角形式出现。因此，劈裂算子方法只需存储 $\mathcal{O}(N)$ 个元素，其计算成本由 FFT 决定，在可应用时能提供很高的加速比。通过考虑式 (3.36) 展开中的更高阶项可以提升其精度，但这会增加算法的额外计算开销。劈裂算子方法可以直接推广到多维空间，然而，与谱方法类似，它很难应用于维度超过二维的系统。此外，劈裂算子也可用于由一般非对易项构成的不同哈密顿量。但此时不同基之间的变换将不再是傅里叶变换，因此无法利用 FFT 算法的优异缩放特性。

3.4.3 偏微分方程求解器#

当前述两种方法均无法使用时，就需应用薛定谔方程的直接积分。式 (3.34) 给出的薛定谔方程形式解的主要问题在于，要计算哈密顿量算子的指数就必须先将其对角化。然而，正如我们在谱方法中看到的，对整个哈密顿量进行对角化会极大地限制方法的适用范围。因此，其思路是计算式 (3.34) 右侧的关系，而无需显式地计算矩阵 $H$ 的指数。利用 $\Delta t$ 可以选得比系统任何能量尺度都小的性质，该问题的一个现成解法当然是将指数按 $\Delta t$ 展开成级数：一阶展开得到欧拉方法（Euler method）

然而，这种粗糙的近似存在一个主要问题：它不守恒态 (state) 的范数。事实上，近似的时间演化算符并非幺正的。通过在 $\Delta t$ 的二阶下计算一个时间步长后时间演化态的范数，可以很容易看出这一点：$\mid \langle \psi ( t + \Delta t ) | \psi ( t + \Delta t ) \rangle \mid^{2} = | 1 - \Delta t^{2} \langle \Delta \hat{H} \rangle^{2} | ^{2}$，其中 $\langle \Delta \hat{H} \rangle^{2} = \langle \hat{H}^{2} \rangle - \langle \hat{H} \rangle^{2}$ 是在时间 t 计算得到的哈密顿量算符的（正定）方差。

为了修正这种展开带来的不便之处，可以通过将时间演化算符表达为两个伴随算符乘积的函数，来强制截断的时间演化算符满足幺正性。事实上，以下关系成立：

式(3.43)的右边部分在构造上就是幺正的，因此解决了范数守恒问题，然而这也引入了计算矩阵逆的需要。这个问题的解决方案被称为克兰克-尼科尔森方法 (Crank-Nicolson method)，它利用了第2章所述，我们能够高效求解线性方程组的能力。事实上，由式(3.34)和(3.43)可推出关系：

因此，要计算演化后的态 $| \psi ( t + \Delta t ) \rangle$，必须求解由上述表达式定义的线性方程组。

最后，我们将展示上述方法仅是偏微分方程数值积分的一个特定示例：正如我们在第3章所见，根据所使用的插值阶数，方法的精度可以大幅提升。事实上，将薛定谔方程写为积分形式：

其中 ${\mathcal{L}} = - i {\hat{H}}$ ，显然求解此问题等价于计算方程右边的积分。因此，我们得到：

其中 $w_{i}$ 是可自由选择的权重。上述关系可以重写为：

这里 $\hat{A} = \mathbb{1} - w_{2} \mathcal{L} \Delta t$ 且 $\hat{B} = \mathbb{1} - w_{1} \mathcal{L} \Delta t$，时间传播算符定义为 $\hat{A}^{- 1} \hat{B}$。梯形法则 (Trapezoidal rule) 规定 $w_{i} = 1 / 2$，这正好对应于克兰克-尼科尔森方法。然而，原则上可以进行不同的选择（例如，用一个高度等于函数一个端点值的矩形来近似积分），选择 $w_{1} = 1 , w_{2} = 0$，就得到了欧拉方法，正如我们所看到的，其精度不如前者。

既然我们已经将通过式(3.46)对李维利算符 (Liouvillian operator) $\mathcal{L}$ 的积分与含时薛定谔方程的求解对应起来，那么现在应该清楚如何利用第3.2节中介绍的高斯积分方案 (Gaussian integration schemes) 来超越克兰克-尼科尔森方法，提高精度了。实际上，我们可以引入 K 个中间函数评估点，并将式(3.46)中的积分计算为：

其中 $\tau_{0} = 0$，$\tau_{K} = \Delta t$。我们在第3.2节中看到，通过高斯积分，原则上可以达到 $\mathcal{O} \big ( \Delta t^{2K} \big )$ 的精度。然而，这里出现了一个额外的复杂性，因为在方程 (3.48) 中引入了 K+1 个额外变量：中间波函数 $| \psi ( t + \tau_{k} ) \rangle , k = 1, ..., K-1$ 是未知的（与之前对已知函数 $f(x)$ 进行积分的情况不同），并且需要与 $| \psi ( t + \Delta t ) \rangle$ 一同计算。解决此问题的方法在于利用被积函数 $\mathcal{L}_{k} \equiv \mathcal{L} | \psi ( t + \tau_{k} ) \rangle$ 的结构，即，利用每个中间时间步的方程 (3.48) 中的关系、第3.4节以及波函数在初始时间 t 附近的泰勒展开。事实上，我们可以将每个时间步 $\tau_{k}$ 的方程 (3.48) 关系写为

其中 $T_{i , j}$ 是一个传递矩阵，包含待确定的系数，这些系数将前 $k-1$ 个中间 $L_{k}$ 与第 k 个联系起来。这样就可以为系数 $w_{i}$（将 $L_{k}$ 与最终状态联系起来）、传递矩阵系数 $T_{i , j}$ 以及中间时间间隔 $\tau_{k}$ 建立一个方程组，以实现高斯积分。下面我们给出 $K = 2$ 时的解，即四阶龙格-库塔方法（Runge-Kutta method）。对于 $K = 2$，方程 (3.48) 变为

回忆 $\mathcal{L} = \mathcal{L} ( \psi , t )$，我们可以将其展开至一阶，得到 $\mathcal{L}_{1} = \mathcal{L}_{0} + \frac{\partial{\mathcal{L}}}{\partial \psi} \Delta \psi + \frac{\partial \mathcal{L}}{\partial t} \tau_{1}$。将方程 (3.49) 代入，并回忆 $\tau_{1} = c_{1} \Delta t$（其中 $c_{1}$ 是某个系数），可将方程 (3.52) 重写为

将上式逐项与 $\psi ( t + \Delta t )$ 的二阶泰勒展开式相等，可得到方程组

其解提供了实现 $\Delta t$ 四阶龙格-库塔求解器所需的系数。在大多数通用平台和编程语言中，都可以找到高阶近似和自适应龙格-库塔方法（adaptive Runge-Kutta methods）的实现。

3.5 习题#

1 考虑由哈密顿量定义的一维量子谐振子

及其本征值 $E_{k}$ 和本征矢 $\vert \psi_{k} \rangle ( \hbar = m = \omega = 1 )$

(a) 给定精确基态 $| \psi_{0} \rangle$，数值计算基态能量期望值 $\langle \psi_{0} | H | \psi_{0} \rangle$，并与其精确值 $E_{0}$ 进行比较

(b) 误差来源于波函数离散化还是积分近似？提示：在固定波函数离散化的情况下提高积分近似精度（梯形法则、辛普森法则等），反之亦然。

2 考虑习题1哈密顿量定义的一维量子谐振子。

(a) 编写一个 FORTRAN 程序，计算前 k 个本征值 $E_{k}$ 和本征矢 $| \psi_{k} \rangle$

(b) 根据课堂上引入的优秀科学软件开发优先级（正确性、稳定性、精确离散化、灵活性、效率）（见附录B），你如何评价你的程序？

3 给定由哈密顿量定义的时间依赖一维量子谐振子

其中 $q_{0} ( t ) = t / T$ 且 $t ~ \in ~ [ 0 ~ : ~ T ]$。给定 $| \psi_{0} \rangle = | n = 0 \rangle$（谐振子的基态），计算不同 $T$ 值下的 $| \psi ( t ) \rangle$