迎风格式

迎风格式（upwind scheme）根据信息传播方向选择空间差分的一类离散格式。对速度为正的输运方程，当前点受到左侧上游的影响，因而用后向差分；速度为负时改用前向差分。它的稳定性依赖传播方向和 Courant 数，不能把“向某一侧作差”当成不分方程的固定规则。

从平移方程确定上游

常速度输运方程为 $u_t+a{u}_x=0$，精确解是 $u(x,t)=u_0(x-at)$。当 $a>0$ 时，$(x,t+\Delta t)$ 的值来自左侧 $x-a\Delta t$。令 $x_i=i\Delta x$，用 $U_i^n$ 近似节点值，显式 Euler 时间差分配合上游空间差分给出 $$U_i^{n+1}=U_i^n-\frac{a\Delta t}{\Delta x}(U_i^n-U_{i-1}^n)=(1-C)U_i^n+C{U}_{i-1}^n,C=\frac{a\Delta t}{\Delta x}.$$ 若 $0\le C\le 1$，足点位于相邻两个节点之间，这恰好等于在足点处作线性插值。$C=1$ 时直接平移一格，因此在恒定速度、周期边界、网格时刻上可精确输运节点数据。

取四个周期节点、$C=1/2$、初值 $(0,1,0,-1)$，一步得到 $(-1/2,1/2,1/2,-1/2)$。例如新零号节点平均的是旧零号和绕回来的旧三号，故为 $-1/2$。若 $a<0$，传播来自右侧，正确更新改为 $U_i^{n+1}=(1-|C|)U_i^n+|C|U_{i+1}^n$。

稳定性与 Courant 条件

对两个数值解之差 $e_i$，凸组合给出 $|e_i^{n+1}|\le (1-C)|e_i^n|+C|e_{i-1}^n|$。因此最大误差不增；周期求和又得 $\sum _i|e_i^{n+1}|\le \sum _i|e_i^n|$。这证明了 $0\le C\le 1$ 时的最大范数与离散一范数稳定性，并解释非负性和无新极值性质。

Fourier 分析还能检验条件的必要性。设模式 $U_i^n=G^n{e}^{i\theta i}$，代入得到 $$G=1-C+C{e}^{-i\theta },|G{|}^2=1-4C(1-C){\sin }^2(\theta /2).$$ 当 $0\le C\le 1$ 时所有频率均不增长；若 $C>1$ 或 $C<0$，任何非零且不为整周的频率都会放大。$C>1$ 还意味着真实信息在一步内跨过一个以上网格，而差分模板只看相邻上游节点；时间步长与数值信息域不匹配。

有解变式用 $C=1.2$ 输运初值 $(1,0,0,0)$，首步成为 $(-0.2,1.2,0,0)$。总和仍为一，然而负值和新峰已经出现。这区分了守恒与稳定；前者不能补救后者。

数值扩散的大小

对光滑真解作 Taylor 展开，时间差分为 $u_t+\Delta t{u}_{tt}/2+\cdots$，后向空间差分为 $u_x-\Delta x{u}_{xx}/2+\cdots$。用输运方程推出 $u_{tt}=a^2{u}_{xx}$，整理得到修正方程 $$u_t+a{u}_x=\frac{a\Delta x}{2}(1-C)u_{xx}+O(\Delta x^2+\Delta t^2),a>0.$$ 稳定范围内右边的主导扩散系数非负。网格越粗或 $C$ 离一越远，波峰更容易被抹平。这个展开适用于光滑区域；不能将不连续处不存在的高阶导数代入当作逐点证明。

周期长度一、40 个节点、C=0.5，80 步后经过一个完整周期。相位在这个特例中回到原位，振幅却乘上 cos(π/40) 的 80 次方，体现数值扩散。

图中的衰减可直接复算。取 $u_0(x)=\sin (2\pi x)$、$\Delta x=1/40$、$C=1/2$。此时 $G=e^{-i\theta /2}\cos (\theta /2)$，$\theta =2\pi /40$。八十步的相移为一整周，振幅则变为 ${\cos }^{80}(\pi /40)$。时间空间联合加密时误差才消失，单独强调稳定并不能说明当前网格足够准确。

可运行的周期节点格式

函数输入节点值、带符号速度、网格间距、时间步长和步数，返回新列表。每步只读取旧列表，不原地覆盖尚未使用的上游值。$M$ 个节点、$N$ 步需 $O(MN)$ 时间和 $O(M)$ 工作空间。它和常速度一阶有限体积法可以有同样代数形式，但这里保存的是节点值，后者保存单元平均值，解释和误差比较对象不同。

Python 3

import math


def upwind(initial, velocity, dx, dt, steps):
    u = list(map(float, initial))
    velocity, dx, dt = map(float, (velocity, dx, dt))
    if not u or not all(math.isfinite(v) for v in u + [velocity, dx, dt]):
        raise ValueError("nonempty finite input required")
    if dx <= 0 or dt <= 0 or not isinstance(steps, int) or steps < 0:
        raise ValueError("positive spacing and nonnegative integer steps required")
    c = abs(velocity) * dt / dx
    if not math.isfinite(c) or c > 1:
        raise ValueError("Courant number exceeds one")
    direction = -1 if velocity >= 0 else 1
    for _ in range(steps):
        new = [(1-c)*u[i] + c*u[(i+direction) % len(u)] for i in range(len(u))]
        if not all(math.isfinite(v) for v in new):
            raise OverflowError("nonfinite state")
        u = new
    return u


print(upwind([0, 1, 0, -1], 1, 1, 0.5, 1))  # [-0.5, 0.5, 0.5, -0.5]
N = 40
values = [math.sin(2*math.pi*i/N) for i in range(N)]
answer = upwind(values, 1, 1/N, 0.5/N, 80)
print(answer[10], math.cos(math.pi/40)**80)  # equal up to roundoff

变速度问题需区分非守恒输运 $u_t+a(x)u_x=0$ 和守恒形式 $u_t+(a(x)u{)}_x=0$。二者相差 $a^{\prime }(x)u$，应按模型选择节点差分或界面通量；不能只把常速度替换成节点数组就认为解的是同一方程。