2.2　常微分方程

海洋数值模型总是包含着对某种形式的偏微分方程的求解，许多模拟问题需要求解常微分方程（ODEs）。1.15节中对温盐环流的区块模拟的讨论和无平流和扩散效应的生物生态系统模拟是典型的例子。这些问题由一组常微分方程表示，这些方程通常是非线性的和耦合的。无论其中所包含方程的阶数是多少，该组方程总是能够简化到N个因变量的N个一阶常微分方程。

在偏微分方程的求解过程中，利用有限差分法、有限元素法或者半光谱法进行空间离散化后，得到的方程都是常微分方程。下文中将要讨论的所有形式都与这个问题有关，要对期望精度、内存需求和CPU需求几方面进行一致考虑。由于在大多数海洋数值模型问题中，用一个已知的初始值对物理系统进行的改进受到很大关注，这些问题称为初始值问题。初始值问题包括对系统的未来状态进行求解，也就是N个变量在有限时间t内的值，给定初始状态下，这些变量在t=0时的值。求解过程包括将时间分成许多小的步长，然后使每一步的解得到改进。所选形式的精度、效率以及稳定性自然与所选择的步长有关系。

对耦合的非线性一阶常微分方程组的初始值问题进行求解的久负盛名的方法是四阶龙格库塔方法（RK）。然而，还有很多高效的方法，例如高效自适应步长嵌入式五阶RK方法和BS方法，这两种方法在所采用的步长下对截断误差大小进行计算，然后在当前步长下用该信息对步长进行降低，从而达到预定的精度，并在下一步中使步长最优化，这些在Press（1992）的观点中有介绍和描述。我们将对这些先进方法的原理进行简单描述，这里只对传统的RK方法和Adama-Bashforth预测校正方法进行描述，因为它们对大多数非重复性应用来说已经足够。

2.2.1　龙格库塔方法

首先，对单一因变量y（t）的单一一阶常微分方程进行考虑，y（t）的形式为

在f（y，t）函数形式下，方程可以是线性的，也可以是非线性的。y（0）值是已知的，要对y（T）进行求取。时间区间t=0到t=T被分成N个子区间（n=1…N+1），那么问题就成为：yn已知的情况下，怎样高效而精确地确定yn+1，时间步长为Δt=tn+1-tn。经典的显式（或前向）欧拉方法只在初始点tn对式（2.2.1）的导数进行计算，这一一阶精度方法的精度和稳定性都不是很高。四阶RK方法除了要计算初始点处的导数外，还要计算两次导数，一次是中点处的导数，还有一次是终点tn+1处的导数，对4个导数值进行加权平均，用它来求tn到tn+1时的导数值为

需要注意的是，每一次新的导数估计都是在之前的最佳估计基础上进行的。如果函数f只包含独立变量t的话，那么该方法就等价于对函数f（t）积分的经典辛普森公式。每一步的截断误差是o（Δt5），但整个区间上的累计截断误差是o（Δt4）。尽管这种方法不是最高效的，但是准确度高，而且鲁棒性很高。扩展为N方程系

这是很简单的（在温盐灾难问题中的应用可以参考本章2.1节，它包含两个耦合的一阶非线性常微分方程构成的式子）。

如果取代使用任意预定的固定步长的话，可以在求解过程中对步长进行自适应控制，RK解的效率可以得到极大改善。自适应步长控制要求每一时间步长下的求解中对截断误差进行估计，这可以在每一步获得两个解得到，其中一个解是在满步长下得到，另一解是在半步长下得到。因为截断误差与步长的关系是已知的，两个解（一个的分辨率为h，而另一个的分辨率为h/2）之间的差异能够方便地对截断误差进行估计，从而使得我们可以确定用什么样的步长来获得达到预定精度的解。如果在当前的一步中没有得到期望精度，便可用所确定的最优步长进行再求解，如果超过了，则可以在下一步中将步长增加到最优值。这一局部外推法的简单原理使我们能够根据需要使用较小的步长或者在步长较小时达不到期望精度的情况下使用较大步长，从而极大地改善计算效率。

另一种估计截断误差的方法是使用嵌入的RK公式，同样可以用两个函数来计算精度不同的两个解。Fehlberg的方法（Press等，1992）利用一个五阶RK公式和一个嵌入的四阶RK公式的解之差来估计截断误差，然后用它自适应地改变步长以获得期望精度。Press等（1992）对这个高效的RK方法及与其相关的代码进行了描述，这种方法用于RK解的效率是一个关键问题的情况。

当解的性质未知的情况下，RK方法具有鲁棒性，它是这种情形下的最佳方法。对于在积分区间上没有任何奇异点的良态函数来说，存在一种更为有效的积分方法，这些方法被称为Bulirsch-Stoer方法（Press等，1992），它在每一时间步长中都包含对零步长的极限进行外推，这是通过在每一步中用不同的步长求取多个解而实现的，把这些解看成步长的解析函数的话，就可以用一个多项式函数外推或一个有理函数外推来确定步长无限小的情况下的解。外推最终得到函数的解及与之相关的误差估计，然后用它来在每一步中增加解的个数而对解进行改进。很显然，每一步得到的解的个数、所采用的步长以及外推方法都决定着这种方法的效率。所采用的步长可以用预先确定的序列得到，但是每一步得到的解的个数则由一个迭代过程确定，起始时为两个，在迭代过程中对外推得到的解的期望误差特征进行检查。如果对得到的解不满意，那么用下一个（减小的）步长对另一个解进行求取，求解的方法使用理查德森外推法，并对误差进行重新检查，一直重复这个过程，直到满足期望误差，或者终止，整个过程用整体上减小的步长不断重复。如上述的嵌入RK方法那样，自适应步长控制方法被用来控制整个步长。Press等（1992）也对此方法进行了非常详细的描述，并介绍了与之相关的代码，在对整个域上解的奇异点的存在性产生怀疑并且对高效性有贡献。

一种被成功应用到地球物理问题的时间离散化中的方法是三阶Adams-Bashforth公式（Durran，1991；Haidvogel等，1997）

该方法克服了海洋和大气模型中对偏微分方程进行时间差分的二阶RK方法和其他常用方法的一些局限性，它还可以用来作为预测校正公式的一部分，预测步长下基于式（2.2.4）得到yn+1的第一个估计值，^表示第一个估计值，然后用它来在校正步骤中对yn+1的值进行更新

这避免了对隐式方程进行求解。在每一步中用该方程进行迭代求解以获得期望精度。对于包含时间导数的复杂表达式的问题来说，例如与时间相关的多维地球物理问题，这些预测校正方法通常比Bulirsch-Stoer方法更好（Press等，1992）。然而，这种方法代价较高，因此也很少使用。

上面讨论的方法是“显式”形式的。当因变量为两个或更多时，就会产生问题中所包含的时间尺度不止一个的可能性，在这种情况下，尽管所考虑的精度可能会允许一个较大的时间尺度或者这种形式可能会变得不稳定，但为了在时间尺度最小的情况下进行求解，时间步长必须要足够小。这些方程被称为刚性方程，当要保持稳定性时，它需要允许使用较大步长（只与期望精度相一致）的“隐式”形式。隐式的RK和BS形式可参考Press等（1992）。在此用一个简单的线性一阶方程对此进行说明

当Δt＞2/α时是不稳定的，因为|yn|随时间步长n单调增加，这种方法就是一个隐式的形式，它不计算步长n处的导数，而是计算步长n+1处的导数（后向欧拉法）

yn随n单调递减，因为它要与原始方程一致，无论步长为多少。然而，对普通的方程组执行隐式形式是很困难的，但利用对步长n+1进行求导的泰勒级数展开式的第一项能使这种执行变得简单。例如，可以用这种半隐式方法对上述方程进行求解（虽然该情形下是不需要的）

在当前讨论的情形中该方程所得结果与全隐式方法所得结果一致。

有时会存在初始点处因变量未知的情况，但是这种未知量很少，在终点（或中间点）的其余变量已知，这种问题被称为两点边界值问题。例如，我们可能知道初始点的浮游植物浓度，而不知道初始点浮游动物的浓度，只可能知道它在终点的浓度，利用初始值问题解法估计未知量在初始点的值并获得终点处的解，便很容易对这类问题进行求解。对估计值进行迭代调整，直到这些未知量的解与其在终点的已知值相匹配时停止迭代，从而完成求解过程。然而，这不属于本书的范畴，对于两点边界值问题的有效解及软件的描述，读者可以参考Press等（1992）。