第7章 算法与数值方法

原著：Xinfeng Zhou，A Practical Guide to Quantitative Finance Interviews（2008）。

本页由 GPT 视觉模型直接依据扫描页识别、翻译和公式排版。

查看本页原始扫描图

虽然量化分析师在编程上花费的时间比例因职位（例如，量化分析师/研究员与量化开发人员）和公司文化而异，但典型的量化分析师通常会花费部分时间通过编程来实现模型。因此，编程技能测试通常是量化面试的一个固有部分。 在很大程度上，量化面试中提出的编程问题与技术面试中提出的问题相似。毫不奇怪，其中许多问题是平台或语言特定的。尽管 C++ 和 Java 仍然主导市场，但我们已经看到向其他编程语言（如 Matlab、SAS、S-Plus 和 R）的日益多样化。由于有许多现有的书籍和网站专门介绍技术面试，本章将不提供对编程问题的全面回顾。相反，它讨论了一些量化面试中偏爱的算法问题和数值方法。

7.1. 算法

在编程中，算法复杂度分析通常使用渐近分析，它忽略了机器相关的常数，并研究运行时间 \(T(n)\)——即原始操作的数量，如加法、乘法和比较——随着输入数量 \(n \to \infty\) 的变化趋势。1脚注 1如果您想复习基础算法，我强烈推荐 Thomas H. Cormen、Charles E. Leiserson、Ronald L. Rivest 和 Clifford Stein 的《算法导论》。它涵盖了本节讨论的所有理论，并包含了许多面试中经常出现的算法。 算法复杂度中最重要的三种记法是“大 O 记法”（big-O notation）、“Ω 记法”（Omega notation）和“Θ 记法”（Theta notation）： \(O(g(n)) = \{ f(n) : \text{存在正的常数 } c \text{ 和 } n_0 \text{，使得对于所有 } n \ge n_0 \text{，都有 } 0 \le f(n) \le cg(n) \}\)。

这是 \(f(n)\) 的渐近上界。\(\Omega(g(n)) = \{ f(n) : \text{存在正的常数 } c \text{ 和 } n_0 \text{，使得对于所有 } n \ge n_0 \text{，都有 } 0 \le cg(n) \le f(n) \}\)。这是 \(f(n)\) 的渐近下界。\(\Theta(g(n)) = \{ f(n) : \text{存在正的常数 } c_1, c_2, \text{ 和 } n_0 \text{，使得对于所有 } n \ge n_0 \text{，都有 } c_1g(n) \le f(n) \le c_2g(n) \}\)。这是 \(f(n)\) 的渐近紧确界。除了记法之外，解释算法复杂度中的两个概念也很重要：

查看本页原始扫描图

最坏情况运行时间 $W(n)$：任意 $n$ 个输入的运行时间的上界。平均情况运行时间 $A(n)$：如果 $n$ 个输入是随机选择的，则期望运行时间。对于许多算法，$W(n)$ 和 $A(n)$ 具有相同的 $O(g(n))$。但正如我们将讨论的，在某些问题中，它们可能大不相同，其相对重要性通常取决于具体问题。一个有 $n$ 个输入的问题可以经常被划分为 $a$ 个子问题，每个子问题有 $n/b$ 个输入。这种范式通常称为分治法。如果需要 $f(n)$ 个原始操作来将问题划分为子问题并合并子问题的解，则运行时间可以表示为递推方程： \[ T(n) = aT(n/b) + f(n) \] 其中 $a > 1$, $b > 1$, 且 $f(n) \ge 0$。

主定理是求解递推方程 $T(n) = aT(n/b) + f(n)$ 的紧确界的一个有价值的工具： 如果 $f(n) = O(n^{\log_b a - \varepsilon})$ 对于某个常数 $\varepsilon > 0$，则 $T(n) = \Theta(n^{\log_b a})$，因为 $f(n)$ 的增长速度慢于 $n^{\log_b a}$。如果 $f(n) = \Theta(n^{\log_b a} \log^k n)$ 对于某个 $k \ge 0$，则 $T(n) = \Theta(n^{\log_b a} \log^{k+1} n)$，因为 $f(n)$ 和 $n^{\log_b a}$ 的增长速率相似。

如果 $f(n) = \Omega(n^{\log_b a + \varepsilon})$ 对于某个常数 $\varepsilon > 0$，并且 $af(n/b) \le cf(n)$ 对于某个常数 $c < 1$，则 $T(n) = \Theta(f(n))$，因为 $f(n)$ 的增长速度快于 $n^{\log_b a}$。

补充说明： 主定理可以直观理解为比较“分解与合并的成本”$f(n)$ 与“叶子节点数量”$n^{\log_b a}$ 的增长速度。如果分解合并成本小，总时间由叶子节点主导；如果成本大，总时间由分解合并主导；如果两者相当，则总时间等于叶子节点数乘以对数因子。

让我们用二分查找来展示主定理的应用。要在数组中查找一个元素，如果数组中的数字是排序的（$a_1 \le a_2 \le \dots \le a_n$），我们可以使用二分查找：算法从 $a_{\lfloor n/2 \rfloor}$ 开始。如果 $a_{\lfloor n/2 \rfloor} = x$，则搜索停止。如果 $a_{\lfloor n/2 \rfloor} > x$，我们只需要搜索 $a_1, \dots, a_{\lfloor n/2 \rfloor - 1}$。如果 $a_{\lfloor n/2 \rfloor} < x$，我们只需要搜索 $a_{\lfloor n/2 \rfloor + 1}, \dots, a_n$。每次比较后，我们可以将需要搜索的元素数量减半。所以我们有 $a = 1$, $b = 2$, 且 $f(n) = 1$。

因此， \[ f(n) = \Theta(n^{\log_2 1} \log^0 n) \] 二分查找的复杂度为 $\Theta(\log n)$。

数字交换

如何在不使用额外存储空间的情况下交换两个整数 $i$ 和 $j$？ 解答： 比较和交换是许多算法的基本操作。最常见的交换技术是使用一个临时变量，但在这个问题中，临时变量需要额外的存储空间，因此是被禁止的。一个简单的

查看本页原始扫描图

数学方法是先存储 i 和 j 的和，然后提取 i 的值并将其赋给 j，最后将 j 的值赋给 i。实现如以下代码所示：

void swap(int &i, int &j) {
    i = i + j; // 存储 i 和 j 的和
    j = i - j; // 将 j 变为 i 的值
    i = i - j; // 将 i 变为 j 的值
}

另一种解决方案利用按位异或（^）函数，其原理是 \(x \oplus x = 0\) 且 \(0 \oplus x = x\):

void swap(int &i, int &j) {
    i = i ^ j;
    j = j ^ i; // j = j ^ (i ^ j) = (j ^ j) ^ i = 0 ^ i = i
    i = i ^ j; // i = (i ^ j) ^ i = (i ^ i) ^ j = 0 ^ j = j
}

唯一元素

如果你有一个已排序的数组，你能写一些代码来提取数组中的唯一元素吗？例如，如果数组是 [1, 1, 3, 3, 3, 5, 5, 5, 9, 9, 9, 9]，那么唯一元素应该是 [1, 3, 5, 9]。 解法： 设 a 是一个 n 个元素的已排序数组，其元素满足 \(a_0 \le a_1 \le \dots \le a_{n-1}\)。当我们遇到已排序数组中的一个新元素 \(a_i\) 时，它的值与前一个元素不同 (\(a_i \ne a_{i-1}\))。利用这个性质，我们可以轻松地提取唯一元素。C++ 的一种实现如下所示的函数：

template <class T> vector<T> unique(T a[], int n) {
    vector<T> vec; // 使用 vector 来避免调整大小的问题
    vec.reserve(n); // 预留空间以避免重新分配
    vec.push_back(a[0]);
    for(int i=1; i<n; ++i) {
        // 在这里添加提取唯一元素的逻辑
    }
    return vec;
}

2脚注 2本章使用 C++ 来演示一些实现。对于其他问题，算法使用伪代码描述。 下面是一个交换两个整数的单行等效函数。但不推荐使用，因为它缺乏清晰度。

void swap(int &i, int &j) { i-=j=(i+=j)-j; };

3脚注 3我应该指出 C++ STL 有通用的算法来执行这个基本操作：unique 和 unique_copy。

查看本页原始扫描图

if(a[i] != a[i-1])
    vec.push_back(a[i]);
}
return vec;

霍纳算法

编写一个算法来计算 \(y = A_0 + A_1x + A_2x^2 + A_3x^3 + \dots + A_nx^n\)。 解答： 一种朴素的方法是分别计算多项式的每一项然后相加，这需要 \(O(n^2)\) 次乘法。我们可以使用霍纳算法将乘法次数减少到 \(O(n)\)。该算法将原始多项式表示为 \(y = ((((A_nx + A_{n-1})x + A_{n-2})x + \dots + A_2)x + A_1)x + A_0\)，并依次计算 \(B_n = A_n, B_{n-1} = B_nx + A_{n-1}, \dots, B_0 = B_1x + A_0\)。我们得到 \(y = B_0\)，最多进行 \(n\) 次乘法。

补充说明： 霍纳算法本质上是通过“嵌套乘法”来减少计算量。例如，对于 \(y = 2 + 3x + 4x^2 + 5x^3\)，可以计算为 \(y = ((5x + 4)x + 3)x + 2\)。这样，计算 \(x^3\) 不需要单独做两次乘法（\(x*x*x\)），而是融入到连续的乘法中。

移动平均

给定一个长度为 \(m\) 的大数组 \(A\)，你能开发一个高效的算法来构建另一个数组，该数组包含原始数组的 \(n\) 元素移动平均值吗？（即 \(B_1, \dots, B_{n-1} = NA, B_i = (A_{i-n+1} + A_{i-n+2} + \dots + A_i)/n, \forall i = n, \dots, m\)） 解答： 当我们计算接下来 \(n\) 个连续数字的移动平均值时，可以重用之前计算过的移动平均值。只需将之前的平均值乘以 \(n\)，减去该移动平均值中的第一个数字，然后加上新的数字，就得到了新的总和。将新的总和除以 \(n\) 就得到新的移动平均值。以下是计算移动平均值的伪代码：

S = A[1] + \dots + A[n]; B[n] = S/n;
for (i=n+1 to m) { S = S - A[i-n] + A[i]; B[i] = S/n; }

补充说明： 这种方法的复杂度是 \(O(m)\)，远优于对每个 \(B_i\) 都重新计算 \(n\) 个元素之和的朴素方法（复杂度 \(O(mn)\)）。它利用了相邻窗口之间大部分元素重叠的特性。

排序算法

你能解释三种对 \(n$ 个不同值 $A_1, \dots, A_n\) 进行排序的算法，并分析每种算法的复杂度吗？ 解答： 排序是一个基本过程，直接或间接地在许多程序中实现。因此，针对不同的

查看本页原始扫描图

目的，已经开发了多种排序算法。这里我们讨论三种这样的算法：插入排序、归并排序和快速排序。

插入排序

插入排序采用一种增量方法。假设我们已经排好了子数组 A[1, ..., i-1]。我们将元素 \(A_i\) 插入到 A[1, ..., i-1] 的适当位置，从而得到排好序的子数组 A[1, ..., i]。从 i=1 开始，逐步增加 i 到 n，我们将得到一个完全排序的数组。对于每一步，期望的比较次数是 \(i/2\)，最坏情况的比较次数是 \(i\)。所以我们有 \[ A(n) = \Theta\left(\sum_{i=1}^{n} i/2\right) = \Theta(n^2) \text{ 且 } W(n) = \Theta\left(\sum_{i=1}^{n} i\right) = \Theta(n^2). \]

补充说明： 插入排序类似于整理手中的扑克牌。每次拿到一张新牌（新元素），就将其插入到手中有序牌序列的正确位置。当数组基本有序时，插入排序效率很高；当数组完全逆序时，效率最低。

归并排序

归并排序采用了分治策略。它将数组分成两个子数组，每个子数组各有 \(n/2\) 个元素，然后对每个子数组分别排序。除非子数组足够小（元素个数不超过几个），否则它会再次被分割以便排序。最后，将排好序的子数组合并成一个有序的整个数组。 该算法可以用如下伪代码表示：

mergesort(A, beginindex, endindex)
if beginindex < endindex
then centerindex ← (beginindex + endindex)/2
merge1 ← mergesort(A, beginindex, centerindex)
merge2 ← mergesort(A, centerindex + 1, endindex)
merge(merge1, merge2)

将两个已排序的、各有 \(n/2\) 个元素的数组合并成一个数组，需要进行 \(\Theta(n)\) 次基本操作。其运行时间 \(T(n)\) 满足以下递归函数： \[ T(n) = \begin{cases} 2T(n/2) + \Theta(n), & \text{if } n > 1 \\ 1, & \text{if } n = 1 \end{cases} \] 对 \(T(n)\) 应用主定理，其中 \(a=2, b=2\)，且 \(f(n)=\Theta(n)\)，我们得到 \(f(n) = \Theta(n^{\log_b a} \log^0 n)\)。所以 \(T(n) = \Theta(n \log n)\)。对于归并排序，平均情况下的运行时间 \(A(n)\) 和最坏情况下的运行时间 \(W(n)\) 都与 \(T(n)\) 相同。

快速排序

快速排序是另一种递归排序方法。它从序列中选取一个元素 \(A_i\)，并将所有其他元素与它进行比较。那些小于 \(A_i\) 的元素被放到 \(A_i\) 左侧的一个子数组中；那些大于 \(A_i\) 的元素被放到 \(A_i\) 右侧的一个子数组中。然后，在两个子数组（以及它们所产生的任何子数组）上重复这个过程，直到所有值都排好序。

查看本页原始扫描图

在最坏情况下，快速排序所需的比较次数与插入排序相同。例如，如果我们总是选择数组（或子数组）的第一个元素作为基准，并将所有其他元素与它进行比较，那么当 \(A_1, \dots, A_n\) 已经是排好序的时候就会出现最坏情况。在这种情况下，一个子数组为空，另一个子数组有 \(n-1\) 个元素。每一步都只是将子数组的大小减一。因此，\(W(n) = \Theta(\sum_{i=1}^{n} i) = \Theta(n^2)\)。 为了估计平均情况下的运行时间，我们假设初始序列是随机排列的，这样每次比较都相当于在 \(A_1, \dots, A_n\) 中随机选择两个元素进行比较。

补充说明： 这里“随机排列”的意思是说，原始数据中元素的大小顺序完全是随机的，没有特定的升序或降序模式。这样，我们选任何元素作为基准，比它大和比它小的元素数量大致各一半的概率更高。

如果我们怀疑原始元素序列具有某种模式，我们总是可以先随机置换该序列，其复杂度为 \(\Theta(n)\)，正如下一个问题中所解释的。设 \(\tilde{A}_p\) 和 \(\tilde{A}_q\) 是最终排序数组中的第 \(p\) 个和第 \(q\) 个元素（\(1 \le p < q \le n\)）。在 \(\tilde{A}_p\) 和 \(\tilde{A}_q\) 之间有 \(q-p+1\) 个数字。\(\tilde{A}_p\) 和 \(\tilde{A}_q\) 被比较的概率是：\(\tilde{A}_p\) 在 \(\tilde{A}_q\) 之前与 \(\tilde{A}_{p+1}, \dots, \tilde{A}_{q-1}\) 中的任何一个被选为基准进行比较的概率。

如果这个基准介于它们之间，则 \(\tilde{A}_p\) 和 \(\tilde{A}_q\) 会被分到不同的子数组中，从而不会再被比较了。这个概率是： \[ P(p,q) = \frac{2}{q-p+1} \] （你可以再次使用对称性论证来推导出这个概率）。总的期望比较次数为 \(A(n) = \sum_{q=2}^{n} \sum_{p=1}^{q-1} P(p,q) = \sum_{q=2}^{n} \sum_{p=1}^{q-1} \left( \frac{2}{q-p+1} \right) = \Theta(n \lg n)\)。尽管理论上快速排序在最坏情况下的性能可能比归并排序慢，但它通常与归并排序一样快，甚至更快。

随机排列

A. 如果你有一个能生成离散或连续均匀分布随机数的随机数生成器，你如何洗一副 52 张牌，使得每种排列的可能性都相等？ 解法：一种对 \(n\) 个元素进行置换的简单算法是基于排序的随机置换。它为每张牌分配一个随机数，然后根据分配的随机数对牌进行排序。 根据对称性，每种可能的顺序（在 \(n!\) 种可能的有序序列中）都是等可能发生的。其复杂度由排序步骤决定，因此 如果我们使用连续均匀分布，理论上任何两个随机数相等的概率都为零。

查看本页原始扫描图

运行时间为 \(\Theta(n \log n)\)。对于较小的 \(n\)，比如一副牌中的 \(n=52\)，其复杂度 \(\Theta(n \log n)\) 是可以接受的。对于较大的 \(n\)，我们可能希望使用一种被称为 Knuth 洗牌法的更快的算法。对于 \(n\) 个元素 \(A[1], \dots, A[n]\)，Knuth 洗牌法使用以下循环来生成一个随机排列：

for (i=1 to n) swap(A[i], A[Random(i, n)]),

其中 Random(\(i, n\)) 是一个从 \(i\) 到 \(n\) 之间离散均匀分布中抽取的随机数。Knuth 洗牌法的复杂度为 \(\Theta(n)\)，并且有一个直观的解释。在第一步中，由于牌号是从 1 到 \(n\) 的离散均匀分布中选取的，所以 \(n\) 张牌中的每一张被选为第一张牌的概率都相等；在第二步中，剩下的 \(n-1\) 张牌中的每一张被选为第二张牌的概率也相等；以此类推。所以很自然地，每个有序序列出现的概率都是 \(1/n!\)。B. 你有一个由字符组成的文件。文件中的字符可以顺序读取，但文件的长度未知。如何从中选取一个字符，使得文件中的每个字符被选中的概率都相等？解法：让我们从选择第一个字符开始。如果有第二个字符，我们以 1/2 的概率保留第一个字符，并以 1/2 的概率将选择替换为第二个字符。

如果有第三个字符，我们以 2/3 的概率保留当前的选择（来自前两个字符），并以 1/3 的概率将选择替换为第三个字符。这个过程一直持续到最后一个字符。换句话说，设 \(C_n\) 是我们在扫描了 \(n\) 个字符并且存在第 \(n+1\) 个字符时所选择的字符，那么保留当前选择的概率是 \(\frac{n}{n+1}\)，而切换到第 \(n+1\) 个字符的概率是 \(\frac{1}{n+1}\)。使用简单的归纳法，我们可以很容易地证明，如果总共有 \(m\) 个字符，每个字符被选中的概率都是 \(1/m\)。

搜索算法

A. 设计一个算法，用不超过 \(3n/2\) 次比较来找到 \(n\) 个数中的最小值和最大值。解法：对于一个未排序的 \(n\) 个数的数组，要找出数组的最小值或最大值需要进行 \(n-1\) 次比较。然而，要同时找出最小值和最大值，最多只需要进行 \(3n/2\) 次比较。如果我们将元素分成 \(n/2\) 对，比较每一对中的元素，并将较小的那个放入 A 组，较大的那个放入 B 组。这一步需要进行 \(n/2\) 次比较。由于整个数组的最小值必定在 A 组中，最大值必定在 B 组中，我们只需要在 A 组中找出最小值，在 B 组中找出最大值即可，这两步各需要 \(n/2 - 1\) 次比较。所以总的比较次数最多为 \(3n/2\)。

5脚注 5如果 \(n\) 是奇数，需要稍作调整，但上界 \(3n/2\) 仍然适用。 B. 给你一个数字数组。从数组开头到某个位置，所有元素都是零；在这个位置之后，所有元素都是非零的。如果你不知道数组的大小，如何找到第一个非零元素的位置？解法：我们可以从第一个元素开始检查；如果它是零，我们检查第二个元素；如果第二个元素是零，我们检查第四个元素……这个过程重复进行，直到第 \(i\) 步时，第 \(2^i\) 个元素是非零的。然后，我们检查第 \(\frac{2^i + 2^{i-1}}{2}\) 个元素。如果它是零，搜索范围就限制在第 \(\frac{2^i + 2^{i-1}}{2}\) 个元素到第 \(2^i\) 个元素之间；

否则，搜索范围就限制在第 \(2^{i-1}\) 个元素到第 \(\frac{2^i + 2^{i-1}}{2}\) 个元素之间……每次，我们将范围减半。这种方法本质上是一种二分查找。如果第一个非零元素位于位置 \(n\)，则该算法的复杂度为 \(\Theta(\log n)\)。C. 你有一个数字组成的方形网格。每行中的数字从左到右递增。每列中的数字从上到下递增。设计一个算法在网格中查找给定的数字。你的算法复杂度是多少？解答： 设 A 是一个 \(n \times n\) 矩阵，代表数字网格，\(x\) 是我们要在网格中查找的数字。搜索从最后一列从上到下开始：\(A_{1,n}, A_{2,n}, \dots, A_{n,n}\)。如果找到该数字，则停止搜索。如果 \(A_{n,n} < x\)，则 \(x\) 不在网格中，搜索也停止。

如果 \(A_{i,n} < x < A_{i+1,n}\)，那么我们知道第 \(1, \dots, i\) 行中的所有数字都小于 \(x\)，因此也被排除。6脚注 6\(i\) 可以为 0，这意味着 \(x < A_{1,n}\)，在这种情况下，我们可以从右向左搜索第一行。 然后我们从右向左搜索第 \((i+1)\) 行。如果数字在第 \((i+1)\) 行中找到，搜索停止。如果 \(A_{i+1,j+1} > x\)，则 \(x\) 不在网格中，因为第 \(i+1\) 行及以上所有行的数字都大于 \(x\)。

如果 \(A_{i+1,j+1} > x > A_{i+1,j}\)，我们排除第 \(j+1, \dots, n\) 列中的所有数字。然后我们可以沿着列从 \(A_{i+1,j}\) 向 \(A_{n,j}\) 搜索，直到找到 \(x\)（或者 \(x\) 不在网格中）。

查看本页原始扫描图

（网格）或一个 $k$ 使得 $A_{k,j} < x < A_{k+1,j}$，然后我们沿着第 $k+1$ 行从 $A_{k+1,1}$ 向 $A_{k+1,j}$ 搜索。使用该算法，搜索最多需要 $2n$ 步。因此，其复杂度为 $O(n)$。

补充说明： 这个算法的核心思想是“逐步缩小搜索范围”。每次比较都能排除掉矩阵的一部分（一行或一列），从而将搜索空间线性减少。最坏情况下，我们需要检查大约 $2n$ 个元素（例如，从右上角开始，先向下走 $n$ 步，再向左走 $n$ 步），因此时间复杂度是 $O(n)$。

Fibonacci 数列

考虑以下计算 Fibonacci 数列的 C++ 程序：

int Fibonacci(int n)
{
if (n <= 0)
return 0;
else if (n==1)
return 1;
else
return Fibonacci (n-1)+Fibonacci (n-2) ;
}

如果对于某个较大的 $n$，计算 Fibonacci($n$) 需要 100 秒，那么计算 Fibonacci($n+1$) 需要多长时间，精确到秒？该算法有效率吗？你会如何计算 Fibonacci 数列？解答： 这个 C++ 函数使用一种相当低效的递归方法来计算 Fibonacci 数列。Fibonacci 数列定义如下递推关系： \[ F_0 = 0, F_1 = 1, F_n = F_{n-1} + F_{n-2}, \forall n \ge 2 \] $F_n$ 有一个闭式解 $F_n = \frac{(1+\sqrt{5})^n - (1-\sqrt{5})^n}{2^n\sqrt{5}}$，这可以通过归纳法轻松证明。从函数来看，可以清楚地看出： \[ T(0) = 1, T(1) = 1, T(n) = T(n-1) + T(n-2) + 1 \] 因此，运行时间与 Fibonacci 数列的序列成正比。

对于一个较大的 $n$，$(1-\sqrt{5})^n \to 0$，所以 \[ \frac{T(n+1)}{T(n)} \approx \frac{\sqrt{5}+1}{2} \] 如果计算 Fibonacci($n$) 需要 100 秒，那么计算 Fibonacci($n+1$) 的时间为 $T(n+1) \approx \frac{\sqrt{5}+1}{2} T(n) \approx 162$ 秒。\[ \phi = \frac{\sqrt{5}+1}{2} \] 被称为黄金分割率。

查看本页原始扫描图

递归算法具有指数复杂度 \( \Theta\left(\left(\frac{\sqrt{5}+1}{2}\right)^n\right) \)，这显然是低效的。原因在于它未能有效利用 Fibonacci 数列中较小 $n$ 对应的 Fibonacci 数的信息。如果我们使用定义按顺序计算 $F_0, F_1, \dots, F_n$，运行时间的复杂度为 $ \Theta(n) $。一种称为递归平方的算法可以进一步将复杂度降低到 $ \Theta(\log n) $。

由于 \[ \begin{bmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} F_n & F_{n-1} \\ F_{n-1} & F_{n-2} \end{bmatrix} \] 且 \[ \begin{bmatrix} F_2 & F_1 \\ F_1 & F_0 \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} , \] 我们可以用归纳法证明 \[ \begin{bmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^n \] 令 $ A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} $，我们可以再次应用分治范式来计算 $ A^n $： \[ A^n = \begin{cases} A^{n/2} \times A^{n/2}, & \text{如果 } n \text{ 是偶数} \\ A^{(n-1)/2} \times A^{(n-1)/2} \times A, & \text{如果 } n \text{ 是奇数} \end{cases} \] 两个 $ 2 \times 2 $ 矩阵的乘法复杂度为 $ \Theta(1) $。

所以 $ T(n) = T(n/2) + \Theta(1) $。应用主定理，我们得到 $ T(n) = \Theta(\log n) $。

补充说明： 递归平方法的关键在于利用矩阵乘法的结合律，将计算 $A^n$ 的问题分解为计算 $A^{n/2}$ 的问题，从而将线性次数的递归调用（如原始递归）减少为对数次数。因为每次递归调用时，问题规模（指数 $n$）减半，所以总的时间复杂度是对数级的。

最大连续子数组

假设你有一个长度为 $n$ 的一维数组 $A$，其中包含正数和负数。设计一个算法来找到 $A$ 的任意连续子数组 $A[i, j]$ 的最大和：$ V(i, j) = \sum_{x=i}^j A[x], 1 \le i \le j \le n $。 解答： 几乎所有的交易系统都需要这样的算法来计算实际交易账簿或模拟策略的最大上涨或最大回撤。因此，这是面试官非常喜欢问的算法问题，尤其是对冲基金和交易部门的面试官。 最显而易见的算法是 $O(n^2)$ 算法，它使用以下方程从头开始顺序计算 $V(i, j)$： 当 $j=i$ 时，$V(i, i) = A[i]$；当 $j>i$ 时，$V(i, j) = \sum_{x=i}^j A[x] = V(i, j-1) + A[j]$。 在计算 $V(i, j)$ 的同时，我们也跟踪 $V(i, j)$ 的最大值以及对应的子数组索引 $i$ 和 $j$。

查看本页原始扫描图

一种更有效的方法使用了分治范式。我们定义 \[ T(i) = \sum_{x=1}^{i} A[x] \quad \text{且} \quad T(0) = 0, \] 则 $ V(i, j) = T(j) - T(i-1), \forall 1 \le i \le j \le n $。显然，对于任意固定的 $ j $，当 $ T(i-1) $ 最小时，$ V(i, j) $ 最大。因此，以 $ j $ 结尾的最大子数组和为 $ V_{\max} = T(j) - T_{\min} $，其中 $ T_{\min} = \min(T(1), \dots, T(j-1)) $。如果我们随着 $ j $ 的增加跟踪并更新 $ V_{\max} $ 和 $ T_{\min} $，我们可以开发出以下 $ O(n) $ 算法：

T = A[1]; V_max = A[1]; T_min = min(0, T)
For j = 2 to n
{
  T = T + A[j];
  If T - T_min > V_max then V_max = T - T_min;
  If T < T_min then T_min = T;
}
Return V_max;

下面是一个相应的 C++ 函数，它在给定数组及其长度的情况下返回 $ V_{\max} $ 以及索引 $ i $ 和 $ j $：

double maxSubarray(double A[], int len, int &i, int &j)
{
  double T = A[0], Vmax = A[0];
  double Tmin = min(0.0, T);
  for (int k = 1; k < len; ++k)
  {
    T += A[k];
    if (T - Tmin > Vmax) { Vmax = T - Tmin; j = k; }
    if (T < Tmin) { Tmin = T; i = (k + 1 < j) ? (k + 1) : j; }
  }
  return Vmax;
}

将它应用于以下数组 A，

double A[] = {1.0, 2.0, -5.0, 4.0, -3.0, 2.0, 6.0, -5.0, -1.0};
int i = 0, j = 0;

查看本页原始扫描图

double Vmax = maxSubarray(A, sizeof(a)/sizeof(A[1]), i, j);

将得到 $V_{\max} = 9$，$i=3$ 和 $j=6$。所以子数组是 [4.0, -3.0, 2.0, 6.0]。

补充说明： 这个 $O(n)$ 算法通常被称为 Kadane 算法 或其变体。其核心思想是动态规划：我们不需要为每个可能的子数组计算和，而是维护两个关键变量：T 代表到当前位置的累积和，Tmin 代表到目前为止遇到的最小累积和。当前最大子数组和就是当前累积和减去之前的最小累积和。这样只需遍历数组一次。

7.2. 2 的幂

世界上只有 10 种人——懂二进制的和不懂的。如果你碰巧听懂了这句笑话，你可能知道计算机使用二进制（基数为 2）数字系统运行。与十进制数字 0-9 不同，每个比特（二进制数字）只有两个可能的值：0 和 1。数字的二进制表示具有一些有趣的特性，这些特性在实践中被广泛探索，也使其成为面试中一个有趣的测试主题。

2 的幂

你如何判断一个整数是否是 2 的幂？ 解：任何形如 \(x = 2^n\)（其中 \(n \ge 0\)）的整数，其二进制表示中只有一位是 1（从右数第 \(n+1\) 位）。例如，8（即 \(2^3\)）的二进制表示为 0...01000。同样容易看出，\(2^n - 1\) 的二进制表示中，从右数的 \(n\) 位全部为 1。例如，7 的二进制表示为 0...00111。因此，\(2^n\) 和 \(2^n - 1\) 没有任何共同的二进制位。于是，通过检查 \(x \& (x-1) == 0\)（其中 \(\&\) 是按位与运算符），就可以简单判断整数 \(x\) 是否是 2 的幂。

乘以 7

给出一种快速方法，在不使用乘法（\(*\)）运算符的情况下，将一个整数乘以 7？ 解：\((x << 3) - x\)，其中 \(<<\) 是左移位运算符。\(x << 3\) 等价于 \(x * 8\)。因此 \((x << 3) - x\) 就是 \(x * 7\)。8脚注 8如果 \(<<\) 导致溢出，结果可能出错。

概率模拟

给定一枚公平的硬币。你能设计一个使用这枚公平硬币的简单游戏，使得你的获胜概率为 \(p\)（其中 \(0 < p < 1\)）吗？9脚注 9提示：计算机存储的是二进制数而非十进制数；二进制数中的每一位都可以用一枚公平硬币来模拟。

查看本页原始扫描图

解： 解决这个问题的关键在于认识到，\(p \in (0,1)\) 也可以表示为一个二进制数，而二进制数的每一位都可以用一枚公平硬币来模拟。首先，我们可以将概率 \(p\) 表示为二进制数： \[ p = 0.p_1 p_2 \cdots p_n = p_1 2^{-1} + p_2 2^{-2} + \cdots + p_n 2^{-n}, \quad p_i \in \{0,1\}, \forall i = 1, 2, \dots, n. \] 然后，我们开始抛掷公平硬币，并将正面记为 1，反面记为 0。设 \(s_i \in \{0,1\}\) 为第 \(i\) 次抛掷的结果，从 \(i=1\) 开始。每次抛掷后，我们将 \(p_i\) 与 \(s_i\) 进行比较。如果 \(s_i < p_i\)，我们就获胜，并且停止抛掷硬币。如果 \(s_i > p_i\)，我们就失败，并且停止抛掷硬币。

如果 \(s_i = p_i\)，我们继续抛掷更多硬币。某些 \(p\) 值（例如，1/3）表示为二进制数时是无穷级数（1/3 \(\to\) 0.010101...）。在这些情况下，随着 \(i\) 增加，最终达到 \(s_i \neq p_i\) 的概率为 1。如果序列是有限的（例如，1/4 = 0.01），并且我们到达最后一步时 \(s_n = p_n\)，那么我们就失败（例如，对于 1/4，只有序列 00 被归类为获胜；其他三个序列 01、10 和 11 都被归类为失败）。这样的模拟将使我们以概率 \(p\) 获胜。

毒酒问题

你为生日派对准备了 1000 瓶葡萄酒。在派对开始前 20 小时，酒庄紧急通知你，其中一瓶酒被下毒了。你恰好有 10 只实验室老鼠可以用来测试一瓶酒是否有毒。这种毒药非常强烈，任何剂量都会在恰好 18 小时内毒死一只老鼠。但在第 18 小时死亡之前，没有任何其他症状。在派对开始前，你能用这 10 只老鼠有把握地找出那瓶有毒的葡萄酒吗？解： 如果老鼠可以依次进行测试，每次排除一半的瓶子，那么这个问题就变成了一个简单的二分查找问题。10 只老鼠最多可以从 1024 瓶葡萄酒中找出有毒的那一瓶。但不幸的是，由于症状要在 18 小时后才会出现，而我们只有 20 小时的时间，所以不能依次测试老鼠。尽管如此，二分查找的思路仍然适用。1 到 1000 之间的所有整数都可以用 10 位二进制格式表示。

例如，瓶 1000 可以标记为 1111101000，因为 \(1000 = 2^9 + 2^8 + 2^7 + 2^6 + 2^5 + 2^3\)。现在，让老鼠 1 从每个第一位（最右边的低位）为 1 的瓶子里喝一小口；让老鼠 2 从每个第二位为 1 的瓶子里喝一小口；……；最后，让老鼠 10 从每个第 10 位（最高位）为 1 的瓶子里喝一小口。18 小时后，如果我们把老鼠从最高位到最低位排成一排，并将活着的老鼠记为 0，死去的老鼠记为 1，那么就可以很容易地反推出有毒瓶子的标签。例如，如果第 6、第 7 和第 9 只老鼠死了，其他的都活着，那么排列得到的序列是 0101100000，有毒瓶子的标签是 \(2^8 + 2^6 + 2^5 = 352\)。

查看本页原始扫描图

7.3 数值方法

许多金融工具的价格没有封闭形式的解析解。这些金融工具的估值依赖于多种数值方法。在本节中，我们将讨论蒙特卡洛模拟和有限差分法的应用。

蒙特卡洛模拟

蒙特卡洛模拟是一种通过使用具有适当概率的随机数作为输入来迭代评估确定性模型的方法。对于衍生品定价，它模拟标的资产的大量价格路径，其概率对应于标的随机过程（通常在风险中性测度下），计算每条路径的衍生品贴现收益，并对贴现收益取平均值以得出衍生品价格。蒙特卡洛模拟的有效性依赖于大数定律。 蒙特卡洛模拟可用于估计衍生品价格，如果其收益仅取决于标的资产的最终值，并且可以改编用于估计价格，如果收益也依赖于路径。尽管如此，它不能直接应用于美式期权或任何其他具有提前行权选项的衍生品。

A. 解释如何使用蒙特卡洛模拟为欧式看涨期权定价？

解答： 如果我们假设股票价格遵循几何布朗运动，我们可以模拟可能的股票价格路径。我们可以将时间 \(t\) 和 \(T\) 之间的时间分割成 \(N\) 个等距的时间步长。10脚注 10对于欧式期权，我们可以简单地设置 \(N=1\)。但对于更一般的期权，特别是路径依赖的期权，我们希望有较小的时间步长，因此 \(N\) 应该较大。 所以 \(\Delta t = \frac{T-t}{N}\) 并且 \(t_i = t + \Delta t \times i\)，对于 \(i=0,1,2,\dots,N\)。

然后我们使用方程 \[ S_i = S_{i-1} e^{(r-\sigma^2/2)(\Delta t) + \sigma\sqrt{\Delta t}\epsilon_i} \] 来模拟风险中性概率下的股票价格路径，其中 \(\epsilon_i\) 是来自标准正态分布的独立同分布随机变量。假设我们模拟了 \(M\) 条路径，每条路径在到期日 \(T\) 产生一个股票价格 \(S_{T,k}\)，其中 \(k=1,2,\dots,M\)。

查看本页原始扫描图

欧式看涨期权的估计价格是期望收益的现值， \[ \sum_{k=1}^{M} \max(S_{T,k} - K, 0) \] 其计算公式为 \[ C = e^{-r(T-t)} \frac{\sum_{k=1}^{M} \max(S_{T,k} - K, 0)}{M} \] B. 如果你的计算机只能生成在 0 到 1 之间服从连续均匀分布的随机变量，你如何生成服从 \(N(\mu, \sigma^2)\)（均值为 \(\mu\)、方差为 \(\sigma^2\) 的正态分布）的随机变量？ 解答：这是一个很好的问题，用于测试对随机数生成（蒙特卡洛模拟的基础）的基本知识。这个问题的解答可以分为两步：

1. 使用逆变换法或拒绝法，从均匀随机数生成器中生成服从 \(x \sim N(0,1)\) 的随机变量。
2. 将 \(x\) 缩放为 \(\mu + \sigma x\)，以生成最终服从 \(N(\mu, \sigma^2)\) 的随机变量。

第二步很简单；第一步需要一些解释。一种流行的生成随机变量的方法是逆变换法：对于任意连续随机变量 \(X\)，其累积分布函数为 \(F\)（\(U = F(X)\)），则随机变量 \(X\) 可以定义为 \(U\) 的反函数：\(X = F^{-1}(U)\)，其中 \(0 \le U \le 1\)。显然 \(X = F^{-1}(U)\) 是一个一一对应的函数，其定义域为 \(0 \le U \le 1\)。因此，任何连续随机变量都可以通过以下过程生成：

• 从标准均匀分布中生成一个随机数 \(u\)。
• 计算满足 \(u = F(x)\) 的值 \(x\)，作为由分布 \(F\) 描述的随机数。

要使该模型有效，必须能够计算 $F^{-1}(U)$。对于标准正态分布， \[ U = F(X) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-t^2/2} dt \] 其反函数没有解析解。理论上，我们可以通过数值积分方法（如欧拉方法）求解常微分方程 $F'(x) = f(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$，从而得到 $X$ 到 $U$ 的一一映射作为数值解。

11脚注 11为了对 $y = F(x)$ 进行积分，已知其一阶导数 $y' = f(x)$ 和一个已知的初始值 $y_0 = F(x_0)$，欧拉方法选择一个小的步长 $h$（$h$ 可以为正或负）来顺序逼近 $y$ 值： \[ y_{i+1} = y_i + h f(x_i) \] 其中 $x_{i+1} = x_i + h$。 然而，这种方法比拒绝抽样法效率低：

补充说明： 欧拉方法的基本思想是“以直代曲”。从已知点 $(x_0, y_0)$ 出发，沿着该点切线的方向（斜率由 $f(x_0)$ 给出）前进一小步 $h$，得到下一个点的近似值 $(x_1, y_1)$，如此反复。对于标准正态分布的累积分布函数，可以从 $F(0)=0.5$ 这个初始值开始计算。

查看本页原始扫描图

某些随机变量的概率密度函数为 $f(x)$，但其累积分布函数的反函数 $F^{-1}(U)$ 没有解析解。在这些情况下，我们可以使用一个概率密度函数为 $g(y)$ 的随机变量，并通过 $Y = G^{-1}(U)$ 来生成具有概率密度函数 $f(x)$ 的随机变量。假设存在一个常数 $M$，使得对于所有的 $y$，都有 $\frac{f(y)}{g(y)} \le M$。 我们可以实现以下接受-拒绝方法：

• 采样步骤： 从 $g(y)$ 中生成随机变量 $y$，并从标准均匀分布 $[0,1]$ 中生成随机变量 $v$。
• 接受/拒绝步骤： 如果 $v \le \frac{f(y)}{Mg(y)}$，则接受 $x = y$；否则，重复采样步骤。12脚注 12$P(X \le x) \propto \int_{-\infty}^x g(y) \frac{f(y)}{Mg(y)} dy = M \int_{-\infty}^x f(y) dy \Rightarrow F(x) = \frac{P(X \le x)}{P(X < \infty)} = \int_{-\infty}^x f(y) dy$

一个参数为 $\lambda=1$ 的指数分布随机变量（$g(x) = \lambda e^{-\lambda x}$）的累积分布函数为 $u = G(x) = 1 - e^{-x}$。因此，其反函数有解析解 $x = -\log(1-u)$，可以方便地模拟指数分布的随机变量。对于标准正态分布，其概率密度函数为 $f(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$。

\[ \frac{f(x)}{g(x)} = \frac{\frac{1}{\sqrt{2\pi}}e^{-x^2/2}}{\frac{2}{\sqrt{\pi}}e^{-x^2/2}} < \frac{\frac{2}{\sqrt{\pi}}e^{-(x-1)^2/2+1/2}}{\frac{2}{\sqrt{\pi}}e^{-x^2/2}} \approx 1.32, \quad \forall 0 < x < \infty \] 因此，我们可以选择 $M=1.32$，并使用接受-拒绝方法生成 $x \sim N(0,1)$ 的随机变量，然后将其缩放到 $N(\mu, \sigma^2)$ 的随机变量。C. 您能否解释几种方差缩减技术，以提高蒙特卡洛模拟的效率？

解答： 蒙特卡洛模拟的基本形式是独立同分布随机变量的均值 \[ Y_1, Y_2, \dots, Y_M: \quad \overline{Y} = \frac{1}{M} \sum_{i=1}^M Y_i \] 由于每个 $Y_i$ 的期望值是无偏的，因此估计量 $\overline{Y}$ 也是无偏的。如果 $Var(Y) = \sigma$ 并且我们生成独立同分布的 $Y_i$，那么 $Var(\overline{Y}) = \sigma / \sqrt{M}$，其中 $M$ 是模拟次数。不出所料，蒙特卡洛 $F(x_0+h) = F(x_0) + f(x_0) \times h, F(x_0+2h) = F(x_0+h) + f(x_0+h) \times h, \dots$。标准正态分布累积分布函数的初始值可以是 $F(0) = 0.5$。

查看本页原始扫描图

如果 $\sigma$ 很大，模拟在计算上会非常密集。通常需要成千上万甚至数百万次模拟才能获得所需的精度。根据具体问题，已经应用了各种方法来减少方差。对偶变量 (Antithetic variable): 对于每一系列 $\epsilon_i$'s，计算其对应的支付 $Y(\epsilon_1, \dots, \epsilon_N)$。然后，将所有 $\epsilon_i$'s 的符号反转，计算对应的支付 $Y(-\epsilon_1, \dots, -\epsilon_N)$。当 $Y(\epsilon_1, \dots, \epsilon_N)$ 和 $Y(-\epsilon_1, \dots, -\epsilon_N)$ 负相关时，方差会减小。矩匹配 (Moment matching): 随机变量的特定样本可能与总体分布不匹配。

我们可以先抽取大量样本，然后重新缩放样本，使样本的矩（均值和方差是最常用的）与期望的总体矩相匹配。控制变量 (Control variate): 如果我们想为衍生品 $X$ 定价，并且存在一个具有解析解的密切相关的衍生品 $Y$，我们可以生成一系列随机数，并使用相同的随机序列为 $X$ 和 $Y$ 定价，得到 $\hat{X}$ 和 $\hat{Y}$。然后，可以通过 $\hat{X} + (Y - \hat{Y})$ 来估计 $X$。本质上，我们使用 $(Y - \hat{Y})$ 来纠正 $\hat{X}$ 的估计误差。

重要性采样 (Importance sampling): 为了从分布 $f(x)$ 估计 $h(x)$ 的期望值，我们不直接从分布 $f(x)$ 中抽取 $x$，而是从分布 $g(x)$ 中抽取 $x$，并使用蒙特卡洛模拟来估计期望值 \[ E_{f(x)}[h(x)] = \int h(x) f(x) dx = \int \frac{h(x) f(x)}{g(x)} g(x) dx = E_{g(x)} \left[ \frac{h(x) f(x)}{g(x)} \right] \] 如果 $\frac{h(x) f(x)}{g(x)}$ 的方差比 $h(x)$ 小，那么重要性采样可以得到一个更有效的估计量。这个方法最好通过一个深度价外期权作为例子来解释。如果我们直接使用风险中性 $f(S_T)$ 作为分布，大多数模拟路径将得到 $h(S_T)=0$，结果是估计方差会很大。

如果我们引入一个分布 $g(S_T)$，它具有更宽的跨度（$S_T$ 的尾部更肥），那么更多的模拟路径将具有正的 $h(S_T)$。缩放因子 $\frac{f(x)}{g(x)}$ 将保持估计量无偏，但该方法将具有较低的方差。

13脚注 13重要性采样本质上是使用测度变换的方差缩减方法。

查看本页原始扫描图

低差异序列 (Low-discrepancy sequence)

不使用随机样本，我们可以生成一个确定性的“随机变量”序列来代表分布。这种低差异序列可能使收敛速度达到 $1/M$。

D. 如果一个期权没有闭式定价公式，你如何估计它的 delta 和 gamma？

解答： 正如我们在问题 A 中讨论过的，无论是否有闭式定价公式，期权的价格都可以通过蒙特卡洛模拟推导出来。同样的方法也可以用来估计 delta 和 gamma，只需将当前标的资产价格从 $S$ 略微改变到 $S \pm \delta S$，其中 $\delta S$ 是一个小的正值。对三个起始价格 $S-\delta S$、$S$ 和 $S+\delta S$ 分别运行蒙特卡洛模拟，我们将得到它们对应的期权价格 $f(S-\delta S)$、$f(S)$ 和 $f(S+\delta S)$。

估计的 delta：$\Delta = \frac{\delta f}{\delta S} = \frac{f(S+\delta S) - f(S-\delta S)}{2\delta S}$ 估计的 gamma：$\Gamma = \frac{\left(f(S+\delta S) - f(S)\right) - \left(f(S) - f(S-\delta S)\right)}{\delta S^2}$ 为了减小方差，通常最好使用相同的随机数序列来估计 $f(S-\delta S)$、$f(S)$ 和 $f(S+\delta S)$。14脚注 14如果支付函数不连续，此方法可能效果不佳。

补充说明： 这种方法本质上是数值微分中的中心差分法。使用相同的随机数序列可以确保对三个不同价格的模拟是在完全相同的“随机场景”下进行的，从而消除随机性带来的额外差异，使得价格差主要反映参数变化的影响，提高了希腊字母估计的精度。

E. 如何使用蒙特卡洛模拟来估计 $\pi$？

解答： 估计 $\pi$ 是蒙特卡洛模拟的一个经典例子。一种标准的估计 $\pi$ 的方法是：在单位正方形内随机选取点（x 和 y 是 0 到 1 之间的独立均匀随机变量），并确定落在圆 $x^2 + y^2 \leq 1$ 内的点的比例。为简化起见，我们关注第一象限。如图 7.1 所示，任何落在圆内的点都满足方程 $x_i^2 + y_i^2 \leq 1$。

圆内点的比例与其面积成正比： \[ \hat{p} = \frac{\text{满足 } x_i^2 + y_i^2 \leq 1 \text{ 的点 }(x_i, y_i) \text{ 的数量}}{\text{正方形内点 }(x_i, y_i) \text{ 的总数}} = \frac{1/4\pi}{1 \times 1} = \frac{1}{4}\pi \Rightarrow \hat{\pi} = 4\hat{p}. \] 因此，我们生成大量独立的 (x, y) 点，估计落在圆内的点占正方形内点的比例，然后将该比例乘以 4，即可得到 $\pi$ 的估计值。图 7.1 仅使用了 1000 个点进行说明。以今天的

查看本页原始扫描图

计算能力，我们可以轻松生成数百万个 (x, y) 点对，以良好的精度估计 $\pi$。在笔记本电脑上，使用 Matlab 进行 1000 次模拟，每次模拟使用 1,000,000 个 (x, y) 点，耗时不到 1 分钟，得到的 $\pi$ 平均估计值为 3.1416，标准差为 0.0015。

有限差分法 (Finite difference method)

有限差分法是另一种流行的衍生品定价数值技术。它通过离散化时间和标的资产价格，数值求解微分方程来估计衍生品价格。我们可以将布莱克-斯科尔斯-默顿方程（一个二阶非线性偏微分方程）转化为热扩散方程（如我们在第6章中所做）。这个新方程表示为 \(\tau\)（到期时间）和 \(x\)（标的资产价格的函数）的函数，是衍生品的一般微分方程。不同衍生品之间的差异在于边界条件。通过构建 \(x\) 和 \(\tau\) 的网格，并利用边界条件，我们可以使用有限差分法递归地计算每个 \(x\) 和 \(\tau\) 处的 \(u\)。

A. 你能简要解释有限差分法吗？

解答：实践中使用的有限差分法有几种版本。让我们简要介绍显式差分法、隐式差分法和

查看本页原始扫描图

克兰克-尼科尔森法。如图7.2所示，如果我们将 \(\tau\) 的范围 \([0, T]\) 划分为 \(N\) 个离散区间，步长为 \(\Delta \tau = T/N\)，并将 \(x\) 的范围 \([x_0, x_J]\) 划分为 \(J\) 个离散区间，步长为 \(\Delta x = (x_J - x_0)/J\)，那么时间 \(\tau\) 和空间 \(x\) 可以表示为网格 \(\tau_n\)（\(n=1, \dots, N\)）和 \(x_j\)（\(j=1, \dots, J\)）。

显式差分法在时间 \(\tau_n\) 处使用前向差分，在 \(x_j\) 处使用二阶中心差分： \[ \frac{\partial u}{\partial \tau} \approx \frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta \tau} = \frac{u_{j+1}^{n} - 2u_{j}^{n} + u_{j-1}^{n}}{(\Delta x)^2} \approx \frac{\partial^2 u}{\partial x^2}. \] 整理各项，我们可以将 \(u_j^{n+1}\) 表示为 \(u_{j-1}^n\)、\(u_j^n\) 和 \(u_{j+1}^n\) 的线性组合： \[ u_j^{n+1} = \alpha u_{j-1}^n + (1-2\alpha)u_j^n + \alpha u_{j+1}^n, \quad \text{其中 } \alpha = \Delta \tau / (\Delta x)^2. \] 此外，我们通常有边界条件 \(u_0^n\)、\(u_J^n\) 和 \(u_j^0\)，适用于所有 \(n=1, \dots, N\)；

\(j=0, \dots, J\)。结合边界

查看本页原始扫描图

条件和方程 \(u_{j}^{n+1} = \alpha u_{j-1}^{n} + (1-2\alpha)u_{j}^{n} + \alpha u_{j+1}^{n}\)，我们可以估计网格上所有 \(u_{j}^{n}\) 的值。隐式差分法在时间 \(t_{n+1}\) 处使用后向差分，在 \(x_j\) 处使用二阶中心差分： \[ \frac{\partial u}{\partial \tau} \approx \frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta \tau} = \frac{u_{j+1}^{n+1} - 2u_{j}^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} \approx \frac{\partial^2 u}{\partial x^2}. \] 克兰克-尼科尔森法在时间 \((t_n + t_{n+1})/2\) 处使用中心差分，在 \(x_j\) 处使用二阶中心差分： \[ \frac{\partial u}{\partial \tau} \approx \frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta \tau} = \frac{1}{2} \left( \frac{u_{j+1}^{n} - 2u_{j}^{n} + u_{j-1}^{n}}{(\Delta x)^2} + \frac{u_{j+1}^{n+1} - 2u_{j}^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} \right) \approx \frac{\partial^2 u}{\partial x^2}. \]

B. 如果你使用显式有限差分法求解抛物线型偏微分方程，时间维度步数过多和空间维度步数过多，哪个更糟糕？

解答：显式有限差分法中 \(u_{j}^{n+1}\) 的方程是 \(u_{j}^{n+1} = \alpha u_{j-1}^{n} + (1-2\alpha)u_{j}^{n} + \alpha u_{j+1}^{n}\)，其中 \(\alpha = \Delta t/(\Delta x)^2\)。为了使显式有限差分法稳定，我们需要满足 \(1-2\alpha > 0 \Rightarrow \Delta t/(\Delta x)^2 < 1/2\)。因此，较小的 \(\Delta t\)（即较多时间步数）是可取的，但较小的 \(\Delta x\)（即过多空间步数）可能导致 \(\Delta t/(\Delta x)^2 > 1/2\)，从而使结果不稳定。从这个意义上说，空间维度步数过多更糟糕。相比之下，隐式差分法始终是稳定且收敛的。