Question

　备忘　

本节要对我们将在本章解决的等式问题给予一个深刻的理解。本节的理解对于本章剩余部分的阅读并不是严格必须的。如果你想要跳过本节，那么请一定要看一下例6-1和例6-2中的算法逻辑，理解我们需要优化的代码。

另外，如果你阅读了本节之后还需要更多的解释，请阅读剑桥大学出版社出版的William Press等人的著作《数值分析方法库》第三版第17章。

为了在本章探索矩阵和矢量计算，我们将重复使用流体扩散的例子。扩散是流体移动并试图均匀混合的一种方式。

本节我们将探索扩散等式背后的数学问题。它看上去可能很复杂，但不要担心！我们会快速简化并让它更容易理解。另外，重要的是要记住，虽然对于我们所需要解决的最终等式的基本理解有助于本章的阅读，但这并不是完全必须的。后续的章节会主要集中于代码中的各种方程式，而不是扩散等式本身。不过，对等式的理解会有助于你的直觉来对代码进行优化。总而言之——理解你代码背后的动机以及算法的难点会让你对可能的优化方案有一个更深刻的理解。

扩散的一个简单例子是水中的染料：如果你在室温的水里滴入几滴染料，它会慢慢扩散直到完全跟水混在一起。由于我们并不搅拌水，水也没有热到出现热对流的程度，扩散将是两种液体混合的主要过程。为了在数学上解决这些等式，我们会特地选取我们所需要的初始条件，之后我们会看看对初始条件的改变会对计算发生怎样的影响（图6-2）。

说了这么多，对于我们的目标来说其实最主要的还是扩散方程式本身。扩散等式可以被写成一个1阶偏微分方程：

这个方程中，u是表示扩散质量的矢量。比如，我们会有一个矢量，里面的值0表示纯水，1表示纯染料（0和1之间的值表示混合液体）。一般来说，这个矢量应该是一个2阶或3阶的矩阵表示液体的实际区域或体积。这样，为了表示一个玻璃杯中的液体，我们就可以设u为一个3阶矩阵，计算所有的轴，而不只是在x方向求导。除此以外，D是一个物理值表示模拟实验中的液体的属性。D越大则表明液体越容易扩散。为了简化，我们在代码中设D = 1，但依然将它包含在计算中。

　备忘　

扩散方程也被称为热方程。此时，u表示一个区域的温度而D描述了材料的热传导能力。解开方程可以告诉我们热如何传导。这样，我们就能够了解CPU产生的热量如何扩散到散热片上而不是水中染料的扩散。

我们会将时空上连续的扩散微分方程改成小块的离散时空。我们用欧拉法来做到这一点。欧拉法将微分方程写成差分形式，如下：

此时dt是一个固定的值，表示一个时间的步长，或者说我们为了解决方程而定的时间的精确度。它可以被想象成我们正在制作的电影的帧率。当帧率提高（或者说dt变小）时，我们就能够得到一个更清晰的图像。实际上，dt越接近0，欧拉近似就越精确（但是注意，这个精确只能在理论上达到，因为计算机只拥有有限的精确度，而精度误差会迅速传播到所有的计算结果）。这样我们就能将方程改写成根据u(x, t)求u(x, t+dt)。这意味着我们可以以某个初始状态u(x, 0)为起点，表示一杯刚加了一滴染料的水，然后通过对初始状态进行“演化”来了解这杯水在未来某个时刻的样子（u(x, dt)）。这种类型的问题被称为“初值问题”或“柯西问题”。

以同样的有限差分近似对x进行求导，我们就能获得最终的方程式：

这里，和dt表示帧率一样，dx表示图像的分辨率——dx越小，则我们矩阵的每一个单元格表示的区域就越小。为了简化，我们设D = 1且dx = 1。这两个值在物理模拟时非常重要，但是因为我们只是以扩散方程为例子来说明问题，它们对我们并不重要。

我们能用这个等式解决几乎所有扩散问题。但是，这个等式有些需要考虑的地方。首先，我们之前说过u方向的空间索引（参数x）代表矩阵中的索引。那么当x处于矩阵的起始索引位置时x-dx的值代表什么呢？这个问题被称为边界条件。你可以这样定义固定的边界条件：“任何超出矩阵边界的值都为0”（或任意其他值）。或者，你也可以定义边界值环绕的周期性边界条件（比如，如果你矩阵的某个维度的长度为N，那么索引为-1的值就等于索引为N-1的值，而索引为N的值则等于索引为0的值。换句话说，当你试图访问索引i上的值时，你得到的就是位于索引(i%N)上的值）。

另一个需要考虑的地方是我们如何存储u的多个时间分量。我们可以让一个矩阵来保存计算所需的所有时间值。看起来我们需要至少两个矩阵：一个保存液体的当前状态，一个保存液体的下一状态。我们会看到这里面有很重要的性能考量。

那么我们如何解决实际的问题？例6-1包含了一些伪代码来说明我们是如何使用这个等式来解决问题的。

例6-1　1阶扩散方程伪代码

# Create the initial conditions
u = vector of length N
for i in range(N):
  u = 0 if there is water, 1 if there is dye

# Evolve the initial conditions
D = 1
t = 0
dt = 0.0001
while True:
  print "Current time is: %f" % t
  unew = vector of size N

  # Update step for every cell
  for i in range(N):
    unew[i] = u[i] + D * dt * (u[(i+1)%N] + u[(i-1)%N] - 2 * u[i])
  # Move the updated solution into u
  u = unew

  visualize(u)

这段代码会根据水中染料的初始条件来告诉我们每一个0.0001秒之后的未来系统是什么样的。结果见图6-1，图中显示的是最浓的染料液滴（高帽函数）随时间的演化。我们可以看到染料是如何均匀混合，直到每处染料的浓度都相等的。

为了本章的目的，我们还将解决方程的2阶版本。这意味着我们将要操作一个2阶矩阵，而不是一个矢量（或者说一个只有1阶的矩阵）。对方程唯一的改动（以及相应的代码改动）是现在我们必须计算第二个y方向上的导数。这意味着原始的方程会变成：

图6-1　1阶扩散的例子

这个2阶扩散方程翻译成伪码如例6-2所示，用的是跟之前一样的方法。

例6-2　计算2阶差分的算法

for i in range(N):
  for j in range(M):
    unew[i][j] = u[i][j] + dt * (          \
       (u[(i+1)%N][j] + u[(i-1)%N][j] - 2 * u[i][j]) + \ # d^2 u / dx^2
       (u[i][(j+1)%M] + u[j][(j-1)%M] - 2 * u[i][j]) \ # d^2 u / dy^2
    )

现在我们可以将所有的代码片段组装到一起，并写出完整的Python2阶扩散代码来作为我们本章接下来进行性能参照的基准。代码虽然看上去更复杂，结果却近似1阶扩散（见图6-2）。

如果你希望阅读更多本节的主题，可以参考扩散方程的维基百科和S.V.Gurevich写的第7章“复杂系统的数值方法”。

图6-2　两组初始条件扩散的例子