# Finite Difference Error Estimate for an Elliptic PDE with an Oscillatory Coefficient

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.

• 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. – Wolfgang Bangerth Dec 1 '14 at 3:31
• @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}$? – Kirill Feb 11 '15 at 2:00
• Yes, I think it's a typo (it's been a while). I suppose I meant to say $\nabla u ={\cal O}(\varepsilon^{-1})$. – Wolfgang Bangerth Feb 11 '15 at 2:18

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

• It would be helpful if you could explain a little bit why you think this formula is helpful in this case. – Paul Dec 16 '14 at 20:43