概述

在数值求解ODE和PDE时，我们通常需要进行时空离散：

  • 空间离散看起来就很麻烦：处理方式与具体的求解格式相关，但是都会涉及到网格剖分，网格加密，单元之间的数据交换，边界处理等；
  • 时间离散看起来就很简单：只需要加上一层for或者while循环，最多加上时间步长的自适应调整即可。

但是就算是看起来非常简单的时间离散，在具体编程中也可能存在一些潜在的问题，下面先回顾一下几种最简单的时间离散情况，然后给出一个由格式和浮点数模型共同造成的最后时间步问题。

固定时间步数

固定时间步数是最简单的情况，在固定时间步数为 \(N\) 时，我们可以直接生成 \(N\) 次循环即可，例如 

1
2
3
4
5
6
7
8
9
const int N = 20;
const double Tend = 1.0;

double tnow = 0;
double dt = Tend / N;
for (int i = 0; i < N; ++i) {
    update(&data, tnow, dt);  // from tnow to tnow+dt
    tnow = tnow + dt;
}

固定时间步长

固定时间步长和固定时间步数的情况并不是等价的，因为在固定时间步长为 \(\Delta t\) 时， 我们不能假定终止时刻 \(T\) 是时间步长 \(\Delta t\) 的整数倍，在最后一步必须使用较小的 \(\Delta \widetilde{t}\)， 在编程中通常有两类做法：

  • 第一类做法是获取正常时间步长对应的时间步数 \(N\)，然后计算最后一步的时间步长 \(\Delta \widetilde{t}\)，对两者分别处理；
  • 第二类做法是使用while循环来处理，在临近终止时刻时修正时间步长。

第一类做法例如 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
const double Tend = 1.0;
const double dt = 0.03;

double tnow = 0;
int N = floor(Tend / dt);
double dt_end = Tend - dt * N;
for (int i = 0; i < N; ++i) {
    update(&data, tnow, dt);  // from tnow to tnow+dt
    tnow = tnow + dt;
}

// from tnow to Tend
update(&data, tnow, dt_end);
tnow = tnow + dt_end;

第二类做法例如 

1
2
3
4
5
6
7
8
9
10
const double Tend = 1.0;

double tnow = 0.0;
double dt = 0.03;
while (tnow < Tend) {
    if (tnow + dt >= Tend) { dt = Tend - tnow; }

    update(&data, tnow, dt);
    tnow = tnow + dt;
}

自适应时间步长

在一些更复杂的情况下，我们需要根据数据的变化动态调整时间步长，此时我们无法预先得知时间步数，必须使用while循环来实现，例如 

1
2
3
4
5
6
7
8
9
10
11
12
const double Tend = 1.0;

double tnow = 0.0;
double dt = 0.03;
while (tnow < Tend) {
    dt = get_dt(&data);

    if (tnow + dt >= Tend) { dt = Tend - tnow; }

    update(&data, tnow, dt);
    tnow = tnow + dt;
}

最后时间步的问题

一个不易察觉的事实是，上面除了固定时间步数之外的时间离散做法都可能产生问题，这个问题由两个原因共同触发。

第一个原因是浮点数累加是存在误差的，在某些情况下会导致最后一个时间步 \(\Delta \widetilde{t}\) 特别的小。考虑下面的例子 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include <stdio.h>

void update(double tnow, double dt) {
    printf("tnow = %.15f dt = %.12f(%e)\n", tnow, dt, dt);
}

int main(int argc, char *argv[]) {
    const double Tend = 1.0;

    double tnow = 0.0;
    double dt = 0.0125;
    while (tnow < Tend) {
        if (tnow + dt >= Tend) { dt = Tend - tnow; }

        update(tnow, dt);
        tnow = tnow + dt;
    }

    return 0;
}

程序运行结果为 

1
2
3
4
5
6
7
8
9
tnow = 0.000000000000000 dt = 0.012500000000(1.250000e-02)
tnow = 0.012500000000000 dt = 0.012500000000(1.250000e-02)
tnow = 0.025000000000000 dt = 0.012500000000(1.250000e-02)

...

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

虽然这里终止时刻恰好是时间步长的整数倍，但是实际的浮点数累和产生了误差，这导致最后一个时间步长只有机器精度量级。

第二个原因是某些数值格式需要均匀的时间步长才能保证精度，在非均匀时间步长下会产生各种问题，下面举两个例子说明。

例如在通常的资料中，多步法的典型代表——BDF2格式的具体公式为 \[ y_{n+2} - \frac43 y_{n+1} + \frac13 y_n = \frac23 \Delta t f(t_{n+2},y_{n+2}) \] 上述公式在保持等距结点时具有二阶精度，但是将其直接应用到不等距结点上只有一阶精度。 在不等距的情况下（\(h_n = t_{n+1} - t_n\)，\(h_{n+1} = t_{n+2} - t_{n+1}\)），正确的BDF2格式应该为 \[ y_{n+2} - \frac{(h_n+h_{n+1})^2}{h_n(h_n+2h_{n+1})} y_{n+1} + \frac{h_{n+1}^2}{h_n(h_n+2h_{n+1})} y_n = \frac{h_{n+1}(h_n+h_{n+1})}{h_n+2h_{n+1}} f(t_{n+2},y_{n+2}) \]

再例如线性常系数对流方程 \(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格式的精度。

上面的两个因素结合到一起：浮点数的累和误差导致可能存在一个仅有机器精度量级的最后时间步， 而某些数值格式本身不允许使用非均匀或剧烈变化的时间步，否则会影响精度，而只有机器精度量级的时间步导致的影响就更大了。

虽然上面提供了两个例子，但是实际上只有少部分数值格式会对小时间步长敏感，大部分数值格式对于极小的时间步长都是无所谓的，它们都满足 \(u_{n+1} = u_n + O(\Delta t)\) 的形式，此时只有机器精度的小时间步长显然不会影响误差。

我们可以通过处理累和误差来避免这种问题，即在终止时刻的判断和最后时间步的计算中，考虑浮点数累和误差的存在，这样就可以避免出现最后时间步为机器精度量级的情况，例如 

1
2
3
4
5
6
7
8
9
10
11
const double Tend = 1.0;
const double ep = 1e-13;

double tnow = 0.0;
double dt = 0.0125;
while (tnow < Tend - ep) {
    if (tnow + dt >= Tend - ep) { dt = Tend - tnow; }

    update(&data, tnow, dt);
    tnow = tnow + dt;
}