微分方程的数值处理

在这一章中，我们将研究常微分方程（ODEs）和偏微分方程（PDEs）的数值解。这些方程在物理学中常用来描述一些现象，如飞机周围的空气流动，或桥梁在各种压力下的弯曲。虽然这些方程通常相当简单，但从它们中得到具体的数字（”如果有一百辆汽车在上面，这座桥会下垂多少”）则比较复杂，往往需要大型计算机来产生所需的结果。这里我们将描述将ODEs和PDEs转化为可计算问题的技术。

我们将首先介绍「初值问题」（IVPs），它描述了随着时间变换的过程。此处我们仅考虑常微分方程：只依赖于时间的标量函数。接下来，我们将研究描述空间过程的「边界值问题」（BVPs）。边界值问题通常涉及到多个空间变量，因此我们可以得到偏微分方程。

最后，我们将考虑 “热方程”，这是一个「初始边界值问题」（IBVP），它同时具有IVP和BVP的特点：它描述了热量在一个物理物体（如杆）上的传播。初始值描述了初始温度，而边界值给出了杆两端的规定温度。

我们在这一章的目的是展示一类重要的计算问题的起源。因此，我们不会去讨论解的存在性、唯一性或条件的理论问题。关于这一点，见[99]或任何专门讨论ODE或PDE的书。为了便于分析，我们还将假设所有涉及的函数都存在任意高阶导数，并且处处光滑。

初值问题

许多物理现象随着时间的推移而变化，物理学定律给出了对变化的描述，而非数值本身的描述。例如，牛顿第二定律 $$ F=ma $$ 是一个关于点质量的位置变化的公式：表示为 $$ a(t)=\frac{d^{2}}{d t^{2}} x(t)=F / m $$ 它指出，加速度线性地取决于施加在质量上的力。对于质量的位置，可以通过分析得出一个封闭式的描述$𝑥(𝑡)=…$，但在许多情况下需要某种形式的近似或数值计算。这也被称为 “数值积分“。

牛顿第二定律是一个常微分方程，因为它表述了变量关于时间变化的函数。此外，牛顿第二定律也是一个初值问题，因为它给定了初始条件随着时间变化的情况。该方程是二阶的，如果我们引入向量，则可以通过引入双分量向量$u$，结合位置$x$和速度$x’$将上述方程降为一阶微分方程： $$ u(t)=\left(x(t), x^{\prime}(t)\right)^{t} $$ 以$u$表示的牛顿方程就变成了 $$ u^{\prime}=A u+B, \quad A=\left(\begin{array}{ll} 0 & 1 \ 0 & 0 \end{array}\right), \quad B=\left(\begin{array}{c} 0 \ F / a \end{array}\right) $$ 为了简单起见，在本课程中我们将只考虑标量方程；那么我们的参考方程是 $$ u^{\prime}(t)=f(t, u(t)), \quad u(0)=u_{0}, \quad t>0 $$ 方程允许过程有明确的时间依赖性，但一般来说，我们只考虑没有这种明确依赖性的方程，即所谓的 “自治（autonomous） “ODE，其形式为 $$ u’(t) = f(u(t)) $$ 其中右手边并不明确地依赖于$t$。

注释 13 非自治ODE可以转化为自治ODE，所以这并不是什么限制。如果$u=u(t)$是一个标量函数，并且$f=f(t, u)$，我们定义$u2(t)= t$，并考虑等同的自治系统$\left(\begin{array}{l} u^{\prime} \ u{2}^{\prime} \end{array}\right)=\left(\begin{array}{c} f\left(u_{2}, u\right) \ 1 \end{array}\right)$

通常情况下，在某个起点（通常选择$𝑡=0$）的初始值是给定的：$𝑢(0)=𝑢_0$，对于某个值$𝑢_0$，我们对$𝑢$随着$𝑡\rightarrow \infty$的行为感兴趣。举例来说，$f(x)=x$给出的方程是：$𝑢′(𝑡)= 𝑢(𝑡)$。这是一个简单的人口增长模型：该方程指出，增长速度等于人口规模。对于某些$f$的选择，方程（4.2）可以用分析法解决，但我们不会考虑这个问题。相反，我们只考虑数值解和这个过程的准确性。

在数值方法中，我们考虑用离散大小的时间步长来近似解决连续的时间依赖过程。由于这引入了一定量的误差，我们将分析在每个时间步长中引入的误差，以及这些误差是如何叠加成一个整体误差的。在某些情况下，限制全局误差的需要会对数值方案施加限制。

误差和稳定性

由于机器运算会产生不精确性，我们希望尽量避免因为初始值的小扰动而造成的干扰。因此，如果与不同初值$𝑢_0$相关的解在$t\rightarrow\infty$时相互收敛，我们将称微分方程为 「稳定」的。

稳定性的一个充分标准是 $$ \frac{\partial}{\partial u} f(u)=\left{\begin{array}{ll}

0 & \text { 不稳定 } \ =0 & \text { 中立稳定 } \

<0 & \text { 稳定 } \end{array}\right. $$ 证明：设 $u^$为$f$ 零点，即：$f(u^)=0$，则常函数$u(t)\equiv u^$为 $u’=f(u)$的一个解，即所谓的“平衡解”。平衡状态下，微小的扰动是不会影响系统的稳定性的。例如：设 $u$ 是PDE的一个解，写作 $u(t)= u^∗ + \eta(t)$，那么我们有 $$ \begin{aligned} \eta^{\prime} &=u^{\prime}=f(u)=f\left(u^{}+\eta\right)=f\left(u^{}\right)+\eta f^{\prime}\left(u^{}\right)+O\left(\eta^{2}\right) \ &=\eta f^{\prime}\left(u^{*}\right)+O\left(\eta^{2}\right) \end{aligned} $$ 忽略二阶项，存在一个解： $$ \eta(t)=e^{f^{\prime}\left(x^{\star}\right) t} $$ 这意味着，如果$f^′(x) < 0$，扰动将被阻尼，如果$f^′(x) > 0$，扰动将被放大。

我们经常会提到简单的例子$f(u)=-\lambda u$，其解$u(t)=u_{0} e^{-\lambda t}$。如果$\lambda>0$，这个问题是稳定的。

有限差分近似法：欧拉显式和隐式方法

为了数值解决这个问题，我们通过研究有限时间/空间步长，将连续问题变成离散问题。假设所有的函数都足够平滑，一个简单的泰勒级数展开就可以得到。 $$ u(t+\Delta t)=u(t)+u^{\prime}(t) \Delta t+u^{\prime \prime}(t) \frac{\Delta t^{2}}{2 !}+u^{\prime \prime \prime}(t) \frac{\Delta t^{3}}{3 !}+\cdots $$ 这就得到了$u’$ $$ \begin{aligned} u^{\prime}(t) &=\frac{u(t+\Delta t)-u(t)}{\Delta t}+\frac{1}{\Delta t}\left(u^{\prime \prime}(t) \frac{\Delta t^{2}}{2 !}+u^{\prime \prime \prime}(t) \frac{\Delta t^{3}}{3 !}+\cdots\right) \ &=\frac{u(t+\Delta t)-u(t)}{\Delta t}+\frac{1}{\Delta t} O\left(\Delta t^{2}\right) \ &=\frac{u(t+\Delta t)-u(t)}{\Delta t}+O(\Delta t) \end{aligned} $$ 如果所有的导数都是有界的，我们可以用一个$O(\Delta t^2)$来近似高阶导数的无限之和。或者，你可以证明这个和等于$\Delta t^{2} u^{\prime \prime}(t+\alpha \Delta t)$ 当$0<\alpha<1$ 。我们看到，我们可以通过有限差分来近似微分算子，其误差是已知的，其数量级是时间步长的函数。

将其代入 $u′ = f(t, u)$，得到 $$ \frac{u(t+\Delta t)-u(t)}{\Delta t}=f(t, u(t))+O(\Delta t) \ $$ 或者 $$ u(t+\Delta t)=u(t)+\Delta t f(t, u(t))+O\left(\Delta t^{2}\right) $$ 注释 14 前面的两个方程是数学上的等式，不应该被理解为对一个给定的函数$u′$进行计算的方法。回顾前面的讨论，你可以看到这样的公式对于小$\Delta t$来说很快就会被取消。关于数值微分的进一步讨论超出了本书的范围，请参见任意一本标准的数值分析教科书。

我们现在使用上述方程来推导一个数值方案：在$t0=0，t{k+1}=tk+\Delta t=\cdots=(k+1)\Delta t$的情况下，我们得到一个差分方程 $$ u{k+1}=u_k+\Delta t f(t_k, u_k) $$ 为$u_k$量，我们希望$u_k$将是对$u(t_k)$的良好近似。这就是所谓的 “显式欧拉 “或 “欧拉正向 “方法。

从微分方程到差分方程的过程通常被称为离散化（discretization），因为我们只在离散的点集中计算函数值。计算的数值本身仍然是实值的。另一种说法是：如果我们计算$𝑘$个时间步长，就可以在有限的二维空间$\mathbb{R}^k$中找到数值解。原问题的解是在$\mathbb{R}\rightarrow \mathbb{R}$的函数空间中找到的。

在上面，我们用一个算子逼近另一个算子，这样做的时候，当$\Delta t \downarrow 0$时，截断误差为$O(\Delta t)$（见附录14对这个数量级的符号的更正式介绍）。这并不意味着差分方程计算出的解就接近于真解。为此，还需要进行一些分析。

我们从分析 “局部误差 “开始：如果假设计算出的解在$𝑘$步骤是精确的，即$𝑢𝑘=𝑢(𝑡𝑘)$，那么$𝑘+1$步时将会出现怎样的错误? 我们有 $$ \begin{aligned} u\left(t{k+1}\right) &=u\left(t{k}\right)+u^{\prime}\left(t{k}\right) \Delta t+u^{\prime \prime}\left(t{k}\right) \frac{\Delta t^{2}}{2 !}+\cdots \ &=u\left(t{k}\right)+f\left(t{k}, u\left(t{k}\right)\right) \Delta t+u^{\prime \prime}\left(t{k}\right) \frac{\Delta t^{2}}{2 !}+\cdots \end{aligned} $$ 和 $$ u{k+1}=u_k+f(t_ku_k)\Delta t $$ 于是 $$ \begin{aligned} L{k+1} &=u{k+1}-u\left(t{k+1}\right)=u{k}-u\left(t{k}\right)+f\left(t{k}, u{k}\right)-f\left(t{k}, u\left(t{k}\right)\right)-u^{\prime \prime}\left(t{k}\right) \frac{\Delta t^{2}}{2 !}+\cdots \ &=-u^{\prime \prime}\left(t{k}\right) \frac{\Delta t^{2}}{2 !}+\cdots \end{aligned} $$ 这表明，在每一步中，我们的误差为$O(\Delta t^2)$。如果我们假设这些误差可以相加，我们发现全局误差为 $$ E{k} \approx \Sigma{k} L_{k}=k \Delta t \frac{\Delta t^{2}}{2 !}=O(\Delta t) $$ 由于全局误差在$\Delta t$中是一阶的，我们称之为 “一阶方法”。需要注意的是，这个误差（衡量真实解和计算解之间的距离）与截距误差（即算子的近似误差）是同阶的$O(\Delta t)$。

欧拉显式方法的稳定性

考虑IVP $u′ = f(t,u)$ 对于 $t \geqslant 0$, 其中$f(t,u)= -\lambda u$，给出初始值$u(0) = u_0$。存在一个精确的解，即$u(t)=u_0e^{-\lambda t}$。从上面的讨论中可知这个问题是稳定的，也就是说，如果$\lambda >0$，解的小扰动最终会被抑制。现在我们将研究数值解的表现是否与精确解相同，也就是说，数值解是否也收敛于零。

这个问题的欧拉正向或显式欧拉方案是 $$ u{k+1}=u{k}-\Delta t \lambda u{k}=(1-\lambda \Delta t) u{k} \Rightarrow u{k}=(1-\lambda \Delta t)^{k} u{0} $$ 为了稳定，我们要求$uk\rightarrow 0$，因为$k\rightarrow \infty$。这就相当于 $$ \begin{aligned} u{k} \downarrow 0 & \Leftrightarrow|1-\lambda \Delta t|<1 \\ & \Leftrightarrow-1<1-\lambda \Delta t<1 \\ & \Leftrightarrow-2<-\lambda \Delta t<0 \\ & \Leftrightarrow 0<\lambda \Delta t<2 \\ & \Leftrightarrow \Delta t<2 / \lambda \end{aligned} $$ 我们看到，数值求解方案的稳定性取决于$\Delta t$的值：只有当$\Delta t$足够小时，该方案才是稳定的。为此，我们称显式欧拉方法为条件稳定的。请注意，微分方程的稳定性和数值方案的稳定性是两个不同的问题。如果$\lambda >0$，连续问题是稳定的；数值问题有一个额外的条件，取决于所用的离散化方案。

请注意，我们刚刚进行的稳定性分析是专门针对微分方程$u′=-\lambda u$的。如果你要处理的是不同的IVP，你必须进行单独的分析。一般情况下，显式方法通常会给出条件稳定性。

欧拉隐式方法

刚才的显式方法很容易计算，但条件稳定性是一个潜在的问题。它可能意味着时间步骤的数量将是一个限制性因素。有一个替代显式方法的方法，不会受到同样的限制。

与其扩展$u(t + \Delta t)$，不如考虑以下对$u(t - \Delta t)$的扩展 $$ u(t-\Delta t)=u(t)-u^{\prime}(t) \Delta t+u^{\prime \prime}(t) \frac{\Delta t^{2}}{2 !}+\cdots $$ 这意味着 $$ u^{\prime}(t)=\frac{u(t)-u(t-\Delta t)}{\Delta t}+u^{\prime \prime}(t) \Delta t / 2+\cdots $$ 如前所述，我们取方程$u’(t)=f(t, u(t))$，用差分公式近似计算$u′(t)$。 $$ \frac{u(t)-u(t-\Delta t)}{\Delta t}=f(t, u(t))+O(\Delta t) \Rightarrow u(t)=u(t-\Delta t)+\Delta t f(t, u(t))+O\left(\Delta t^{2}\right) $$ 我们再次定义固定点$𝑡𝑘=𝑘𝑡$，并定义一个数值方案： $$ u{k+1}=uk+\Delta t f(t{k+1, u_{k+1}}) $$ 其中$u_k$是$u(t_k)$的近似值。

与显式方案的一个重要区别是，$u{k+1}$现在也出现在方程的右侧。也就是说，$u{k+1}$的计算现在是隐含的。例如，让$f(t,u)=-u^3$，那么$u{k+1}=u_k-\Delta tu^3{k+1}$。换句话说，$u_{k+1}$是方程$\Delta tx=u_k$的解。这是一个非线性方程，通常可以用牛顿法求解。

隐式欧拉方法的稳定性

让我们再看看这个例子$f(t, u(t))=-\lambda u$。用隐式方法计算，可以得到 $$ u{k+1}=u{k}-\lambda \Delta t u{k+1} \Leftrightarrow(1+\lambda \Delta t) u{k+1}=u{k} $$ 那么 $$ u{k+1}=\left(\frac{1}{1+\lambda \Delta t}\right) u{k} \Rightarrow u{k}=\left(\frac{1}{1+\lambda \Delta t}\right)^{k} u_{0} $$ 如果$\lambda>0$，这是一个稳定方程的条件，我们发现，对于所有的$u_k$和$\Delta t$的值，$\lambda \rightarrow 0$。这种方法被称为无条件稳定。与显式方法相比，隐式方法的一个优势显然是稳定性：可以采取较大的时间步长而不用担心非物理行为。当然，大的时间步长会使收敛到稳定状态（见附录15.4）变慢，但至少不会出现发散。

另一方面，隐式方法更加复杂。正如你在上面看到的，它们可能涉及到在每个时间步长中需要解决的非线性系统。在$u$是矢量值的情况下，比如下面讨论的热方程，我们会发现隐式方法需要解决一个方程组。

练习 4.1 分析以下IVP $𝑢′(x) = f(x)$ 方案的准确性和计算性。 $$ u{i+1}=u{i}+h\left(f\left(x{i}\right)+f\left(x{i+1}\right)\right) / 2 $$ 相当于把欧拉显式和隐式方案加在一起。我们不需要分析这个方案的稳定性。

练习4.2 考虑初值问题 $y′(t) = y(t)(1-y)$。请注意，$y\equiv 0$和$y\equiv 1$是解决方案。这些被称为 “平衡解”。

1. 一个解决方案是稳定的，如果扰动 “收敛到解决方案”，意味着对于$\varepsilon$足够小。 $$ \text{if}\quad y(t) = \varepsilon \quad \text{for some }𝑡, \text{then} \lim{t\rightarrow \infty} y(t) = 0 $$ 并且 $$ \text{if}\quad y(t) = 1 + \varepsilon \quad \text{for some }𝑡, \text{then} \lim{t\rightarrow \infty} y(t) = 1 $$ 这就要求，例如 $$ y(t)=\varepsilon \Rightarrow y′(t)<0. $$ 0是一个稳定的解决方案吗？1是吗？

2. 考虑显式方法 $$ y{k+1}=y_k+\Delta ty_k(1-y_k) $$ 用于计算微分方程的数值解。说明 $$ y_k \in (0,1)\Rightarrow y{k+1}>yk,\quad y_k>1\Rightarrow y{k+1}<y_k $$

3. 编写一个小程序来研究在不同的$\Delta t$的选择下数值解的行为。在你提交的作业中包括程序清单和一些运行的情况。

4. 通过运行你的程序，你发现数值解会出现振荡。请对$\Delta t$提出一个条件，使数值解呈单调性。只要说明$yk<1\Rightarrow y{k+1}<1$，以及$y_k>1\Rightarrow y_{k+1}>1$即可。

5. 现在考虑隐式方法 $$ y{k+1}-\Delta ty{k+1}(1-y{k+1})=y_k $$ 并说明$y{k+1}$可以从$y_k$计算出来。编写一个程序，并研究在不同的$\Delta t$选择下的数值解的行为。

6. 显示出对所有$\Delta t$的选择，隐式方案的数值解都是单调的。

边界值问题

在上一节中，我们看到了初值问题，它模拟的是随时间变化的现象。现在我们将转向边界值问题，一般来说，边界值问题在时间上是静止的，但它描述的是与位置有关的现象。例如，桥梁在负载下的形状，或窗玻璃中的热量分布，因为外面的温度与里面的不同。

二阶一维BVP的一般形式是 $$ u^{\prime \prime}(x)=f\left(x, u, u^{\prime}\right) \text { 对于 } x \in[a, b] \text { 其中 } u(a)=u{a}, u(b)=u{b} $$ 但这里我们只考虑简单的形式 $$ -u^{\prime \prime}(x)=f(x) \text { for } x \in[0,1] \text { with } u(0)=u{0}, u(1)=u{1} $$ 在一个空间维度上，或 $$ -u{x x}(\bar{x})-u{y y}(\bar{x})=f(\bar{x}) \text { for } \bar{x} \in \Omega=[0,1]^{2} \text { with } u(\bar{x})=u_{0} \text { on } \delta \Omega $$ 在两个空间维度上。这里，$\delta \Omega$是域$\Omega$的边界。由于我们在边界上规定了$u$的值，这样的问题被称为边界值问题（Boundary Value Problem，BVP）。

注释 15 边界条件可以更普遍，涉及区间端点上的导数。这里我们只看Dirichlet边界条件，它规定了域的边界上的函数值。

一般性PDE理论

有几种类型的PDE，每种都有不同的数学特性。最重要的属性是影响区域：如果我们对问题进行修补，使解决方案在一个点上发生变化，那么还有哪些点会受到影响。

双曲方程

PDEs的形式为 $$ A u{x x}+B u{y y}+\text { lower order terms }=0 $$ $A，B$的符号相反。这样的方程描述的是波，或更一般的对流现象，是保守的，不倾向于稳定状态。

直观地说，在任何一点上改变波浪方程的解，只会改变未来的某些点，因为波有一个传播速度，使得一个点不可能影响到空间上太远的近期的点。这种类型的PDE将不会在本书中讨论。

抛物线方程

PDEs的形式为 $$ A u{x}+B u{y y}+\text { no highter terms in }x=0 $$ 并且它们描述了类似于扩散的现象；这些现象往往趋于稳定状态。描述它们的最好方法是考虑在空间和时间的每一点上的解决方案受到空间的每一点上的某个有限区域的影响。

注释 16 这导致了一个限制IBVP时间的条件，即所谓的Courant-Friedrichs-Lewy条件http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition。它描述了这样一个概念：在精确问题中，$u(x, t)$取决于$ u(x′, t-\Delta t)$的数值范围；数值方法的时间步长必须小到足以让数值解考虑到所有这些点。

热方程是抛物线类型的标准例子。

椭圆方程

PDEs的形式为 $$ A u{x x}+B u{y y}+\text { lower order terms }=0 $$ 其中$A，B>0$；它们通常描述已经达到稳定状态的过程，例如抛物线问题中的$𝑡 \rightarrow \infty$。它们的特点是所有的点都相互影响。这些方程经常描述结构力学中的现象，如梁或膜。直观地讲，压下膜上的任何一点都会改变其他每一点的高度，无论多么微小。泊松方程（4.2.2节）是这种类型的标准例子。

一维空间的泊松方程

算子$\Delta$ $$ \Delta u = u{xx}+u{yy}, $$ 是二阶微分算子，方程（4.7）是二阶PDE。具体来说，问题是 $$ -\Delta u=-u{x x}(\bar{x})-u{y y}(\bar{x})=f(\bar{x}) \text { for } \bar{x} \in \Omega=[0,1]^{2} \text { with } u(\bar{x})=u_{0} \text { on } \delta \Omega $$ 被称为泊松方程（Poisson equation），定义在单位平方上。二阶PDEs相当常见，它们描述了流体和热流以及结构力学中的许多现象。

首先，为了简单起见，我们考虑一维泊松方程 $$ -u_{xx}=f(x). $$ 下面我们考虑的是二维的情况；然后扩展到三维的情况。

为了找到一个数值方案，我们像以前一样使用泰勒级数，用$𝑢(𝑥+h)$和$u(x-h)$来表示。

$u$及其在$x$的导数。设$h>0$，则 $$ u(x+h)=u(x)+u^{\prime}(x) h+u^{\prime \prime}(x) \frac{h^{2}}{2 !}+u^{\prime \prime \prime}(x) \frac{h^{3}}{3 !}+u^{(4)}(x) \frac{h^{4}}{4 !}+u^{(5)}(x) \frac{h^{5}}{5 !}+\cdots $$ 以及 $$ u(x-h)=u(x)-u^{\prime}(x) h+u^{\prime \prime}(x) \frac{h^{2}}{2 !}-u^{\prime \prime \prime}(x) \frac{h^{3}}{3 !}+u^{(4)}(x) \frac{h^{4}}{4 !}-u^{(5)}(x) \frac{h^{5}}{5 !}+\cdots $$ 我们现在的目标是对$u’’(x)$进行近似。我们看到，这些方程中的$u’$项在加法下会被抵消，剩下$2u(x)$。 $$ u(x+h)=u(x)+u^{\prime}(x) h+u^{\prime \prime}(x) \frac{h^{2}}{2 !}+u^{\prime \prime \prime}(x) \frac{h^{3}}{3 !}+u^{(4)}(x) \frac{h^{4}}{4 !}+u^{(5)}(x) \frac{h^{5}}{5 !}+\cdots $$ 于是 $$ -u^{\prime \prime}(x)=\frac{2 u(x)-u(x+h)-u(x-h)}{h^{2}}+u^{(4)}(x) \frac{h^{2}}{12}+\cdots $$ 那么，上述数值方案的基础是观察 $$ \frac{2 u(x)-u(x+h)-u(x-h)}{h^{2}}=f\left(x, u(x), u^{\prime}(x)\right)+O\left(h^{2}\right) $$ 这表明我们可以用一个差分算子来近似微分算子，在$h\downarrow 0$时有一个$O(h^2)$。截断误差为$h\downarrow 0$。

为了得出一种数值方法，我们将区间[0, 1]划分为等距的点。$xk=k_h$ 其中$h=1/(n+1)，k=0 … n + 1$. 有了这些，有限差分（FD）公式（4.9）导致了一个形成方程组的数值方案。 $$ -u{k+1}+2uk-u{k-}=h^2f(x_k)\quad \text{for } k=1,…,n $$ 这种使用FD公式近似解决PDE的过程被称为有限差分法（Finite Difference Method，FDM）。

对于大多数的$k$值，这个方程将$uk$未知数与$u{k-1}$和$u{k+1}$的未知数联系起来。例外情况是$k=1$和$k=n$。在这种情况下，我们记得$u_0$和$u{n+1}$是已知的边界条件，我们把左边是未知数，右边是已知量的方程写为 $$ \left{\begin{array}{l} -u{i-1}+2 u{i}-u{i+1}=h^{2} f\left(x{i}\right) \ 2 u{1}-u{2}=h^{2} f\left(x{1}\right)+u{0} \ 2 u{n}-u{n-1}=h^{2} f\left(x{n}\right)+u{n+1} \end{array}\right. $$ 我们现在可以将这些方程总结为$uk ，k = 1 … n- 1$的矩阵方程。 $$ \left(\begin{array}{cccc} 2 & -1 & & \ -1 & 2 & -1 & \ & \ddots & \ddots & \ddots \end{array}\right)\left(\begin{array}{c} u{1} \ u{2} \ \vdots \end{array}\right)=\left(\begin{array}{c} h^{2} f{1}+u{0} \ h^{2} f{2} \ \vdots \end{array}\right) $$

其形式为$Au=f$，$A$为完全已知矩阵，$f$为完全已知向量，$u$为未知向量。请注意，右边的向量在第一个和最后一个位置有问题的边界值。这意味着，如果你想用不同的边界条件解决同一个微分方程，只有向量$f$会发生变化。

练习 4.3 $u(0)= u_0$这种类型的条件被称为狄利克雷边界条件。在物理学上，这相当于，例如，知道一个棒端点的温度。其他边界条件也存在。如果我们对流体流动进行建模，并且已知$x=0$时的流出率，那么为导数指定一个值，$u’(0)=u_0’$，而不是为函数值指定一个值，是合适的。这就是所谓的诺伊曼边界条件。诺伊曼边界条件 $u’(0) = u_0’$ 可以通过以下方式来模拟 $$ \frac{u_0-u_1}{h}=u_0’. $$ 证明与狄利克雷边界条件的情况不同，这影响了线性系统的矩阵。

证明在两端都有诺伊曼边界条件会产生一个奇异矩阵，因此线性系统没有唯一的解。(提示：猜测特征值为零的向量）。

在物理学上这是有意义的。例如，在一个弹性问题中，狄利克雷边界条件说明杆子被夹在一定的高度；诺伊曼边界条件只说明它在端点的角度，这使得它的高度无法确定。

让我们列举一些$𝐴$的属性，你会发现这些属性与解决此类方程组有关。

• 矩阵非常稀疏：非零元素的百分比很低，零元素不是随机分布的，而是位于主对角线周围的一个带状结构中。在一般情况下，我们称之为带状矩阵（banded matrix），而在这个特定情况下称之为三对角矩阵（tridiagonal matrix）。带状结构是典型的PDEs，但不同应用中的稀疏矩阵可能不太规则。

• 矩阵非对称。这一特性并不涉及由BVP分解而来的所有矩阵，但如果没有奇数阶（指第一、第三、第五……）导数，例如$ux，u{xxx}，u_{xy}$。

• 矩阵元素在每个对角线上是恒定的，也就是说，在每一组点${(𝑖,𝑗)∶𝑖 - 𝑗 = 𝑐}$，对于某个𝑐。这只对非常简单的问题是正确的。如果微分方程$\frac{d}{d x}\left(a(x) \frac{d}{d x} u(x)\right)$。如果我们在区间内设置h变量，它也不再成立，例如，因为我们想更详细地模拟左端点周围的行为。

• 矩阵元素符合以下符号模式：对角线元素是正的，而非对角线元素是非正的。这一属性取决于所使用的数字方案，但它通常是真实的。连同下面的确定性属性，这被称为$M-$矩阵。关于这些矩阵有一整套数学理论[12]。

• 矩阵是正定的：$x^tAx>0$，适用于所有非零向量$x$。如果仔细选择数值方案，这一属性将从原始连续问题中继承下来。虽然这一点的用途目前看来并不明确，但以后你会看到取决于它的线性系统的求解方法。

严格地说，方程的解很简单：$u=A^{-1}f$。然而，计算$A^{-1}$并不是找到$u$的最好方法。正如刚才所观察到的，矩阵$A$只有$3N$个非零元素可以存储。另一方面，它的逆矩阵则没有一个非零元素。虽然我们不会证明这一点，但这种说法对大多数稀疏矩阵都是成立的。因此，我们希望以一种不需要储存$O(n^2)$的方式来解决$Au=f$。

练习 4.4 你将如何解决这个三对角方程组？证明系数矩阵的LU因子化给出了双对角矩阵形式的因子：它们有一个非零对角线和正好一个非零子对角线或超对角线。解决三对角方程组的总操作数是多少？一个向量与这样的矩阵相乘的运算次数是多少？这种关系并不典型！

二维空间的泊松方程

上面的一维BVP在很多方面是不典型的，特别是与由此产生的线性代数问题有关。在本节中，我们将研究二维泊松问题。你会看到，它构成了一维问题的非微观概括。三维的情况与二维的情况非常相似，所以我们将不讨论它。

上面的一维问题有一个函数$u=u(x)$，现在变成了二维的$u=u(x,y)$。那么我们感兴趣的二维问题就是 $$ -u{xx}-u{yy}=f, \quad (x,y)\in [0,1]^2, $$ 其中边界上的数值是给定的。我们通过在$x$和$y$方向上应用方程得到我们的离散方程 $$ \begin{array}{l} -u{x x}(x, y)=\frac{2 u(x, y)-u(x+h, y)-u(x-h, y)}{h^{2}}+u^{(4)}(x, y) \frac{h^{2}}{12}+\cdots \ -u{y y}(x, y)=\frac{2 u(x, y)-u(x, y+h)-u(x, y-h)}{h^{2}}+u^{(4)}(x, y) \frac{h^{2}}{12}+\cdots \end{array} $$ 或者说，合起来看 $$ 4 u(x, y)-u(x+h, y)-u(x-h, y)-u(x, y+h)-u(x, y-h)=1 / h^{2} f(x, y)+O\left(h^{2}\right) $$ 令 $h=1/(n+1)$，定义$xi=ih$，$y_j=jh$；让$u{ij}$是对$u(x_i,y_j)$的近似，那么我们的离散方程为

$$ 4 u{i j}-u{i+1, j}-u{i-1, j}-u{i, j+1}-u{i, j-1}=h^{-2} f{i j} $$ 我们现在有$n\times n$未知数$u{ij}$。为了像以前一样将其转化为一个线性系统，我们需要将它们放在一个线性排序中，我们通过定义$I = I{ij} = j+ i\times n$来实现。这被称为字典序（lexicographic ordering），因为它将坐标$(i, j)$当作字符串来排序。

使用这个排序，我们可以得到$N=n^2$的方程 $$ 4 u{I}-u{I+1}-u{I-1}-u{I+n}-u{I-n}=h^{-2} f{I}, \quad I=1, \ldots, N $$

而线性系统看起来像 $$ A=\left(\begin{array}{ccccc|ccccc|c} 4 & -1 & & & \varnothing & -1 & & & & \varnothing & \ -1 & 4 & -1 & & & & -1 & & & & \ & \ddots & \ddots & \ddots & & & & \ddots & & & \ & & \ddots & \ddots & -1 & & & & \ddots & & \ \varnothing & & & -1 & 4 & \varnothing & & & & -1 & \ \hline-1 & & & & \varnothing & 4 & -1 & & & -1 & \ & -1 & & & & -1 & 4 & -1 & & & -1 \ & \uparrow & \ddots & & & \uparrow & \uparrow & \uparrow & & & \uparrow \ & k-n & & & & k-1 & k & k+1 & & -1 & k+n \ & & & & -1 & & & & -1 & 4 & & \ \hline & & & & & \ddots & & & & \ddots & \end{array}\right) $$ 矩阵的大小为$N\times N$，其中$N=n^2$。与一维的情况一样，我们看到BVP会产生一个稀疏的矩阵。

以矩阵形式考虑这个线性方程组似乎是很自然的事情。然而，以明确未知数的二维联系的方式呈现这些方程可能更有洞察力。为此，图4.1展示了领域中的变量，以及方程如何通过有限差分模板将它们联系起来。从现在开始，在制作这样的领域图片时，我们将只使用变量的索引，而省略 “$u$ “标识符。

方程矩阵和一维的情况一样是带状的，但是和一维的情况不同，带状的内部有零点。因为该矩阵有五个非零对角线，所以被称为五对角线结构。

你也可以在矩阵上放一个块状结构，把在域的一行中的未知数组合起来。这被称为分块矩阵（block matrix），在块的层面上，它有一个三对角的矩阵结构，所以我们称之为分块三对角矩阵（block tridiagonal matrix）。请注意，对角线块本身是三对角的；非对角线块是减去单位矩阵的。

这个矩阵和上面的一维例子一样，有恒定的对角线，但这又是由于问题的简单性质造成的。在实际问题中，这不会是真的。也就是说，这样的 “恒定系数 “问题是有的，当它们在矩形域上时，有非常有效的方法来解决线性系统，其时间复杂性为$N\log N$。

练习 4.5 矩阵的块状结构，所有的对角线块都有相同的大小，这是因为我们在方形域上定义了我们的BVP。画出方程离散化所产生的矩阵结构，同样是中心差分，但这次是定义在一个三角形域上；见图4.2。说明同样有一个块状三对角矩阵结构，但现在块的大小不一。提示：先画出一个小例子。对于$n=4$，你应该得到一个$10\times 10$的矩阵，其块结构为$4\times 4$。

对于更加不规则的域，矩阵结构也将是不规则的。

有规律的块状结构也是由我们决定将未知数按行和列排序造成的。这被称为自然排序或词法排序；其他各种排序也是可能的。一种常见的未知数排序方式是红黑排序或棋盘排序，这对并行计算有好处。这将在第6.7节讨论。

关于BVP的分析方面还有更多要说的（例如，解有多平滑，它是如何取决于边界条件的），但这些问题超出了本课程的范围。这里我们只关注矩阵的数值方面。在线性代数一章中，特别是第5.4和5.5节，我们将讨论从BVP中求解线性系统。

stencils 差分

离散化通常被表述为应用stencils差分 $$ \begin{array}{ccc} .& -1 &, \ -1 & 4 & -1 \ .& -1 &. \end{array} $$ 到函数$u$。给定一个物理域，我们将stencils应用于该域中的每一个点，得出该点的方程。图 4.1 说明了一个由$n\times n$点组成的正方形域的情况。将此图与上述方程联系起来，你会发现同一条线上的连接产生了主对角线和第一条上下对角线；与下一条线和上一条线的连接则成为非对角线块的非零点。

这种特殊的模版通常被称为 “五点星 “或五点模版。还有其他不同的网板；其中一些网板的结构在图4.3中得到描述。只有水平或垂直方向连接的网板被称为 “星形网板”，而有交叉连接的网板（如图4.3中的第二种）则被称为 “星形网板”。

练习 4.6 考虑图4.3中的第三个模版，用于正方形域上的BVP。如果我们再次将变量按行和列排序，所得到的矩阵的稀疏性结构是什么样子的？

除了五点星形之外，还可以使用其他stencils来达到更高的精度，例如，给出一个$O(h^4)$的截断误差。它们还可以用于上面讨论的微分方程以外的其他微分方程。例如，不难看出，对于方程$ux+u{yyy}=f$，我们需要一个同时包含$x$ ，$𝑦\pm h$和$x$ ，$𝑦\pm 2h$连接的钢网，如图中的第三个钢网。相反，使用5点stencils，没有系数值的情况下，四阶问题的离散化小于$O(1)$截断误差。

虽然到目前为止讨论的是二维问题，但对于诸如$-u_x- u_y - u_z= f$这样的方程，它可以被推广到更高维。例如，5点stencils的直接泛化，在三维空间中变成了7点stencils。

其他离散化技术

在上面，我们用有限差分法来寻找微分方程的数值解。还有其他各种技术，事实上，在边界值问题的情况下，它们通常比有限差分更受欢迎。最流行的方法是有限元法和有限体积法。尤其是有限元方法很有吸引力，因为它比有限差分更容易处理不规则的形状，而且它更适合于近似误差分析。然而，在这里讨论的简单问题上，它给出的线性系统与FD方法相似甚至相同，所以我们将讨论限制在有限差分上，因为我们主要关注的是线性系统的计算方面。

初始边界值问题

现在我们将继续讨论初始边界值问题（IBVP），正如你可能从名字中推断的那样，它结合了IVP和BVP的各个方面。在这里，我们将把自己限制在一个空间维度。

我们考虑的问题是棒材中的热传导问题，其中$T(x,t)$描述了在时间上$x$的温度，对于$x\in [a,b]，t>0$。所谓的热方程（见附录15对一般的PDE，特别是热方程的快速介绍）是 $$ \frac{\partial}{\partial t} T(x, t)-\alpha \frac{\partial^{2}}{\partial x^{2}} T(x, t)=q(x, t) $$ 其中

• 初始条件$T(x,0) = T_0(x)$ 描述初始温度分布。

• 边界条件$T(a,t)=T_a(t)$，$T(b,t)=T_b(t)$描述棒的末端，例如可以固定在一个已知温度的物体上。

• 杆的材料由一个参数$\alpha >0$来模拟，即热扩散率，它描述了热量在材料中扩散的速度。

• 强迫函数$q(x, t)$描述了外部应用的加热，作为时间和地点的函数。

IBVP和BVP之间有一个简单的联系：如果边界函数$T_a$和$T_b$是常数，并且$q$不依赖于时间，只依赖于位置，那么直观地说，$T$将收敛到一个稳定状态。这方面的方程式是$-\alpha u’’(x)=q$。

离散化

现在我们将空间和时间都离散化，即$x{j+1}=x_j+\Delta x{k+1}=tk+\Delta t$，边界条件$x_0=a$，$x_n=b$，并且$t_0=0$。 我们写出$T{jk}$的数值解，$x=x_j$，$t=t_k$；运气好的话，这将近似于精确解$T(x_j，t_k)$。

对于空间离散化，我们使用中心差分公式 $$ \left.\frac{\partial^{2}}{\partial x^{2}} T\left(x, t{k}\right)\right|{x=x{j}} \Rightarrow \frac{T{j-1}^{k}-2 T{j}^{k}+T{j+1}^{k}}{\Delta x^{2}} $$ 对于时间离散化，我们可以使用前面中的任何一种方案。我们将再次研究显式和隐式方案，对所产生的稳定性有类似的结论。

显式方案

通过明确的时间步进，我们将时间导数近似为 $$ \left.\frac{\partial}{\partial t} T\left(x{j}, t\right)\right|{t=t{k}} \Rightarrow \frac{T{j}^{k+1}-T{j}^{k}}{\Delta t} $$ 将此与空间的中心差异结合起来，我们现在有 $$ \frac{T{j}^{k+1}-T{j}^{k}}{\Delta t}-\alpha \frac{T{j-1}^{k}-2 T{j}^{k}+T{j+1}^{k}}{\Delta x^{2}}=q{j}^{k} $$ 我们将其改写为 $$ T{j}^{k+1}=T{j}^{k}+\frac{\alpha \Delta t}{\Delta x^{2}}\left(T{j-1}^{k}-2 T{j}^{k}+T{j+1}^{k}\right)+\Delta t q_{j}^{k} $$

在图4.4中，我们将其呈现为一个差分stencils。这表示每个点的函数值是由前一个时间层次上的点的组合决定的。

将给定的$k$和所有$j$值的方程组用矢量形式概括为 $$ T^{k+1}=\left(I-\frac{\alpha \Delta t}{\Delta x^{2}} K\right) T^{k}+\Delta t q^{k} $$ 其中 $$ K=\left(\begin{array}{cccc} 2 & -1 & & \ -1 & 2 & -1 & \ & \ddots & \ddots & \ddots \end{array}\right), \quad T^{k}=\left(\begin{array}{c} T{1}^{k} \ \vdots \ T{n}^{k} \end{array}\right) $$ 这里的重要观察是，从$T_{k+1}$得出向量$T_k$的主要计算是一个简单的矩阵-向量乘法。 $$ T^{k+1} \leftarrow A T^{k}+\Delta t q^{k} $$ 其中$A= I- \frac{\alpha \Delta t}{\Delta x^2}K$。这是第一个迹象，表明稀疏矩阵-向量乘积是一个重要的操作。使用显式方法的实际计算机程序通常不形成矩阵，而是评估方程。然而，为了分析的目的，线性代数公式是更有见地的。

在后面的章节中，我们将考虑操作的并行执行。现在，我们注意到显式方案是琐碎的并行的：每个点都可以只用周围几个点的信息进行更新。

隐式方案

在上述方程中，我们让$T{k+1}$从$T_k$定义。我们可以通过从$T{k-1}$定义$Tk$来扭转这一局面，正如我们在第4.1.2.2节中对IVP所做的那样。对于时间离散化，这就得到 $$ \left.\frac{\partial}{\partial t} T(x, t)\right|{t=t{k}} \Rightarrow \frac{T{j}^{k}-T{j}^{k-1}}{\Delta t} \quad \text { or }\left.\quad \frac{\partial}{\partial t} T(x, t)\right|{t=t{k+1}} \Rightarrow \frac{T{j}^{k+1}-T{j}^{k}}{\Delta t} $$ 整个热力方程的隐式时间步长离散化，在$t{k+1}$中进行评估，现在变成了。 $$ \frac{T{j}^{k+1}-T{j}^{k}}{\Delta t}-\alpha \frac{T{j-1}^{k+1}-2 T{j}^{k+1}+T{j+1}^{k+1}}{\Delta x^{2}}=q{j}^{k+1} $$ 或者 $$ T{j}^{k+1}-\frac{\alpha \Delta t}{\Delta x^{2}}\left(T{j-1}^{k+1}-2 T{j}^{k+1}+T{j+1}^{k+1}\right)=T{j}^{k}+\Delta t q{j}^{k+1} $$ 图4.5将其渲染成一个模版；这表达了当前时间层上的每一个点都会影响到下一层的点的组合。我们再一次用矢量的形式来写这个。

$$ \left(I+\frac{\alpha \Delta t}{\Delta x^{2}} K\right) T^{k+1}=T^{k}+\Delta t q^{k+1} $$ 与显式方法相比，在显式方法中，矩阵-向量乘法就足够了，从$𝑇^{𝑘+1}$推导出的向量$𝑇^𝑘$现在涉及一个线性系统解决方案。 $$ T^{k+1} \leftarrow A^{-1}\left(T^{k}+\Delta t q^{k+1}\right) $$ 其中$A= I+ \frac{\alpha \Delta t}{\Delta x^2}K$ 一个比矩阵-向量乘法更难的操作。在这种情况下，不可能像上面那样，直接评估方程（4.21）。使用隐式方法的代码实际上形成了系数矩阵，并以此来解决线性方程组。解决线性系统将是第5章和第6章的重点。

与显式方案相比，我们现在没有明显的并行化策略。线性系统的并行求解将在第6.6节及以后的章节中占据我们的位置。

练习4.7 证明隐式方法的一个时间步骤的flop数与显式方法的一个时间步骤的flop数是相同的。(这只适用于一个空间维度的问题)。至少给出一个论据，说明为什么我们认为隐式方法在计算上 “更难”。

我们在这里使用的数值方案是时间上的一阶和空间上的二阶：截断误差为$O(\Delta t+ \Delta x^2)$。也可以通过使用时间上的中心差分来使用一个时间上的二阶方案。另外，见练习4.8。

稳定性分析

现在我们在一个简单的情况下分析显式和隐式方案的稳定性。让$q\equiv 0$，并假设$Tj^k=\beta^ke^{i\ell x_j}$，对于某些$\ell$。这一假设在直觉上是站得住脚的：由于微分方程没有 “混合 “$x$和$t$的坐标，我们推测解决方案将是$T(x,t) = v(x) ⋅ \omega (t)$ 的单独解决方案的乘积。 $$ \left{\begin{array}{ll} v{x x}=c{1} v, & v(0)=0, v(1)=0 \ w{t}=c{2} w & w(0)=1 \end{array}\right. $$ 唯一有意义的解发生在$c_1, c_2 < 0$，在这种情况下我们发现。 $$ v{x x}=-c^{2} v \Rightarrow v(x)=e^{-i c x}=e^{-i \ell \pi x} $$ 其中我们用$c=\ell \pi$代替，以考虑到边界条件。 $$ w(t)=e^{-c t}=e^{-c k \Delta t}=\beta^{k}, \quad \beta=e^{-c k} $$ 如果对这种形式的解决方案的假设成立，我们需要$|\beta |<1$的稳定性。将推测的$T{j}^k$形式代入显式方案，可以得到 $$ \begin{aligned} T{j}^{k+1} &=T{j}^{k}+\frac{\alpha \Delta t}{\Delta x^{2}}\left(T{j-1}^{k}-2 T{j}^{k}+T{j+1}^{k}\right) \ \Rightarrow \beta^{k+1} e^{i \ell x{j}} &=\beta^{k} e^{i \ell x{j}}+\frac{\alpha \Delta t}{\Delta x^{2}}\left(\beta^{k} e^{i \ell x{j-1}}-2 \beta^{k} e^{i \ell x{j}}+\beta^{k} e^{i \ell x{j+1}}\right) \ &=\beta^{k} e^{i \ell x{j}}\left[1+\frac{\alpha \Delta t}{\Delta x^{2}}\left[e^{-i \ell \Delta x}-2+e^{i \ell \Delta x}\right]\right] \ \Rightarrow \beta &=1+2 \frac{\alpha \Delta t}{\Delta x^{2}}\left[\frac{1}{2}\left(e^{i \ell \Delta x}+e^{-\ell \Delta x}\right)-1\right] \ &=1+2 \frac{\alpha \Delta t}{\Delta x^{2}}(\cos (\ell \Delta x)-1) \end{aligned} $$ 为了保持稳定，我们需要$|\beta|<1$:

• $\beta<1 \Leftrightarrow 2 \frac{\alpha \Delta t}{\Delta x^{2}}(\cos (\ell \Delta x)-1)<0$ : this is true for any $\ell$ and any choice of $\Delta x, \Delta t$.
• $\beta>-1 \Leftrightarrow 2 \frac{\alpha \Delta t}{\Delta x^{2}}(\cos (\ell \Delta x)-1)>-2:$ this is true for all $\ell$ only if $2 \frac{\alpha \Delta t}{\Delta x^{2}}<1$, that is $\Delta t<\frac{\Delta x^{2}}{2 \alpha}$

后一个条件对允许的时间步长提出了很大的限制：时间步长必须足够小，方法才能稳定。这与IVP的显式方法的稳定性分析类似；然而，现在时间步长也与空间离散化有关。这意味着，如果我们决定需要更多的空间精度，并将空间离散化$\Delta x$减半，时间步数将乘以4。

现在让我们考虑隐式方案的稳定性。将解的形式$Tj^k=\beta^ke^{i\ell x_j}$替换到数值方案中，可以得到 $$ \begin{aligned} T{j}^{k+1}-T{j}^{k} &=\frac{\alpha \Delta t}{\Delta x^{2}}\left(T{j{1}}^{k+1}-2 T{j}^{k+1}+T{j+1}^{k+1}\right) \ \Rightarrow \beta^{k+1} e^{i \ell \Delta x}-\beta^{k} e^{i \ell x{j}} &=\frac{\alpha \Delta t}{\Delta x^{2}}\left(\beta^{k+1} e^{i \ell x{j-1}}-2 \beta^{k+1} e^{i \ell x{j}}+\beta^{k+1} e^{i \ell x_{j+1}}\right) \end{aligned} $$ 除去$e^{i\ell x_j}\beta^{k+1}$，可得 $$ \begin{array}{l} 1=\beta^{-1}+2 \alpha \frac{\Delta t}{\Delta x^{2}}(\cos \ell \Delta x-1) \ \beta=\frac{1}{1+2 \alpha \frac{\Delta t}{\Delta x^{2}}(1-\cos \ell \Delta x)} \end{array} $$ 由于$1-\cos l\Delta x\in (0,2)$，分母严格>1。因此，无论l的值如何，也无论$\Delta x$和$\Delta t$的选择如何，条件$|\beta|<1$总是被满足的：该方法总是稳定的。

练习 4.8 我们在这里考虑的方案是时间上的一阶和空间上的二阶：它们的离散化顺序是$O(\Delta t)+ O(\Delta x^2)$。推导出由显式和隐式方案平均化得到的Crank-Nicolson方法，说明它是不稳定的，并且在时间上是二阶的。