有限体积法

有限体积法（finite volume method）通过控制体内的守恒量及其边界通量离散偏微分方程。它首先积分守恒定律，再用数值通量近似相邻控制体之间的交换。相同的界面通量以相反符号进入两个相邻单元，因此局部计算能够保持整个区域的离散守恒。

四个储格之间的输运

设一条环形通道分成四个宽为 $\Delta x=1$ 的单元，某物质以速度 $a=1$ 向右流动。单元平均浓度初值为 $(1,0,0,0)$。时间步长取 $\Delta t=0.5$，则半个单元的量在一步中流向右邻。更新为 $$(1,0,0,0)⟼(0.5,0.5,0,0)⟼(0.25,0.5,0.25,0).$$ 浓度逐渐铺开，但每一步总量仍为一。这个四格例子使用周期边界，最右单元流出的量会进入最左单元；若边界有流入或流出，区域总量应相应变化。

Courant 数为 0.5。每种颜色的柱高之和保持为一；物理单元宽均为一，图中窄柱仅为并列比较。数值扩散使尖脉冲逐步变宽。

从积分守恒律到数值通量

一维守恒律写为 $$u_t+f(u{)}_x=0.$$ 令单元 $I_i=[x_{i-1/2},x_{i+1/2}]$ 的宽度为 $\Delta x_i$，平均值为 ${\overline{u}}_i=\Delta x_i^{-1}{\int }_{I_i}u\mathrm{d}x$。对单元积分得 $$\frac{\mathrm{d}{\overline{u}}_i}{\mathrm{d}t}=-\frac{f(u(x_{i+1/2},t))-f(u(x_{i-1/2},t))}{\Delta x_i}.$$ 这一式子只在边界处需要通量。将真实边界通量改为单值的数值通量 ${\hat{f}}_{i+1/2}$，再用显式 Euler 推进： $${\overline{u}}_i^{n+1}={\overline{u}}_i^n-\frac{\Delta t}{\Delta x_i}({\hat{f}}_{i+1/2}^n-{\hat{f}}_{i-1/2}^n).$$ 所谓单值，是界面左侧和右侧共用同一个数，不能各算一份不相等的通量。Clawpack 的波传播与 Riemann 问题文档采用这一守恒观点。

对 $f(u)=au$，信号从上游到达界面。令 $a^{+}=max(a,0),a^{-}=min(a,0)$，迎风数值通量为 $$\hat{f}(u_L,u_R)=a^{+}{u}_L+a^{-}{u}_R.$$ 它满足一致性 $\hat{f}(v,v)=f(v)$。若 $a>0$、网格均匀，格式简化为 $${\overline{u}}_i^{n+1}=(1-C){\overline{u}}_i^n+C{\overline{u}}_{i-1}^n,C=\frac{a\Delta t}{\Delta x}.$$ 上文四格算例就是 $C=1/2$。若速度改为 $a=-1$，上游变成右邻，第一步为 $(0.5,0,0,0.5)$，不能只给原公式换一个负系数。

守恒与不产生新极值的证明

把一般更新式乘以单元宽度再求和，内部界面出现两次：左单元将其减去，右单元将其加上。它们逐项抵消，得到 $$\sum _i\Delta x_i{\overline{u}}_i^{n+1}=\sum _i\Delta x_i{\overline{u}}_i^n-\Delta t({\hat{f}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}-{\hat{f}}_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}).$$ 周期边界下最后两项也抵消，所以离散总量精确保留。这一证明与解是否光滑无关，因此有限体积法适合含激波的守恒问题；不过守恒本身还不能挑选物理正确的熵解。

在线性迎风例中，若 $0\le C\le 1$，新平均值是两个旧平均值的凸组合。因此它处于原最小值与最大值之间，非负数据保持非负，常数状态也保持不变。相反，取 $C=1.2$ 和初值 $(1,0,0,0)$，首步便出现 $(-0.2,1.2,0,0)$：总量虽仍是一，却已产生负浓度。守恒和稳定性是两个需要分别验证的性质。

重构、精度与数值扩散

最简单的有限体积法把每个单元内的值视为常数，只保存平均值。它的空间精度通常为一阶。较高阶方法在单元内重构线性或更高次多项式，将界面左右迹送入 Riemann 求解器，并配合高阶时间推进。重构必须保持给定的单元平均值；在不连续附近还需要限制器等机制控制振荡。

对光滑线性输运，把上面的迎风差分作 Taylor 展开，可得到主导修正项 $$u_t+a{u}_x=\frac{a\Delta x}{2}(1-C)u_{xx}+\text{更高阶项},a>0.$$ 其中用 $u_{tt}=a^2{u}_{xx}$ 消去了时间二阶导数。$C<1$ 时正的二阶项像扩散，使脉冲变宽；它解释图形，但该展开要求光滑，不能直接当作激波处的逐点误差公式。$C=1$ 时单元值恰好平移一格，对常速度的周期平均值更新可精确对齐网格。

周期网格的 Python 实现

函数输入有限的单元平均值、常速度、正网格宽度、正时间步长和非负步数，返回新的列表。它不修改输入，检查 $|a|\Delta t/\Delta x\le 1$。每步只访问相邻单元，$M$ 个单元、$N$ 步的时间为 $O(MN)$，工作空间为 $O(M)$。

Python 3

import math


def finite_volume_advection(initial, a, dx, dt, steps):
    u = [float(value) for value in initial]
    a, dx, dt = map(float, (a, dx, dt))
    if not u or not all(math.isfinite(v) for v in u + [a, dx, dt]):
        raise ValueError("nonempty finite data 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(a) * dt / dx
    if not math.isfinite(c) or c > 1:
        raise ValueError("Courant number must not exceed one")
    count = len(u)
    for _ in range(steps):
        flux = [max(a, 0) * u[i] + min(a, 0) * u[(i+1) % count]
                for i in range(count)]
        new = [u[i] - dt / dx * (flux[i] - flux[i-1]) for i in range(count)]
        if not all(math.isfinite(v) for v in new):
            raise OverflowError("nonfinite flux or state")
        u = new
    return u


print(finite_volume_advection([1, 0, 0, 0], 1, 1, 0.5, 2))
# [0.25, 0.5, 0.25, 0.0]
print(finite_volume_advection([1, 0, 0, 0], -1, 1, 0.5, 1))
# [0.5, 0.0, 0.0, 0.5]

周期条件只是此实现的边界选择。入口浓度、固壁和开边界应分别根据物理信息设置通量；将周期代码直接用于有入口出口的管道，会改变模型本身。非线性系统还需要相应波速、通量与稳定性分析。