On May 15, 9:04 pm, Ondrej Certik <ond...@certik.cz> wrote: > <vinzent.steinb...@googlemail.com> wrote: > > Please note that you also could use the brand-new nsolve() to do this > > without calling mpmath directly and using a simpler syntax. (If besselj > > () was implemented in sympy.) > > > By the way, nsolve() should support scipy too. Do you think it should > > be implemented in nsolve() or in mpmath's findroot()? > > If there is a way to take advantage of symbolic stuff, then in sympy, > if it's just pure numerics, then mpmath is fine too.
Yeah, but in this case it would simplify the code, because you would not need all the lambdas. > > Btw, how is it possible that findroot returns complex numbers > sometimes? E.g. I had to do the following hacks: > > + if method == "sympy": > + # The findroot() is fragile, it sometimes returns complex > numbers, > + # so we chop all complex parts (that are small anyway). Also we > + # need to set the tolerance, as it sometimes fail without it. > + def f_real(x): > + return f(complex(x).real) > + root = findroot(f_real, x, solver="muller", tol=1e-9) > + root = complex(root).real Why are you using method='muller'? This method is known for converging towards complex roots even with real starting points! It's not "fragile", it's actually an expected behavior. See the docstring [1]. > > Which is not nice. Also I had to setup tol=1e-9, otherwise the > doctests run when executed directly, but fail when executed using > ./setup test (e.g. after running the regular testsuite), so it looks > like there is some global state in findroot. This is an general issue of mpmath, you need to set explicitly mp.dps = 15 at the beginning of every doctest, because if a doctest that changes dps is executed before, you might get a precision you did not expect. Currently mp is being made context-dependent rather than global. > > Compare this to scipy: > > + elif method == "scipy": > + root = newton(f, x) > > it just works. It can be a problem in evaluating the jn() function, > since in the sympy case, I am using the sin/cos representation, while > in the scipy case I am just calling sph_jn() from scipy. > > Ondrej root = findroot(f, x) should just work too. Vinzent [1] http://mpmath.googlecode.com/svn/trunk/doc/build/calculus/optimization.html --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~----------~----~----~----~------~----~------~--~---
