Introduction to Linear Homotopy

Béla  Paláncz

Abstract

Solving nonlinear systems, especially algebraic polynomial systems, is a fundamental problem which occurs frequently in various models of science and engineering like robotics, computer vision, computational geometry, signal processing. The application of  symbolic elimination methods, like reduced Gröbner basis and Dixon resultant is  restricted by the size and complexity of the system to be solved, therefore mostly numerical methods should be employed. Frequently, we have no proper information about the locations of the roots, which is important for using local numerical methods. Local numerical methods, like the very popular Newton-Raphson method,  require initial guess, sufficiently close to the roots being sought. Therefore we need global numerical methods, which are very robust concerning the start values of the iteration.  Homotopy continuation proved to be  as a reliable and efficient method to solve nonlinear systems as well as  polynomial systems for the last two decades, is a globally convergent method. It does not need accurate initial guess and can be used to locate all  geometrically isolated solutions of a square polynomial system. In this  tutorial type study linear  homotopy continuation methods are introduced. The definitions and  basic ideas  are considered and illustrated with many examples in Mathematica. The general algorithm of  linear homotopy method in form of  repeating solution of homotopy equation via Newton - Raphson method  as well as in form of  an initial value problem of a system of ordinary differential equations, are implemented in Mathematica. In addition, a function providing  automatic construction of a start system for polynomial systems is also developed. The efficiency of  these functions have been illustrated  by solving  nonlinear as well as polynomial equation systems.

Introduction

Solving nonlinear systems, especially algebraic polynomial systems, is a fundamental problem which occurs frequently in various models of science and engineering like robotics, computer vision, computational geometry, signal processing. Using numerical algorithms to solve polynomial systems with tools from algebraic geometry is the main activity in the so called Numerical Algebraic Geometry, see Sommese et al (2003). This is a new developing field on the crossroads of algebraic geometry, numerical analysis, computer science and engineering. Homotopy continuation method is a global numerical method to solve nonlinear system in general, and also polynomial systems, see Gritton at al (2001). This method is used to locate all geometrically isolated roots as well as decompose positive dimensional solution sets into irreducible components,  Verschelde (2007). Homotopy continuation has been established as a reliable and efficient method to solve polynomial systems for the last two decades, originating from the works of Drexler (1977) and Garcia and Zangwill (1979).

Some computer codes have been developed for homotopy methods. The most efficient is written in Maple, Verschelde (2007) and also in Matlab, Decarolis et al. Recently, codes for fix point and Newton homotopy were developed in Mathematica, too, see Binous (2007a) and Binous (2007b).

In this  tutorial type study linear  homotopy continuation methods are introduced, see Kotsireas (2001). The definitions and  basic ideas  are considered and illustrated with many examples in Mathematica. The general algorithm of  linear homotopy method in form of  repeating solution of homotopy equation via Newton - Raphson method  as well as in form of  an initial value problem of a system of ordinary differential equations, are implemented in Mathematica. In addition, a function providing  automatic construction of a start system for polynomial systems is also developed. The efficiency of  these functions have been illustrated  by solving  nonlinear as well as polynomial equation systems.

Definition of homotopy

The continuous deformation of an object to an other object is called as homotopy. Consider two univariate polynomials (Subscript[p, 1]) and (Subscript[q, 1]) of the same degree,

Example 1

(Subscript[p, 1])[x_]:=(x^2)-3

(Subscript[q, 1])[x_]:=-(x^2)-x+1

Define the linear convex function, H  in variables x and  λ,  called as homotopy function as,

H(x,λ)=(1-λ) (Subscript[p, 1])(x)+λ (Subscript[q, 1])(x)

H[x_,λ_]:=(1-λ) (Subscript[p, 1])[x]+λ (Subscript[q, 1])[x]

In geometric terms, the homotopy H provides us a continuous deformation from (Subscript[p, 1]) -which is obtained for λ = 0 by H(x, 0) - to (Subscript[q, 1]) - which is obtained for λ = 1 by H(x, 1). We call it linear homotopy because H is a linear function of the variable λ.

Show[{Plot[Table[H[x,λ],{λ,0,1,0.1}],{x,-2,1.5},AxesLabel->{"x"}],Graphics[Text["(Subscript[p, 1])",{-0.5,-3.}]],Graphics[Text["(Subscript[q, 1])",{0.2,1.}]]}]

Now let us consider two univariate polynomials (Subscript[p, 2]) and (Subscript[q, 2]) of the different degrees,

Example 2

(Subscript[p, 2])[x_]:=(x^2)-1

(Subscript[q, 2])[x_]:=((x^2)-(1/4))((x^2)-4)

Define the linear convex homotopy again with,

H[x_,λ_]:=(1-λ) (Subscript[p, 2])[x]+λ (Subscript[q, 2])[x]

Show[{Plot[Table[H[x,λ],{λ,0,1,0.1}],{x,-2.5,2.5},AxesLabel->{"x"}],Graphics[Text["(Subscript[p, 2])",{-1.5,2.}]],Graphics[Text["(Subscript[q, 2])",{0.9,-2.4}]]}]

The (Subscript[p, i])'s are called start systems and (Subscript[q, i])'s are called target systems.

Solving nonlinear equation via homotopy

Homotopy continuation method deforms  the known roots of the start system into the roots of the target system. Now, let us look at how homotopy can be used to solve a simple polynomial equation. Consider the polynomial equation of degree two,

q(x)=(x^2)+8 x-9=0

Example 3

q[x_]:=(x^2)+8 x-9

By deleting the middle term, we can get a simpler equation, which can be solved easily by inspection,

p[x_]:=(x^2)-9

This equation also has two roots and will be considered as start system for the target system (2). The linear homotopy can be defined as it follows

H[x_,λ_]:=(1-λ) p[x]+λ q[x]

or

H[x,λ]//Expand

Show[{Plot[Table[H[x,λ],{λ,0,1,0.1}],{x,-10,10},AxesLabel->{"x"}],Graphics[Text["p(x)",{8,40.}]],Graphics[Text["q(x)",{8,140}]]}]

Homotopy continuation method deforms p(x) = 0, the known roots of the start system,  into q(x) = 0, the roots of the target system. Let us solve the equation  H(x, λ) = 0 for different values of λ. Consider   (Subscript[x, 0]) = 3 one of the solutions of  p(x) = 0 as initial guess value, and solve H(x, (Subscript[λ, 1])) = 0, where let (Subscript[λ, 1])= 0.2, employing Newton-Raphson method, repeatingly

(Subscript[x, 0]) = 3;

(Subscript[λ, 1])=0.2;

(Subscript[x, 1])=x/.FindRoot[H[x,(Subscript[λ, 1])]==0,{x, (Subscript[x, 0])}]

(Subscript[λ, 2])=0.4;

using the result as guess value for the next solution

(Subscript[x, 2])=x/.FindRoot[H[x,(Subscript[λ, 2])]==0,{x,(Subscript[x, 1])}]

(Subscript[λ, 3])=0.6;

(Subscript[x, 3])=x/.FindRoot[H[x,(Subscript[λ, 3])]==0,{x,(Subscript[x, 2])}]

(Subscript[λ, 4])=0.8;

(Subscript[x, 4])=x/.FindRoot[H[x,(Subscript[λ, 4])]==0,{x,(Subscript[x, 3])}]

(Subscript[λ, 5])=1;

(Subscript[x, 5])=x/.FindRoot[H[x,(Subscript[λ, 5])]==0,{x,(Subscript[x, 4])}]

Show[{Plot[Table[H[x,λ],{λ,0,1,0.2}],{x,0,5},AxesLabel->{"x"}],ListPlot[{{3,0},{(Subscript[x, 1]),0},{(Subscript[x, 2]),0},{(Subscript[x, 3]),0},{(Subscript[x, 4]),0},{(Subscript[x, 5]),0}},PlotStyle->{PointSize[0.016],RGBColor[1,0,0]}],Graphics[Text["p(x)",{5,10.}]],Graphics[Text["q(x)",{4,45}]]}]

(Subscript[λ, 0])=0.;

Show[{ListPlot[Table[{(Subscript[λ, i]),(Subscript[x, i])},{i,0,5}],Joined->True,AxesLabel->{"λ","x(λ)"}],ListPlot[Table[{(Subscript[λ, i]),(Subscript[x, i])},{i,0,5}],PlotStyle->{PointSize[0.016],RGBColor[1,0,0]}],Graphics[Text["p(3)= 0",{0.1,3.}]],Graphics[Text["q(1)=0",{0.95,1.25}]]}]

Tracing homotopy curves as initial value problem

Comparing homotopy solution with the traditional Newton-Raphson solution,  if  Δλ is small enough, the convergence may be ensured in every step.

One can consider this root tracing problem as an initial value problem of an ordinary differential equation. While H (x, λ) = 0 for every  λ ∈ [0, 1], therefore

dH (x,λ)=((∂H)/(∂x))dx+((∂H)/(∂λ))dλ ≡ 0               λ∈[0,1]

Then the initial value problem is

(Subscript[H, x])((dx(λ))/(dλ))+(Subscript[H , λ])= 0

with

x(0) = (Subscript[x, 0])

Here (Subscript[H, x]) is the Jacobian of H with respect to (Subscript[x, i]), i = 1,..., n, in case of n nonlinear equations with n variables.

In our single variable case, the two partial derivates of the homotopy function are

dHdλ=D[H[x,λ],λ]

dHdx=D[H[x,λ],x]

Then the right hand side of the differential equation to be solved is

deqrhs=-((dHdλ)/dHdx)/.x->x[λ]

The equation itself

deq=D[x[λ],λ]==deqrhs

The initial value is

x0=(Subscript[x, 0])

The numerical solution

sol=NDSolve[{deq,x[0]==x0},{x[λ]},{λ,0,1}]

The trajectory is

Plot[x[λ]/.sol,{λ,0,1},AxesLabel->{"λ","x(λ)"},Epilog->{PointSize[0.016],Blue,Point[{0,x0}],PointSize[0.016],Red,Point[{1,First[(x[λ]/.sol)/.λ->1]}]}]

Then the  correspondig root of q(x) is

First[x[λ] /. sol/.λ->1]

Types of linear homotopy

General linear homotopy

As we have seen in Example 3, the start system can be constructed intuitively, reducing the original system (target system) to a more simple system, which roots can be easily computed. In order to get all of the roots of the target system, the start system should have so many roots as the target system has. 
The start system can be constructed in many ways, however there two typical type of  the start sytems, which are usually employed.

Fixedpoint homotopy

In (1), the start system can be considered as

p(x) = x-(Subscript[x, 0])

where (Subscript[x, 0]) is a guess value for the root of the target system q(x) = 0.

In that case the homotopy function is

H(x,λ) =(1-λ)(x-(Subscript[x, 0]))+λ q(x)

Newton homotopy

An other construction for the start system, when we consider p (x) as

p(x) = q(x)-q((Subscript[x, 0]))

in that case the homotopy function is

H(x,λ) =(1-λ)(q(x)-q((Subscript[x, 0])))+λ q(x)

or

H(x,λ) = q(x)-(1-λ)q((Subscript[x, 0]))

There are other methods to  construct start system for linear homotopy. In a section, following later, we shall see how one can define start system for polynomial systems automatically.

Newton-Raphson vs. homotopy method

It is known that local methods, like Newton-Raphson method, require initial guess of a root, which is close to the intended root for the particular application. Global methods, like homotopy continuation method can find solution from a starting guess, which is far from the solution. Now, we shall illustrate some situations when Newton-Raphson method fails, but homotopy method can be successfull. We shall consider two usual problems with  the  Newton-Raphson method, namely when the Jacobi matrix is singular and when the convergency is too slow.

Singular Jacobi matrix

Example 4

Let us consider the following monomial of degree three,

q[x_]:=27 (x^3)+108 (x^2)+108 x+40

This polynomial has one real and  two conjugate complex roots.

Plot[q[x],{x,-3,1},AxesLabel->{"x","p(x)"}]

Let us try to find the real root using Newton Raphson method, with an initial guess value (Subscript[x, 0])= -2,

FindRoot[q[x],{x,-2}]

This method also fails when (Subscript[x, 0])= -1.5

FindRoot[q[x],{x,-1.5}]

Now we use homotopy method. Let us employ fixpoint homotopy, see (7)

x0=-2;

H[x_,λ_]:=(1-λ) (x-x0)+λ q[x]

Show[{Plot[Table[H[x,λ],{λ,0,1,0.1}],{x,-3,-2},AxesLabel->{"x","p(x), q(x)"},PlotRange->{-10,40}],Graphics[Text["p(x) = x - (Subscript[x, 0])",{-2.15,2}]],Graphics[Text["q(x)",{-2.45,33}]]}]

The right hand  side of the differential equation is,

deqrhs=-((D[H[x,λ],λ])/(D[H[x,λ],x]))/.x->x[λ]

and

deq=D[x[λ],λ]==deqrhs

The solution is

sol=NDSolve[{deq,x[0]==x0},{x[λ]},{λ,0,1}]

Plot[x[λ]/.sol,{λ,0,1},PlotRange->{-1.9,-2.9},AxesLabel->{"λ","x(λ)"},Epilog->{PointSize[0.016],Blue,Point[{0,x0}],PointSize[0.016],Red,Point[{1,First[(x[λ]/.sol)/.λ->1]}]}]

A value of the real root is

First[x[λ]/.sol/.λ->1]

This fixpoint homotopy also works with (Subscript[x, 0]) = - 1.5.

Slow convergency

Let us consider now a two- dimensional problem, Binous (2007a). The target sytem is

Example 5

f1[x_,y_]:=x^2-x y-1

f2[x_,y_]:=Sqrt[x+1]+y

Show[{ContourPlot[{f1[x,y]==0,f2[x,y]==0},{x,-1.1,2},{y,-2,1},FrameLabel->{"x","y"}],Graphics[Text["(Subscript[f, 1])(x,y)= 0",{1.38,0}]],Graphics[Text["(Subscript[f, 2])(x,y)=0",{-0.5,-1}]],Graphics[Text["(Subscript[f, 1])(x,y)=0",{-0.5,0.3}]]}]

We try  Newton-Raphson method with (Subscript[x, 0])= - 0.1 and (Subscript[y, 0])=  0.5,

x0=-0.1;y0=0.5;

FindRoot[{f1[x,y],f2[x,y]},{x,x0},{y,y0}]

Now let us employ Newton homotopy, see (10),

H[x_,y_,λ_]:=Flatten[-(1-λ) (({{f1[x0,y0]}, {f2[x0,y0]}}))+ (({{f1[x,y]}, {f2[x,y]}}))]

H[x,y,λ]

Show[{ContourPlot[{H[x, y, 0][[1]] == 0, H[x, y, 1][[1]] == 0}, {x, -1, 2}, {y, -2, 1},FrameLabel->{"x","y"}],Graphics[Text["H(x,y,1)= 0",{-0.7,0.}]],Graphics[Text["H(x,y,0)=0",{0.25,0.5}]],Graphics[Text["H(x,y,1)=0",{0.9,-1.2}]],ContourPlot[{H[x, y, 0][[2]] == 0, H[x, y, 1][[2]] == 0}, {x, -1, 2}, {y, -2, 1}]}]

Let us define the Jacobi computation as it follows

Jacobi[F_,X_]:=Outer[D,F,X]

Then the total Jacobian, (((∂H)/(∂x)) ((∂H)/(∂y))((∂H)/(∂λ)))

Jxλ=Jacobi[H[x,y,λ],{x,y,λ}];MatrixForm[Jxλ]

The Jacobian respecting to x and y

Jx=Take[Jxλ,{1,2},{1,2}];MatrixForm[Jx]

The right hand side of the linear system

b=-Take[Jxλ,{1,2},{2+1}];MatrixForm[b]

Solving the linear system,

rhs=(LinearSolve[Jx,b]//Simplify)/.Map[#->#[λ]&,{x,y}]//Flatten

Then the left hand side of the differential equation system is,

lhs=Map[D[#[λ],{λ,1}]&,{x,y}]

The differential equation system itself

deqs=MapThread[#1==#2&,{lhs,rhs}]

The numerical solution

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1}]

Show[{Plot[{x[λ],y[λ]}/.sol,{λ,0,1},FrameLabel->{"λ","x(λ), y(λ)"},Frame->True,Epilog->{PointSize[0.016],Blue,Point[{{0,x0},{0,y0}}],PointSize[0.016],Red,Point[{{1,First[(x[λ]/.sol)/.λ->1]},{1,First[(y[λ]/.sol)/.λ->1]}}]}],Graphics[Text["x(λ)",{0.6,-0.6}]],Graphics[Text["y(λ)",{0.4,0.4}]]}]

The solution of the system,

({x[λ],y[λ]}/.sol)/.λ->1

Regularization of  homotopy

Sometimes this form of homotopy methods may also fail, because of singularties resulting diverging path. Let us consider the following system,

Example 6

f1[x_,y_]:=x^2+ y-3

f2[x_,y_]:=x+1/8 y^2-1

Show[{ContourPlot[{f1[x,y]==0,f2[x,y]==0},{x,-10,10},{y,-10,10},FrameLabel->{"x","y"}],Graphics[Text["(Subscript[f, 1])(x,y)= 0",{0,-9}]],Graphics[Text["(Subscript[f, 2])(x,y)=0",{-8,6.5}]]}]

First, let us  try to solve the problem using Newton Rapson method, starting with  (Subscript[x, 0])= - 1 and (Subscript[y, 0])=  -1,

FindRoot[{f1[x,y],f2[x,y]},{x,-1},{y,-1}]

This result is not correct

{f1[x, y], f2[x, y]} /.%

or with (Subscript[x, 0])=  1 and (Subscript[y, 0])=  -1,

FindRoot[{f1[x,y],f2[x,y]},{x,1},{y,-1}]

again, this is not a solution

{f1[x, y],f2[x, y]}/.%

These mean that none of them is a correct solution, Newton - Rapson fails.

Let us try to solve the problem with homotopy method! We consider the following obvious functions as start system, deleting the low  order terms, we get

g1[x_,y_]:=x^2-1

g2[x_,y_]:=y^2-1

The zeros of this system are (1, 1), (-1, 1), (-1, -1) and (1, -1). Let us consider

x0=1;y0=-1;

The homotopy function is

H[x_,y_,λ_]:=Flatten[(1-λ) (({{g1[x,y]}, {g2[x,y]}}))+ λ(({{f1[x,y]}, {f2[x,y]}}))]

H[x,y,λ]

Show[{ContourPlot[{H[x,y,0][[1]]==0,H[x,y,0.2][[1]]==0,H[x,y,1][[1]]==0},{x,-10,10},{y,-10,10},FrameLabel->{"x","y"}],ContourPlot[{H[x,y,0][[2]]==0,H[x,y,0.2][[2]]==0,H[x,y,1][[2]]==0},{x,-10,10},{y,-10,10}],Graphics[Text["H(x, y, 1)= 0",{-5,5.5}]],Graphics[Text["H(x, y, 0)= 0",{7,1.5}]],Graphics[Text["H(x, y, 1)= 0",{5,-8}]],Graphics[Text["H(x, y, 0.5)= 0",{-8,2.5}]],Graphics[Text["H(x, y, 0)= 0",{4,8}]],Graphics[Text["H(x, y ,0.5)= 0",{-3.5,-9}]]}]

The start system (vertical and horizontal paralell lines,  with blue color) are deforming into the target system, with yellow color.

The total Jacobian,

Jxλ=Jacobi[H[x,y,λ],{x,y,λ}];MatrixForm[Jxλ]

The Jacobian respecting to x and y

Jx=Take[Jxλ,{1,2},{1,2}];MatrixForm[Jx]

The right hand side of the linear system

b=-Take[Jxλ,{1,2},{2+1}];MatrixForm[b]

Solving the linear system,

rhs=(LinearSolve[Jx,b]//Simplify)/.Map[#->#[λ]&,{x,y}]//Flatten

Then the differential equation system is,

lhs=Map[D[#[λ],{λ,1}]&,{x,y}]

deqs=MapThread[#1==#2&,{lhs,rhs}]

The numerical solution

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1},Method->StiffnessSwitching]

The integration has failed because of singularity of the Jacobi matrix of  H(x, y,  λ) occured.

The denominator of the right hand side of the differential equations

(2 (λ^2)+(-8+7 λ) x[λ] y[λ])/.sol/.λ->0.6638182503426503

The determinant of  the Jacobi matrix of  H (x, y,  λ)

Jx/.{x->x[λ],y->y[λ]}/.sol/.λ->0.6638182503426503

Det[First[%]]

In order to avoid singulatrity in real field, we consider a modified complex homotopy function,

H(x,λ)=γ(1-λ) (Subscript[p, 1])(x)+λ (Subscript[q, 1])(x)

where γ is a random complex number. For almost all choices of a complex constant γ, all solution paths defined by the homotopy above are regular, i.e.: for all λ∈ [0, 1], the Jacobian matrix of  H (x, λ) is regular and no path diverges. Let us now consider,

γ=1+I;

Then the homotopy

H[x_,y_,λ_]:=Flatten[γ (1-λ) (({{g1[x,y]}, {g2[x,y]}}))+ λ(({{f1[x,y]}, {f2[x,y]}}))]

H[x,y,λ]

The total Jacobian,

Jxλ=Jacobi[H[x,y,λ],{x,y,λ}];MatrixForm[Jxλ]

The Jacobian respecting to x and y

Jx=Take[Jxλ,{1,2},{1,2}];MatrixForm[Jx]

The right hand side of the linear system

b=-Take[Jxλ,{1,2},{2+1}];MatrixForm[b]

Solving the linear system,

rhs=(LinearSolve[Jx,b]//Simplify)/.Map[#->#[λ]&,{x,y}]//Flatten

Then the differential equation system is,

lhs=Map[D[#[λ],{λ,1}]&,{x,y}]

deqs=MapThread[#1==#2&,{lhs,rhs}]

The numerical solution

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1}]

The solution of the system,

({x[λ],y[λ]}/.sol)/.λ->1

Now we have got a  complex  root,

Show[{ParametricPlot[{Re[y[λ]],Im[y[λ]]}/.sol,{λ,0,1},PlotRange->All,FrameLabel->{"Re","Im"},Frame->True,Epilog->{PointSize[0.02],Blue,Point[{Re[y0],Im[y0]}],PointSize[0.02],Red,Point[{First[(Re[y[λ]]/.sol)/.λ->1],First[(Im[y[λ]]/.sol)/.λ->1]}],PointSize[0.02],Blue,Point[{Re[x0],Im[x0]}],PointSize[0.02],Red,Point[{First[(Re[x[λ]]/.sol)/.λ->1],First[(Im[x[λ]]/.sol)/.λ->1]}]}],ParametricPlot[{Re[x[λ]],Im[x[λ]]}/.sol,{λ,0,1},PlotRange->All],Graphics[Text["x(λ)",{1.5,0.5}]],Graphics[Text["y(λ)",{0.5,-2}]]}]

Basically, our original  system, the target system has four solutions,

NSolve[{f1[x, y]==0, f2[x, y]==0},{x,y}]

The other roots also can be computed using different starting values, namely

x0=1;y0=1;

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1}];

({x[λ],y[λ]}/.sol)/.λ->1

and the real roots

x0=-1;y0=1;

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1}];

({x[λ],y[λ]}/.sol)/.λ->1

x0=-1;y0=-1;

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1}];

({x[λ],y[λ]}/.sol)/.λ->1

These complex "tails" can be eliminated by using high precision computation. For example

sol=NDSolve[Join[deqs,{x[0]==x0,y[0]==y0}],{x[λ],y[λ]},{λ,0,1},AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->30];

({x[λ],y[λ]}/.sol)/.λ->1

Chop[%]//N

Starting system for polynomial systems

How can we find the proper start system, which will provide all of the solutions of the target system? This problem can be solved if the nonlinear system is specially a system of polynomial equations.

Let us consider the case when we are looking for  the homotopy solution of  f(x) = 0, where f(x) is a polynomial system,   f(x): (ℝ^n)->(ℝ^n). To get all of the solutions,  one should find out a proper polynomial system, as start system,  g(x) = 0, where g(x): (ℝ^n)->(ℝ^n) with known and easily computable solutions.

An appropriate start system can be generated in the following way, Hazaveh et al. (2003)

Let (Subscript[f, i])((Subscript[x, 1]),...,(Subscript[x, n])), i = 1,...,n be a system of n polynomials. We are interested of the common
 zeros of the system, namely  f = ((Subscript[f, 1])(x),...,(Subscript[f, n])(x)) = 0.

Let (Subscript[d, j]) denote the  degree of the jth polynomial - that is the degree of the highest order monomial
in the equation. Then such a starting system is,

(Subscript[g, j])(x)= (e^(i (Subscript[ϕ, j ]) )) (((Subscript[x, j])^(Subscript[d, j]))-((( (e^(i (Subscript[θ, j ]) ))))^(Subscript[d, j])) = 0,      j = 1,...,n

where (Subscript[ϕ, j ]) and (Subscript[θ, j ])are random real numbers in the interval [0, 2 π]. The equation above has the obvious particular solution (Subscript[x, j])= (e^(i (Subscript[θ, j ]) ))and the complete set of the starting solutions for  j = 1, ..., n is given by

(e^((i (Subscript[θ, j ])+((2 π i k)/(Subscript[d, j]))))),  k=0,1,...,(Subscript[d, j])-1

Bézout's theorem, see Kotsireas (2001), states that the number of isolated roots of such a system is bounded by the total degree of the system,  (Subscript[d, 1])(Subscript[d, 2])...(Subscript[d, n.])

Let us consider the following system

Example 7

f1[x_,y_]:=x^2+y^2-1

f2[x_,y_]:=x^3+y^3-1

The degrees of the polynomials are

d1=2; d2=3;

Indeed, this system as it is expected, has the following six roots ((Subscript[d, 1])(Subscript[d, 2]) = 2 * 3 = 6)

NSolve[{f1[x,y]==0,f2[x,y]==0},{x,y}]

Now we compute the start system considering  (12) and (13). We generate random real numbers in the interval [0, 2π] as it follows

ϕ1=Random[Real,{0,2 π}]

ϕ2=Random[Real,{0,2 π}]

θ1=Random[Real,{0,2 π}]

θ2=Random[Real,{0,2 π}]

The starting system is, see (12)

g1[x_,y_]:=(E^(I ϕ1)) ((x^d1)-((((E^(I θ1))))^d1))

g2[x_,y_]:=(E^(I ϕ2)) ((y^d2)-((((E^(I θ2))))^d2))

Then the complete set of the starting solutions, see (13)

Xi=Table[Exp[I θ1+2 π I k/d1],{k,0,d1-1}]

and

Yi=Table[Exp[I θ2+2 π I k/d2],{k,0,d2-1}]

We need all of the combination of the initial values {(Subscript[X, i]),(Subscript[Y, j])}, i=1, 2,  j=1, 2

X0=Tuples[{Xi, Yi}]

So we have six  initial values for solving our homotopy equations. These initial values will provide the start points of the six homotopy paths. The end points of these paths are the six solutions.

Before carring out the computation, it is high time to implement four Mathematica functions to carry out linear homotopy computations in Mathematica environment.

Implementation of homotopy methods  in Mathematica

Function for computing start system for polynomials

This function will compute the start system as well as the initial values for the computation of the homotopy paths in case of polynomial system.

Input variables
F   - list of functions of the polynomial system to be solved, F = {(Subscript[f, 1])(x), (Subscript[f, 2])(x),...,(Subscript[f, n])(x)}
X  - list of the independent variables X =  {(Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, n])}

Output variables

G    -  list of  start  system, G = {(Subscript[g, 1])(x), (Subscript[g, 2])(x),...,(Subscript[g, n])(x)}
X0  -  list of initial values, X0 = {{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_1,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_2,...,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_m}
         where m is the number of the roots of  F

StartingSystem[F_,X_]:=Module[{n,d,ϕ,θ,G,Xi,X0},
n=Length[X];
d=Map[Max[Exponent[#,X]]&,F];
ϕ=Table[Random[Complex,{0,2 π}],{n}];
θ=Table[Random[Complex,{0,2 π}],{n}];
G=MapThread[(E^(I #1)) (((#4)^(#2))-((((E^(I #3))))^(#2)))&,{ϕ,d,θ,X}];
Xi=MapThread[Table[Exp[I #1+2 π I k/#2],{k,0,#2-1}]&,{θ,d}];
X0=Tuples[Xi];
{G,X0}]

Now, let us  create the start system with this function. 
The list of the functions of the system we considered above, see Example 7

Example 8

F8={f1[x,y],f2[x,y]}

The list of variables

X8={x,y}

Now, employing our function

sol8=StartingSystem[F8,X8];

The start system is

G8=sol8[[1]]

and the initial values are

X08=sol8[[2]]

Function for direct path tracing

This function will compute the homotopy paths using successive Newton-Raphson method, see Example 3.

Input variables
F   - list of functions of the polynomial system to be solved, F = {(Subscript[f, 1])(x), (Subscript[f, 2])(x),...,(Subscript[f, n])(x)}
G  - list of the start system, G = {(Subscript[g, 1])(x),(Subscript[g, 2])(x),...,(Subscript[g, n])(x)}, 
X  - list of the independent variables X =  {(Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, n])}
X0  -  list of initial values, X0 = {{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_1,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_2,...,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_m}
         where m is the number of the roots of  F 
 γ   -  list of complex numbers, {(Subscript[γ, 1]), (Subscript[γ, 2]),...,(Subscript[γ, n])}, it means (Subscript[γ, i]) can be different for every (Subscript[g, i])(x), 
          see (11). If start system is generated for polynomials all  (Subscript[γ, i]) = 1.
 n  -   number of the subintervalls  in [0, 1].

Output variables

sol[1]    -  list of the i-th solutions correspondig to the i-th initial values,  {(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_i, i = 1...m
sol[2]    -  list of the path of  i-th solutions in form of interpolating functions of the variables  
                 correspondig to the i-th  initial values,  {(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_i,   i = 1...m,
                 {{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_1,{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_2,...,{(Subscript[(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }, m])},  (Subscript[φ, i])=(Subscript[φ, i])(λ)

LinearHomotopyFR[F_,G_,X_,X0_,γ_,n_]:=Module[{H,X0L,λ0,i,β,m,R,RR,j,k,sol},
Off[FindRoot::"lstol"];
λ0=Table[i (1/n),{i,0,n}];
H=Flatten[(1-β) Thread[G  γ]+β F];
m=Length[X];
k=Length[X0];
RR={};
Do[
X0L={X0[[j]]};
Do[AppendTo[X0L,Map[#[[2]]&,FindRoot[H/.β->λ0[[i+1]],MapThread[{#1,#2}&,{X,X0L[[i]]}]]]],{i,1,n}];
R={};
Do[AppendTo[R,Interpolation[MapThread[{#1,#2[[i]]}&,{λ0,X0L}]]],{i,1,m}];AppendTo[RR,{Map[Chop[N[#[1]]]&,R],R}],
{j,1,k}];
sol=Transpose[RR];{sol[[1]],Table[MapThread[#1[λ]->#2[λ]&,{X,sol[[2,i]]}],{i,1,Length[X0]}]}]

Remark
It is reasonable to carry out computation with n = 100, however sometimes one may need more  subintervalls, but it also possible that less will enough to get satisfactory result.

Let us compute the roots of the system above,

γ={1,1};

sol= LinearHomotopyFR[F8,G8,X8,X08,γ,100];

The roots are

sol1FR=Chop[sol[[1]],(10^(-8))]

NSolve (#) provides the same result

NSolve[F8,X8]

The interpolation functions

sol2FR=sol[[2]]

Function for  path tracing by integration

Computing inverse

This function will compute the homotopy paths using numerical integration of the differencial equation system, see (4) and (5) using Runge-Kutta method. In this implementation we compute the of the Jacobian inverse, namely the initial value problem  (4), (5) will be expressed as

((dx(λ))/(dλ)) = - ((Subscript[H, x])^(-1)) (Subscript[H , λ])

with

x(0) = (Subscript[x, 0])

Here ((Subscript[H, x])^(-1)) is the inverse of the  Jacobian of (Subscript[H, x]) with respect to (Subscript[x, i]), i = 1,..., n, in case of n nonlinear equations with n variables. The inverse will be computed in symbolic way therefore the number of equations may be limited by n = 6 or 8!

Input variables
F   - list of functions of the polynomial system to be solved, F = {(Subscript[f, 1])(x), (Subscript[f, 2])(x),...,(Subscript[f, n])(x)}
G  - list of the start system, G = {(Subscript[g, 1])(x),(Subscript[g, 2])(x),...,(Subscript[g, n])(x)}, 
X  - list of the independent variables X =  {(Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, n])}
X0  -  list of initial values, X0 = {{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_1,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_2,...,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_m}
         where m is the number of the roots of  F 
 γ   -  list of complex numbers, {(Subscript[γ, 1]), (Subscript[γ, 2]),...,(Subscript[γ, n])}, it means (Subscript[γ, i]) can be different for every (Subscript[g, i])(x), 
          see Eq.(11). If start system is generated for polynomials all  (Subscript[γ, i]) = 1.
 p  -   indicating numerical precision., p = 0 standard precison, p =1 high precision computation

Output variables

sol[1]    -  list of the i-th solutions correspondig to the i-th initial values,  {(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_i, i = 1...m
sol[2]    -  list of the path of  i-th solutions in form of interpolating functions of the variables  
                 correspondig to the i-th  initial values,  {(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_i,   i = 1...m,
                 {{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_1,{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_2,...,{(Subscript[(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }, m])},  (Subscript[φ, i])=(Subscript[φ, i])(λ)

LinearHomotopyNDS01[X_,F_,G_,X0_,γ_,p_]:=Module[{H,n,Jλ,J,u,i,j,b,A,a,B,rhs,lhs,deqs,sol,pa,pp,wp},

Off[NDSolve::"precw"];H=Flatten[(1-λ) Thread[G  γ]+λ F];
n=Length[X];
Jλ=Outer[D,H,Join[X,{λ}]];
J=Take[Jλ,{1,n},{1,n}];
u=-Take[Jλ,{1,n},{n+1}];
A=Table[(Subscript[a, i,j]),{i,1,n},{j,1,n}];
B=Table[(Subscript[b, i]),{i,1,n}];
rhs=(Simplify[LinearSolve[A,B]]/.(Join[Table[(Subscript[b, i])->u[[i]],{i,1,n}],Flatten[Table[(Subscript[a, i,j])->J[[i,j]],{i,1,n},{j,1,n}]]])/.Map[#->#[λ]&,X])//Flatten;
lhs=Map[D[#[λ],{λ,1}]&,X];
deqs=MapThread[#1==#2&,{lhs,rhs}];
{pa,pp,wp}=If[p==1,{AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->30},{{},{},{}}];
sol=Map[Flatten[NDSolve[Join[deqs,MapThread[#1[0]==#2&,{X,#}]],Map[#[λ]&,X],{λ,0,1},pa,pp,wp]]&,X0];
{Map[Chop[N[(Map[#[λ]&,X]/.#/.λ->1)]]&,sol],sol}]

Now we solve the same problem. First, we use standard precision,

sol= LinearHomotopyNDS01[X8,F8,G8,X08,γ,0];

sol[[1]]

Then with high precision

sol= LinearHomotopyNDS01[X8,F8,G8,X08,γ,1];

sol1NDS=sol[[1]]

or

Chop[%,(10^(-7))]

The interpolation functions for the paths

sol2NDS=sol[[2]]

Remarks :
   1) To get explicit form for the differential equation system, the linear equation system is solved in symbolic way. Therefore, the maximum number of equation can be about  n = 6 - 8
   2) High precision solution (p = 1) would take much longer time than  standard precision one ( p = 0)

Without  computation of inverse

In this implementation we consider an new parameter in (3), namely

dH (x(t),λ(t))=((∂H)/(∂x))dx(t)+((∂H)/(∂λ))dλ(t) ≡ 0               λ(t)∈[0,1]

this means with  s(t)  = {x (t), λ(t)} we get

(Subscript[H, s]) (s(t)) (ds/dt) =0

where (Subscript[H, s]) is the Jacobian of the homotopy function respect to vector s  having n+1 dimensions. Then the i-th derivate function can be expressed as

((Subscript[ds, i])/dt)=(((-1))^(i+1)) det(((Subscript[H, s])^((1))),...,((Subscript[H, s])^((i-1))),((Subscript[H, s])^((i+1))),...((Subscript[H, s])^((n+1))))

where ((Subscript[H, s])^((i))) is the i-th column of the Jacobian matrix (Subscript[H, s]) (s(t)). Garcia and Zangwill (1981).

However,  in that case one need to check the upper bound of the integration parameter, while integration should be carried out up to (SuperStar[t]), where λ ((SuperStar[t])) = 0. see Binous (2007a).

Therefore it is reasonable to eliminate the last derivate of  vector s, namely

((Subscript[ds, n+1])/dt)=((dλ)/dt)

Considering that

((Subscript[ds, i])/dt)=((d (Subscript[x, i]))/dt)= ((d (Subscript[x, i]))/(dλ)) ((dλ)/dt)=((d (Subscript[x, i]))/(dλ))((Subscript[ds, n+1])/dt)

for i = 1, 2, ..., n.

((d (Subscript[x, i]))/(dλ))=((d (Subscript[x, i]))/dt)/((dλ)/dt) =((Subscript[ds, i])/dt)/((Subscript[ds, n+1])/dt)

and now the integration can be carried out with independent variable λ on [0, 1].

LinearHomotopyNDS02[X_,F_,G_,X0_,γ_,p_]:=Module[{H,n,Jλ,rhs,rhs0,i,lhs,deqs,sol,pa,pp,wp},
Off[NDSolve::"precw"];H=Flatten[(1-λ) Thread[G  γ]+λ F];
n=Length[X];
Jλ=Outer[D,H,Join[X,{λ}]];
rhs0=Table[(((-1))^(i+1))Det[Transpose[Drop[Transpose[Jλ],{i,i}]]],{i,1,n+1}];
rhs=Drop[Map[#/Last[rhs0]&,rhs0],{n+1,n+1}]/.Map[#->#[λ]&,X]//Flatten;
lhs=Map[D[#[λ],{λ,1}]&,X];
deqs=MapThread[#1==#2&,{lhs,rhs}];
{pa,pp,wp}=If[p==1,{AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->30},{{},{},{}}];
sol=Map[Flatten[NDSolve[Join[deqs,MapThread[#1[0]==#2&,{X,#}]],Map[#[λ]&,X],{λ,0,1},pa,pp,wp]]&,X0];
{Map[Chop[N[(Map[#[λ]&,X]/.#/.λ->1)]]&,sol],sol}]

The input and out are the same as  in case of  the function LinearHomotopyNDS01(#).

Example 9

F={x^2+y^2-1,Sin[x]-y}

FindRoot[F,{{x,0},{y,0}}]

Let us employ fixpoint homotopy

G={x-0,y-0}

X={x,y}

X0={{0,0}}

γ={1,1}

sol=LinearHomotopyNDS02[X,F,G,X0,γ,1];

sol[[1]]

Function for visualization of paths

This function plots the homotopy paths computed by either direct method or by integration.

Input variables

X   -  the list of the independent variables X = {(Subscript[x, 1]), (Subscript[x, 2]),...,(Subscript[x, n])}
sol  - list of the interpolation functions of the variables  
         {{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_1,{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_2,...,{(Subscript[(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }, m])}
X0  -  list of initial values, X0 = {{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_1,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_2,...,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_m}
         where m is the number of the roots of  F

Paths[X_,sol_,X0_]:=GraphicsArray[Table[ParametricPlot[{Re[X[[i]][λ]],Im[X[[i]][λ]]}/.sol[[j]],{λ,0,1},PlotRange->All,BaseStyle->{FontSize->10,FontFamily->"Times"},Axes->None,FrameLabel->{"Re","Im"},Frame->True,AspectRatio->0.6,PlotLabel->StringJoin[ToString[X[[i]]],"(λ) "],Epilog->{PointSize[0.02],Blue,Point[{Re[X0[[j,i]]],Im[X0[[j,i]]]}],PointSize[0.02],Red,Point[{(Re[X[[i]][λ]]/.sol[[j]])/.λ->1,(Im[X[[i]][λ]]/.sol[[j]])/.λ->1}]}],{j,1,Length[X0]},{i,1,Length[X]}]];

Let us plot the result of the solution of Example 8. above

pNDS= Paths[X8, sol[[2]], X08]

Examples

Nonlinear systems

In this section we solve  some test examples using different methods.

Nonlinear equation systems

Example 10

F={(x1^2)+(x2^2)-2,(E^((1.216x1-1)))+(x2^3)-2}

GraphicsArray[{ContourPlot[{F[[1]]==0,F[[2]]==0},{x1,-1.5,1.5},{x2,-1.5,1.5}],ContourPlot[{F[[1]]==0,F[[2]]==0},{x1,0.8,1.5},{x2,-0.5,1.}]}]

Try to find the root, which is roughly (Subscript[x, 1]) =1.4, (Subscript[x, 2]) =  -0.2. Use Newton-Raphson with initial values (Subscript[x, 1]) =1.4, (Subscript[x, 2]) =  0.

FindRoot[F,{x1,1.4},{x2,0}]

It fails. However starting from closer a point

FindRoot[F,{x1,1.35},{x2,-0.2}]

Let us employ fixpoint homotopy

G={x1-1.4,x2-0}

with

X0={{1.4,0}}

X={x1,x2}

and

γ={1+I,1-2 I};

We get the solution even with the original initial values

sol=LinearHomotopyNDS01[X,F,G,X0,γ,1];

sol[[1]]

This is a solution indeed

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

The paths of the solution

Paths[X, sol[[2]], X0]

Example 11

Let us consider the following system,

f1[x_,y_]:=x-Sin[2x+3 y]+Cos[3 x-5 y]

f2[x_,y_]:=y-Sin[x-2 y]+Cos[x+3 y]

ContourPlot[{f1[x,y]==0,f2[x,y]==0},{x,-2,2.5},{y,-2,2}]

F={f1[x,y],f2[x,y]}

FindRoot[F,{{x,1.5},{y,1.5}}]

But it is not a solution

F/.%

Let us employ fixpont homotopy

G={x-1.5,y-1.5}

X={x,y};

X0={{1.5,1.5}};

γ={1,1};

sol= LinearHomotopyFR[F,G,X,X0,γ,400];

sol[[1]]

It is a solution

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

Paths[X, sol[[2]], X0]

(x[λ]/.sol[[2]])/.λ->1

Show[{ContourPlot[{f1[x,y]==0,f2[x,y]==0},{x,0,2.1},{y,-0.5,1.8},Frame->True,FrameLabel->{"x", "y"}],ParametricPlot[{x[λ],y[λ]}/.sol[[2]]/.λ->u,{u,0,1}]},Epilog->{PointSize[0.016],Blue,Point[{1.5,1.5}],PointSize[0.016],Red,Point[{First[(x[λ]/.sol[[2]])/.λ->1],First[(y[λ]/.sol[[2]])/.λ->1]}]}]

Using an other starting point

FindRoot[F,{{x,-1},{y,-1}}]

G={x+1,y+1}

X={x,y};

X0={{-1,-1}};

γ={1+I,1-2 I};

sol= LinearHomotopyFR[F,G,X,X0,γ,200];

sol[[1]]

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

Paths[X, sol[[2]], X0]

Polynomial systems

Example 12

let us consider the Davenport's  polynomial system, in case of  n = 4.

F={x1+x2+x3+x4,x1 x2+x2 x3+x3 x4+x4 x1, x1 x2 x3+x2 x3 x4+x3 x4 x1+x4 x1 x2,x1 x2 x3 x4-1}

Try to solve with Newton - Raphson method

FindRoot[F,{{x1,1},{x2,1},{x3,1},{x4,1}}]

In order to solve the problem with homotopy, we construct a start system

Now we try to solve a not well formed polynomial system,

X={x1,x2,x3,x4};

{G,X0}=StartingSystem[F,X];

where is the start system

G

and X0 is the set of start values

X0

Let us employ the direct method first,

γ={1,1,1,1};

sol= LinearHomotopyFR[F,G,X,X0,γ,100];

sol[[1]]

Chop[%,(10^(-7))]

this a complex solution

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

Paths[X, sol[[2]], X0]

However we can get an other,  real solution, too,

sol[[1]](I)

Chop[%,(10^(-7))]

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

or in more general form

sol[[1]] (I^n)//ComplexExpand

Chop[%,(10^(-7))]

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

It means we have infinite solutions. Indeed

NSolve[F,X]

Remark

This problem can be solved in symbolic way with Dixon resultant too,

DixonResultant[F,{x2,x3,x4},{α,β,δ}]

Solve[%==0,x1]

However Gröbner basis does not work

GroebnerBasis[F,X,{x2,x3,x4}]

In addition NSolve can detect other solutions , too

NSolve[F,X]

for example

F/.{x1->-0.7531888931634305`,x2->-1.3276881922673471`,x3->0.7531888931634305`,x4->1.3276881922673471`}

Let us consider an other problem.

Example 13

This  is a not well formed polynomial system, a system of radicals

F={(Sqrt[x])+y z-33,(x^2)+(y^2)+(z^2)-81,((x y)^(1/3))+(Sqrt[z])-4};

Let the start system as

G={(Sqrt[x])-33,(y^2)-81,(Sqrt[z])-4};

The solutions of this start system are {1089, - 9, 16} and {1089, 9, 16}. Therefore we have two start value sets,

X0 ={ {1089, -9, 16},{1089, 9, 16}};

The list of variables is

X={x,y,z};

The list of complex numbers is,

γ={(1+I),(1-I),(1+ 2I)};

First let us compute the zeros with standard precision, using integration method

sol=LinearHomotopyNDS01[X,F,G,X0,γ,0];

sol[[1]]

With high precision the computation takes much longer,

sol=LinearHomotopyNDS01[X,F,G,X0,γ,1];

sol[[1]]

So we have got a real and one complex root.

Let us display the homotopy paths for the first and second initial values

Paths[X,sol[[2]],X0]

Remarks

Reduced Gröbner Basis gives but a wrong result,

grbx=GroebnerBasis[F,X,{y,z}];

while we have (Sqrt[z]) in this univariate polynomial of  x

Coefficient[grbx,(Sqrt[z])]//N

Dixon Resultant does not work at all

DixonResultant[F,{y,z},{α,β}]

NSolve works properly and it is fast

NSolve[F,{x,y,z}]

We also may visualize the two real solutions

Show[{ContourPlot3D[F[[1]]==0,{x,0,3},{y,3,10},{z,3,10},AxesLabel->Automatic],ContourPlot3D[F[[2]]==0,{x,0,3},{y,3,10},{z,3,10},AxesLabel->Automatic],ContourPlot3D[F[[3]]==0,{x,0,3},{y,3,10},{z,3,10},AxesLabel->Automatic],Graphics3D[{Red,Sphere[{1,8,4},0.1]}],Graphics3D[{Red,Sphere[{0.39,4.03,8.04},0.1]}]},ViewPoint->{2.8766,1.25936,1.260619},ViewVertical->{0.25881,0.96541459,0.031488}]

Even Newton - Raphson can work and provides all of the roots with the proper start values,

FindRoot[F,{{x,1089},{y,9},{z,16}}]

FindRoot[F,{{x,1089},{y,-9},{z,16}}]//Chop

FindRoot[F,{{x,1},{y,2},{z,3}}]//Chop

FindRoot[F,{{x,1089 -I},{y,9+I},{z,16-I}}]//Chop

However using these start values, employing Newton homotopy with direct solution, we can get also the other roots, too,

f1[x_,y_,z_]:=F[[1]];

f2[x_,y_,z_]:=F[[2]];

f3[x_,y_,z_]:=F[[3]];

The other real solution is

G={f1[x,y,z]-f1[1,2,3],f2[x,y,z]-f2[1,2,3],f3[x,y,z]-f3[1,2,3]};

X0={{1,2,3}};

sol=LinearHomotopyFR[F,G,X,X0,γ,100];

sol[[1]]

The other complex solution is

G={f1[x,y,z]-f1[1089-I,9+I,16-I],f2[x,y,z]-f2[1089-I,9+I,16-I],f3[x,y,z]-f3[1089-I,9+I,16-I]};

X0={{1089-I,9+I,16-I}};

sol=LinearHomotopyFR[F,G,X,X0,γ,100];

sol[[1]]

Example 14

Our last polynomial example is the  Reimer 5 system

f1[x_,y_,z_,t_,u_]:=2 (x^2) -2 (y^2) +2 (z^2)-2 (t^2)+2 (u^2)-1

f2[x_,y_,z_,t_,u_]:=2 (x^3) -2 (y^3) +2 (z^3)-2 (t^3)+2 (u^3)-1

f3[x_,y_,z_,t_,u_]:=2 (x^4) -2 (y^4) +2 (z^4)-2 (t^4)+2 (u^4)-1

f4[x_,y_,z_,t_,u_]:=2 (x^5) -2 (y^5) +2 (z^5)-2 (t^5)+2 (u^5)-1

f5[x_,y_,z_,t_,u_]:=2 (x^6) -2 (y^6) +2 (z^6)-2 (t^6)+2 (u^6)-1

F={f1[x,y,z,t,u],f2[x,y,z,t,u],f3[x,y,z,t,u],f4[x,y,z,t,u],f5[x,y,z,t,u]}

X={x,y,z,t,u};

Let us consider the initial values for the Newton Raphson method as

{x0,y0,z0,t0,u0}={9.,-2.,-5.,1,-1.}

FindRoot[F,{x,x0},{y,y0},{z,z0},{t,t0},{u,u0}]

Let us employ fix point homotopy

G={x-x0,y-y0,z-z0,t-t0,u-u0}

X0={{x0,y0,z0,t0,u0}}

γ={1-I, 1+I,1-2 I,-1-2I,-1-I};

sol= LinearHomotopyFR[F,G,X,X0,γ,500];

sol[[1]]

this a complex solution

F/.MapThread[#1->#2&,{X,Flatten[sol[[1]]]}]//Chop

Paths[X, sol[[2]], X0]

Fig .26 A solution of the polynomial system Reimer 5

NSolve (#) gives all of the roots

NSolve[F,X]//TableForm

References

Binous H. (2007a) Homotopy Continuation Method to Solve A   System of Nonlinear Algebraic Equations.
    http://library.wolfram.com/infocenter/MathSource/6710/
    
    
Binous H. (2007b) Homotopy Continuation Method  to Find All Root of a Polynomial Equation I - II.
   http://library.wolfram.com/infocenter/MathSource/6717/
   
   
Decarolis F., Mayer R. and Santamaria M.  Homotopy Continuation Methods. http://home.uchicago.edu/~fdc/H-topy.pdf


Drexler, F.J. (1977) Eine methode zur Berechnung sämtlicher Lösungen von Polynomgleichungssytemen, Numer. Math. 29. pp.45-58.


Garcia, C.B. and Zangwill W.I. (1979) Determining all solutions to certain systems of nonlinear equations, Mathematics of Operation Research 4. pp. 1-14.

Garcia C.B. and W.I. Zangwill (1981) Pathways to Solutions, Fixed Points and Equilibrai. Prentice Hall, Englewood Cliffs. NJ.  

   
Gritton K.S., Seader J.D. and Wen-Jing Lin (2001) Global homotopy continuation procedures for seeking all roots of a nonlinear equation.
   Computers and Chem. Engn. 25. pp.1003-1019
   
   
Hazaveh K, Jeffrey D.J, Reid G.J, Watt S.M and Wittkopf A.D (2003) An exploration of homotopy solving in Maple. 	
    http://www.apmath.uwo.ca/~djeffreyOffprints/ascm2003pdf   
    
Kotsireas, I. S. (2001) Homotopies and polynomial systems solving I: Basic Principles, ACM SIGSAM Bulletin , 35.,No. 1. pp.19-32
      
      
Sommese A.J., Verschelde J. and Wampler C.W. (2003) Introduction to Numerical Algebraic Geometry. 
  http://mate.dm.uba.ar/~visita16/cimpa/notes/Verschelde.pdf     
   
  
Verschelde J. (2007) Homotopy Methods for Solving Polynomial Systems tutorial at ISSAC'05, Beijing, China.
   http://www.math.nic.edu/~jan/tutorial.pdf