A very rough implementation of Newton-Raphson method to find polynomial roots. Method taken from an interesting research paper..
The algorithm is recursively defined as $$x_{n+1} = x_n + \frac{f(x_n)}{f^{\prime} (x_n)}$$
where \(x_1\) is an initial guess. Notably, convergence point depends crucially on inital guess. Finding a set of initial guesses which allow us to find each root is a non-trivial task, mainly because a Newton map(which shows us where a particular starting point converges) is a fractal.
An initial point \(z_1\) is in the basin of a root \(r\) if the sequence \(z_{n+1} = z_{n} + \frac{p(z_n)}{p^{\prime} (z_n)}, n \geq 1\) converges to \(r\).
Newton map for \(p(z) = z^3 -2z +2\). Colours indicate the basins of each root. Red points don't converge at all.
This paper gives a recipe to construct a finite set of points such that, for every root, at least one of these points will converge to this root. Their main breakthrough is that only \(1.11 d \log^2 d\) inital points are required to guarantee the above condition, where \(d\) is degree of polynomial.The actual method
Because I am too dumb to understand the whole paper(hopefully this changes one day), I just followed the recipe given in section 9, aptly named as A recipe for dessert. Here is a small explanation of what I did shamelessly based on the recipe:
- Get a polynomial \(p(x) = a_n x^n + a_{n-1} x^{n-1} + \dots + a_1 x + a_0\)
- Find a radius \(R\) such that all roots lie within it. We can compute all bounds given on Wikipedia and choose the smallest one.
- Construct \(v\) circles and \(j\) points on each of them placed equidistant to each other. These are the initial guesses. The specifics of this step can be read in section 9
- Run newton's algorithm till we get \(p(x) < 0.000001\) or till we reach a certain number of cycles without convergence.
- Because coefficients of \(p(x)\) are real, we are guaranteed that \(\overline{z}\) is a root if \(z\) is a root.