SymPy - Kronecker Delta Function Evaluation-CodePudding

I am using SymPy for the numerical analysis of large sets of equations. Part of my equation contains a Kronecker Delta function acting as an impulse such that when q = 0 -> dirac_delta = 1, otherwise dirac_delta = 0. I need to perform this for values of q = - 10 -> 10 in integer steps of 1.

A simplified example of my code is:

import sympy as sp
import numpy as np

modules = ["numpy", "sympy"]

# Create symbols
q = sp.symbols('q', integer=True)

# Create P_q symbol as a function of q
P_q = sp.symbols('P_q', cls=sp.Function)
P_q = P_q(q)

# Define Equation
# Simplified example KroneckerDelta - when q = 0, P_q = 1, otherwise P_q = 0
P_q_eq = sp.Eq(P_q, sp.KroneckerDelta(0,q))
P_q = sp.KroneckerDelta(0,q)
display(P_q_eq)

# Create a lambda function for fast numerical calculation 
lam_P_q = sp.lambdify(q, P_q, modules)

# Define the values of q
num_points = 21
data = np.linspace(-10, 10, num_points, dtype=int)
#print(data)

ans = lam_P_q(data)
print(ans)

On run I receive an error:

ValueError Traceback (most recent call last) in 36 #print(data) 37 ---> 38 ans = lam_P_q(data) 39 print(ans)

in _lambdifygenerated(q) 1 def _lambdifygenerated(q): ----> 2 return ((1 if 0 == q else 0))

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

My understanding is that I need to add the .any() or .all() because the lambdify is comparing the array q, to a single value 0. So when I modify the input data with .any() or .all() it then returns a single value.

However, I require a response of 0 or 1 for each value of q - such that it is an impulse response depending on the value of q.

print(q)

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

I have tried to provide the "0" comparison value as an array of zeroes with equal length to q, but this did not work.

I know that this can be generated with the scipy signal function "signal.unit_impulse(21, 'mid')", but I am unsure how to implement this in the SymPy format for lambdify to output q as outlined just above. I have tried to create a custom module to replace the sp.KroneckerDelta function to perform this but couldn't get a working solution (likely due to a poor implementation):

def impulse(num, p):
    val = signal.unit_impulse(num, 'mid')
    j = p   (num-1)/2
    kron_p = int(val[j])
        
    return kron_p

kronecker = {'KroneckerDelta': impulse(21,q)}
modules = [kronecker, "numpy", "sympy"]

Do I need to substitute the values of q into the lambdify function instead - if so, how do I specify a range of values to substitute?

I feel like I am doing something fundamentally wrong in my approach and would appreciate help getting what I thought would be a relatively simple thing to do in SymPy to work. (I am quite new to SymPy and definitely still trying to get my head around it). Thanks.

CodePudding user response：

This is a bug in lambdify. Please open a sympy issue:

https://github.com/sympy/sympy/issues

You can work around it by rewriting the KroneckerDelta to a Piecewise:

In [12]: P_q
Out[12]: 
δ   
 0,q

In [13]: P_q.rewrite(Piecewise)
Out[13]: 
⎧0  for q ≠ 0
⎨            
⎩1  otherwise

In [14]: f = lambdify(q, P_q.rewrite(Piecewise), modules)

In [15]: f(data)
Out[15]: 
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0.])

This is the generated code:

In [16]: import inspect

In [18]: print(inspect.getsource(f))
def _lambdifygenerated(q):
    return select([not_equal(q, 0),True], [0,1], default=nan)

CodePudding user response：

Ur created variable data is a list type looking like this:

data = [-10  -9  -8  -7  -6  ... 8   9  10]

lam_P_q(data) is not meant to take a variable with type list. It is meant to take a single number. In your case for example: data = -10

To feed lam_P_q(data) with all single numbers from the data variable, a for loop can be used. Below is your modified example:

import sympy as sp
import numpy as np
from IPython.display import display


modules = ["numpy", "sympy"]

# Create symbols
q = sp.symbols('q', integer=True)

# Create P_q symbol as a function of q
P_q = sp.symbols('P_q', cls=sp.Function)
P_q = P_q(q)
print(P_q)

# Define Equation
# Simplified example KroneckerDelta - when q = 0, P_q = 1, otherwise P_q = 0
P_q_eq = sp.Eq(P_q, sp.KroneckerDelta(0,q))
P_q = sp.KroneckerDelta(0,q)
display(P_q_eq)

# Create a lambda function for fast numerical calculation 
lam_P_q = sp.lambdify(q, P_q, modules)

# Define the values of q
num_points = 21
data = np.linspace(-10, 10, num_points, dtype=int)
print(data)
#feeding lam_P_q with numbers and saving answers into the list results
results=[]
for num in data:
    ans = lam_P_q(num)
    results.append(ans)

print(results)