0

The GSL polynomial solver gives me an incorrect result. I have attempted to solve the following polynomial using GSL cubic polynomial solver:

x^3-1.96848e20 x^2+9.07605e28 x+9.07605e28 = 0

On Wolframalpha the results are:

  • -1
  • 4.61069e+8
  • 1.96848e+20

In my program the results are:

  • 2.30531e+08
  • 2.30531e+08
  • 1.96848e+20

The code I used:

   #include <iostream>
    #include <gsl/gsl_poly.h>

    using namespace std;

    int main() {
      double result[3]={0,0,0};
      gsl_poly_solve_cubic(-1.96848e20,9.07605e28,9.07605e28, &result[0], &result[1], &result[2]);
      cout << result[0] << endl;
      cout << result[1] << endl;
      cout << result[2] << endl;
      return 0;
    }

What could be wrong?

Edit: Actually, Wolframalpha is also giving weird answers

Update: Both Wolframalpha is wrong (a bit closer to the right answer) and my code is wrong.

Here's the same solution in Python (Numpy library)

In [11]: np.roots(p)
Out[11]: array([  1.96848000e+20,   4.61068948e+08,  -9.99999998e-01])
Oha
  • 11
  • 2
  • It seems that jenkins-traub polynomial roots algorithm gives more accurate results, whereas the GSL cubic polynomial root gives rounded answers. In this case the order of the answers are around -1 ~ 2e20, which might give hint why the GSL algorithm fails. The Jenkins-Traub solver seems to give more accurate answers in this case. – Oha Sep 10 '15 at 08:44
  • No, there is no rounding involved in the GSL routine, only a numerically unstable case detection method that falsely concludes for a double root. Avoiding that branch would have resulted in the correct answer. – Lutz Lehmann Sep 13 '15 at 06:43

1 Answers1

1

The numerical method used in GSL is not accurate enough in this case.

Use e.g. Jenkins-Traub algorithm (code here)

Oha
  • 11
  • 2