jverzani / SymPyCore.jl

Package to help connect Julia with the SymPy library of Python
https://jverzani.github.io/SymPyCore.jl/
MIT License
4 stars 5 forks source link

issues with sinc function #31

Closed Djuleeo closed 9 months ago

Djuleeo commented 11 months ago

i have an issue with lambdification of the sinc function, due to the different definitions of the sinc


using SymPy
x = symbols("x", real=true)

fun      = SymPy.sinc(x/pi)
f        = lambdify(fun)
x_test   = 1.0
body     = convert(Expr, fun)                                           # correctly translated  :(x ^ -1 * sin(x))
f(x_test) - sin(x_test)/x_test                                          # is equal to zero, as expected
  #
#### BUT if i use the library to rewrite it :

fun2    = sinc(x/pi).rewrite(sinc)      
body2   = convert(Expr, fun2)               # :(sinc(x))
f2  = lambdify(fun2)                            # becomes sin(π x)/ πx   !!!!

f2(x_test) - sin(x_test)/x_test                             # different than 0 !
jverzani commented 11 months ago

So, I think this is a different in defintions:

In Julia, sinc is sin(pi x)/(pi x) and defined by

Base.sinc(x::Sym) = iszero(x) ? one(x) : sinpi(x)/(PI*x)

In Sympy (https://docs.sympy.org/latest/modules/functions/elementary.html#sympy.functions.elementary.trigonometric.sinc) it is sin(x)/x.

You are basically finding the difference between sympy.sinc (sinc(x)) and SymPy.sinc (sinpi(x)/(PI*x)). To smooth this out, we could modify the walk_expression function to handle the case of sinc (the first one above), though perhaps just documenting it better would be the best case forward.

Djuleeo commented 11 months ago

maybe it could be changed such that:

code                              symbolic
sinc(x)              ->    sinc( pi x ) (= sin(pi x) /(pi x))
SymPy.sinc(x)        ->    sinc(x) (= sin(x) /x)                              using the sympy.sinc definition
as currently  there is no difference between the two:
sinc(x)              ->    sinc( pi x ) 
SymPy.sinc(x)        ->    sinc( pi x)                  

provided that the difference is well documented.

in any case, in the code above the issue is that i get the wrong version of sinc when i lambdify:

code                                           symbolic
SymPy.sinc(x/pi)                       ->     sin(x)/x                    # uses the sinpi(x)/(pi x) definition
SymPy.sinc(x/pi).rewrite(sinc)         ->     sinc(x)                     # symbolic manipulation, sin(x)/x definition
everything is good so far, but if i lambdify
symbolic                                                  numeric
sinc(x) (=sin(x)/x)                      ->     sinc(x)  (=sin(pi x)/(pi x))  !!!!
instead it should be:
 sinc(x) (=sin(x)/x)                     ->     sinc(x/pi)  ( =sin(pi x/pi )/(pi x/pi) = sin(x)/x   )

in the current way, the net result is that using .rewrite(sinc) or a .simplify() that leads to a (symbolic) sinc and then lambdifying gives back a different function, when these methods are meant to change the form, but not modify the math.

jverzani commented 11 months ago

I think the issue is the distinction between sinc(x), SymPy.sinc(x), and sympy.sinc(x). The first two are the same, the latter different (it being the python libraries version sin(x)/x.) However, lambdify rewrites sympy.sinc(x) to Julia's sinc and so we get:

sympy.sinc(x)(1) --> sin(1)/1 lambdify(sympy.sinc(x))(1) --> sin(pi*1)/(pi*1)

It seems like the proper thing to do would be to modify lambdify to reexpress sinc(x) using the Python definition, not the Julia one. This would be the same as this:

SINC(x) = sinc(x/pi)
lambdify(sympy.sinc(x), fns=Dict("sinc"=>SINC))(1)