search

JuliaPy / SymPy.jl

Julia interface to SymPy via PyCall
http://juliapy.github.io/SymPy.jl/
MIT License
268 stars 62 forks source link

Error calculating the inverse Laplace transform using floats #497

Closed runjaj closed 1 year ago

runjaj commented 1 year ago

Description of the problem

I'm trying to calculate the inverse Laplace transform of a function that has floating coefficients:

using SymPy
@syms s t::real
sympy.inverse_laplace_transform(0.2/(20s^2 +3.52s+1)*1/s, s, t)

I get a long error message, you can find it bellow.

Thankfully the problem is easy to solve:

julia> sympy.inverse_laplace_transform(sympy.nsimplify(0.2/(20s^2 +3.52s+1)*1/s, rational=true), s, t)
⎛      11⋅t                                            ⎞  -11⋅t      
⎜      ────                                            ⎟  ──────     
⎜      125                ⎛√2641⋅t⎞           ⎛√2641⋅t⎞⎟   125       
⎜2641⋅ℯ     - 22⋅√2641⋅sin⎜───────⎟ - 2641⋅cos⎜───────⎟⎟⋅ℯ      ⋅θ(t)
⎝                         ⎝  250  ⎠           ⎝  250  ⎠⎠             
─────────────────────────────────────────────────────────────────────
                                13205

I have two questions:

  1. Why does this error happen? I'm just curious about it
  2. Is there a way to avoid writing sympy.nsimplify(...)? I know that Maxima makes this change by default. I have no problem writing it, but I'm afraid it will confuse my students.

Anyway, thanks for your amazing work creating and mataining SymPy.jl!

Error message

ERROR: PyError ($(Expr(:escape, :(ccall(#= /Users/javier/.julia/packages/PyCall/twYvK/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'sympy.polys.polyerrors.PolynomialError'>
PolynomialError(PolynomialError('RisingFactorial(_t + 1, 1.0) contains an element of the set of generators.'))
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 2097, in inverse_laplace_transform
    return InverseLaplaceTransform(F, s, t, plane).doit(**hints)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 171, in doit
    fn, T = self._try_directly(**hints)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 135, in _try_directly
    T = self._compute_transform(self.function,
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 2044, in _compute_transform
    return _inverse_laplace_transform(F, s, t, self.fundamental_plane, **hints)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 255, in wrapper
    res = func(*args, **kwargs)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 1962, in _inverse_laplace_transform
    f = Add(*[_inverse_laplace_transform(X, s, t, plane, simplify)\
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 1962, in <listcomp>
    f = Add(*[_inverse_laplace_transform(X, s, t, plane, simplify)\
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 255, in wrapper
    res = func(*args, **kwargs)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 1967, in _inverse_laplace_transform
    f, cond = inverse_mellin_transform(F, s, exp(-t), (None, S.Infinity),
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 936, in inverse_mellin_transform
    return InverseMellinTransform(F, s, x, strip[0], strip[1]).doit(**hints)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 171, in doit
    fn, T = self._try_directly(**hints)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 135, in _try_directly
    T = self._compute_transform(self.function,
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 875, in _compute_transform
    return _inverse_mellin_transform(F, s, x, strip, **hints)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 255, in wrapper
    res = func(*args, **kwargs)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/integrals/transforms.py", line 803, in _inverse_mellin_transform
    h = hyperexpand(G)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/simplify/hyperexpand.py", line 2494, in hyperexpand
    return f.replace(hyper, do_replace).replace(meijerg, do_meijer)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/core/basic.py", line 1581, in replace
    rv = walk(self, rec_replace)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/core/basic.py", line 1566, in walk
    rv = F(rv)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/core/basic.py", line 1574, in rec_replace
    v = _value(expr, result)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/core/basic.py", line 1503, in <lambda>
    _value = lambda expr, result: value(*expr.args)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/simplify/hyperexpand.py", line 2490, in do_meijer
    r = _meijergexpand(G_Function(ap[0], ap[1], bq[0], bq[1]), z,
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/simplify/hyperexpand.py", line 2375, in _meijergexpand
    slater2, cond2 = do_slater(tr(func.bm), tr(func.an), tr(func.bq), tr(func.ap),
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/simplify/hyperexpand.py", line 2301, in do_slater
    hyp = _hyperexpand(Hyper_Function(nap, nbq), harg, ops,
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/simplify/hyperexpand.py", line 2042, in _hyperexpand
    formula = try_lerchphi(func)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/simplify/hyperexpand.py", line 1779, in try_lerchphi
    part = apart(numer/denom, t)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/utilities/decorator.py", line 35, in threaded_func
    return func(expr, *args, **kwargs)
  File "/Users/javier/.julia/conda/3/x86_64/lib/python3.10/site-packages/sympy/polys/partfrac.py", line 84, in apart
    raise PolynomialError(msg)

Stacktrace:
  [1] pyerr_check
    @ ~/.julia/packages/PyCall/twYvK/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/.julia/packages/PyCall/twYvK/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/exception.jl:96
  [4] macro expansion
    @ ~/.julia/packages/PyCall/twYvK/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:473 [inlined]
  [7] __pycall!
    @ ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{Sym, Sym, Sym}, nargs::Int64, kw::Ptr{Nothing})
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{Sym, Sym, Sym}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:11
 [10] (::PyCall.PyObject)(::Sym, ::Vararg{Sym}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:86
 [11] (::PyCall.PyObject)(::Sym, ::Vararg{Sym})
    @ PyCall ~/.julia/packages/PyCall/twYvK/src/pyfncall.jl:86
 [12] top-level scope
    @ REPL[11]:1

Description of the system

My system is:

Julia Version 1.9.0-rc2 Commit 72aec423c2a (2023-04-01 10:41 UTC) Platform Info: OS: macOS (x86_64-apple-darwin21.4.0) CPU: 8 × Intel(R) Core(TM) i7-1060NG7 CPU @ 1.20GHz WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-14.0.6 (ORCJIT, icelake-client) Threads: 1 on 8 virtual cores

And this is my environment:

[f6369f11] ForwardDiff v0.10.35 [b964fa9f] LaTeXStrings v1.3.0 [91a5bcdd] Plots v1.38.10 [438e738f] PyCall v1.95.1 [f2b01f46] Roots v2.0.12 [24249f21] SymPy v1.1.8

jverzani commented 1 year ago

For this particular problem, all SymPy.jl does is pass along arguments to the underlying sympy functions. So, the issue sits there. Simplification can often help when using SymPy, as can assigning the variable type. You have already done both.

jverzani commented 1 year ago

Another response, as I didn't realize I hadn't closed this. All nsimplify does in this case is convert floats to rationals. I'd suggest having your students do that using //, as in sympy.inverse_laplace_transform((2//10)/(20s^2 +(352//100)s+1)*1/s, s, t). There are a few conversions where you lose the exactness that students expect (pi is converted to sympy.pi, but 2pi is a float in Julia, so a float in SymPy on conversion). So you have to explain that aspect anyways, and rationals are, perhaps, the easiest for them to be mindful of.

runjaj commented 1 year ago

Thanks for the answer.