使用Sympy求解含DiracDelta积分时的变量处理问题

你需要求解的含狄拉克δ函数的积分形式如下：

式中参数定义：

• v 为运动速度
• C_i 为时间步t₀对应的浓度值
• 积分变量为t₀（代码中简写为to），积分区间为[0, t_s]

• 接口混用错误：scipy.integrate.quad是纯数值积分函数，无法直接解析SymPy定义的DiracDelta符号对象，不能直接传入SymPy符号表达式做数值计算
• 传参逻辑错误：quad要求被积函数的第一个入参必须是积分变量，其余入参通过args按顺序传入，你原有代码的函数参数顺序、args传参顺序不匹配，且把观测时刻t错误固定为x/v，实际t就是你遍历的积分上限ts
• 方法选型冗余：狄拉克δ函数有明确的积分筛选性质，完全可以先通过符号推导得到解析解，再转数值计算，效率和精度远高于通用数值积分

利用狄拉克δ函数的积分筛选性质：
$$\int_a^b f(\tau)\delta(g(\tau))d\tau = \sum_{\tau_0 \in [a,b], g(\tau_0)=0} \frac{f(\tau_0)}{|g'(\tau_0)|}$$
针对你的被积函数，$g(t_0) = x/v - t + t_0$，唯一零点为$t_0 = t - x/v$，且$g'(t_0)=1$，因此积分结果可以直接判定：

• 若零点$t_0 = t - x/v$落在积分区间[0, t_s]内，积分结果等于该点对应的浓度c
• 若零点不在积分区间内，积分结果为0

你可以先用SymPy做符号积分验证结论：

import sympy as smp
t0, x, c, v, t, ts = smp.symbols('t0 x c v t ts', real=True)
integrand = c * smp.DiracDelta(x/v - t + t0)
# 对t0在[0, ts]区间做符号积分
res = smp.integrate(integrand, (t0, 0, ts))
print(res)
# 输出：c*Heaviside(ts - t + x/v) - c*Heaviside(-t + x/v)

输出结果中的Heaviside为单位阶跃函数，和我们手动推导的判定逻辑完全一致。

基于解析结果写数值计算逻辑，天然支持c为常量或与t0对齐的数值数组：

import numpy as np

def calc_delta_integral(c, v, x, t_arr):
    """
    计算各观测时刻t对应的积分结果
    参数：
        c: 浓度值，支持标量常量，或与t0采样点一一对应的数组
        v: 速度
        x: 位置
        t_arr: 观测时刻数组，即每个积分的上限值
    """
    results = np.zeros_like(t_arr, dtype=np.float64)
    # 计算每个观测时刻对应的δ函数零点t0
    t0_zero = t_arr - x / v
    # 筛选零点落在积分区间[0, t]内的有效位置
    valid_mask = (t0_zero >= 0) & (t0_zero <= t_arr)
    if np.isscalar(c):
        results[valid_mask] = c
    else:
        # 若c是t0对应的数组，按t0索引取对应浓度即可，需保证c的索引和t0采样对齐
        # 示例默认t0步长为1，和你代码中t_arr的步长一致
        valid_t0_idx = t0_zero[valid_mask].astype(int)
        results[valid_mask] = c[valid_t0_idx]
    return results

# 测试用例
c_const = 0.5
v = 0.1
x = 10
t_arr = np.arange(0, 1000, 1)
res = calc_delta_integral(c_const, v, x, t_arr)
# 结果符合物理预期：t < 100（即x/v）时结果为0，t>=100时结果为0.5

如果后续扩展被积函数形式、必须用数值积分，不要混用SymPy符号对象，直接用窄峰近似狄拉克δ函数即可：

from scipy.integrate import quad
import numpy as np

def f_num(t0, c, v, x, t):
    # 用半宽eps的矩形近似δ函数，eps取时间步长的1/2即可
    eps = 0.5
    t0_center = t - x/v
    return c/(2*eps) if abs(t0 - t0_center) < eps else 0

c_const = 0.5
v = 0.1
x = 10
t_arr = np.arange(0, 1000, 1)
integrals = [quad(f_num, 0, ts, args=(c_const, v, x, ts))[0] for ts in t_arr]

注意：该方法精度和计算效率都远低于解析解法，非必要不使用。

内容的提问来源于stack exchange，提问作者Luna Paluna