1

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.

Luna Paluna
  • 63
  • 5
  • `quad` is a numerical integration routine and simpy a symbolic algebra package and they just can't be mixed freely with a numerical integration routine like quad. Also x is missing from args, quad takes single variable functions. the limit of integration should be ts too (if one could just mix things). Finally, if C is a constant you just have the integral of a dirac delta which is trivial no? – GuillemB Jul 07 '22 at 15:18
  • Yes, you are right in this case it will be trivial but after I will need to calculate the mean of cf which sort of complicates things. I'll fix the x in args and the ts in the limits. What I expected to obtain from this is the ci value at the given tt. – Luna Paluna Jul 07 '22 at 15:30

1 Answers1

2

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.

Oscar Benjamin
  • 12,649
  • 1
  • 12
  • 14