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]$.

  • 0
    Do you really mean a $3$-dimensional polynomial, or a polynomial of degree $3$? If by a $3$-dimensional polynomial you mean a polynomial containing $3$ variables, then it doesn't make sense to talk of one given interval. That article doesn't seem to be freely available to non-members.2011-07-10
  • 0
    You might want to look up Sturm sequences.2011-07-10
  • 0
    I'm talking about a 3-dimensional polynomial (means 3 variables) and a 3-dimensional interval which is a volume (http://en.wikipedia.org/wiki/Interval_%28mathematics%29#Multi-dimensional_intervals). The degree is about 6.2011-07-10
  • 0
    @man: "3-dimensional interval" is not a standard use of the word "interval." If you mean a product of intervals, say so.2011-07-10
  • 0
    I hopefully clarified it now in the question.2011-07-10
  • 0
    Sturm sequences looks promising, I'll do some research...2011-07-10
  • 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
    Wow, thank you. But what still confuses me is that you are using RootFinding[Isolate] which only works, if there is a finite number of solutions.2011-07-11
  • 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