How can you compute all the roots of an arbitrary polynomial if you are only allowed to use an ode solver ?
Let us consider two polynomials p_1 and p_2 of degree n. A way of finding the the roots of p_1 by homotopy from the (known) roots of p_2 is to consider the two-variable function F:\mathbb{C}\times [0,1] \longrightarrow \mathbb{C} defined by
where \gamma\in\mathbb{C}^* is arbitrary (its role will be clarified later). Knowing that
if the roots of p_2 are known and distinct, we can for each of these roots (z_k)_{k=0\dots n-1} calculate the graph of the application t\longrightarrow z(t) where F(z(t),t)=0, by proceeding by homotopy from z(0)=z_k to z(1), a root of p_2. If everything works as expected, each path starting from z_k will lead to a given (and different) root of p_2.
There are several methods of doing this, for example by taking an arbitrary sequence (t_k) of t values and using the Newton’s method for solving F(z,t_k)=0 starting from the solution of F(z,t_{k-1})=0. In fact the (t_k) sequence cannot be that arbitrary and tuning it is not easy (backtracking is necessary in case of non-convergence).
Davidenko method
In my humble opinion a smarter way is to use the method proposed by a Russian mathematician (D.F. Davidenko, On a new method of numerical solution of systems of nonlinear equations (in Russian), Dokl. Akad. Nauk. SSSR, 88, 601-602, (1953)). His idea of taking the derivative of both sides of F(z(t),t)=0 with respect to t (which results in an ode) is not new, but he was convinced that ode solvers will gain automatic step size control mechanisms in the future, which would make following the path from t=0 to t=1 an easy task. However, as the actual implementation of such solvers started more than ten years later (e.g. the Runge-Kutta-Fehlberg method in 1969) the Davidenko method almost fell into oblivion and other path following methods became more popular. Hence, consider the following analysis and Scilab code as a modest tribute to his work.
Let us consider the method in the general case, i.e. when z(t)\in\mathbb{R^m}. We have
and its total derivative of both sides with respect to t gives
where \frac{\partial F}{\partial z} denotes the Jacobian matrix of F w.r.t. to z. If \frac{\partial F}{\partial z}(z(t),t) has full rank for all t\in[0,1], then t\rightarrow z(t) is well defined as the solution of the ode
with the initial condition z(0)=z_0 such that F(z_0,0)=0. Without any analysis, one cannot know if a solution exists on the whole domain [0,1], as the right hand side can present different types of singularity.
Application to polymomial root finding
In our application to polynomial roots (scalar case) we have F(z,t)=t\gamma p_1(z)+(1-t)p_2(z) and the ode writes
where we will take \gamma=1, p1(z)=\prod_{k=1}^n (z-k-1/2) and p_2(z)=z^n-1. We already know the roots of p_2
and p_1 roots are k+1/2 for k=1\dots n.
Scilab code
Below we take n=5 and use Scilab polynomials handling features:
n = 5;
p1 = poly((1:n)+0.5,"z");
p2 = %z^n-1;
d1 = derivat(p1);
d2 = derivat(p2);
We then define the Scilab function computing the rhs and the initial conditions:
function dzdt=f(t,z,gam)
dzdt = horner(p2-gam*p1,z)./horner(t*gam*d1+(1-t)*d2,z);
end
z0 = exp(2*%pi*%i*(0:n-1)/n);
We then call the arkode solver (with a high order method, Runge Kutta 8th order) and pass \gamma=1 to the rhs:
[t,y,info] = arkode(list(f,1),[0 1],z0,method="ERK_8",refine=4);
We were lazy by giving the whole vector z0 as initial condition as the n odes are not coupled, but this will result in a more compact code for plotting the trajectories below:
u = linspace(0,2*%pi,512);
h = plot(cos(u),sin(u),":k",real(z0),imag(z0),"^",(1:n)+1/2,zeros(1,n),'o');
plot(real(y.'),imag(y.'));
xlabel("Re(z(t))")
ylabel("Im(z(t))")
isoview on
hl = legend(h(2:3),"$p_2$ roots","$p_1$ roots");
hl.interpreter = "latex";
We miss 2 roots because two pairs of unit roots have paths that lead to the same root of p_1. The points where the paths merge are bifurcation points. Although the ode solver spent a lot of very small steps there (as the trajectory is not differentiable), it did not help to follow different paths. Studying and handling the singularity of these points is another subject, and we will solve the problem differently. In fact, the problem is the canonical (1,1-t) combination we used and the value \gamma=1 is a very bad choice: if we take \gamma=1+i/10 instead we recover all the expected roots:
[t,y,info] = arkode(list(f,1+%i/10),[0 1],z0,method="ERK_8",refine=4);
clf
h = plot(cos(u),sin(u),":k",real(z0),imag(z0),"^",(1:n)+1/2,zeros(1,n),'o');
plot(real(y.'),imag(y.'));
xlabel("Re(z(t))")
ylabel("Im(z(t))")
isoview on
hl = legend(h(2:3),"$p_2$ roots","$p_1$ roots");
hl.interpreter = "latex";
In the litterature this way of fixing the bifurcation problem has been called the gamma trick (see e.g. the book A. J. Sommese and C W. Wamplersee, Numerical Solution of Polynomial Systems Arising in Engineering and Science, World Scientific, 2005) and \gamma is usually taken as a random uniform variable on the unit circle.