Rays and the ray equation: behaviour in optical systems, valid solutions to Maxwell's equations In my studies of geometrical optics, I have encountered variations of the claim that, from what I recall, when the amplitude of the wave is changing slowly, the ray equation is a valid solution to Maxwell's equations, and can therefore be used to evaluate optical systems, or something of the sort. However, this is never elaborated upon, and I'm unsatisfied with my understanding (or lack thereof) of what this really means.
For instance, often accompanying the aforementioned claim is a diagram such as the following:

If my understanding is correct, the ray equation is not valid along the boundaries of the wavefronts in the above diagram, since that's where refraction occurs. But I don't actually understand why this is the case, or, put another way, I don't really appreciate what this means. 
Furthermore, I'm not completely sure what qualifies as a boundary of a wavefront in the above diagram. This is made more confusing by the fact that the author of the diagram has drawn lines depicting rays along (what I'm presuming are) the boundaries of the wavefronts, despite the fact that, as I just said, the ray equation is invalid along the boundaries of the wavefronts.
And lastly, what is happening at the focus in such a diagram? Are foci actually boundaries, and, therefore, the ray equation is invalid at foci? And, tying this back into the aforementioned claim, what happens to the amplitude of the wave at such a point?
I would greatly appreciate it if people could please take the time to explain all of this.

EDIT:
I just found the following on page 388 of Electromagnetic Fields, by Jean G. Van Bladel:

An important example of failure occurs when the rays converge to a focus, where the theory predicts an infinite value for the power density of $I$ in (8.101).

This still isn't an explanation, but at least it confirms what I suspected: The ray equation is invalid at foci.
 A: Consider a packet composed of a roughly uniform train of waves
spread out over a region that is substantially  longer and
wider  than their
mean  wavelength.
Any particular
point of space and  time, ${\bf x}$ and $t$,  it has a definite  phase
$\Theta({\bf x},t)$. Once we know this phase, we can define the
local frequency $\omega$ and wave-vector ${\bf k}$ by
$$
\omega= -\left(\frac{\partial\Theta}{\partial
t}\right)_x,\qquad  k_i = \left(\frac{\partial\Theta}{\partial
x_i }\right)_t. 
$$
These definitions are motivated by the idea that 
$$
\Theta({\bf x},t)\sim {\bf k}\cdot {\bf x} -\omega t,
$$
at least locally.
We wish to understand  how ${\bf k}$ changes as the wave
propagates through a slowly varying  medium. Assume that the dispersion equation
$\omega=\omega ({\bf k})$, which is initially derived for a
uniform medium, can  be extended to $\omega =\omega({\bf
k}, {\bf x})$, where the $\bf x$ dependence arises, for
example, as a result of a position-dependent  refractive
index. This assumption is only an approximation, but  it is
a good approximation when the  distance over which the
medium changes is much larger than the distance between
wavecrests.
Applying the equality of  mixed partials to the definitions
of ${\bf k}$ and $\omega$ gives us
$$
\left(\frac{\partial \omega}{\partial
x_i }\right)_t= -\left(\frac{\partial k_i}{\partial
t }\right)_{\bf x},\quad \left(\frac{\partial k_i}{\partial
x_j }\right)_{x_i}= \left(\frac{\partial k_j}{\partial
x_i }\right)_{x_j}. \quad (\star)
$$
The subscripts indicate  what is being left fixed when we
differentiate. We must be  careful about this,  because we
want to use the dispersion equation to express  $\omega$ as
a function of ${\bf k}$ and  ${\bf x}$, and the wave-vector
${\bf k}$ will itself be a function of ${\bf x}$ and $t$. 
Taking this dependence into account, we write
$$
\left(\frac{\partial \omega}{\partial
x_i }\right)_t= \left(\frac{\partial \omega}{\partial
x_i }\right)_{\bf k} + \left(\frac{\partial \omega}{\partial
k_j }\right)_{\bf x}\left(\frac{\partial k_j}{\partial
x_i }\right)_t.
$$
We now  use ($\star$)  to rewrite this as
$$
\left(\frac{\partial k_i}{\partial
t }\right)_{\bf x}+ \left(\frac{\partial \omega}{\partial
k_j }\right)_{\bf x}\left(\frac{\partial k_i}{\partial
x_j }\right)_t = -\left(\frac{\partial \omega}{\partial
x_i }\right)_{\bf k}.
$$
Interpreting the left hand side as a convective derivative 
$$
\frac{dk_i}{dt}=\left(\frac{\partial k_i}{\partial
t }\right)_{\bf x}+({\bf V}_g\cdot\nabla) k_i,
$$
we read off that
$$
\frac {dk_i}{dt}= -\left(\frac{\partial \omega}{\partial
x_i }\right)_{\bf k}
$$
provided we are moving at velocity 
$$
\frac {dx_i}{dt}=({\bf V}_g)_i= \left(\frac{\partial \omega}{\partial
k_i }\right)_{\bf x}.
$$
Since this is the group velocity, the  packet of
waves is actually travelling at this speed.  The last two  equations 
therefore  tell us how the orientation and wavelength of
the  wave train evolve if we ride along with the packet as
it is refracted by the inhomogeneity.
The formulae
$$
\dot {\bf k} = -\frac{\partial \omega}{\partial
{\bf x} },\\
\dot {\bf x}= \phantom -\frac{\partial \omega}{\partial
{\bf k} },
$$
are Hamilton's ray equations and are the basis for geometric optics. 
These Hamilton equations are 
identical in form to Hamilton's equations for classical
mechanics
except that  ${\bf k}$ is playing the role of the canonical
momentum, ${\bf p}$, and  $\omega({\bf k}, {\bf
x})$  replaces  the Hamiltonian, $H({\bf p}, {\bf
x})$.
