Hessian gives a worse approximation of a multivariate function I have a real, smooth, multivariate (with 10 variables or many more) function, for which I have the exact Jacobian and Hessian. It turns out that unless the norm of the increment of the function is very small the linear approximation (just with the Jacobian) is usually better than the quadratic one (with also the Hessian).
I guess this is due to the radius of convergence of the Taylor series being very small (although in several variables it is a domain rather than a radius of convergence). 
I'm using this for trust-region optimization. In the books, however, they don't mention this possible issue. On the other hand, one rarely has the exact Hessian. I do but cannot benefit from it. Anybody could provide some thoughts or links? For instance I'd like to estimate a good error bound for the quadratic approximation. 
Thank you very much in advance for your time.

If $f$ is $C^2$ and the Hessian eigenvalues are contained in $[\underline{\lambda},\overline{\lambda}]$ then you have the bounds
$$f(x) +\frac{\partial f(x)}{\partial x}h+ \frac{1}{2} \underline{\lambda} \|h|^2 \leq f(x+h) \leq f(x) +\frac{\partial f(x)}{\partial x}h+ \frac{1}{2} \overline{\lambda} \|h|^2$$
It is hard to make broad generalizations about gradient vs. Hessian unless one knows more about the function. Many optimization schemes use gradients when 'far' from a solution (measured in some way) and transition to incorporating Hessian information (or approximations such as, eg, BFGS) when approaching a presumed solution. 
Using Hessian information appropriately can drastically reduce the number of iterations required. Also, Hessians (or approximations) can provide very good scaling.
To some extent, one of the points of using trust region methods was to take advantage of the Hessian information, particularly along directions corresponding to negative eigenvalues. The trust region size is adjusted to satisfy some 'sufficient decrease' criterion.
