I briefly took a look at this. The GSL cubic routine is failing because 'R' and 'Q' are not precisely zero. If the CR2==CQ3 test succeeded, the routine might also give the correct answer, except that Q is actually negative in this context, so that sqrt(Q) is not finite.
The real solution is to either use gsl_poly_complex_solve_cubic, or to tailor the gsl_poly_solve_cubic routine to your specific physical problem. I do not see a way to improve the GSL code at present, as most modifications will solve this problem by creating others. I also examined the accuracy of gsl_poly_complex_solve_cubic, and compared it to my implementation of CERNLIB's cubic solver (as reimplemented in C at http://o2scl.sourceforge.net) and found that GSL actually performed slightly better for your particular example. Take care, Andrew On Dec 20, 2007 8:37 AM, nick kama <[EMAIL PROTECTED]> wrote: > Hello, > Thank you for your answer, i think this bad behaviour would need some more > attention. > This is an unphysical situation, it rarely happens to have a diagonal > matrix with two or more equal double, but it could happen and i think this > is an unwanted behaviour and should be fixed or documented in some-way. > > Niccolo' > > > > On Dec 20, 2007 3:15 AM, Ralph Silva <[EMAIL PROTECTED]> wrote: > > > Hi Niccolo, > > > > it looks like a kind of round error of the *gsl_poly_solve_cubic* . I have > > tested your program with other GSL functions and they gave quite the same > > results. Take a look at the *gsl_poly_complex_solve_cubic *and * > > gsl_poly_complex_solve*. You will see nearly the expected result. > > > > Cheers, > > Ralph. > > > > gcc -static test.c -lgsl -lm > > GSL 1.10 > > Fedora 8 - Kernel 2.6.23.8-63.fc8 > > gcc 4.1.2 20070925 (Red Hat 4.1.2-33) > > Intel Core 2 Quad > > > > > > On Dec 20, 2007 11:34 AM, nick kama <[EMAIL PROTECTED] > wrote: > > > > > Hello, > > > > > > I have this problem with the solver of the cubic equations in gsl: > > > > > > if i take an equation of the form > > > > > > (x-1)^3=0 > > > > > > solver finds 3 roots > > > 1 1 1 > > > > > > if i took the equation > > > (x-1.1)^3 > > > > > > it finds only one root equal to 1.1,this is a little of code that can > > > show > > > this behaviour. > > > Documentation says :"the case of coincident roots is not considered > > > special. > > > For example, the equation (x-1)^3=0 will have three roots with exactly > > > equal > > > values " so this should'nt happen, pleas let me know where i'm wronging. > > > > > > Many Thanks > > > > > > Niccolo' > > > > > > > > > #include <iostream> > > > #include <gsl/gsl_matrix.h> > > > #include <gsl/gsl_eigen.h> > > > #include <gsl/gsl_vector.h> > > > #include <gsl/gsl_poly.h> > > > > > > > > > > > > main() > > > { > > > > > > double landa1=0; > > > double landa2=0; > > > double landa3=0; > > > //(x-1.1)^3 > > > int numero = gsl_poly_solve_cubic ( - 3.3,3.63, -1.331, &landa1, > > > &landa2, > > > &landa3); > > > std::cout << numero<< std::endl; > > > std::cout << landa1<< std::endl; > > > std::cout << landa2<< std::endl; > > > std::cout << landa3<< std::endl; > > > > > > // (x-1)^3 > > > numero = gsl_poly_solve_cubic ( -3,3, -1, &landa1, &landa2, &landa3); > > > std::cout << numero<< std::endl; > > > std::cout << landa1<< std::endl; > > > std::cout << landa2<< std::endl; > > > std::cout << landa3<< std::endl; > > > > > > } > > > _______________________________________________ > > > Help-gsl mailing list > > > [email protected] > > > http://lists.gnu.org/mailman/listinfo/help-gsl > > > > > > > > _______________________________________________ > Help-gsl mailing list > [email protected] > http://lists.gnu.org/mailman/listinfo/help-gsl > _______________________________________________ Help-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/help-gsl
