I did use the nsolve before (that was the topic of my first post :)), but I need to find all the solutions in a range, and solutions can be very close to each other depending on the constants. I am trying to optimize a problem, so I would like to check the final outcome for all valid solutions. I thought about symbolically solving the equation and then .evalf() all solutions, to make sure I had them all, but ended up with the extrange imaginary parts. Speed would not be that much of a issue in this case... yet :)
At the end, I added a loop to the script to search solutions using the nsolve function and initial points that range from the minimum to maximum possible (phisical) solutions, with a fine stepsize. I append the solutions to a list as they come if they have not been already found with a different initial point. Regards, anartz On Aug 11, 10:58 pm, Vinzent Steinberg <[email protected]> wrote: > There are three variants to solve this using sympy, and you chose the > slowest one. ;) > > Why do you want to solve it symbolically? It's a polynomial of 6th > degree, so there are probably very complex symbolic solution, if any. > > The best way to get the solutions is using mpmath.polyroots, see [1]. > > Just run the following script: > > from sympy import * > a = 5.75 > e = 8e-6 > c = 0.25 > x = 20. > o = 0.02 > v = 0.005 > b = 0.12 > s = 0.23 > w = 0.23 > d = 0. > f = 20. > r = Symbol('r') > > myEqConstants = (-4*a**2*r**2*e**2*(0.002*c*x*(o + v) + 2*a*b*r*e + > 2*a*s*r*e > + 2*a*w*r*e + a*e*r**2)**2 -4*a**2*r**2*e**2*(0.002*c*x*(o + v) > - (2*a*b*r*e + 2*a*s*r*e + 2*a*w*r*e + a*e*r**2))**2 > + 16*a**4*d**2*r**4*e**4 + 16*a**4*r**4*f**2*e**4) > pretty_print(myEqConstants.expand()) > print > print 'first variant:' > # we can solve it numerically, using a starting point > print nsolve(myEqConstants, r, -30) > print nsolve(myEqConstants, r, 0) > print nsolve(myEqConstants, r, 30) > print > print 'second variant:' > # mpmath has a solver for polynomials, but we have to convert it to a > list of > # coefficients (please not that the results are not very accurate, you > can refine > # them using an iterative solver) > mypoly = Poly(myEqConstants.expand(), r) > print mpmath.polyroots(mypoly.coeffs) > print > print 'third variant:' > # solve it symbolically > pretty_print([e.evalf(50) for e in solve(myEqConstants, r)]) # this is > slow > > I get this output: > > 2 4 > 5 6 > - 1.058e-15⋅r + 2.86075194816512e-14⋅r - 8.310158336e-17⋅r - > 3.5819648e-17⋅r > > first variant: > -29.4436689512432 > 0.0 > 27.1235615040264 > > second variant: > [mpf('-29.461988293519347'), mpf('0.036987320809019009'), mpf > ('27.105000972710329')] > > third variant: > [27.123561504028409502282500262935396379311391281644 - > 2.74926057288475386433396561255 > 12261631e-13⋅ⅈ, 0.19236859708818558635612456726563255832985519147261 + > 2.7526569045186 > 83537962579904305246576705e-13⋅ⅈ, > -0.1922611498753081213912536102857789787529633972320 > 1 + 2.527748870415793485879011339093847593146e-13⋅ⅈ, > -29.44366895124128636328604582383 > 0091968432655634478 - 2.5311452020497231595076256308478680067e-13⋅ⅈ, > 0] > > Indeed it's strange that the imaginary parts don't vanish, even if you > use higher precision for evaluating. This smells like a bug (assuming > the roots are really real). Fredrik, what do you think? > Please note that the second variant is somewhat inaccurate, you can > however use it as starting points for the first variant (see > documentation). > > Vinzent > > [1]http://mpmath.googlecode.com/svn/trunk/doc/build/calculus/polynomials... > > On 11 Aug., 14:14, New2Sympy <[email protected]> wrote: > > > I think I am in agreement with Vinzent, this would probably solve the > > issue I am having now... > > > I am trying to obtain the solutions using the linex below. When I > > solve for the constants assigned to the variables, I get complex > > solutions with very small imaginary components, that should be real > > solutions. I assume this happens due to some precision settings. And > > when I try to solve it symbolically, it hangs.... > > > The equation to solve is the numerator of function I got using the > > _.as_numer_denom() > > > >>> from sympy import var,solve > > >>> #where 'a', 'e', 'c', 'x', 'o', 'v', 'b', 's', 'w', 'd', 'f' are > > >>> constants, and want to solve it for 'r' > > >>> a=5.75 > > >>> e=8e-6 > > >>> c=0.25 > > >>> x=20. > > >>> o=0.02 > > >>> v=0.005 > > >>> b=0.12 > > >>> s=0.23 > > >>> w=0.23 > > >>> d=0. > > >>> f=20. > > > >>> var('r') > > r > > > >>> myEqConstants=-4*a**2*r**2*e**2*(0.002*c*x*(o + v) + 2*a*b*r*e + > > >>> 2*a*s*r*e + 2*a*w*r*e + a*e*r**2)**2 -4*a**2*r**2*e**2*(0.002*c*x*(o + > > >>> v) - (2*a*b*r*e + 2*a*s*r*e + 2*a*w*r*e + a*e*r**2))**2 > > >>> +16*a**4*d**2*r**4*e**4 + 16*a**4*r**4*f**2*e**4 > > > >>> mySols=solve(myEqConstants,r) > > >>> print [aa.evalf() for aa in mySols] > > > [27.1235615040284 - 2.75e-13*I, -29.4436689512413 - 2.53e-13*I, > > 0.192368597088186 + 2.7527e-13*I, 0, -0.192261149875308 + > > 2.5278e-13*I] > > > >>> var('a r e c x o v b s w d f') > > > (a, r, e, c, x, o, v, b, s, w, d, f) > > > >>> myEqSymbolic=-4*a**2*r**2*e**2*(0.002*c*x*(o + v) + 2*a*b*r*e + > > >>> 2*a*s*r*e + 2*a*w*r*e + a*e*r**2)**2-4*a**2*r**2*e**2*(0.002*c*x*(o + > > >>> v) - (2*a*b*r*e + 2*a*s*r*e + 2*a*w*r*e + > > >>> a*e*r**2))**2+16*a**4*d**2*r**4*e**4 + 16*a**4*r**4*f**2*e**4 > > > >>> mySols=solve(myEqSymbolic,r) > > > Thanks for your help! > > > regards, anartz > > > On Aug 6, 6:26 pm, Vinzent Steinberg > > > <[email protected]> wrote: > > > On 6 Aug., 19:57, smichr <[email protected]> wrote: > > > > > On Aug 5, 9:29 pm, Vinzent Steinberg > > > > > <[email protected]> wrote: > > > > > On Aug 2, 12:14 am, smichr <[email protected]> wrote: > > > > > > > >>> num,den = (1/(.001+a)**3-6/(.9-a)**3).as_numer_denom() > > > > > > >>>nsolve(num,a,.3) # no need for sympy now > > > > > > Maybe solve()/nsolve() should do this for you. What do you think? > > > > > > Vinzent > > > > > I was thinking the same thing, but you would have to watch out that > > > > you didn't introduce spurious roots: e.g. > > > > > >>> (x/(x-1)).diff(x) > > > > > -1/(1 - x) - x/(1 - x)**2>>> n,d=_.as_numer_denom();n/d > > > > > (-x*(1 - x) - (1 - x)**2)/(1 - x)**3 > > > > > Unless we cancel out the common factors (something that isn't working > > > > in general yet but will hopefully work with the new polys module) we > > > > will get a root of x=1. > > > > It should be easy to check for false roots numerically, but yeah, it > > > has to be implemented. > > > > > I was also wondering if the routine shouldn't > > > > try to solve the expression symbolically to see if it's possible > > > > before trying to solve it numerically...or should sympy just trust the > > > > user and not waste the time checking symbolic solutions if a numeric > > > > solution has been requested? > > > > It would do it this way: The user wants to solve() an equation, if we > > > can't do it symbolically, we return an instance of RootOf(), which can > > > be evalf()'ed to get a numerical approximation.nsolve() should solve > > > it always numerically, possibly doing some symbolic optimizations, > > > which can be disabled. > > > > BTW: It's now fixed to find the root of 1/(0.001+a)**3-6/(0.9-a)**3. > > > > Vinzent > > > > Vinzent > > --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected] To unsubscribe from this group, send email to [email protected] For more options, visit this group at http://groups.google.com/group/sympy?hl=en -~----------~----~----~----~------~----~------~--~---
