Euler、Heun与可分离的采样误差

任务：把已学速度场与数值误差分开验收

上一章说明了应学习什么边缘速度，本章处理如何积分。令状态 $x(t)$ 满足 $dx/dt=v(x,t)$，初值来自标准正态，目标区间为零到一。对固定初值与固定场，ODE轨迹是确定的；样本多样性来自不同初始噪声，不需要每一步额外注入随机数。项目中的Euler和Heun都是ODE方法，不是反向SDE离散器。往其中随手加入噪声会改变密度演化，不能称为同一个求解器的小修。

设把区间分为 $N$ 个等长步，$h=1/N$，当前网格时刻为 $t_n=nh$，数值状态为 $x_n$。本课只用固定步长，明确记录速度场调用次数。固定步数、固定调用数和固定真实耗时是三种不同预算；网络推理、数据搬运与求解器额外操作都可能影响时间。本章会用有精确解的简单方程检查程序，再将同样代码用于学习场，避免用看起来像目标的样本反推实现一定正确。

一、Euler从积分恒等式得到

精确状态满足 $x(t+h)=x(t)+\int_t^{t+h}v(x(s),s)ds$。Euler用左端点速度近似整个小区间的积分，得到当前状态加步长乘当前速度。若场足够光滑，沿真实解求导得到第二阶项包含时间偏导与空间雅可比乘速度；Euler丢掉这一项，所以单步从精确状态出发的局部误差为步长平方量级。这个量级不意味着执行许多步后的总误差仍为平方，因为每步都会产生误差并传播。

假设速度对空间满足Lipschitz界，即状态相差多少，速度差至多是一个常数 $L$ 乘该差，且局部截断常数有界。误差递推可被上界为前一步误差乘 $(1+hL)$ 再加一个二阶局部项。把这条递推展开，约有一步数倒数那么多次局部误差累加，最终得到一阶全局误差。这个结论有条件，不适用于场不连续、数值已经爆炸或没有落在误差界成立区域的任意运行。

二、Heun用一次预测修正速度

Heun先做一个Euler预测点 $x_{\rm pred}=x_n+h v(x_n,t_n)$，再在预测点与下一时刻重新计算速度，取起点速度和预测终点速度的平均作更新。第二次调用必须使用下一时刻，而不是仍使用原时间；也必须使用预测状态，而不是直接拿原状态代入。一个名字叫Heun但这两项有误的函数，并不具有我们将推导的二阶性质。

把第二次速度在起点展开，得到 $k_2=v+h(\partial_t v+D_xv\,v)+O(h^2)$。代回平均更新后，恰好恢复精确解的二阶项，单步局部误差成为三阶；在类似的稳定与光滑条件下，全局误差为二阶。这是方法的阶数，不是说任意给定步长下Heun误差必然是Euler误差的固定比例，更不是说同样步数代表同样成本。

还有一个容易忽略的缓存错误：Heun第二次速度是在预测点算的，而最终状态是校正点，二者通常不同。不能不加推导就把本步第二次速度当成下一步的第一速度，以减少调用。某些专门方法支持端点导数复用，本项目的这套Heun实现没有该约定。调用次数测试会逐步累计，要求Euler为N次、Heun为二N次。

算例一：有精确解才能独立检查阶数

取标量方程 $dx/dt=x$，初值一，终点精确解是自然常数e。用步长一半做两步，Euler每步乘一点五，结果二点二五；Heun每步乘一点六二五，结果二点六四零六二五。相对精确值约二点七一八二八二，二者绝对误差分别约零点四六八三和零点零七七七。这里Heun用四次调用，Euler只有两次，不能直接把这次差异叫同预算优势。

若把Euler也给四次调用，步长四分之一，结果是一点二五的四次方，约二点四四一四零六，误差约零点二七六九。在这一个光滑例子和这一个预算上，Heun仍较准确；结论的范围就是这个例子。项目单元测试还比较十步与二十步误差比，Euler应趋近二、Heun应趋近四，从而核对阶数，而非只检查某个偶然的最终小数。

三、稳定性与学习误差是不同问题

对收缩方程 $dx/dt=-\lambda x$、$\lambda>0$，精确解随时间减小。Euler放大因子为 $1-\lambda h$，要求其绝对值小于一才能保持渐近稳定，因此步长必须小于二除以lambda。Heun的放大因子为 $1-\lambda h+(\lambda h)^2/2$，在负实轴上也受同样的区间限制。二阶方法没有自动解决所有刚性问题；若速度对位置变化特别剧烈，两种显式方法都可能需要很细步长。

即使对一个网络场积分得非常精确，若网络本身与正确边缘场有偏差，最后分布仍可能不对。反过来，拥有精确解析教师也不能弥补一步巨大Euler更新带来的离散化错误。项目同时提供学习场和解析场，并对相同初始噪声使用多种步数：解析场帮助定位求解器误差，网络对解析场的验证误差帮助定位学习误差，最终统计指标才检查输出分布。三类证据不可拼成一个模糊的“质量提升”。

算例二：精确教师加一步Euler也能彻底失败

本课独立线性桥在零时刻满足边缘速度 $v_*(x,0)=\mathbb E[X_1]-x=-x$，因为目标总体均值为零，初始点与目标独立。如果Euler只走一步，步长为一，则终点为 $x+(-x)=0$，每个不同噪声都被映到原点。目标却有四个簇，方差四点零九。这是完全可预测的数值失败，不是教师公式错误，也不是网络没学会；它直接暴露了把瞬时速度当整个区间恒定速度的错误。

用Heun走一步时，预测点先落到零，下一时刻教师在零处仍为零，校正结果为原噪声的一半，也不是四簇目标。因此“多看一次速度”不是任何步长下的万能补救。随着网格细化，两者才可能在正则条件下逼近正确流。我们会保留这两个粗网格输出，不因它们难看而删掉，因为它们是检查评价程序能否发现明显错误的阳性对照。

四、对齐比较与保存输入

比较两种求解器时，先冻结同一个模型，再固定相同的一批初始噪声。若还同时换网络、训练步数或随机种子，最终差异便不再只来自积分法。细网格解可作为近似数值参照，但不是数学真值；应继续加密网格，检查参照本身是否稳定。逐点误差检查需要配对初始样本，分布指标则不要求两种方法逐点对应某个真实目标点。

项目评价命令支持显式给出步数一、十、二十、四十，这样既能看粗细变化，也能比较四十次Euler调用与二十步Heun的四十次调用。记录中的耗时是该次运行的观察，不是通用性能定律。三组采样种子共享一个训练模型，只反映采样随机性；它们不是三次独立训练，不足以据此声称训练稳定性已经被重复验证。

可以把两类误差关系再写得更清楚。假设在所有相关位置上，学习场与真实场的差都不超过某个统一界，真实场又有Lipschitz常数，那么两条精确ODE轨迹的距离满足“距离导数不超过常数乘距离加场误差”的不等式。积分后得到随时间放大的场偏差上界，再加数值求解的离散误差。然而验证集的平均平方场误差并不是这个统一界；它不能排除极少数位置出现很大错误，所以不能把一个很小的平均MSE直接代成严格的逐轨迹保证。

项目中的细网格配对还可发现另一类失败：若每次调用积分器都重新抽初始点，那么输出差同时包含输入变化和求解器误差。程序应先抽一次初始张量，再将同一张量传给不同方法，积分器内部只克隆而不改写调用者输入。单元测试会检查原张量未被原地修改。这个对象身份约定看似细小，却决定多步数对照是否真正在比较同一初值问题。

from math import e
def solve(n,method):
    x=1.0;h=1/n;calls=0
    for i in range(n):
        k1=x;calls+=1
        if method=="euler":x+=h*k1
        else:
            k2=x+h*k1;calls+=1;x+=h*(k1+k2)/2
    return x,calls
assert solve(2,"euler")== (2.25,2)
assert solve(2,"heun")== (2.640625,4)
for method in ["euler","heun"]:
    a,_=solve(10,method);b,_=solve(20,method)
    print(method,a,b,abs(a-e)/abs(b-e))

练习与完整解答

# integrate
def integrate(field,initial,grid,method='heun'):
    raise NotImplementedError('ode')

# integrate
def integrate(field,initial,grid,method='heun'):
 if method not in ['euler','heun']or len(grid)<2:raise ValueError('method/grid')
 times=[float(t)for t in grid]
 if not all(torch.isfinite(torch.tensor(t))for t in times)or any(b<=a for a,b in zip(times,times[1:])):raise ValueError('strict finite increasing grid')
 x=initial.clone();path=[x.clone()];calls=0
 for t,end in zip(times,times[1:]):
  h=end-t;k1=field(x,t);calls+=1
  if method=='euler':x=x+h*k1
  else:k2=field(x+h*k1,end);calls+=1;x=x+h*(k1+k2)/2
  path.append(x.clone())
 return torch.stack(path),calls