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.html#polynomial-roots-polyroots
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
-~----------~----~----~----~------~----~------~--~---
