help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] gsl_poly_solve_cubic problem


From: Andrew W. Steiner
Subject: Re: [Help-gsl] gsl_poly_solve_cubic problem
Date: Sat, 5 Jan 2008 11:31:40 -0500

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 <address@hidden> 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 <address@hidden> 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 <address@hidden > 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
> > > address@hidden
> > > http://lists.gnu.org/mailman/listinfo/help-gsl
> > >
> >
> >
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl
>




reply via email to

[Prev in Thread] Current Thread [Next in Thread]