Computing maximal value of a function with constraints using software - how to make the computation rigorous?

I have a fairly messy rational function $$D$$ in variables $$u_1,u_2,u_3,q$$. I am trying to compute the maximum of this function on the compact set $$K$$ defined by $$u_1^2+u_2^2+u_3^2=1, 0 \leq q \leq e^{-2\pi}$$. I would like to be able to compute this maximum rigorously, or at least be able to find a reasonable constant $$C$$ and give a proof that $$D \leq C$$ on $$K$$. As the function is messy, and this is only the $$n=3$$ case, and I would like to generalise this to larger $$n$$, I don't want to do this by hand.

Mathematica has code which claims to find the maximum of the the function exactly, e.g. the function Maximise. This seems to run endlessly in my case. It also has a function NMaximise, which computes a numerical approximation to the maximum, which gives an answer almost immediately, which I am certain is correct.

What I want to know is how to convert this certainty into a rigorous statement, one I could state in a paper on number theory. E.g. can I obtain an error bound on the number which Mathematica returns? Or is there a different software that will do this for me?

Any help would be appreciated. I'm not quite sure what to google, this problem seems to go under the name of deterministic global optimization, but most applications are understandably in industry, not pure mathematics.

If it helps, $$D$$ is a homogenous polynomial of degree 6 in $$u_1,u_2,u_3$$, and its coefficients are rational functions in $$q$$.