关于差分中的库朗数和计算稳定性的一个注记

龚明，中国科学技术大学

我们经常用差分方法计算一些时间演化问题，包括扩散问题、薛定谔方程、电磁学、流体力学方程等。这些方程的共同特点是涉及时间差分$\delta t$和空间差分（假设一个一维系统，其空间差分为 $\delta x$）。这两个间隔不能随便取，并不是说这些差分越小越好。为了保证稳定性，这两个差分$\delta x$和$\delta t$必须满足一定的条件。这个条件叫做Courant数（即库朗数）, 有时也叫Courant–Friedrichs–Lewy（CFL）数。它是一个无量纲的数，在很多方程中普遍存在，而且有明确的解释。本文给出一个简单的推导证明，并解释之。

一个简单的可解模型

首先，我们分析下面的方程\begin{equation*} \partial_t u + a \partial_x u =0. \end{equation*}其中$a$表示速度。这个方程当然很容易处理，甚至可以解析求解，所以便于数值模拟检验。这个方程可能有广泛的意义，因为它和波动方程有关，也和麦克斯韦方程有关。我们定义$u(x, t) \rightarrow u_k^n$, 其中$k$, $n$分别表示空间和时间。最常见的差分法是中心差分\begin{equation*}{u_k^{n+1} - u_k^n \over \delta t} + a {u_{k+1}^{n} - u_{k-1}^n \over 2\delta x} =0. \end{equation*}这个方程可以定义一个无量纲数（库朗数）\begin{equation*}C = {a \delta t \over \delta x}.\end{equation*}同时，我们假设其解$u_{k}^n = A_n e^{iq k \delta x}$, 这里$q$是某个参数，比如它代表平面波的动量。代入这个方程中，我们得到\begin{equation*}A_{n+1}/A_n = 1 - iC \sin(q\delta x). \end{equation*} 显然, $|A_{n+1}/A_n|^2 =1 + C^2 > 1$。它表明，随着时间推移，$A_n$的幅度会越来越大，即$A_n$会发散：这也预示着不稳定。所以这个差分法对任意间隔$\delta t$和$\delta x$都是不稳定的。这个叫做无条件不稳定（unconditional unstable）。

换一个差分法可能可以保证其稳定性，比如采用下面的形式\begin{equation*}{u_k^{n+1} - u_k^n \over \delta t} + a {u_{k}^{n} - u_{k-1}^n \over \delta x} =0. \end{equation*}这是一个很小的变化。采用上面的解，我们得到\begin{equation*}A_{n+1}/A_n = 1-C + C \cos(q\delta x) - i C \sin(q\delta x). \end{equation*} 我们要求对任意$k$都有$|A_{n+1}/A_n|\le 1$，则要求$2C(1-C) \ge 0$，所以\begin{equation*} 0 \le C \le 1. \end{equation*}此时，我们看到，这个差分法可以导致稳定的结果。但是下面的向前差分却总是不稳定的\begin{equation*}{u_k^{n+1} - u_k^n \over \delta t} + a {u_{k+1}^{n} - u_{k}^n \over \delta x} =0. \end{equation*}我们省略这个证明细节。这个例子传递几个重要的信息，比如，它可以很容易解释库朗数的存在性，同时，它也说明，我们通过选择合适的差分法是有可能让它稳定的。在实际计算中，我们需要在$\delta x$和$\delta t$上取一个平衡，它不能无限靠近这个上限。这是因为它是一个粗糙的估计，真实的情况可能不会如此理想。通常，我们会选择一个合适的$C$，比如在这个上限的一半或者三分之一，从而保证在计算过程中，在计算速度和精度上选一个合适的平衡。不同的数值方法，需要在提高数值精度的同时，保证其数值稳定性。所以，在一些模拟软件中，它是一个很重要的参数。

扩散方程

那么对于扩散方程呢？扩散方程可以写成\begin{equation*}u_t = D \partial_x^2 u. \end{equation*}其中$D$为扩散系数。我们采用下面的差分形式\begin{equation*} {u_k^{n+1} - u_k^n \over \delta t} = D {u_{k+1}^n + u_{k-1}^n -2 u_k^n \over \delta x^2}. \end{equation*} 同样可以定义无量纲的库朗数以及定义它的解\begin{equation*} C = {D \delta t \over \delta x^2}, \quad u_{k}^n = A_n e^{i q k \delta x}. \end{equation*}代入这个方程得到\begin{equation*}A_{n+1} = A_n + C A_n (2 \cos(q \delta x) -2). \end{equation*}所以利用稳定条件$|A_{n+1}/A_n| \le 1$可以给出\begin{equation*} 0 \le C \le {1\over 2}. \end{equation*}可见，对扩散方程，我们也得到\begin{equation*} C \sim \mathcal{O}(1). \end{equation*}这个关系可能是普适的。很多不同的模型，它们的库朗数都是在这个范围，而且一般都有$C < 1$。我们定义这个最大值为$C_\text{max}$，其具体数值，还需要看它的具体定义，不能简单看它的数值。但上述条件，$C_\text{max} = \mathcal{O}(1)$，看起来却几乎是普适的。

库朗数的物理意义

（图片来自网络；这图片中的CFL为Courant–Friedrichs–Lewy缩写； 图片来自Courant number in CFD simulations，https://www.idealsimulations.com/resources/courant-number-cfd/）

这个库朗数应该有明确的物理意义。比如上图的解释就非常好，也很直观。上图的这句话很关键

If the Courant number is greater than one, it means that the information propagates through more than one grid cell at each time step, making the solution inaccurate and potentially leading to nonphysical results or divergence of the solution in certain integration schemes.

或者，在其它地方，有下面的描述，表达相同的意思

The Courant number is a dimensionless value representing the time a particle stays in one cell of the mesh. It must be below 1 and should ideally be below 0.7. If the Courant number exceeds 1, the time step is too large to see the particle in one cell, it “skips” the cell.

也就是说，如果$C > 1$, 它表明粒子传递的速度太快（比如超过系统的最大速度），一个单位时间跨越了多个格子，从而导致非物理的结果，即数值解会发散。所以，库朗数受制约于它的最快传播速度。而任何一个系统，从声学系统到光学系统，都存在这样的速度上限。理解到了上图的图像，便可以理解库朗数了。

总结

总结一下，很多运动方程的差分法，其稳定性由一个无量纲的库朗数决定。这个数普遍存在。所以时间和空间差分间隔不能随便取，它们不是完全独立的。不同的差分法有可能给出不同的稳定性，有些差分法甚至是完全不稳定的。所以它是一个很重要的概念，在具体数值计算中，我们需要调节好这个参数，从而在数值精度和计算速度上取一个平衡。

本文没有讨论的是，在有限元方法中（比如FETD，有限元时域方法），也存在这个数。有不少文章专门讨论这个临界库朗数。

2022年10月15日《计算物理》简单介绍到这个概念。但我主要介绍的是电磁学中的库朗数（参考文献【1】）, 此时$a = c$（光速）, 一般取$\delta x = \lambda/4$，$\lambda$为波长，那么$\delta t$就存在一个上限，它应该是周期的1/4（或者类似的值，但都要求$\delta t < T$，$T$为周期）。这个结论本质上和第一个模型是一致的，它们都是一阶偏微分方程；只是电磁学方程稍微复杂一点而已，因为涉及电场和磁场分量。扩散方程的讨论是今天想到而临时写出来的。后来我找到一个网页提到这个库朗数的物理解释（上图），而这个解释和我的理解是一样的；但它的图像更直观更清晰。上课的时候，我对这个问题介绍得不是很详细，算是一个补充，是为记。

参考资料：

1. FDTD方法简介：https://my.ece.utah.edu/~ece6340/LECTURES/lecture%2014/FDTD.pdf 

注：这便文章对FDTD介绍得很清楚，也提到了库朗数，但是没有分析其稳定性。它完全可以用本文提到的方法分析。电磁学中的FDTD方法，即Kane Yee的Leapfrog技巧，也不能完全消除这种不稳定性。

2. Numerical methods for the advanced equations (200年 ICTP School on Astrophysical Fluid Dynamics PPT), https://indico.ictp.it/event/a06220/session/18/contribution/10/material/0/2.pdf

3. An Adaptive, Courant-numberdependent implicit scheme for vertical advection in oceanic models, Alexander F. Shchepetkin, UCLA, Honolulu, HI, February 25, 2014, 

http://people.atmos.ucla.edu/alex/ROMS/OS_Honolulu2014Talk.pdf

4. Time-stepping, numerical stability, dispersion and conservation Dr Hilary Weller, 2016. 

http://www.met.reading.ac.uk/~sws02hs/teaching/dcmip2016/Weller_2_lec.pdf

注：欢迎大家转发这篇文章，无需征求我的同意。