discriminant.html

MCS 595 Thursday 4 September 2003

> restart;

This very simple worksheet illustrates the main theorem of elimination theory applied to justify numerical homotopy continuation methods. For a simple homotopy, we set up the discriminant system which defines all singular points. From this system we eliminate the orginal independent variables to arrive at one polynomial in the continuation parameter t. This shows there are only finitely many singular values for the continuation parameter t. By random choice of complex coefficients in the start system we can ensure that all roots of this discriminant polynomial will miss the interval [0,1].

Consider the intersection of two simple quadratic equations:

> f := [x^2 + 4*y^2 - 4,2*y^2- x];

By the "method of degeneration", we replace each equation in f by an equation of the same degree, but much simpler:

> g := [x^2 - 1, y^2 - 1];

So g is the start system in the homotopy:

> h := t*f + (1-t)*g;

The continuation parameter t will move from 0 to 1. At t=0, we have the start system, at t=1, we have the target system (the system we want to solve).

For further manipulations in Maple, we will expand the homotopy:

> eh := expand(h);

We are interested in the singular points of the homotopy. Therefore, we compute the Jacobian matrix of the homotopy:

> jh := matrix(2,2,[[diff(eh[1],x),diff(eh[1],y)],[diff(eh[2],x),diff(eh[2],y)]]);

Then all singular points of the homotopy are defined by the following polynomial system:

> sys := [eh[1],eh[2],linalg[det](jh)];

With a pure lexicographical Groebner basis we can eliminate the variables x and y, to arrive at a polynomial in t only:

> gb := grobner[gbasis](sys,[x,y,t],plex);

We find this polynomial in t as the last polynomial of our lexicographical Groebner basis:

> gb[nops(gb)];

As the degree of this "discriminant polynomial" is seven, we can expect seven complex roots:

> fsolve(gb[nops(gb)],t,complex);

We should be particularily troubled by the root around 0.4. This means that as t moves from 0 to 1, we will encounter a singularity.

We will choose a random constant, which we call gamma. The instruction which is commented out shows how we should really do it, but in view of Maple's "grobner", we just choose a much simpler gamma.

> unprotect(gamma):

> #gamma := exp(stats[random,uniform[0,2*Pi]](1)*I);

> gamma := 1+I;