sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.94k stars 4.42k forks source link

use numerical solver automatically? #3571

Open certik opened 16 years ago

certik commented 16 years ago

bc.. This is taken from: http://trac.sagemath.org/sage_trac/ticket/1235 In [1]: e = exp(1)

In [2]: t = Symbol("t")

In [3]: a = 0.004_(8_e(-300_t) - 8_e(-1200t))(720000_e(-300_t)-11520000_e(-1200t)) + 0.004(9600_e(-1200_t) - 2400_e(-300t))*2

In [4]: a Out[4]: 2
⎛ -300_t -1200t⎞ ⎛
0.004000000000000000083266726847
⎝- 2400_ℯ + 9600*ℯ ⎠ + ⎝- 115200

-1200*t           -300*t⎞ ⎛                                   -1200*t     

00_ℯ + 720000ℯ ⎠⎝- 0.03200000000000000066613381478_ℯ + 0.

                           -300*t⎞

03200000000000000066613381478*ℯ ⎠

In [5]: solve(a==0, t)

Not a polynomial equation. Can't solve it, yet.Traceback (most recent call last)

/home/ondra/sympy/

/home/ondra/sympy/sympy/solvers/solvers.py in solve(eq, syms, simplified) 73 solutions = list(set(roots(equ, syms[0]))) 74 except PolynomialException: ---> 75 raise "Not a polynomial equation. Can't solve it, yet." 76 77 if simplified == True:

Not a polynomial equation. Can't solve it, yet.: None

But it should call the iterative solver to return at least an approximative answer. Or it should return some kind of RootOf() or something, so that the user can manipulate with the answer further.

p. Original issue for "#3571":https://github.com/sympy/sympy/issues/3571: "http://code.google.com/p/sympy/issues/detail?id=472":http://code.google.com/p/sympy/issues/detail?id=472

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

smichr commented 12 years ago
This shouldn't be necessary since the solver can solve this equation. But probably a hint is needed to make it work:

>>> from sympy.solvers.solvers import _tsolve as ts
>>> e=exp(1)
>>> a = 0.004*(8*e**(-300*t) -8*e**(-1200*t))*(720000*e**(-300*t)-11520000*e**(-
1200*t)) +0.004*(9600*e**(-1200*t) - 2400*e**(-300*t))**2
>>> ta=terms_gcd(a)
>>> ta
(46080.0 - 576000.0*exp(-900*t) + 737280.0*exp(-1800*t))*exp(-600*t)
>>> _.subs(e**(-600*t),u)
u*(-576000.0*u**(3/2) + 737280.0*u**3 + 46080.0)
>>> solve(_,u)
[0, 0.201541162404989, 0.781429110349124]
>>> [ts(w-e**(-600*t),t) for w in _]
[[zoo], [0.00266960273088699], [0.000411051404934985]]
>>> ans = _[1:]
>>> [a.subs(t,ai[0]).n(2) for ai in ans]
[-2.8e-11, 8.0e-11]

_tsolve has to be used or else several hundred results are going to come back:

>>> solve(e**(-t)-x,t)
[log(1/x)]
>>> solve(e**(-2*t)-x,t)
[log(-sqrt(1/x)), log(sqrt(1/x))]
>>> solve(e**(-3*t)-x,t)
[log(-(1/x)**(1/3)/2 - sqrt(3)*I*(1/x)**(1/3)/2), log(-(1/x)**(1/3)/2 + sqrt(3)*
I*(1/x)**(1/3)/2), log((1/x)**(1/3))]

Should a hint like 'principle' or 'real' be used here to just get a simple inversion of the e**(-600*t) - u -> t = log(u)/-600 ?

**Labels:** Solvers  

Original comment: http://code.google.com/p/sympy/issues/detail?id=472#c1 Original author: https://code.google.com/u/117933771799683895267/

asmeurer commented 12 years ago
**Status:** Valid  

Original comment: http://code.google.com/p/sympy/issues/detail?id=472#c2 Original author: https://code.google.com/u/asmeurer@gmail.com/

oscarbenjamin commented 3 years ago

Still valid:

In [58]: a = (-2400*exp(-300*t) + 9600*exp(-1200*t))**2/250 + (4*exp(-300*t)/125 - 4*exp(-1200*t)/125)*(720000*exp(-300*t) - 11520000*exp(-1200*t))

In [59]: a
Out[59]: 
                                2                                                                
⎛        -300⋅t         -1200⋅t⎞    ⎛   -300⋅t      -1200⋅t⎞                                     
⎝- 2400⋅ℯ       + 9600⋅ℯ       ⎠    ⎜4⋅ℯ         4⋅ℯ       ⎟ ⎛        -300⋅t             -1200⋅t⎞
───────────────────────────────── + ⎜───────── - ──────────⎟⋅⎝720000⋅ℯ       - 11520000⋅ℯ       ⎠
               250                  ⎝   125         125    ⎠                                     

In [60]: solve(a, t)    # <----- hangs
^C---------------------------------------------------------------------------
KeyboardInterrupt

In [61]: solve(a.subs(t, z/300), z)
Out[61]: 
⎡   ⎛3 ___ 3 ____________⎞             ⎛3 ___ 3 ____________⎞             ⎛3 ___ 3 ____________⎞             ⎛3 ___ 3 __________
⎢   ⎜╲╱ 2 ⋅╲╱ 25 - 3⋅√41 ⎟   2⋅ⅈ⋅π     ⎜╲╱ 2 ⋅╲╱ 25 - 3⋅√41 ⎟   2⋅ⅈ⋅π     ⎜╲╱ 2 ⋅╲╱ 3⋅√41 + 25 ⎟   2⋅ⅈ⋅π     ⎜╲╱ 2 ⋅╲╱ 3⋅√41 + 2
⎢log⎜────────────────────⎟ - ─────, log⎜────────────────────⎟ + ─────, log⎜────────────────────⎟ - ─────, log⎜──────────────────
⎣   ⎝         2          ⎠     3       ⎝         2          ⎠     3       ⎝         2          ⎠     3       ⎝         2        

__⎞             ⎛    ____________⎞     ⎛    ____________⎞⎤
5 ⎟   2⋅ⅈ⋅π     ⎜   ╱ 25   3⋅√41 ⎟     ⎜   ╱ 3⋅√41   25 ⎟⎥
──⎟ + ─────, log⎜3 ╱  ── - ───── ⎟, log⎜3 ╱  ───── + ── ⎟⎥
  ⎠     3       ⎝╲╱   4      4   ⎠     ⎝╲╱     4     4  ⎠⎦