# Polynomial roots, Davidenko method (1953) and the gamma trick **URL:** \ **Category:** Modeling and simulation\ **Created:** [May 21, 2025, 2:11pm UTC](https://scilab.discourse.group/t/polynomial-roots-davidenko-method-1953-and-the-gamma-trick/950 "2025-05-21T14:11:45Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [May 21, 2025, 2:11pm UTC](https://scilab.discourse.group/t/polynomial-roots-davidenko-method-1953-and-the-gamma-trick/950/1 "2025-05-21T14:11:45Z")
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 F(z,t)=t\gamma p\_1(z)+(1-t)p\_2(z), where \gamma\in\mathbb{C}^\* is arbitrary (its role will be clarified later). Knowing that F(z,0)=p\_2(z), \,F(z,1)=\gamma p\_1(z), 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 F(z(t),t)=0, \forall t\in[0,1], and its total derivative of both sides with respect to t gives \frac{d}{dt}F(z(t),t) = \frac{\partial F}{\partial z}(z(t),t)z'(t) + \frac{\partial F}{\partial t}(z(t),t)=0, 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 z'= -\left(\frac{\partial F}{\partial z}(z,t)\right)^{-1}\frac{\partial F}{\partial t}(z,t),\,t\in]0,1], 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 z'=-\cfrac{\gamma p\_1(z)-p\_2(z)}{t\gamma p\_1'(z)+(1-t)p\_2'(z)}, 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 z\_k=\exp\left( \frac{2k\pi i}{n}\right),\,k=0\dots n-1, 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: ```scilab 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: ```scilab 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: ```auto [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: ```scilab 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"; ``` ![graph2](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/4dcaeba23ff596416893f1388c303d034bbb9e19.svg) 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: ```auto [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"; ``` ![graph2](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/1f2fa832969077fd8526dccab94c60a716f84dc0.svg) 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. ---
**Author:** ![arctica](https://avatars.discourse-cdn.com/v4/letter/a/b9e5f3/32.png) [@arctica](https://scilab.discourse.group/u/arctica)\ **Post date:** [May 23, 2025, 2:55pm UTC](https://scilab.discourse.group/t/polynomial-roots-davidenko-method-1953-and-the-gamma-trick/950/2 "2025-05-23T14:55:07Z")
A very handy explanation and code implementation. For some reason the line: `legend(h(2:3),"$p_2$ roots","$p_1$ roots").interpreter="latex";` causes Scilab to crash; tested on 2025.0.0 and 2025.1.0 Lester ---
**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\ **Post date:** [May 23, 2025, 3:30pm UTC](https://scilab.discourse.group/t/polynomial-roots-davidenko-method-1953-and-the-gamma-trick/950/3 "2025-05-23T15:30:27Z")
Hello Lester, Thanks for testing. In fact the recursive assignment does not crash on macOS (I also tested and it crashes under Windows, an issue will be created on gitlab), I have modified the above code in order to make it pass for all platforms. S.