I don't think you can ensure that with gradients though: if f and f' are computed in machine arithmetic, cancellation errors might pile up.
One thing to bear in mind, though: there are numerical algorithms which are stable and unimpeachably useful but which are unstable in interval arithmetic owing to the natural interval inclusion behaving poorly. A great example of this is the Clenshaw recursion for evaluating Chebyshev series. It is the method of choice for evaluating a Chebyshev series at a single point, but the corresponding natural interval inclusion is actually unstable itself (the width tends to blow up exponentially fast in the degree of the series, regardless of the coefficients).
I tend to view "classical numerical methods" and interval arithmetic as somewhat different camps. There's definitely a significant overlap and they can both profit from each other, but you need to put on a different hat when you're doing each. Sometimes interval arithmetic is just too expensive, in which case you need to select a good, stable numerical algorithm. On the other hand, if you're doing interval arithmetic because you must (you need a certified, correct answer), then you will often need to use a completely different set of algorithms and spend some time fussing over what the most appropriate interval inclusion is.
Again, another good example here relates to Chebyshev series. If you want to find all of the roots of an analytic function on an interval, expand it in a Chebyshev series and solve the colleague matrix eigenvalue problem. The result could be subtly wrong due to numerical roundoff, and you'd have a hard time recovering from this position; but this method is undeniably extremely useful. On the other hand, if you want to compute rigorous enclosures of all zeros of that same analytic function, you'd need to develop a suitable interval inclusion for it, and then run some type of interval rootfinding algorithm based on a combination of interval contractors (Krawczyk) and subdivision. Basically, a branch and bound algorithm.
Does modal interval analysis help with these cases? I still don't understand it, but there are a lot of things I don't understand in numerical methods.
Yes, you can. You need to do automatic differentiation with interval arithmetic. This gives you a mathematically correct result for the derivative.
Always keep in mind with interval arithmetic that (-inf, inf) is also a mathematically correct result. Keeping intervals small is the challenge with interval arithmetic and means that every algorithm not specifically developed for it is likely going to lead to useless results.
1) compute f with interval arithmetic
2) compute f normally and f' with interval arithmetic
3) compute f rounding towards zero, compute f' from f rounded towards infinity, and round f' up (if f positive) or round f' down (if f negative).
In all 3 cases you can use what you computed to figure out bounds on f, (1) is direct, the other two need extra work.
But that kinda defeats the point of replacing interval arithmetic with gradients, though.