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
17  3437 
In article <ma***************************************@python. org>,
Nick Coghlan <nc******@iinet.net.au> wrote: 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/S...thonManual/Sci
entific_9.html
[2]
http://starship.python.net/~hinsen/S...thonManual/Sci
entific_13.html
(Hm, I had the impression that scipy != Konrad Hinsen's Scientific
module.)
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?
Just
In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> 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.
In article <cv**********@e250.ripco.com>, jg*****@ripco.com (John M. Gamble) wrote: In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> 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
In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> 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.
In article <cv**********@e250.ripco.com>, jg*****@ripco.com (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
Just wrote:
<snip>
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
The GNU Scientific Library has a nice root finder for polynomials with
real coefficients. I have wrapped this with Pyrex to work with my
ratfun module see: http://calcrpnpy.sourceforge.net/ratfun.html
If this will suit your needs, I can send you an alpha release of the
package with the root finder. It is not pure Python. I requires Pyrex
and a C compiler to install. My guess is that it will work on OSX as
well as it does on Linux. This functionality will be included in the
next release of the ratfun package but I still have to unit test a
number of components and update the documentation. Consequently, an
official release will not happen soon.
Ray
In article <I6***************@twister.rdc-kc.rr.com>,
"Raymond L. Buvel" <le******@wi.rr.com> wrote: Just wrote:
<snip>
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
The GNU Scientific Library has a nice root finder for polynomials with
real coefficients. I have wrapped this with Pyrex to work with my
ratfun module see:
http://calcrpnpy.sourceforge.net/ratfun.html
If this will suit your needs, I can send you an alpha release of the
package with the root finder. It is not pure Python. I requires Pyrex
and a C compiler to install. My guess is that it will work on OSX as
well as it does on Linux. This functionality will be included in the
next release of the ratfun package but I still have to unit test a
number of components and update the documentation. Consequently, an
official release will not happen soon.
Thank you, I'll check this out. I had come across GSL, but not Python
bindings. (GPL is probably a problem for my project, but it's very good
to know anyway.)
On the other hand, I just finished translating the relevant portions of
Math::Polynomial::Solve to Python, so I'm probably all set, at least for
now. Thanks everyone for the responses, especially to John Gamble!
Just
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
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
| 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
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.
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
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
| 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 ....
Thanks again for the link to the Python GNU Scientific Library ....
I managed to build it under Debian GNU/Linux
and the solution you posted above works as is ....
The PyGSL results agree with the numarray eigenvalue
solution posted by Alex Renelt and a few timing tests
for small problems seem to indicate that it runs a bit quicker ....
I'm looking forward to further use of this package
in the future ....
--
Stanley C. Kitching
Human Being
Phoenix, Arizona
| ....
| 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
This thread has been closed and replies have been disabled. Please start a new discussion. Similar topics
by: mike420 |
last post by:
I think everyone who used Python will agree that its syntax is
the best thing going for it. It is very readable and easy
for everyone to learn. But, Python does not a have very good
macro...
|
by: Glen Wheeler |
last post by:
Hello all,
My program uses many millions of integers, and Python is allocating
way too much memory for these. I can't have the performance hit by
using disk, so I figured I'd write a C...
|
by: Kyler Laird |
last post by:
I'm trying to do some georeferencing - using points of known
location (ground control points, GCPs) on an image to develop
a polynomial that can be used to approximate the locations of
other...
|
by: Joseph Aldred |
last post by:
I am new to C and I am writing a program for a class I am taking. I need to
use a exponent but I can not figure how to do it. I figure I need to use
exp() or pow() but I can not figure out how to...
|
by: sekitoleko |
last post by:
c program for bisection method to find roots of a polynomial
|
by: MathNewbie |
last post by:
Hi,
I'm trying to do get my head around some, probably basic, math
calculations. I checked some .net math libs, but the only thing I learn
from studying those is humility ;-).
Anyway, i have...
|
by: Dilip |
last post by:
Recently in our code, I ran into a situation where were stuffing a
float inside a double. The precision was extended automatically
because of that. To make a long story short, this caused...
|
by: Jordan |
last post by:
Hi everyone,
I'm a big Python fan who used to be involved semi regularly in
comp.lang.python (lots of lurking, occasional posting) but kind of
trailed off a bit. I just wrote a frustration...
|
by: curiously enough |
last post by:
I need a practical method to find the eigenvalues of a matrix in C++ because the one I know(the only one I know) is to subtract the elements of the diagonal by the eigenvalue and then find the...
|
by: isladogs |
last post by:
The next Access Europe meeting will be on Wednesday 2 August 2023 starting at 18:00 UK time (6PM UTC+1) and finishing at about 19:15 (7.15PM)
The start time is equivalent to 19:00 (7PM) in Central...
|
by: linyimin |
last post by:
Spring Startup Analyzer generates an interactive Spring application startup report that lets you understand what contributes to the application startup time and helps to optimize it. Support for...
|
by: kcodez |
last post by:
As a H5 game development enthusiast, I recently wrote a very interesting little game - Toy Claw ((http://claw.kjeek.com/))。Here I will summarize and share the development experience here, and hope it...
|
by: isladogs |
last post by:
The next Access Europe meeting will be on Wednesday 6 Sept 2023 starting at 18:00 UK time (6PM UTC+1) and finishing at about 19:15 (7.15PM)
The start time is equivalent to 19:00 (7PM) in Central...
|
by: Rina0 |
last post by:
I am looking for a Python code to find the longest common subsequence of two strings. I found this blog post that describes the length of longest common subsequence problem and provides a solution in...
|
by: DJRhino |
last post by:
Private Sub CboDrawingID_BeforeUpdate(Cancel As Integer)
If = 310029923 Or 310030138 Or 310030152 Or 310030346 Or 310030348 Or _
310030356 Or 310030359 Or 310030362 Or...
|
by: lllomh |
last post by:
How does React native implement an English player?
|
by: Mushico |
last post by:
How to calculate date of retirement from date of birth
|
by: DJRhino |
last post by:
Was curious if anyone else was having this same issue or not....
I was just Up/Down graded to windows 11 and now my access combo boxes are not acting right. With win 10 I could start typing...
| |
