[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] Polynomial root finding
From: |
Will Leckie |
Subject: |
[Help-gsl] Polynomial root finding |
Date: |
Tue, 10 Jan 2006 11:02:16 -0500 |
Hi,
Great work on the GSL, I love using it in my research.
I'm using the GSL in my research to solve arbitrary polynomials of order up to
4. I've written a little wrapper around the GSL to use
gsl_poly_complex_solve_quadratic() for quadratics and
gsl_poly_complex_solve_cubic() for cubics, and I use gsl_poly_complex_solve()
to solve quartics.
I also implemented my own analytical quartic solver based on Ferrari's method
(http://mathworld.wolfram.com/QuarticEquation.html) using the GSL complex math
operations such as gsl_complex_mul etc. to manipulate the complex numbers
involved in Ferrari's solution.
When comparing the speed of the GSL gsl_poly_complex_solve() solution of a
quartic to my analytical version, I find that the analytic routine runs just a
touch faster:
*********** timing GSL version...
GSL version took 1.078 seconds for 100000 runs; averaging 1.078e-05 seconds per
run.
*********** Analytic version...
Analytic version took 0.797 seconds for 100000 runs; averaging 7.97e-06 seconds
per run.
(I'm running a 3.2 GHz Pentium IV with 1GB RAM, Cygwin on Windows XP for the
timing tests, using the cstdlib rand() function to generate "random"
coefficients between 0.0 and 10.0 for the quartics for each run)
When comparing the "precision" of the two quartic solvers, I find a very small
(machine-level?) difference between the roots that they find, typically 10^-30,
using gsl_complex_sub() to calculate the difference between roots.
I also find that sometimes the GSL gsl_poly_complex_solve() routine returns a
non-zero status int meaning that it didn't find roots or it found
unstable/overflow/underflow conditions, I presume.
So, my question is, what is my advantage in using the GSL
gsl_poly_complex_solve() routine to solve my quartics? I have a feeling the
GSL gsl_poly_complex_solve() routine is more robust to large and small
coefficients in the polynomial, because I haven't taken too much care in
looking out for those in my analytic code. I'm sure the GSL developers looked
at Ferrari's method when considering whether or not to implement a
gsl_poly_complex_solve_quartic() routine, since an analytic solution is
acheivable, and I'm wondering why they didn't implement one? Is there some
obvious flaw when implementing Ferrari's method computationally? What are the
other advantages of solving quartics using the iterative method in
gsl_poly_complex_solve() instead?
If you'd like to see any of the code I've used, either for Ferrari's method, or
for testing, I'd be happy to share it.
Thanks again for the great library.
Cheers,
Will
- [Help-gsl] Polynomial root finding,
Will Leckie <=