使用 sympy 求解狄拉克 δ 积分

Question

我正在尝试使用 sympy 求解以下积分：

v 是速度，C_i 是时间步 t_0 处的浓度。这是我到目前为止所拥有的：

import sympy as smp
from scipy.integrate import quad
to,x = smp.symbols(('to','x'), real=True)

def f(to,x,c,v,t):
    return c*smp.DiracDelta((x/v) - t + to)

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

我不确定如何处理 dt0 和变量。任何见解都值得赞赏。在本例中，我使用 c_arr 作为常量，以使其更简单，否则，它将是一个值数组。

Answer 1

这是使用 SymPy 符号计算积分的方法：

In [53]: C_i = Function('C_i')

In [54]: t, t0, x, v = symbols('t, t0, x, v', positive=True)

In [55]: g = lambda x, t: DiracDelta(x/v - t + t0)

In [56]: C_f = Integral(C_i(t0)*g(x,t-t0), (t0, 0, t))

In [57]: C_f
Out[57]: 
t                              
⌠                              
⎮         ⎛            x⎞      
⎮ Cᵢ(t₀)⋅δ⎜-t + 2⋅t₀ + ─⎟ d(t₀)
⎮         ⎝            v⎠      
⌡                              
0                              

In [58]: C_f.doit()
Out[58]: 
    ⎛t⋅v - x⎞  ⎛-(t⋅v - x) ⎞     ⎛t⋅v - x⎞  ⎛    t⋅v - x⎞
  Cᵢ⎜───────⎟⋅θ⎜───────────⎟   Cᵢ⎜───────⎟⋅θ⎜t - ───────⎟
    ⎝  2⋅v  ⎠  ⎝    2⋅v    ⎠     ⎝  2⋅v  ⎠  ⎝      2⋅v  ⎠
- ────────────────────────── + ──────────────────────────
              2                            2             

In [59]: C_f.doit().simplify()
Out[59]: 
⎛     ⎛-t⋅v + x⎞⎞   ⎛t⋅v - x⎞
⎜1 - θ⎜────────⎟⎟⋅Cᵢ⎜───────⎟
⎝     ⎝  2⋅v   ⎠⎠   ⎝  2⋅v  ⎠
─────────────────────────────
              2 

这里的 theta 是 Heaviside 函数。现在你只需要对这个积分进行数值计算，你可以用

lambdify

:

来完成

In [73]: C_fs = C_f.doit().simplify()

In [74]: f = lambdify((x, v, t), C_fs, ['scipy', {'C_i': lambda e: 0.5}])

In [75]: f(10, 0.1, t_arr)
Out[75]: 
array([0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,
       0.   , 0.125, 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 ,
       0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 ,
       0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 ,
       0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 ,
...

如果您有

C_i

的值数组，那么您可以使用 interpolate 将其转换为可调用函数。