Bounding a bicubic polynomial My actual situation is working with bicubic polynomials, (that is, polynomials of the form $\sum_{i=0}^3 \sum_{j=0}^3 a_{i,j} x^iy^j$) defined on the unit square $[0,1]^2$ (actually these are splines).  I also have simple "boxes" of the form $[x_1,x_2]\times [y_1,y_2]\times [z_1,z_2]$ for reasonable fixed real numbers there.  The question is, how can I efficiently determine, using numerical methods, whether the graph of the polynomial touches my box?
I don't need a perfect analytic solution (and, I suspect, none exists for general $a_{i,j}$), but I need one that can work "arbitrarily well" or however numerical analysts like to put it.
One method I have considered is to reduce the problem to this: if the graph is bounded by a horizontal plane, say $z=c$ for some fixed $c$, on the range of parameters $[x_1,x_2]\times [y_1,y_2]$.  There are well-known numerical methods to do this in one variable, so checking boundary solutions becomes easy, but checking for interior extrema seems intractable (setting the derivatives to zero, then solving, is not feasible in general unless I'm missing something). 
Another method is to bound derivatives, evaluating the function at a middle point, and using the easy linear bounds you get from that and MVT.  But this is not very accurate, so you get a lot of false positives and zooming in when the coefficients $a_{i,j}$ are large or of different signs, so the algorithm doesn't run very quickly.  One solution would be to get better bounds on the derivatives, but if I could do that, I could probably solve the whole problem.
So, answers to any of the following would solve the problem:


*

*Is there a well-known solution to this specific problem (bounding a bicubic polynomial on a range)?

*Is there a general technique for finding (or even detecting) interior maxima (which is efficient)?

*Is there a general method, which is at all pertinent, beyond what I've already said?


I don't know much about numerical analysis, so I just assume they've sort of solved this problem before.  So just a reference to something classical would be great (if it exists).
Please and thank you :)
 A: Convert the polynomial from "power" form to Bezier-Bernstein form. Then, the convex hull property of Bezier curves and surfaces gives you bounds on the location of the graph. Differentiating gives you some new polynomials in Bezier-Bernstein form, so you can then get bounds on derivatives, too, if you need them.
Or, if your splines are already in b-spline form, you can use the b-spline coefficients (rather than Bezier-Bernsteain coefficients) in the same way.
The bounds computed this way are not very tight, but they are good enough for quick eliminination of many non-touching cases.
If you really want tight bounds, you'll have to use numerical methods to compute maximum and minimum values of your functions. Possibly there are some shortcuts, because your functions are just cubic, so derivatives will be quadratics.
There's also a brute force approach. Just calculate thousands of function values and take the min/max of those values. There are computationally efficient ways to do this using forward differencing. Or, if you know how to write GPU code, the computations can be done very quickly in parallel.
The middle ground is to use the brute force approach on a medium-density grid, and then use the max/mins you find as starting points for iterative numerical methods. Since a bicubic can't oscillate many times, a $5 \times 5$ grid on each bicubic patch is probably about right. Evaluate your function at each of the 25 grid points, and remember the $(x,y)$ value that produced the maximum function value. Use this $(x,y)$ value as a starting point, and iterate from there to a nearby (local) maximum. For the iteration, any decent optimization code will suffice. You can find a large selection of good ones at the NEOS site. Maybe start with this page if you're not familiar with the field. In the past, I've had success with the donlp2 package. But it seems like overkill for such a simple problem, and I'd guess you can find something less sophisticated that will work fine. For example, the "Numerical Recipes" book has some fairly simple codes that will probably work for you.
