在数值求解PDE时，通常需要分别进行时空离散：

  • 空间离散比较复杂：与具体格式相关，涉及到网格剖分与加密，单元之间的数据交换，边界处理等；
  • 时间离散比较简单：只需要加上一层for或者while循环，最多加上时间步长的自适应调整即可。

但是看起来非常简单的时间离散，在具体编程中也可能存在一些潜在的问题，下面用几个例子说明。

例一

考虑下面的时间离散 

1
2
3
4
5
6
t = 0.0
while t < T:
    dt = get_dt(data)

    update(data, t, dt)
    t += dt

这段代码的问题比较明显：最后一个时间迭代可能越界，出现 \(t + \Delta t > T\) 的情况。

加上越界判断即可 

1
2
3
4
5
6
7
8
9
t = 0.0
while t < T:
    dt = get_dt(data)

    if t + dt >= T:
        dt = T - t

    update(data, t, dt)
    t += dt

例二

一个不易察觉的事实是：浮点数的累加是存在误差的，这在一些情况下会导致最后一个时间步特别小，仅有机器精度量级。

考虑下面的例子 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def update(tnow, dt):
    print(f"tnow = {tnow:.15f} dt = {dt:.12f}({dt:.6e})")


T = 1.0
tnow = 0.0
dt = 0.0125

while tnow < T:
    if tnow + dt >= T:
        dt = T - tnow

    update(tnow, dt)
    tnow = tnow + dt

运行结果为 

1
2
3
4
5
6
tnow = 0.000000000000000 dt = 0.012500000000(1.250000e-02)

...

tnow = 0.987499999999998 dt = 0.012500000000(1.250000e-02)
tnow = 0.999999999999998 dt = 0.000000000000(1.554312e-15)

虽然这里的终止时刻恰好是时间步长的整数倍，但是实际的浮点数累和产生了误差，这导致最后一个时间步长只有机器精度量级。这个极小的 \(\Delta t\) 可能有什么影响？

  • 对于多步法，这表明最后几步可能不是时间等距的，而且相邻时间步的差异非常大，多步法的格式通常假定时间步等距（不等距的格式在形式上会更加复杂），直接应用必然会掉阶。
  • 对于单步格式，通常并不会出现问题，因为大部分单步格式都满足如下形式 \[ U_{n+1} = U_n + O(\Delta t) \] 加上一个机器精度量级的时间迭代几乎没什么影响。

但是在单步格式中是存在例外情况的：考虑线性常系数对流方程 \(u_t + a u_x = 0\) 的Lax-Friedrich格式 \[ v_{j}^{n+1} = \frac{v_{j+1}^n+v_{j-1}^n}{2} - \frac{a}2 \frac{\Delta t}{\Delta x} (v_{j+1}^n-v_{j-1}^n) \]

这是一个有条件相容的有限差分格式，局部截断误差为 \(O(\Delta t + \frac{\Delta x^2}{\Delta t})\)， 只有在比例 \(\frac{\Delta t}{\Delta x}\) 固定时才会具有一阶精度。在空间网格不变的情况下，如果存在极小的时间步 \(\Delta t\)，就可能影响 Lax-Friedrich 格式的精度。 （同理其它有条件相容的格式也有一样的问题）

可以加上额外的检查来避免出现只有机器精度量级的最后时间步 

1
2
3
4
5
6
7
8
9
10
eps = 1e-12
t = 0.0
while t < T - eps:
    dt = get_dt(data)

    if t + dt >= T - eps:
        dt = T - t

    update(data, t, dt)
    t += dt