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
>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_poly.h>
int main(void){
double landa1=0;
double landa2=0;
double landa3=0;
int i,numero;
double a[4] = {-1.331,3.63,-3.3,1.0};
double z[6];
gsl_complex ax,ay,az;
gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(4);
/* (x-1)^3 */
numero = gsl_poly_solve_cubic ( -3,3, -1, &landa1, &landa2, &landa3);
printf("--------------------\n");
printf("Numero = \%d\t(x-1.0)^3\tgsl_poly_solve_cubic\n",numero);
printf("Root[1] = %5.10lf\n",landa1);
printf("Root[2] = %5.10lf\n",landa2);
printf("Root[3] = %5.10lf\n",landa3);
/* (x-1.1)^3 */
landa1 = 0;
landa2 = 0;
landa3 = 0;
numero = gsl_poly_solve_cubic( -3.3000,3.63000,-1.331000, &landa1, &landa2,&landa3);
printf("--------------------\n");
printf("Numero = \%d\t(x-1.1)^3\tgsl_poly_solve_cubic\n",numero);
printf("Root[1] = %5.10lf\n",landa1);
printf("Root[2] = %5.10lf\n",landa2);
printf("Root[3] = %5.10lf\n",landa3);
/* (x-1.1)^3 */
numero = gsl_poly_complex_solve_cubic (-3.3000,3.63000,-1.331000, &ax, &ay, &az);
printf("--------------------\n");
printf("Numero = \%d\t(x-1.1)^3\tgsl_poly_complex_solve_cubic\n",numero);
printf("Root[1-Real] = %+.10lf\tRoot[1-Im] = %+.10lf\n",ax.dat[0],ax.dat[1]);
printf("Root[2-Real] = %+.10lf\tRoot[2-Im] = %+.10lf\n",ay.dat[0],ay.dat[1]);
printf("Root[3-Real] = %+.10lf\tRoot[3-Im] = %+.10lf\n",az.dat[0],az.dat[1]);
/* coefficients of P(x) = (x-1.1)^3 */
gsl_poly_complex_solve (a, 4, w, z);
gsl_poly_complex_workspace_free (w);
printf("--------------------\n");
printf("Numero = 3\t(x-1.1)^3\tgsl_poly_complex_solve\n",numero);
for (i = 0; i < 3; i++)
printf ("Root[%d-Real] = %+.10lf\tRoot[%d-Im] = %+.10lf\n",i+1,z[2*i],i+1,z[2*i+1]);
return 0;
}
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
