数值计算中的分析方法

我发现很多同学对数值计算中的一些基本概念的分析方法不是特别扎实，所以特此整理一份学习笔记， 希望能够起到一个抛砖引玉的效果。笔记多来源于自己以前上课的course notes，如有错漏请轻拍。

1.1 函数与泰勒展开

人们用以建模分析物理系统的理论基石是泰勒展开, 在这里, 假设读者具备基础数学知识, 一个动态的物理系统通常由时空域向量函数抽象表示

它在时空中任意一点的邻域内可近似表达为

其中, 上式右边第二项中的 通常叫做系统的Jacobian, 而 是系统的二阶导, 通常叫做Hessian

Jacobian的形状是一个矩阵, 它每一个元素的意思是f ⃗中的某一个维度关于x ⃗的某一个维度的微分, 记为

Hessian更为复杂一些, 它相当于是整个系统沿着单位正交方向的”加速度”. 为了帮助大家理解， 我们可以联系高中时学过的位移公式作类比（针对线性加速度）：

事实上，另外一种看待 的方式是使用方向导数：

物理意义上，这就相当于是 这个系统，沿着 方向的“加速度”。

泰勒展开为我们能够研究动态的物理系统提供了坚实的理论基础，从泰勒级数的推导中我们可以得出一个朴素的方法：即系统的下一个时刻的状态，可以通过当前时刻已知的状态+一定的计算得到。这就形成了“数值积分” 。

第二章 数值积分

大家都知道世界上第一台计算机名字叫埃尼阿克(ENIAC), 可能很少有人知道， 它被发明的目的和主要工作，就是进行数值积分计算。其英文全称是Electronic Numerical Integrator And Computer(电子数值积分计算器)。

数值积分算法与求解器的设计，是计算物理永恒的主题，这其中涉及几个比较重要的概念，我们分几个子章节进行探讨。

2.1 Consistency

Consistency要求所给出的数值积分格式，与函数的泰勒展开表达式相差的是高阶无穷小，不然的话，我们的数值积分不可能逼近原函数。

不失一般性地讲，对于所给出的数值积分方法，

我们需要左边减右边得到的截断误差，是高阶无穷小。

例子 2.1.1

对于微分方程 讨论前向欧拉方法的Consistency

根据泰勒展开，有

左减右， 得

则前向欧拉是一个Consistent的近似。事实上，T时刻的精确值是

我们得到的近似是

根据积分中值定理，我们有

是一个对原方程的解的一阶近似。

至此，我们初步了解了Consistency的概念，那是否只要满足这个条件，我们就能够给出任何微分方程的足够好的近似呢？

答案是否定的， 这关系到我即将引出的第二个概念：Stability

2.2 有关Stability的一些概念

事实上，由于浮点数表达和对于原函数在离散的时空位置进行的近似，我们无法连续地近似真正的积分过程， 而Stability所分析的问题就是在时空离散下的解的行为。

概念Zero-Stability:

一个方法是Zero-Stable的， 当任意时刻产生的计算误差是会被该方法缩小的。

重要定理：

数值积分能够逼近真实解，当且仅当所使用的的方法是Consistent且Zero-Stable的

案例 2.2.1

考虑微分方程 以及前向积分方法

当 时，我们在计算过程中连续两步之间的计算误差是会被缩小的， 则满足此条件的前向积分是Zero-Stable的。

当 , 即 时， 我们的解会以越来越大的幅度做正负向震荡，最终数值爆炸。像这样的需要满足一定条件才能执行的数值积分方法又叫做conditionally stable。

如果采用了这种数值方法进行物理仿真，就有可能会产生我们熟知的名词“蝴蝶效应”。

比如当数值积分器求解的地球地理仿真发散时，便会产生“非洲的蝴蝶拍一拍翅膀， 美国的西海岸发生海啸”这样无稽的错误结果。

案例 2.2.2

考虑与上一个案例同样的微分方程， 以及后向积分方法

移项并约分，可得

则对于任意的 , 我们的数值解都不会发散，且稳定地趋向一致的结果。

这样的方法又叫A-Stable, 或者unconditionally Stable，一个A-Stable的方法, 必定是Zero-Stable的.

注意稳定不意味着更准确， 在这两个例子中， 无论前向还是后向积分， 采用的都是一阶格式，所以收敛的时候对真实值的近似是 .

心得：隐式积分通过在导数侧用下一个时刻的状态来表达导数来得到更稳定的积分格式， 一种形象地理解是：

显式积分：

系统未来状态=系统当前状态+F(系统当前导数,Δt)

隐式积分：

系统未来状态-F(系统未来导数,Δt)=F(系统当前状态,Δt)

通常地，这将把原本简单的函数计算转换为解方程计算， 对于复杂的系统，甚至有可能出现关于变量的非线性函数， 届时，或将对系统进行线性化，或将使用牛顿迭代法来求解非线性系统。这看起来很有道理，毕竟我们为此付出了额外的智慧和计算量。它带来的益处就是我们可以使用百倍于显式积分法的Δt，当我们一次求解系统的代价优于进行更多的显式积分时，我们实际上的积分效率提升了。

隐式积分通常比显式积分更为稳定，但通常都需要付出一些代价，比如大型线性方程求解或者大型非线性系统的牛顿法求解，或者带约束的非线性优化求解， 在很少的情况下， 我们能够得到A-Stable的显式积分方法， 这是极好的， 且这样的算法需要对物理系统更多的洞见和其发现者的灵机一动。

案例 2.2.3

考虑一个ODE系统，其中包含多个物理变量 , 满足方程

如果我们对矩阵A进行特征值分解，即找到 , 同时令 , 则我们有

特征值分解将原来的N维常微分系统转化为了N个各自独立的常微分方程， 最有帮助的是， 我们可以很方便地分析这个系统，我们有：

1. 这个系统可以找到稳定的方法， 1.当且仅当所有的λ_i都是负的。
2. 1.这个系统可以走的时间步长， 是由它最大的 决定的。若采用之前讨论的显式积分，则要求 .
   1. 然而， 越大的特征值意味着对应的物理结构随时间的消逝越快， 越小的特征值意味着对应的物理量具有越显著的运动特性.
   2. 这意味着我们将用很小的时间步长花很长时间来模拟那些转瞬即逝的现象，而只能很慢地积分那些现象显著的部分。
3. 这样的问题， 是一个Stiff Problem，通常来讲对于Stiff Problem的定义没有固定的说法， 但通常这类问题的一个显著的特征是，系统的最大特征值和最小特征值差别巨大
4. 对于更真实的物理问题来说， 矩阵A就是所研究的物理系统的Jacobian。
5. 在求解Stiff Problem时，隐式的方法， 通常比显式的方法好。

2.3 数值积分的分析方法

在之前的章节中， 我们接触了一些计算数值积分过程中的核心概念， 在实际应用中， 针对不同的物理方程和需求， 往往需要不同的算法选择，这就需要相关的工程人员有分析问题和设计算法的能力。这一章中，我会介绍一些比较强大的分析工具和思维导图。

2.3.1 思维范式

通常的，直接对3维的，复杂的， 完整的物理问题的求解方法进行分析， 是不现实的， 也是很难做到的。 我们首先需要做的是化完整为局部， 化高维为低维，化物理问题为数学模型。幸好，由于我们在解决的是物理数学问题， 所以往往在低维问题上分析的结果，是容易扩展到高维的。

对于低维问题的分析， 即使不能告诉我们什么是最好的， 但它能让我们知道什么是一定不行的。

2.3.2 Modified Equation

在我们得到足够简化的问题后， 我们就有可能对问题施加十分准确的分析， 这其中的一种工具就是Modified Equation方法。

在之前的章节中， 我们通过时空离散来用数值积分近似原有的连续函数， 得到近似解。 Modified Equation则反其道而行：它试图找到连续的函数， 使得该函数是对于数值积分方法更高阶的近似。我们将通过案例来认识ME方法。

案例 2.3.1

分析二维常微分方程

使用前向积分和后向积分时解的行为。

易得，这个方程的解是一个绕原点运动的半径为1的圆。因为在每一个路径点处的运动方向垂直于该点向原点的连线。 因为：

该系统的前向积分由下式给出：

为分析该数值方法的行为， 我们不妨假设一个时空连续的函数系统，它在任意时刻（而非离散时刻）都以更高阶的精度满足给定的数值方法，即

泰勒展开，得:

消元，移项，整理，舍去高阶小量后得：

由于 , 我们有

这个新的ODE系统直接揭示了数值方法的行为，即我们实际得到的路径曲线，它任意时刻的速度方向不再是与半径方向垂直的， 而是同时再往法线方向外扩一点螺线，事实上以上ME的解析解就是一条半径逐渐增大的螺线

略去推导过程， 对于后向积分法，我们在极限意义下得到的则是一条半径不断向内缩的螺线：

而这， 与我们的实际算法执行结果也相符

而对于这个简单系统的分析， 告诉了我们一个重要的道理，像这种连走圆都走不好的积分方法，我们不可以拿来对行星轨道，卫星运动， 甚至巨量星体运动等进行分析，得到的结果将与事实千差万别（比如这样的方法可能会认为地球在以快得多地多地速度靠近太阳，造成末日恐慌）

案例 2.3.2

对前一案例所讨论的ODE系统，拟用Adams-Bashforth方法进行积分，试对该方法进行分析。

Adams-Bashforth方法由下式给出

经过泰勒展开， 带入， 化简， 舍去高阶无穷小项后，我们可以得到该方法的ME为（感兴趣的读者可以自己作业）：

这预示着， 采用Adams-Bashforth方法求解得到的曲线，将仍然是一个半径为1的圆（速度方向总是垂直于半径方向），只是越往后，该轨迹运行的角速度会越快！但这个加快幅度是很小的O(Δt^2) ，由ME的解析解可以看到

2.3.3 冯诺依曼分析（Von-Neumann Analysis）

回想一下我们在案例2.2.3中使用特征值来研究系统的刚性的方法， 当时我们假设了任意的ODE系统矩阵，而现实中， A往往来自于我们对于物理问题的空间离散， 并最终和时间离散一起构成了“偏微分方程组”，事实上在经过足够的假设和简化后，对于物理微分算子的离散， 常常是可以被傅里叶变换对角化的，此时，我们有

即通过令 y ̂=FT(x), x=FT^T (y ̂)， 我们将全局耦合的偏微分方程转换为了傅里叶空间的解耦的常微分方程

然后我们就可以通过研究离散微分算子的频率空间性质来研究数值方法的作用，这样的分析方法，就是冯诺依曼分析。

案例 2.3.3

考虑一个1D的对流方程：

已知中心差分格式在任意情况下都会导致系统爆炸， 为什么？迎风格式是如何得到的？

事实上，冯诺依曼分析的出现，正正在当时给与了以上这个系统强大的洞见，回答了中心差分在该系统失效的原因，且引导人们找到了解决这个系统的良药。

为了简化问题，我们假设ϕ这个函数是周期性的， 并对它进行傅里叶变换，考虑到其中心差分格式为：

任何一个 都能被表达为傅里叶级数，即，

则原方法在傅里叶空间中任意一个 上的特征ODE就是

的模， 无关乎Δt, Δx怎么取， 都恒大于+1，则中心差分方法在这个方程上的特征方程不具有Zero-Stability， 于是该方法必然数值爆炸。就算尝试再如何使用更小的Δt来求解这个方程， 苦苦等待， 也不会有结果！

于是科学家根据这个分析的启发，构造了迎风格式，迎风格式通过判断速度u的方向，而总是从上风方向来构造 的差分近似， 比如，

这个差分格式的特征值 ，经过计算后， 为：

随着算法和数学理论的发展， 针对以上的方程，人们已经找到了可以超越CFL Condition的显式积分方法, 也是熟知的Semi-Lagrangian方法。