Re: any Python equivalent of Math::Polynomial::Solve?
|
| Did you perhaps use a list (type(p) == type([])) for p?
|
Alex
Using the coefficients in an array instead of a list
was the key in the solution to my problems
Your other suggestions regarding floating p
and the off-by-one error that I had with the
polynomial degree were also included
The results agree with solutions from PyGSL
as suggested by Pierre Schnizer, but seem
to run just a bit slower on my machine
Thanks again for your assistance
The version that I tested with follows
Stanley C. Kitching
# -
#!/usr/bin/env python
'''
NewsGroup comp.lang.python
Date . 2005-02-27
Subject .. any Python equivalent of Math::Polynomial::Solver
Reply_By . Alex Renelt
Edited_By Stanley c. Kitching
I'm writing a class for polynomial manipulation.
The generalization of the above code
for providing eigenvalues of a polynomial
is
definitions
1.) p = array( [ a_0 , a_i , , a_n ] )
P( x ) = \sum _{ i = 0 } ^n a_i x^i
2.) deg( p ) is its degree
3.) monic( p ) makes P monic
monic( p ) = p / p[ -1 ]
'''
import sys
import time
from numarray import *
import numarray.linear_algebra as LA
print '\n ' , sys.argv[ 0 ] , '\n'
def report( n , this_data ) :
print 'Coefficients ..' , list_coeff[ n ]
print
print 'Roots .'
print
dt , roots = this_data
for this_item in roots :
print '%s' % this_item
print
print 'Elapsed Time %.6f Seconds' % dt
print
print
def roots( p ) :
p = p / float( p[ -1 ] ) # monic( p )
n = len( p ) - 1 # degree of polynomial
z = zeros( ( n , n ) )
M = asarray( z , typecode = 'f8' )# typecode = c16, complex
M[ : -1 , 1 : ] = identity( n - 1 )
M[ -1 , : ] = -p[ : -1 ]
return LA.eigenvalues( M )
list_coeff = [ array( ( 2. , 3. , 1. ) ) ,
array( ( 1. , 3. , 5. , 7. , 9. ) ) ,
array( ( 10. , 8. , 6. , 4. , 2. , 1. , 2. , 4. , 6. ) ) ]
list_roots = [ ]
for these_coeff in list_coeff :
beg = time.time()
rs = roots( these_coeff )
end = time.time()
dt = end - beg
list_roots.append( [ dt , list( rs ) ] )
i = 0
for this_data in list_roots :
report( i , this_data )
i += 1
print
--
http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Cousin Stanley wrote: Alex Thanks for posting your generalized numarray eigenvalue solution It's been almost 30 years since I've looked at any characteristic equation, eigenvalue, eignevector type of processing and at this point I don't recall many of the particulars Not being sure about the nature of the monic( p ) function, I implemented it as an element-wise division of each of the coefficients This is correct. The aim is that p has a leading coefficient that equals 1. Is this anywhere near the correct interpretation for monic( p ) ? Using the version below, Python complained about the line . M[ -1 , : ] = -p[ : -1 ] Works for me. Did you perhaps use a list (type(p) == type([])) for p? Then python does not know what -p means (numeric or numarray does). So, in view of you comments about slicing in you follow-up, I tried without the slicing on the right .. M[ -1 , : ] = -p[ -1 ] That's wrong because you don't set a slice but a single item! Old code should work. We need all coefficients but the leading coeff. So take the slice p[:-1]. The following code will run and produce results, but I'm wondering if I've totally screwed it up since the results I obtain are different from those obtained from the specific 5th order Numeric solution previously posted here . from numarray import * . . import numarray.linear_algebra as LA . . def monic( this_list ) : . . m = [ ] . . last_item = this_list[ -1 ] . . for this_item in this_list : . . m.append( this_item / last_item ) . . return m your function equals the following one: def monic(p): return p / p[-1] But you have to ensure that p is an array object. It does element-wise operations per default. Remember we need future division or take float(p[-1]) or the denominator. . . . def roots( p ) : . . p = monic( p ) . . n = len( p ) # degree of polynomial The degree is len(p) -1 or something smaller (some people are calling len(p) -1 a degree bound instead). It is smaller if p contains leading zeros which should be deleted, e.g. P(x) = x^2 + 4 x + 4 could be entered as p = array([4, 4, 1, 0, 0]) which would produce a deg of 4 instead of 2. . . z = zeros( ( n , n ) ) . . M = asarray( z , typecode = 'f8' ) # typecode = c16, complex . . M[ : -1 , 1 : ] = identity( n - 1 ) . . M[ -1 , : ] = -p[ -1 ]# removed : slicing on the right . . return LA.eigenvalues( M ) . . . coeff = [ 1. , 3. , 5. , 7. , 9. ] remember to use an array or convert it on-the-fly inside your roots function: M[-1,:] = - asarray(p)[:-1] . . print 'Coefficients ..' . print . print '%s' % coeff . print . print 'Eigen Values .. ' . print . . eigen_values = roots( coeff ) . . for this_value in eigen_values : . . print '%s' % this_value . Any clues would be greatly appreciated Hope that helps. Alex -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Raymond L. Buvel wrote:
Alex Renelt wrote:
Alex Renelt wrote:
in addition:
I'm writing a class for polynomial manipulation. The generalization
of the above code is:
definitions:
1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
P(x) = \sum _{i=0} ^n a_i x^i
2.) deg(p) is its degree
3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]
then you get:
from numarray import *
import numarray.linear_algebra as la
def roots(p):
p = monic(p); n = deg(p)
M = asarray(zeros((n,n)), typecode = 'f8')
# or 'c16' if you need complex coefficients
M[:-1,1:] = identity(n-1)
M[-1,:] = -p[:-1]
return la.eigenvalues(M)
Alex
uhh, I made a mistake:
under definitions, 3.)
its monic(p) = p / p[-1] of course
Alex
Alex,
If you want a class for polynomial manipulation, you should check out my
ratfun module.
http://calcrpnpy.sourceforge.net/ratfun.html
Ray
Ray,
thanks a lot for your hint but I'm writing it for a students paper in a
german math class so I believe I should better do some work alone ;-)
In addition I only need a class for polynomials and not for rational
functions and I'm testing different iterative polynomial solvers. So I'm
happy to have my own small class which I understand 100%.
Generally I'm against rediscovering the wheel again and again!
Alex
--
http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
| In case you are still interested pygsl wraps the GSL solver. | | from pygsl import poly | | pc = poly.poly_complex( 3 ) | | tmp, rs = pc.solve( ( 2 , 3 , 1 ) ) | | print rs | | | You get pygsl at http://sourceforge.net/projects/pygsl/ Pierre I am still interested and have downloaded the PyGSL source from SourceForge and will attempt to build under Debian Linux Thanks for the information -- Stanley C. Kitching Human Being Phoenix, Arizona -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Carl Banks wrote: If you don't have a great need for speed, you can accomplish this easily with the linear algebra module of Numeric/numarray. Suppose your quintic polynomial's in the form a + b*x + c*x**2 + d*x**3 + e*x**4 + x**5 The roots of it are equal to the eigenvalues of the companion matrix: 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 -a -b -c -d -e The method zeros() in Scientific.Functions.Polynomial uses exactly that trick for finding the zeros of a general polynomial. If you need to do more with polynomials than just finding the zeros, the Polynomial class is probably better than an on-the-spot solution. Root finding through eigenvalues is not the fastest method, but it's simple and stable, and not terribly bad either. Sorry for not making that comment earlier, I don't have the time to follow this list at the moment (to my great regret), but I was made aware of this thread through PythonURL. Konrad. -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Alex Thanks for posting your generalized numarray eigenvalue solution It's been almost 30 years since I've looked at any characteristic equation, eigenvalue, eignevector type of processing and at this point I don't recall many of the particulars Not being sure about the nature of the monic( p ) function, I implemented it as an element-wise division of each of the coefficients Is this anywhere near the correct interpretation for monic( p ) ? Using the version below, Python complained about the line . M[ -1 , : ] = -p[ : -1 ] So, in view of you comments about slicing in you follow-up, I tried without the slicing on the right .. M[ -1 , : ] = -p[ -1 ] The following code will run and produce results, but I'm wondering if I've totally screwed it up since the results I obtain are different from those obtained from the specific 5th order Numeric solution previously posted here . from numarray import * . . import numarray.linear_algebra as LA . . def monic( this_list ) : . . m = [ ] . . last_item = this_list[ -1 ] . . for this_item in this_list : . . m.append( this_item / last_item ) . . return m . . . def roots( p ) : . . p = monic( p ) . . n = len( p ) # degree of polynomial . . z = zeros( ( n , n ) ) . . M = asarray( z , typecode = 'f8' ) # typecode = c16, complex . . M[ : -1 , 1 : ] = identity( n - 1 ) . . M[ -1 , : ] = -p[ -1 ]# removed : slicing on the right . . return LA.eigenvalues( M ) . . . coeff = [ 1. , 3. , 5. , 7. , 9. ] . . print 'Coefficients ..' . print . print '%s' % coeff . print . print 'Eigen Values .. ' . print . . eigen_values = roots( coeff ) . . for this_value in eigen_values : . . print '%s' % this_value . Any clues would be greatly appreciated -- Stanley C. Kitching Human Being Phoenix, Arizona -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
In case you are still interested pygsl wraps the GSL solver. snip from pygsl import poly pc = poly.poly_complex(3) tmp, rs = pc.solve((2,3,1)) print rs /snip You get pygsl at http://sourceforge.net/projects/pygsl/ Pierre -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Just wrote:
In article [EMAIL PROTECTED],
Carl Banks [EMAIL PROTECTED] wrote:
It should be pretty easy to set up a Numeric matrix and call
LinearAlgebra.eigenvalues. For example, here is a simple quintic
solver:
. from Numeric import *
. from LinearAlgebra import *
.
. def quinticroots(p):
. cm = zeros((5,5),Float32)
. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
. cm[4,0] = -p[0]
. cm[4,1] = -p[1]
. cm[4,2] = -p[2]
. cm[4,3] = -p[3]
. cm[4,4] = -p[4]
. return eigenvalues(cm)
now-you-can-find-all-five-Lagrange-points-ly yr's,
Wow, THANKS. This was the answer I was secretly hoping for... Great
need for speed, no, not really, but this Numeric-based version is about
9 times faster than what I translated from Perl code yesterday, so from
where I'm standing your version is blazingly fast...
Thanks again,
Just
in addition:
I'm writing a class for polynomial manipulation. The generalization of
the above code is:
definitions:
1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
P(x) = \sum _{i=0} ^n a_i x^i
2.) deg(p) is its degree
3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]
then you get:
from numarray import *
import numarray.linear_algebra as la
def roots(p):
p = monic(p); n = deg(p)
M = asarray(zeros((n,n)), typecode = 'f8')
# or 'c16' if you need complex coefficients
M[:-1,1:] = identity(n-1)
M[-1,:] = -p[:-1]
return la.eigenvalues(M)
Alex
--
http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Alex Renelt wrote:
in addition:
I'm writing a class for polynomial manipulation. The generalization of
the above code is:
definitions:
1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
P(x) = \sum _{i=0} ^n a_i x^i
2.) deg(p) is its degree
3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]
then you get:
from numarray import *
import numarray.linear_algebra as la
def roots(p):
p = monic(p); n = deg(p)
M = asarray(zeros((n,n)), typecode = 'f8')
# or 'c16' if you need complex coefficients
M[:-1,1:] = identity(n-1)
M[-1,:] = -p[:-1]
return la.eigenvalues(M)
Alex
uhh, I made a mistake:
under definitions, 3.)
its monic(p) = p / p[-1] of course
Alex
--
http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Alex Renelt wrote:
Alex Renelt wrote:
in addition:
I'm writing a class for polynomial manipulation. The generalization of
the above code is:
definitions:
1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
P(x) = \sum _{i=0} ^n a_i x^i
2.) deg(p) is its degree
3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]
then you get:
from numarray import *
import numarray.linear_algebra as la
def roots(p):
p = monic(p); n = deg(p)
M = asarray(zeros((n,n)), typecode = 'f8')
# or 'c16' if you need complex coefficients
M[:-1,1:] = identity(n-1)
M[-1,:] = -p[:-1]
return la.eigenvalues(M)
Alex
uhh, I made a mistake:
under definitions, 3.)
its monic(p) = p / p[-1] of course
Alex
Alex,
If you want a class for polynomial manipulation, you should check out my
ratfun module.
http://calcrpnpy.sourceforge.net/ratfun.html
Ray
--
http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
On 2005-02-26, Just [EMAIL PROTECTED] wrote: While googling for a non-linear equation solver, I found Math::Polynomial::Solve in CPAN. It seems a great little module, except it's not Python... I'm especially looking for its poly_root() functionality (which solves arbitrary polynomials). Does anyone know of a Python module/package that implements that? Just Although I dont' really work on it any more, the Py-ML module which interfaces Python to Mathematica would almost certain be able to do this, although you'd need an installation of Mathematica. http://sourceforge.net/projects/py-ml/ -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Carl Banks wrote: . from Numeric import * . from LinearAlgebra import * . . def quinticroots(p): . cm = zeros((5,5),Float32) . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0 . cm[4,0] = -p[0] . cm[4,1] = -p[1] . cm[4,2] = -p[2] . cm[4,3] = -p[3] . cm[4,4] = -p[4] . return eigenvalues(cm) Here's an improved version. It uses 64-bit numbers (I had used type Float32 because I often use a float32 type at work, not in Python, unfortunately), and array assignment. . def quinticroots(p): . cm = zeros((5,5),Float) . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0 . cm[4,:] = -array(p) . return eigenvalues(cm) -- CARL BANKS -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Just wrote: While googling for a non-linear equation solver, I found Math::Polynomial::Solve in CPAN. It seems a great little module, except it's not Python... I'm especially looking for its poly_root() functionality (which solves arbitrary polynomials). Does anyone know of a Python module/package that implements that? Just Does SciPy do what you want? Specifically Scientific.Functions.FindRoot [1] Scientific.Functions.Polynomial [2] Regards, Nick. [1] http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Scientific_9.html [2] http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Scientific_13.html -- Nick Coghlan | [EMAIL PROTECTED] | Brisbane, Australia --- http://boredomandlaziness.skystorm.net -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Just [EMAIL PROTECTED] wrote in message news:[EMAIL PROTECTED] Does SciPy do what you want? Specifically Scientific.Functions.FindRoot [1] Scientific.Functions.Polynomial [2] http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci entific_9.html [2] http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci entific_13.html (Hm, I had the impression that scipy != Konrad Hinsen's Scientific module.) www.scipy.org (first hit on Python SciPy google search) Terry J. Reedy -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
In article [EMAIL PROTECTED], Just [EMAIL PROTECTED] wrote: While googling for a non-linear equation solver, I found Math::Polynomial::Solve in CPAN. It seems a great little module, except Thank you. it's not Python... Sorry about that. I'm especially looking for its poly_root() functionality (which solves arbitrary polynomials). Does anyone know of a Python module/package that implements that? Are you looking for that particular algorithm, or for any source that will find the roots of the polynomial? The original source for the algorithm used in the module is from Hiroshi Murakami's Fortran source, and it shouldn't be too difficult to repeat the translation process to python. -- -john February 28 1997: Last day libraries could order catalogue cards from the Library of Congress. -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
In article [EMAIL PROTECTED], [EMAIL PROTECTED] (John M. Gamble) wrote: In article [EMAIL PROTECTED], Just [EMAIL PROTECTED] wrote: While googling for a non-linear equation solver, I found Math::Polynomial::Solve in CPAN. It seems a great little module, except Thank you. it's not Python... Sorry about that. Heh, how big are the odds you find the author of an arbitrary Perl module on c.l.py... I'm especially looking for its poly_root() functionality (which solves arbitrary polynomials). Does anyone know of a Python module/package that implements that? Are you looking for that particular algorithm, or for any source that will find the roots of the polynomial? Any will do. As I wrote in another post, I'm currently only looking for a quintic equation solver, which your module does very nicely. The original source for the algorithm used in the module is from Hiroshi Murakami's Fortran source, and it shouldn't be too difficult to repeat the translation process to python. Ah ok, I'll try to locate that (following the instruction in Solve.pm didn't work for me :( ). Just -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
In article [EMAIL PROTECTED], Just [EMAIL PROTECTED] wrote: Heh, how big are the odds you find the author of an arbitrary Perl module on c.l.py... Hey, that's why it's called lurking. Any will do. As I wrote in another post, I'm currently only looking for a quintic equation solver, which your module does very nicely. The original source for the algorithm used in the module is from Hiroshi Murakami's Fortran source, and it shouldn't be too difficult to repeat the translation process to python. Ah ok, I'll try to locate that (following the instruction in Solve.pm didn't work for me :( ). Ouch. I just did a quick search and found that that site has undergone a few changes, and the code that i reference is missing. A few other links in the docs are stale too. I need to update the documentation. Anyway, doing a search for 'hqr' and Eispack got me a lot of sites. In particular, this one is pretty friendly: http://netlib.enseeiht.fr/eispack/ Look at the source for balanc.f (does the prep-work) and hqr.f (does the solving). Minor annoyance: the real and imaginary parts of the roots are in separate arrays. I combined them into complex types in my perl source, in case you want to make a comparison. Of course, all this may be moot if the other suggestions work out. -- -john February 28 1997: Last day libraries could order catalogue cards from the Library of Congress. -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
In article [EMAIL PROTECTED], [EMAIL PROTECTED] (John M. Gamble) wrote: The original source for the algorithm used in the module is from Hiroshi Murakami's Fortran source, and it shouldn't be too difficult to repeat the translation process to python. Ah ok, I'll try to locate that (following the instruction in Solve.pm didn't work for me :( ). Ouch. I just did a quick search and found that that site has undergone a few changes, and the code that i reference is missing. A few other links in the docs are stale too. I need to update the documentation. Anyway, doing a search for 'hqr' and Eispack got me a lot of sites. In particular, this one is pretty friendly: http://netlib.enseeiht.fr/eispack/ Look at the source for balanc.f (does the prep-work) and hqr.f (does the solving). Minor annoyance: the real and imaginary parts of the roots are in separate arrays. I combined them into complex types in my perl source, in case you want to make a comparison. Thanks! I'll check that out. Of course, all this may be moot if the other suggestions work out. SciPy indeed appear to contain a solver, but I'm currently stuck in trying to _get_ it for my platform (OSX). I'm definitely not going to install a Fortran compiler just to evaluate it (even though my name is not Ilias ;-). Also, SciPy is _huge_, so maybe a Python translation of that Fortran code or your Perl code will turn out to be more attractive after all... Just -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Just wrote: (Hm, I had the impression that scipy != Konrad Hinsen's Scientific module.) You're probably right :) I had played with [1], but it only calculates one root, and I need all roots (specifically, for a quintic equation). [2] doesn't seem to be a solver? Actually, I was curious whether the 'zeros' method in [2] did the right thing. Cheers, Nick. -- Nick Coghlan | [EMAIL PROTECTED] | Brisbane, Australia --- http://boredomandlaziness.skystorm.net -- http://mail.python.org/mailman/listinfo/python-list
Re: any Python equivalent of Math::Polynomial::Solve?
Just wrote: While googling for a non-linear equation solver, I found Math::Polynomial::Solve in CPAN. It seems a great little module, except it's not Python... I'm especially looking for its poly_root() functionality (which solves arbitrary polynomials). Does anyone know of a Python module/package that implements that? If you don't have a great need for speed, you can accomplish this easily with the linear algebra module of Numeric/numarray. Suppose your quintic polynomial's in the form a + b*x + c*x**2 + d*x**3 + e*x**4 + x**5 The roots of it are equal to the eigenvalues of the companion matrix: 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 -a -b -c -d -e It should be pretty easy to set up a Numeric matrix and call LinearAlgebra.eigenvalues. For example, here is a simple quintic solver: . from Numeric import * . from LinearAlgebra import * . . def quinticroots(p): . cm = zeros((5,5),Float32) . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0 . cm[4,0] = -p[0] . cm[4,1] = -p[1] . cm[4,2] = -p[2] . cm[4,3] = -p[3] . cm[4,4] = -p[4] . return eigenvalues(cm) now-you-can-find-all-five-Lagrange-points-ly yr's, -- CARL BANKS -- http://mail.python.org/mailman/listinfo/python-list
