Search code examples
pythonsympynumerical-integration

Solving an integral with dirac delta using sympy

I'm trying to solve the following integral using sympy:

Equation

v is velocity, C_i is the concentration at the time step t_0. This is what I have so far:

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]

I'm not sure how to handle the dt0 and the variables. Any insight is appreciated. I'm using c_arr as a constant in this case to make it simpler, otherwise, it would be an array of values.

Solution

  • This is how you can compute the integral symbolically using 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

    The theta here is the Heaviside function. Now you just need to evaluate that integral numerically which you can do with 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 ,
...

    If you have an array of values for C_i then you can just use interpolate to turn that into a callable function.