Just like in
ordinary analysis of global error, you have a truncation error that
behaves asymptotically as
$$ \tau \sim h^2 \partial^4 u + h^2\epsilon^{-1}\partial^3u, $$
and the global error approximately solves the equation
$$ \nabla\cdot(a\nabla E) = -\tau. $$
Ordinarily (I remember only 1d $a=1$ explicitly), this leads to an estimate
like $\|E\|\sim Ch^2\|\partial^2u\|$.
When $u$ oscillates on a scale of $\epsilon$, its derivatives become
large, behaving as $\partial^ku=O(\epsilon^{-k})$. Therefore, keeping
$h$ constant and $\epsilon$ variable (and small), the truncation error
for a given $\epsilon$ behaves as
$$ \tau_\epsilon \sim \epsilon^{-4}\tau. $$
Now when solving the equation $\nabla\cdot(a\nabla
E_\epsilon)=-\tau_\epsilon$, the only question is to estimate the
magnitude of $E_\epsilon$.
Because $\tau_\epsilon$ oscillates rapidly, one might expect some
amount of "smoothing" by the elliptic equation. For this analysis I
found Introduction to Perturbation Methods by Holmes, Chapter 5 very
helpful. The relevant result is as follows.
In solving $\nabla\cdot(a\nabla u) = f$ we can write $u=u_0+\epsilon u_1+\cdots$ to
leading order, and then the solution is given by the equation
$$ \nabla_x\cdot\big(\langle a \partial_{y_i}w\rangle \cdot
\nabla_xu_0 + \langle a\rangle\partial_{x_i} u_0\big) =
\langle f\rangle, $$
(the expression in the brackets is a vector with components indexed by $i$), where $w$ is a vector of functions satisfying the equations
$$ \nabla_y\cdot(a \nabla_y w_k) = -\partial_{y_k}a, $$
and $\langle\cdots\rangle$ denotes averaging over the periodic
fast-scale variable:
$$ \langle q(x, y)\rangle = \int_0^1 q(x, y)\,dy, $$
So to estimate global error, we need
$\langle\tau(x)\rangle$. The local truncation error $\tau(x)$ consists
of derivatives of the exact solution $v(x, y)$, and scales as
$\epsilon^{-4}$. If I can use the one-dimensional expression, it is
$$ \tau(x) \sim C_1 h^2\partial^4v + C_2h^2\epsilon^{-1}\partial^3v,
\qquad \partial = \epsilon^{-1}\partial_y + \partial_x. $$
So the leading term is $\epsilon^{-4}h^2(C_1v_{yyyy}+C_2 v_{yyy})$.
When this is averaged over the period $y$, periodicity causes the
result to vanish to first order:
$$ \epsilon^{-4}h^2\big( C_1 v_{yyy} + C_2 v_{yy})|_{y=0}^1. $$
But periodicity is only valid to first order, so $$v_{yyy}|_{y=0}^{1} \sim
\epsilon v_{xyyy}.$$ So the averaged local truncation error behaves
as
$$ \langle \tau\rangle \sim \epsilon^{-3}h^2(C_1 v_{xyyy}+C_2v_{xyy}), $$
and the magnitude of the global error scales as $\epsilon^{-3}$.
This slightly ad-hoc and crufty, and doesn't give a precise error estimate, just the scaling power for $\epsilon$. I don't know how to derive from this an estimate on the global error $E_\epsilon$ in terms of some norm of the exact solution. I understand that this is very
hand-wavy, but this is the most intuitive I could do, and the order doesn't depend on problem dimensions, so it
addresses your other question as well.
