Implementation of Poincaré–Miranda theorem

Question

To test whether a continuous function has a root one simple root in a given interval [x0, x1] is relatively easy: according to Intermediate value theorem when the sign of function value at x0 is opposite to the sign at x1, there is (at least) one root.

For example, given a quadratic function:

g(x): a*x**2 + b*x + c = 0

The test looks like:

if sign of g(x0) is opposite of sign of g(x1)
then return true
else return false

For multivariate case there is Poincaré–Miranda theorem but I have a bit of difficulty to implement the test correctly from reading the linked article.

Given two quadratic bivariate functions:

g1(x, y): a1*x**2 + b1*y**2 + c1*x*y + d1*x + e1*y + f1 = 0
g2(x, y): a2*x**2 + b2*y**2 + c2*x*y + d2*x + e2*y + f2 = 0

and a rectangular region [x0, x1] x [y0, y1], how to check if there is at least one root in the region?

I mean, I assume the test should look somewhat like (this doesn't work):

if (sign of g1(x0, y0) is opposite of sign of g1(x1, y0) and
    sign of g1(x0, y1) is opposite of sign of g1(x0, y1) and
    sign of g2(x0, y0) is opposite of sign of g2(x1, y0) and
    sign of g2(x0, y1) is opposite of sign of g2(x0, y1))
then return true
else return false

Please, does anyone know what pairs of functions, interval end-points, and logical operators to check and in what order?

Answer 1

You need to check whether your bivariate functions

g1(x, y): a1*x**2 + b1*y**2 + c1*x*y + d1*x + e1*y + f1 = 0
g2(x, y): a2*x**2 + b2*y**2 + c2*x*y + d2*x + e2*y + f2 = 0

satisfy

I).  g1(x0,y) < 0, for all y in [y0,y1]
II). g2(x,y0) < 0, for all x in [x0,x1]

and

III). g1(x1,y) > 0, for all y in [y0,y1]
IV).  g2(x,y1) > 0, for all x in [x0,x1]

Your functions are quadratic, so this can be done without sampling values along all 4 boundaries for the 4 cases. For example, for the first condition on g1(x0,y), simply plug in x0 for x, obtaining a quadratic equation in y:

G1(y) = b1*y**2 + c1*x0*y + e1*y + (f1 + d1*x0 + a1*x0**2)

We need to check whether G1 is ever positive for y in [y0,y1]. Since G1 is quadratic, its maximum either occurs where {G1' = 0, G1'' < 0} or at the endpoints. So:

a. express G1' analytically, use simple bisection to find a root in [y0,y1]
b. if there is one, say y*, express G1'' analytically and compute G1''(y*)
c. if G1''(y*) is also < 0 then you have your maximum y*
d. if G1(y*) > 0 then the condition are violated, you may break
e. if not, then test the endpoints G1(y0), G1(y1).
f. if any of those are > 0 then break.

If your function passes these tests without breaking, you have satisfied the first of the 4 conditions (I) above.

The conditions (II-IV) can be tested in a similar manner. If all conditions hold, the Miranda test holds and you have a coincident root of the two functions. If not, then you are in the "maybe" case - the functions may still have a common root, but you'll have to use a different method to prove existence.

Answer 2

First of all, your original "intermediate value"-based code doesn't quite do what it is advertised to do:

To test whether a continuous function has a root in a given interval [x0, x1] is relatively easy: according to Intermediate value theorem when the sign of function value at x0 is opposite to the sign at x1, there is (at least) one root.

The test looks like:

if sign of g(x0) is opposite of sign of g(x1)
then return true
else return false

This "test" has one-sided error, as pointed out by David Eisenstat. If the signs are indeed opposite, then return true is okay, but if the signs are not opposite, then return false should perhaps be return maybe or something...

---

Second of all as regards the Poincare Miranda theorem, in higher dimensions comparing the signs of a few points doesn't give you enough information to apply the theorem.

Consider n continuous functions of n variables. Assume that for every variable x_i, the function f_i is constantly negative when x_i = 0 and constantly positive when x_i = 1. Then there is a point in the n-dimensional unit cube in which all functions are simultaneously equal to 0.

There is no black-box test if a continuous function is "constantly negative" on some region.

You would need to assume something more, like, you assume it is actually a low-degree polynomial and you sample it at enough points to discover its coefficients, etc.

---

If we assume like you stated that we have two bivariate quadratics, and we actually have (or deduce) the coefficients... it is possible.

What I would do is simply, substitute in the value for x_i in each function as required, so it reduces to a univariate quadratic, and then solve for its roots (if any) using the quadratic formula like we learned in grade school. Then check whether they occur in the region of interest. Then test a point in-between the roots to determine the sign. Then you'll know if the theorem can be applied.

It's possible that you can solve for a precise condition in closed-form but I'm not sure if this will actually help you write a better (simpler / more efficient) implementation.

Here is some pseudocode:

 def quadratic_positive_in_region(p, x_0, x_1)
   ASSERT(p is univariate)
   ASSERT(x_0 <= x_1)

   // If one of the roots lies in the region then
   // we are zero there, and thus not positive
   def roots = quadratic_formula(p)
   for r in roots:
     if x_0 <= r and r <= x_1 then return false

   // If there are no roots in the region then
   // we are either always positive or always negative,
   // so test a single point to determine.
   if p(x_0 + x_1 / 2) > 0 then return true       
   return false

 def poincare_miranda(g1, g2, x_0, x_1, y_0, y_1)
   return quadratic_positive_in_region(-g1 | y = y_0, x_0, x_1) and
          quadratic_positive_in_region( g1 | y = y_1, x_0, x_1) and
          quadratic_positive_in_region(-g2 | x = x_0, y_0, y_1) and
          quadratic_positive_in_region( g2 | x = x_1, y_0, y_1)

 def generalized_test(g1, g2, x_0, x_1, y_0, y_1)
   return poincare_miranda( g1,  g2, x_0, x_1, y_0, y_1) or
          poincare_miranda(-g1,  g2, x_0, x_1, y_0, y_1) or
          poincare_miranda(-g1, -g2, x_0, x_1, y_0, y_1) or
          poincare_miranda( g1, -g2, x_0, x_1, y_0, y_1)

I'm using some notation here where - operator can be applied to a polynomial, also | notation represents substitution of a value for a variable in a polynomial.