One of the ideas which I’ve written about often is, that when certain Computer Algebra Software needs to compute the root of an equation, such as of a polynomial, an exact Algebraic solution, which is also referred to as the analytical solution, or symbolic Math, may not be at hand, and that therefore, the software uses numerical approximation, in a way that never churned out the Algebraic solution in the first place. And while it might sound disappointing, often, the numerical solution is what Engineers really need.
But one subject which I haven’t analyzed in-depth before, was, how this art might work. This is a subject which some people may study in University, and I never studied that. I can see that in certain cases, an obvious pathway suggests itself. For example, if somebody knows an interval for (x), and if the polynomial function of (x), that being (y), happens to be positive at one end of the interval, and negative at the other end, then it becomes feasible to keep bisecting the interval, so that if (y) is positive at the point of bisection, its value of (x) replaces the ‘positive’ value of (x) for the interval, while if at that new point, (y) is negative, its value for (x) replaces the ‘negative’ value of (x) for the interval. This can be repeated until the interval has become smaller than some amount, by which the root is allowed to be inaccurate.
But there exist certain cases in which the path forward is not as obvious, such as what one should do, if one was given a polynomial of an even degree, that only has complex roots, yet, if these complex roots nevertheless needed to be found. Granted, in practical terms such a problem may never present itself in the lifetime of the reader. But if it does, I just had lots of idle time, and have contemplated an answer.
(Updated 1/30/2019, 13h00 … )
(As of 1/21/2019 : )
I observed that a method which gets used to blend two values together in CGI, is as follows:
x = xa + t (xb – xa)
x = (1 – t) xa + t xb
And that for blend-values outside the range [0 .. +1], this method still works; it only generates an extrapolation instead of an interpolation. Thus, inversely, (t) can be found like so:
t = (x – xa) / (xb – xa)
Or, in some cases, if (x=0):
t = -xa / (xb – xa)
Well I puzzled together a hypothetical algorithm, which would be
far too complicated for me to implement, and I apologize in advance, if this does not conform to some people’s ideas of pseudo-code. But these are the schema in which I think.
(Updated 1/26/2019, 1h10 : )
There exists a question which the work-sheet above does not answer, which would be, What to do about doubled roots, aka multiplicities, in which 2 or more roots occur at the same value of (x). These tend to be important in the real-number domain because, if the curve changes direction before reaching the x-axis, complex roots exist, while, if the curve changes direction after crossing the x-axis, two different, real roots exist. What makes doubled roots problematic, is the fact that they have a derivative of zero, exactly where the root is. And while measured values should not work out that way, Mathematical examples exist that do.
In order to solve the doubled roots, what needs to be done is that the derivative of the polynomial must be computed, which in itself is easy to do. But then, the roots of the derivative might be found, and the original polynomial tested where its derivative is zero, to see whether the original one is sufficiently close to zero at the same value of (x). If so, a doubled root has been found.
The problem with the summary above is, that it needs to be made recursive. I.e., while the original polynomial might have a doubled root at (x), its derivative might also have one, but at a different value of (x)… The main problem with this last situation is, the fact that it could interfere in the ability to compute the roots at all, of the derivative.
If an easier opportunity doesn’t present itself, to find the roots of the first-order derivative, then the next thing to try, is to compute all the derivatives down to the linear equation, and for each root of a derivative found, to factorize that out of the representation of the parent equation, without polishing the division, if the parent equation is sufficiently close to zero at the same value of (x), so that its doubled root, no longer needs to be solved for. (:1)
The doubled root needs to be factorized out twice, from an equation ‘one up’ of the derivative, if the parent equation is close to zero where the derivative is. But from an equation ‘two up’ of the derivative, the root needs to be factorized out three times, and from an equation ‘three up’ it needs to be factorized four times… Assuming that all the parent equations were zero, where the lowest derivative was also zero. (:2)
And then, in the parts of all the equations which remain after factorizing, again, the roots need to be found, and potentially tested for in the respective parent equations. But, this does not require that more derivatives be computed. This is a situation in which the division may be polished.
For such reasons I can expect, that Computer Algebra Systems either
- Don’t check for multiplicities fully, or
- Only offer doing so as an option, which the user must enable.
The latter seems to be the case with “Sage”.
Another approach which could work would be, ‘only’ to compute the derivatives down to a quadratic, since the quadratic can always be solved using its general solution. In that case, whether it has doubled roots or not, can be determined according to whether:
b^2 == 4*a*c
The peculiarity of this approach would be, that in addition to solving the quadratic, the CAS would also need to record the multiplicity of its root(s) as just belonging to the quadratic, since next, the parent equation must be checked at those roots. And then, whether the quadratic produced two different, or one double root, needs to be treated differently when processing the parent equations(s).
(Updated 1/26/2019, 1h10 : )
As another way to organize what this posting is explaining, a subroutine can be devised, which can be given a suggested root of a polynomial as a parameter, and which may factorize that root out of the polynomial as long as the remaining equation is close to zero, at the suggested value of (x), without polishing the division. This subroutine may have as return value, the total number of times this was true.
Once the return-value is zero, the operation can be discontinued, to apply the suggested root to successive parent equations, of the present derivative.
If this is done, then that subroutine must also be given the task which my work-sheet above described, which is, in case a complex root is suggested, to factorize it out of the polynomial along with its conjugate, in the form of a quadratic equation, causing an equation to remain that only has real coefficients.
At some point when writing actual code, the decision would need to be made, which subroutines are responsible for what. And in many high-level languages, subroutines take the form of so-called ‘functions’.
It can be observed however, that:
- In my worksheet I defined a derived value for epsilon, which takes the present polynomial’s derivative into account, so that the original epsilon will define accuracy of (x) as much as possible, rather than accuracy of (y). But, some version of epsilon must be used to compare the absolute of (y) with, when testing the roots of the respective derivative-equation.
- If the roots of the derivatives have already been tested for in the present polynomial, and any occurrence as roots here factorized out, then none of the roots which remain to be found, should be doubled. And in that case, those roots may be polished,
- If the amount of error in the derivative is much smaller than (1), because the present polynomial is its integral, the amount of deviation of (y) from zero, in the present polynomial, that exists in spite of a doubled root, should be even smaller.
Therefore, an approach should still be possible, that uses a derived value of epsilon, actually to find roots.
(Update 1/30/2019, 13h00: )
This algorithm is no longer hypothetical. I have written actual C++ source-code to test it, and it does work.