2D Bisection Method - Roots Finding

Question

I'm trying to use a Bisection Method to solve two highly nonlinear equations. Let us say; f(x,y) = 0 with degree eight and g(x,y) = 0 with degree six;

I need a matlab code for 2D Bisection Method to solve f(x,y) = 0 and g(x,y) = 0 and find all possible roots.

Answer 1

Your idea does not work. With two polynomials in two variables you entered the area of algebra with Gröbner bases, resultants, Bezout theorems and solutions at infinity.

The Bezout theorem tells you that there are at most 6*8=48 solutions, real and complex. I do not think that matlab has a Gröbner package, use Maple, Magma or Singular for that.

With resultants one can reduce the number of variables. The resultant for x of f and g is a univariate polynomial in y that then can be solved by standard methods. It is the determinant of the Sylvester matrix containing polynomials in y, in principle one should be able to evaluate this in matlab, but it's not very practical.

---

The idea of resultants lends itself (over not that obvious detours) to homotopy methods. This was extensively done for fully numerical methods by Verschelde, Wampler et al.

If I remember correctly, the initial solvable problem for the homotopy is constructed by considering

g0(x,y)=(y-a1*x-b1)*...*(y-a6*x-b6) 

with random coefficients a1,..,a6,b1,...,b6. Then the initial solutions can be determined for every linear factor of g0 from the univariate polynomials

0=f1(x)=f(x,a1*x+b1),..., 0=f6(x)=f(x,a6*x+b6)

using Jenkins-Traub or Laguerre, giving roots xjk and yjk=aj*xjk+bj. Bisection or quadrisection of the complex plane is not very helpful, but exists. See the work of Yakoubsohn and Didieu.

Now introduce a homotopy parameter t going from 0 to 1 in a straight line or a curve t=s+c*s*(1-s), s in [0,1], c random small imaginary, in the complex plane and consider the systems

0=f(x(t),y(t)),
0=t*g(x(t),y(t))+(1-t)*g0(x(t),y(t)),

starting with x(0)=xjk and y(0)=yjk for all j=1,...,6, k=1,...,8.

For general values of the coefficients, the paths that are followed will not intersect, thus are regular from start to end, and all roots can be found. One tricky part is to decide that when a path drifts away if it is because of a very large solution or a non-solution at infinity.

---

Another homotopy starts with computing the roots x1,...,x8 of f(x,y0) for a fixed arbitrary y0 and then the roots y11,...,y61, y12,...,y68 of g(x1,y),...,g(x8,y). The homotopy is then given as

0=f( x(t), (1-t)*y0+t*y(t)), 
0=g( x(t), y(t)),

where x(0)=xk and y(0)=yjk.