Re: any Python equivalent of Math::Polynomial::Solve?

[Alex Renelt]

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