Suppose I want to solve the PDE:

$$-\nabla\cdot\left( a_\epsilon\left(x\right)\nabla u\right)=0 \text{ in } \Omega$$ $$u=g(x) \text{ on }\partial\Omega$$

Here, I assume that $a_\epsilon(x)=a(x,\frac{x}{\epsilon})$ for some $a(x,y)$ such that $a$ is 1-periodic in $y$.

I know that if I apply the standard finite difference scheme on this PDE with a uniform grid of stepsize $h$, I can obtain a numerical solution $u_h$ an error estimate of the form

$$||u-u_h||\leq Ch^2$$

where $C$ is a constant independent of $h$

I've found a paper which suggests that the error estimate can be further refined so that the constant $C$ is also independent of $\epsilon$ as well. It is stated (without proof or reference) that for the above problem in 2D, the error estimate can be written as

$$||u-u_h||\leq C*\frac{h^2}{\epsilon^3}$$

where $c^*$ only depends on the properties of the function $a$.

For the life of me, I can't seem to derive this estimate. The above PDE is a textbook example of an asymptotic homogenization problem, but all the sources I've found seem to take this as a given without any proof. Any insight into how this factor of $1/\epsilon^3$ is obtained would be greatly appreciated. I highly suspect that this formula may only apply in the 2D case as stated above. I wonder if there is a generalized formula for arbitrary dimensions.

  • 2
    $\begingroup$ I don't know if the formula is correct, but the general idea is that if $a_\epsilon$ oscillates with frequency $\epsilon$, then $\nabla u_h$ must do so as well (and with the same magnitude). So $\nabla^2 = {\cal O}(\epsilon^{-1})$ and the a priori error estimate $\|u-u_h\| \le C h^2 \|u\|_{H^2}$ already yields a factor of $1/\epsilon$. I don't know where the other two powers would come from, though. $\endgroup$ Dec 1, 2014 at 3:31
  • 1
    $\begingroup$ @WolfgangBangerth I think the extra two powers might come from the idea that the "effective" step size is $h/\epsilon$. (Imagine scaling the whole region by $1/\epsilon$ and using the same number of nodes to discretize the equation, then scaling back.) How do you get $O(\epsilon^{-1})$ for $\nabla^2$? If a function oscillates on a scale $\epsilon$, shouldn't its derivatives behave as $\nabla^k=O(\epsilon^{-k})$, so $\epsilon^{-2}$? $\endgroup$
    – Kirill
    Feb 11, 2015 at 2:00
  • $\begingroup$ Yes, I think it's a typo (it's been a while). I suppose I meant to say $\nabla u ={\cal O}(\varepsilon^{-1})$. $\endgroup$ Feb 11, 2015 at 2:18

2 Answers 2


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.


May be a look to the following article of Efendiev and al. may help you (eq 3.8). http://users.cms.caltech.edu/~hou/papers/MsFEM_nonlinear.pdf

Cheers, Thomas

  • $\begingroup$ It would be helpful if you could explain a little bit why you think this formula is helpful in this case. $\endgroup$
    – Paul
    Dec 16, 2014 at 20:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.