在科学计算与机器学习交叉领域，一个反复出现的尴尬场景是：模型在训练数据上拟合误差极小，预测出的轨迹却明显违反物理定律。例如预测的摆锤运动能量随时间不断增大，或者流体模拟中质量凭空消失。这类问题被称为物理定律违反（physics law violation）。单纯依赖数据驱动的符号回归往往难以避免该问题，而数值求解器本身虽遵守离散格式的守恒性，却缺乏从数据中发现新表达式的能力。将两者结合，是近年来科学发现领域一条行之有效的技术路线。
为什么纯数据驱动的符号回归会违反物理定律
符号回归（Symbolic Regression）的目标是从数据中搜索出一个数学表达式，使得表达式在给定输入上的输出与观测值尽可能接近。经典实现如遗传编程（GP）或近期流行的 PySR，都以拟合误差作为主要适应度指标。问题在于，拟合误差小并不等于物理规律正确。
第一，观测数据通常是离散采样点，采样点之间表达式可以任意波动，只要恰好穿过采样点即可。一条能量不守恒的曲线完全可以在采样点上与真实轨迹重合。第二，适应度函数只衡量逐点误差，不衡量导数结构、守恒量或对称性，搜索算法自然不会偏好满足这些约束的解。第三，数据中往往含有噪声，噪声会引导搜索偏向更复杂但不合物理的表达式。
例如对简谐振子数据做符号回归，可能得到 x(t) = sin(1.05*t) + 0.02*t 这样的表达式。多出来的线性项在短期内降低了拟合误差，但长期演化会导致振幅持续增大，违反能量守恒。这类错误在没有物理验证时几乎无法被自动发现。
数值求解器能提供什么：守恒性与一致性验证
数值求解器（如 Runge-Kutta 方法、辛积分器、有限元求解器）在计算数学中积累了大量保证物理一致性的技术。辛积分器（Symplectic Integrator）在长时间积分哈密顿系统时能保持能量近似守恒；守恒格式能保证流体质量与动量的离散守恒。这些性质正是符号回归所缺失的验证手段。
结合的核心思路是：符号回归负责产生候选方程，数值求解器负责检验候选方程的演化行为。具体做法是将候选表达式视为一个微分方程的右端项，用数值求解器从初始条件出发积分一段时间，然后比较积分结果与真实数据的偏差，以及守恒量的漂移程度。这样评估不再是逐点的，而是轨迹层面的，能暴露表达式在采样点之间的病态行为。
以简谐振子为例，把候选表达式 dx/dt = f(x, v) 用四阶 Runge-Kutta 积分后，计算总能量 E = 0.5*v^2 + 0.5*x^2 随时间的漂移。若漂移超过阈值，即便逐点拟合误差很小，也应淘汰该候选。这种双重评估显著提高了发现表达式的可靠性。
在适应度函数中嵌入物理约束惩罚项
最直接的结合方式是改造符号回归的适应度函数，在数据误差之外加入物理约束惩罚项。常见约束包括三类：守恒量约束、对称性约束和边界条件约束。
守恒量约束要求某个量沿轨迹保持不变，例如能量、质量或电荷。对称性约束要求方程在坐标变换（如平移、旋转）下形式不变。边界条件约束则限制表达式在特定区域的取值范围，比如浓度不能为负。惩罚项的权重需要仔细调节：权重过小约束形同虚设，权重过大则搜索空间被压缩得难以找到拟合解。
import numpy as np
from scipy.integrate import solve_ivp
def fitness_with_physics(expr, t_train, x_train, x0):
# expr: 可调用的候选右端项函数 f(t, x)
# 先做数值积分，得到候选轨迹
sol = solve_ivp(expr, (t_train[0], t_train[-1]), x0,
t_eval=t_train, method='RK45',
rtol=1e-8, atol=1e-10)
if not sol.success:
return 1e6 # 积分失败直接给极大惩罚
x_pred = sol.y[0]
# 数据误差项：均方误差
mse = np.mean((x_pred - x_train) ** 2)
# 物理惩罚项：能量漂移，简谐振子能量 E = 0.5*x^2 + 0.5*v^2
E = 0.5 * sol.y[0] ** 2 + 0.5 * sol.y[1] ** 2
energy_drift = np.max(np.abs(E - E[0])) / max(abs(E[0]), 1e-12)
# 加权组合：lambda 控制物理约束强度
lam = 10.0
return mse + lam * energy_drift
上面的代码展示了惩罚项设计的基本骨架。solve_ivp 承担了数值求解器的角色，把逐点拟合升级为轨迹拟合。energy_drift 度量候选方程积分后的能量漂移比例，作为物理违反程度的量化指标。实际工程中还可以加入导数平滑性惩罚、稀疏性正则（如 SINDy 方法中的 L1 正则）来进一步约束搜索方向。
完整的结合流程与工具选型
一套实用的结合流程分为四个阶段。第一阶段是数据准备，最好使用密集采样的轨迹数据而非散乱的点集，因为轨迹数据天然包含动力学信息。第二阶段是符号回归搜索，推荐使用 PySR 或 gplearn，PySR 性能更好且支持自定义损失函数，便于嵌入物理惩罚。第三阶段是候选验证，对搜索得到的前若干个候选表达式，分别用高精度数值求解器积分并检查守恒量。第四阶段是模型修正，对轻微违反物理的候选可以做符号化简或手工调整。
from pysr import PySRRegressor
model = PySRRegressor(
niterations=200,
binary_operators=["+", "-", "*", "/"],
unary_operators=["sin", "cos", "exp"],
# 自定义带物理约束的损失函数
loss_function="""
function loss(prediction, target)
data_err = (prediction - target) .^ 2 |> mean
# 对预测值做二阶差分，惩罚剧烈波动（近似能量漂移检测）
smooth = diff(diff(prediction)) .^ 2 |> mean
return data_err + 10.0 * smooth
end
""",
population_size=100,
maxsize=25,
)
model.fit(X, y)
print(model.sympy()) # 输出可解释的符号表达式
在工具选型上，除了 PySR，还可以考虑 AI Feynman（内置物理先验的符号回归）和 SINDy（稀疏识别动力学方程，天然适合与数值积分结合验证）。数值求解一侧，SciPy 的 solve_ivp 适合常微分方程，偏微分方程则需要 FEniCS 或 Dedalus 这类专业求解器。对于哈密顿系统，建议使用辛积分器（如 solve_ivp 配合自定义格式或第三方库 diffrax 中的辛方法），它在长时间演化中能更好地保持结构。
实践中的常见陷阱与调优建议
第一个陷阱是数值积分的刚性问题。某些候选表达式会产生刚性右端项，普通显式方法需要极小步长，搜索效率急剧下降。解决办法是自适应切换方法，例如检测到刚性时改用 method='Radau' 或 'BDF'，或者直接给积分失败的候选一个大惩罚快速淘汰。
第二个陷阱是惩罚权重的平衡。实践中推荐课程学习策略：搜索初期物理惩罚权重较小，让算法先找到拟合方向；后期逐步增大权重，把不符合物理的候选挤出种群。PySR 支持在回调中动态修改损失权重，实现这种渐进式约束。
第三个陷阱是评估成本。每个候选都要数值积分，计算量远高于逐点评估。缓解手段包括：先用小子集和粗精度积分做初筛，只对进入前 K 名的候选做高精度验证；利用并行评估，PySR 默认的多进程架构可以直接受益。此外，约束条件尽量写成可微形式，若后续想用梯度方法微调表达式参数，可微的惩罚项可以通过自动微分无缝接入。
综合来看，符号回归提供了发现物理规律的表达能力，数值求解器提供了验证物理一致性的裁判能力，两者结合让科学发现从单纯的拟合走向了受物理定律约束的可靠推断。这条路线在动力学系统发现、材料物性建模、气候参数化等场景中都已展现出实际价值，值得相关方向的工程与科研人员深入掌握。