A possibly new example: a ball bouncing on a parabola
($=$ graph of a quadratic polynomial) is never chaotic.
The associated dynamical system has an extra invariant
that makes it "integrable" (if I have the terminology correctly).
The general orbit is an elliptic curve $\cal E$ equipped with a point $P$
such that going from one bounce to the next corresponds to translation by $P$
on $\cal E$. The general orbit looks like one of these two pictures,
depending on whether the curve has one or two real components:
(source)
$\phantom\infty$
(source)
In the two-component case, $P$ is on the non-identity component.
So in either case the orbit is dense in $\cal E$, unless $P$ is torsion in which case
the orbit is finite. The envelope curves are themselves parts of two
parabolas, which meet at right angles in the one-component case.
In the transition between the two cases the lower parabola narrows to
a vertical line that meets the upper parabola at its center,
as suggested by the next two pictures on either side of the transition:
(source)
$\phantom\infty$
(source)
At the transition point itself the ball quickly approaches
a fixed point where it bounces vertically on the parabola's vertex.
There are two other limiting cases.
In one direction, the ball bounces very close to the surface,
and in the limit rolls on it:
$\phantom{centeringkludge}$
(source)
Here $\cal E$ remains smooth but $P$ approaches $0$.
In the other limiting case, $\cal E$ curve degenerates to two points
related by translation by $P$; the second picture shows the one-parameter
family of degenerate trajectories of this kind:
(source)
$\phantom\infty$
(source)
There's more to say here than is reasonable for a Mathoverflow answer,
but let me at least get to the the derivation of the
extra invariant and the elliptic curve.
The general setting can be described as follows.
Choose units so that the gravitational constant is $1$,
and assume that the ball bounces on the graph $(x, f(x)/2)$
of a differentiable function (the factor of $1/2$ will be convenient later).
We further assume that $f$ is bounded below,
or at least that $f(x) > -cx^2$ for all $c>0$,
so that the ball always keeps hitting the graph.
Then ball's path is completely described by a sequence of vectors
$v_n = (x_n, y_n, z_n)$ where the $n$-th bounce occurs at $(x_n, f(x_n)/2)$
and shoots the ball with velocity $(y_n,z_n)$, with $z_n \geq f'(x_n) y_n/2$
so that $(y_n,z_n)$ is in the half-plane
above the tangent to the graph at $(x_n, f(x_n)/2)$.
This puts the ball on the path
$(x_n + y_n t, \frac12 f(x_n) + z_n t - \frac12 t^2)$ with $t \geq 0$,
so $x_{n+1}$ is the value of $x_n + y_n t$ at the smallest positive solution of
$$
f(x_n) + 2 z_n t - t^2 = f(x_n + y_n t).
$$
At that time the ball has velocity $(y_n, z_n - t)$.
The new velocity vector $(y_{n+1}, z_{n+1})$ is then obtained by reflecting
$(-y_n, t-z_n)$ about the normal to the graph at $(x_{n+1}, f(x_{n+1})/2)$,
which is the line spanned by $(-\frac12 f'(x_{n+1}, 1)$
(orthogonal to the tangent $(1, \frac12 f'(x_{n+1})$).
We're thus iterating a map $\beta: H \rightarrow H$ where
$$
H = \{ (x_n, y_n, z_n) \in {\bf R}^3 : z_n \geq f'(x_n) y(n) / 2 \},
$$
and can factor $\beta$ as a product $\rho \sigma$ of involutions
where $\sigma$ takes $(x_n, y_n, z_n)$ to $(x_{n+1}, -y_n, t-z_n)$
and $\rho$ fixes $x$ and reflects $(y,z)$ about $(-\frac12 f'(x),1)$.
Time reversal corresponds to the inverse map $\beta^{-1} = \sigma \rho$.
Here's a picture:
$\phantom{centeringkludge}$
(source)
Black arrows are normal to the graph; blue tangent arrows $(y,z)$ coming from
$(x,f(x)/2)$ represent the point $(x,y,z)$; the two arrows coming from
the same point are related by the reflection $\rho$; and the two arrows
at either end of a parabolic segment are related by $\sigma$.
By conservation of energy, $\sigma$ preserves
$$
E(x,y,z) := f(x) + y^2 + z^2
$$
(that's why I wanted the factor of $1/2$ in $(x,f(x)/2)$);
since the collisions are perfectly elastic, $\rho$ preserves $E$ as well.
Moreover $\rho$ and $\sigma$, and thus also $\beta$,
are volume-preserving bijections. So we have a non-dissipative
dynamical system on each of the $2$-dimensional fibers of $E$.
Now suppose $f$ is quadratic.
We translate and scale our coordinates (including the time coordinate $t$)
to put the focus of the parabola at $(0,\frac12)$
while keeping the gravitational constant equal $1$.
Then $f(x) = x^2$, and $E(x,y,z)$ is simply $x^2+y^2+z^2$;
and $\frac12 f'(x)=x$, so $\beta$ moves $v = (x,y,z)$ around
the part of the sphere of radius $E^{1/2}$ on which $z \geq xy$.
The equation $f(x_n) + 2 z_n t - t^2 = f(x_n + y_n t)$ for the
time between bounces $n$ and $n+1$ is now the quadratic equation
$x_n^2 + 2 z_n t - t^2 = (x_n + y_n t)^2$. This equation has
the known solution $t=0$, so the other solution must be
a rational function of $x_n,y_n,z_n$.
Thus $\sigma$, like $\rho$, is a rational map.
We find that $t = 2(z_n - x_n y_n) / (1+y_n^2)$, and thus that
$\sigma$ takes $(x,y,z)$ to
$$
\left(
\frac{(1-y^2)x + 2yz}{1+y^2}, \phantom. -y, \phantom. \frac{(1-y^2)z-2yx}{1+y^2}
\right).
$$
This resembles the formula for $\rho$ for good reason: while
$\rho$ fixes $x$ and reflects $(y,z)$ about the line generated by $(-x,1)$,
$\sigma$ fixes $y$ and reflects $(x,z)$ about the line generated by $(-y,1)$,
then multiplies both $x$ and $y$ by $-1$.
It follows that both $\rho$ and $\sigma$ leave $z-xy$ invariant,
as may be seen either by direct computation or by noting that
$z-xy$ can be written as an inner product in two ways
$$
z - xy = \langle (y,z), (-x,1) \rangle = \langle (x,z), (-y,1) \rangle
$$
and that reflecting one vector about the span of another does not change
their inner product. Therefore $\beta$ also preserves $z-yx$.
It follows that the vectors $(x_n,y_n,z_n)$ are all on the intersection
of a sphere $x^2+y^2+z^2 = E$ with another quadric $z-yx = F$,
which for general $E,F$ (and thus for general initial conditions
$(x_0,y_0,z_0)$) gives a genus-$1$ curve ${\cal E} = {\cal E}_{E,F}$.
Finally, since $\rho,\sigma$ are rational involutions on ${\cal E}$
that are not translations, $\beta=\rho\sigma$ is translation on ${\cal E}$
by some point that depends on $E$ and $F$, and we're done.
It seems reasonable to guess that one can prove, for any smooth $f$ with
a local minimum, that the bouncing-ball dynamical system will be non-chaotic
for $E$ close enough to the minimum, by deforming the second invariant
$z-yx$ of the approximating parabola; but I'm too far from my expertise here
to carry out such a proof.
