Hi again,
I am getting confused now. See the script and the results below.
First Variant (using nsolve) finds some extra solutions, that probably
would be gone changed the tolerance, but using the option tol=1e-25
crashes, in fact, I could not go beyond tol=1e-21.
Second Variant (using Poly) kind of works, but there is a small fix
required. "mypoly.coeffs" ommits the zeros, therefore if the
coefficients for x^0, x^1 and x^3 are 0, they would not appear. Is
this a bug or a feature?
Third Variant (using solve) is completely messed now, solutions don't
get even close to the real solutions.
I plotted the function in a spreadsheet, and found that the solutions
provided by the Second Variant (with the manual fix) are the correct
ones.
Any ideas on how to fix all this to get the same results from all
methods? how can I believe any of the results otherwise?
Thanks!
regards, anartz
--------------------------------
-- The script --
--------------------------------
from sympy import *
from numpy import arange
a = 5.75
e = 8e-6
c = 0.5
x = 20.
o = 0.02
v = 0.005
b = 0.12
s = 0.23
w = 0.23
d = 0.
f = 8.
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)
print "Equation to solve is: "
pretty_print(myEqConstants.expand())
print
print 'first variant:'
print "Solutions using nsolve function in the range are:"
myRoot=[]
for aa in arange(-20.,20.,0.1):
NewRoot=nsolve(myEqConstants, r, aa, tol=1e-20)
if ((str(NewRoot) not in myRoot) and (NewRoot>0.0001 or
NewRoot<-0.0001)): # don't get loads of roots near 0, also excludes 0
myRoot.append(str(NewRoot))
for aa in myRoot:
print aa
print
print 'Second Variant:'
mypoly = Poly(myEqConstants.expand(), r)
print "mypoly.coeffs:" + str(mypoly.coeffs)
print "instead of: (-3.58196480000000e-17, -8.31015833600000e-17,
4.53671602565120e-15, 0, -4.23200000000000e-15,0,0)"
print
print "Using correct coefficients:"
for aa in mpmath.polyroots((-3.58196480000000e-17,
-8.31015833600000e-17, 4.53671602565120e-15, 0,
-4.23200000000000e-15,0,0)):
print aa
print
print 'Third variant:'
print "Solutions obained using solve function:"
NewSols=solve(myEqConstants, r)
TheSols=[N(e, n=4, chop=True) for e in NewSols]
for aa in TheSols:
pprint(aa)
--------------------------------
-- End Script --
--------------------------------
--------------------------------
-- Results --
--------------------------------
Equation to solve is:
2 4
5 6
- 4.232e-15⋅r + 4.5367160256512e-15⋅r - 8.310158336e-17⋅r -
3.5819648e-17⋅r
first variant:
Solutions using nsolve function in the range are:
-12.4399173136662
-10.304955620728
-10.151265495134
0.978377565046911
-0.960891724232247
8.71027033335483
10.1024314728516
Second Variant:
mypoly.coeffs:(-3.58196480000000e-17, -8.31015833600000e-17,
4.53671602565120e-15, -4.23200000000000e-15)
instead of: (-3.58196480000000e-17, -8.31015833600000e-17,
4.53671602565120e-15, 0, -4.23200000000000e-15,0,0)
Using correct coefficients:
-12.4399173136662
-0.960891724232247
-9.52592891648496e-19
5.88734784480316e-19
0.978377565046911
10.1024314728516
Third variant:
Solutions obained using solve function:
-7.834e-1
-12.62
9.870
0
1.211
--------------------------------
-- End Results --
--------------------------------
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
