Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

share|improve this question
    
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. –  joriki Jul 10 '11 at 22:34
    
You might want to look up Sturm sequences. –  Yuval Filmus Jul 10 '11 at 22:40
    
I'm talking about a 3-dimensional polynomial (means 3 variables) and a 3-dimensional interval which is a volume (en.wikipedia.org/wiki/…). The degree is about 6. –  man Jul 10 '11 at 22:44
    
@man: "3-dimensional interval" is not a standard use of the word "interval." If you mean a product of intervals, say so. –  Qiaochu Yuan Jul 10 '11 at 22:50
4  
Interesting -- "n-dimensional interval" does turn up about 50,000 Google hits -- I'd never seen that expression before. –  joriki Jul 10 '11 at 23:00

2 Answers 2

up vote 7 down vote accepted

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)) then next end if;
  if `and`(seq(Q[i][1] >= 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

share|improve this answer
    
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. –  man Jul 11 '11 at 20:10
    
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]); –  Robert Israel Jul 12 '11 at 0:31

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.

share|improve this answer
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. –  Steven Stadnicki Aug 24 '11 at 17:08

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.