3
$\begingroup$

I'm writing a program where I have a 3-dimensional polynomial (3 variables) of which I have to check if it is positive in a given product of 3 intervals (a volume).

I found a paper which solved this problem but I didn't really understand it, nor am I able to extract an algorithm out of it: http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1084864

So can this approach be formulated as an algorithm for a computer?

Thanks in advance.

Edit:

To make things more specific.

My polynomial is always:

z^3*(z*(6*z-15)+10)(-y^3(y*(6*y-15)+10)(x^3(x*(6*x-15)+10)*(g3z*z-g2z*z+g3y*(y-1)-g2y*(y-1)-g2x*x+g3x*(x-1))-x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))+g2z*z-g0z*z-g0y*y+g2y*(y-1)+g2x*x-g0x*x)-x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))-g0z*z+g4z*(z-1)+y^3*(y*(6*y-15)+10)(g6z(z-1)-g4z*(z-1)+x^3*(x*(6*x-15)+10)(g7z(z-1)-g6z*(z-1)+g7y*(y-1)-g6y*(y-1)-g6x*x+g7x*(x-1))-x^3*(x*(6*x-15)+10)(g5z(z-1)-g4z*(z-1)+g5y*y-g4y*y-g4x*x+g5x*(x-1))-g4y*y+g6y*(y-1)+g6x*x-g4x*x)+x^3*(x*(6*x-15)+10)(g5z(z-1)-g4z*(z-1)+g5y*y-g4y*y-g4x*x+g5x*(x-1))+g4y*y-g0y*y+g4x*x-g0x*x)+y^3*(y*(6*y-15)+10)(x^3(x*(6*x-15)+10)*(g3z*z-g2z*z+g3y*(y-1)-g2y*(y-1)-g2x*x+g3x*(x-1))-x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))+g2z*z-g0z*z-g0y*y+g2y*(y-1)+g2x*x-g0x*x)+x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))+g0z*z+g0y*y+g0x*x

where all variables starting with $g$ are constant. The given intervals for $x,y,z$ are limited to $[0,1]$.

  • 4
    Interesting -- "n-dimensional interval" does turn up about 50,000 Google hits -- I'd never seen that expression before.2011-07-10

2 Answers 2

7

The paper you cite contains the following theorem:

$ P(x) > 0, \ \ \forall x \in D^n \subseteq R^n $

where $P(x)$ is a real polynomial, iff:

1) All the $n-1$-dimensional polynomials arrived at by substituting the end points of the intervals of $x_i$ $(i=1,2,\cdots,n)$ in $P(x)$ (each one at a time), are positive in $D^{n-1}$,

2) The set of $n$ equations in $n$ real variables $ \eqalign{P(x) &= 0 \cr \frac{\partial P}{\partial x_i} &= 0, \ \ \forall i \in \{1,2,\cdots,n-1\}\cr}$ has no solution in $D^n$.

I would program this in Maple (for the case of finite intervals) as follows:

ispositive:= proc(P::polynom, R::list(name=range)) local x,D,n,i,S,s,Q; x:= map(lhs,R); D:= map(convert,map(rhs,R),list); if D = [] then return is(P > 0) end if; n:= nops(D); for i from 1 to n do   if not (ispositive(subs(x[i]=D[i][1], P),subsop(i=NULL,R))     and ispositive(subs(x[i]=D[i][2], P), subsop(i=NULL,R))) then return false  end if end do; S:= RootFinding[Isolate]({P, seq(diff(P,x[i]),i=1..n-1)},x,output=interval);  for s in S do   Q:= map(rhs,s);   if `or`(seq(Q[i][1] > D[i][2],i=1..n),seq(Q[i][2] = D[i][1],i=1..n),seq(Q[i][2]<=D[i][2],i=1..n))          then return false    else error "Roots returned with too little precision, try increasing Digits"   end if; end do; true end proc; 

For example,

ispositive(16*x^2 + 16*y^2 - 64*x - 16*y + 68,[x=1..10, y=-1..1]);

false

  • 0
    Yes, so this may not work in some "degenerate" cases. For example, ispositive((x^2+y^2-1)^2,[x=-2..2,y=-2..2]);2011-07-12
2

I solved this by converting the polynomial to a bezier volume with regular control points and splitting it to lie in my asked intervals. Then I can use the fact that the curve lies completely inside the convex hull of its control points.

  • 1
    This is a very good approach to the problem. In general you can use Interval Arithmetic (along with some subdivision) to solve this sort of thing, but the right form for the polynomial can make a huge difference in the results. The convex-hull property of Bezier surfaces that you mention is a perfect example of this.2011-08-24