Several months ago, I wrote a C++ program (with C coding-style) that computes numerical approximations, of the roots of an arbitrary polynomial. Even though I changed the code slightly since March, mainly to improve the chances that it will compile on other people’s computers, the basic algorithm hasn’t changed. Yet, because I wrote the algorithm, I know that it has certain weaknesses, and I’m about to divulge one of those…
The main search algorithm isn’t able to determine whether (x) is sufficiently close to a root. It’s really only able to determine whether (y) is sufficiently close to zero. And so an idea that an implementer could have about this would be, to name the maximum allowable error in (x) ‘epsilon’, to compute the derivative of the polynomial for the given value of (x), and then to multiply ‘epsilon’ by this derivative, perhaps naming the derived, allowable error in (y) ‘epsilond’. But I consciously chose not to go that route because:
- Computing the real derivative of the polynomial each time, though certainly doable, would have made my code complex beyond what I was willing to tolerate,
- Even had I done so, certain polynomials will have a uselessly low derivative near one of their roots. This can happen because of a doubled root, or simply because of two roots that are too close together. ‘epsilon’ is as small as what would be practical, given the number of bits the variables have to work with, that were, a 52-bit significand.
So my program, instead, computes a ‘notional derivative’ of the polynomial, near (x), that is really just the derivative of the highest exponent by itself. This little detail could actually make it difficult for some observers to understand my code. At the same time, my program limits the lower extent of this derivative, to an arbitrary, modest fraction, just so that (x) won’t need to be impossibly precise, just to satisfy the computation of how precise (y) would need to be.
But then, a predictable problem becomes, that the real derivative of a given polynomial might be much smaller than just, the derivative of this highest exponent. And if it is, then satisfying … for (y) will no longer satisfy … for (x).
The following output illustrates an example of this behaviour:
"Maxima" Preparation: (%i1) expand((x-100)*(x-101)*(x-102)); (%o1) x^3-303*x^2+30602*x-1030200 (%i2) expand((x-100)*(x-110)*(x-120)); (%o2) x^3-330*x^2+36200*x-1320000 (%i3) expand((x-100)*(x-101)); (%o3) x^2-201*x+10100 Approximations from Program: dirk@Phosphene:~$ poly_solve 1 -303 30602 -1030200 101.00000000013 + 0I 99.999999999938 + 0I 101.99999999993 + 0I dirk@Phosphene:~$ poly_solve 1 -330 36200 -1320000 100 + 0I 110 + 0I 120 + 0I dirk@Phosphene:~$ poly_solve 1 -201 10100 100 + 0I 101 + 0I dirk@Phosphene:~$
(Updated 9/22/2019, 5h40 … )