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?
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
