It is probably a student level question but I can't exactly make it cleat to myself. Why is it more accurate to use non-uniform grids in the numerical methods? I am thinking in the context of some finite-difference method for the PDE of the form $u_t(x,t)=u_{xx}(x,t)$. And assume I am interested in a solution at the point $x^{\ast}$. So, I can see that if I approximate the second derivative, for example on a uniform grid using three point approximation, the error is second order $O(h^2)$. Then I can construct non-uniform grid via a mapping and find coefficients for the three points that are used to approximate the derivative. I can do the Taylor expansions and again obtain a bound for the derivative to be a second order $O(h^2)$, where $h$ is the distance on a uniform grid from which I obtained mapping to a non-uniform grid. Both estimates contain derivatives and it is not clear to me why the solution would be more accurate on non-uniform grid as it depends on the magnitude of the corresponding derivatives in error estimates?

up vote 18 down vote accepted

The rationale for non-uniform meshes goes like this (all equations understood to be qualitative, i.e., in general true but without the pretense to be provably so in all circumstances and for all equations or all possible discretizations):

When solving an equation with, say, linear finite elements, then you typically have an error estimate of the kind $$ \|u-u_h\|_{L^2(\Omega)} \le C h_{\text{max}}^2 \|\nabla^2 u\|_{L^2(\Omega)}, $$ or, equivalently but in a form better suited to the following: $$ \|u-u_h\|_{L^2(\Omega)}^2 \le C h_{\text{max}}^4 \|\nabla^2 u\|_{L^2(\Omega)}^2. $$ However, this is an overestimate. In fact, one can in many instances show that the error is actually of the form $$ \|u-u_h\|_{L^2(\Omega)}^2 \le C \sum_{K \in {\mathbb T}} h_K^4 \|\nabla^2 u\|_{L^2(K)}^2. $$ Here, $K$ are the cells of the triangulation $\mathbb T$. This shows that in order to make the error small, it isn't actually necessary to reduce the maximal mesh size $h_\text{max}$. Rather, the most efficient strategy will be to equilibrate the cellwise error contributions $h_K^4 \|\nabla^2 u\|_{L^2(K)}^2$ -- in other words, you should choose $$ h_K \propto \|\nabla^2 u\|_{L^2(K)}^{-1/2}. $$ In other words, the local mesh size $h_K$ should be small where the solution is rough (has large derivatives) and large where the solution is smooth, and the formula above provides a quantitative measure for this relationship.

  • 1
    I would add that anisotropy is most efficiently represented with an anisotropic ansatz space (i.e., an anisotropic mesh). Since the anisotropy may not be aligned with some initial coarse mesh, an isotropic AMR algorithm may be very inefficient. Anisotropy causes some extra problems because many methods are not uniformly stable with respect to aspect ratio. – Jed Brown Jan 3 '13 at 3:12

Prove it to yourself with this example. What is the maximum error when interpolating sqrt(x) on the interval [0,1] with piecewise linear interpolation on a uniform mesh?

What is the maximum error when interpolating on a mesh in which the ith of n points is given by (i/n)^s, and s is a carefully chosen mesh grading parameter?

  • This is actually simple and intuitive. In fact, if I write the error it depends on some derivative * $h_i$, so when the derivative is big I kill that with a smaller $h_i$. And thus, I use the non uniformity where the function has high gradients, or where has some wiggles. Am I right, that I should put more points in that region, even if it is not the region where I estimate the solution? Because at first I thought to put more points in the region of interest of the answer, but from this discussion, it not there. – Kamil Jan 4 '13 at 2:25

The typical reason why a non-uniform can lead to higher accuracy is that the PDE to solve is not of the form $u_t(x,t)=u_{xx}(x,t)$, but of the form $u_t(x,t)=(D(x)u_x(x,t))_x$. If your non-uniform grid allows you to represent the real $D(x)$ more accurately, you will get a more accurate solution. Because $D(x)$ is normally determined by material properties, it is likely constant within each material, so you normally have a piecewise constant function and should really align the grid accordingly.

A different reason might be that $u(x,0)$ has more variation in certain regions than in others. One can try to compensate this a bit with an adaptively refined non-uniform grid. (However, in my opinion there are other techniques which can handle this situation even better than a non-uniform grid.)

  • could you please specify, what are other techniques you would use in order to have a closer "look" at the regions of discontinuities on the initial data, for example? – Kamil Jan 3 '13 at 1:51
  • @Kamil I have two things in mind here. The first thing is to compute the projection of the initial data into the "representation used on the grid" with sufficient accuracy. (This typically includes things like oversampling or simple analytic computations at jump discontinuities.) I know that this is just good style and too simple to even mentioning it, but in my experience it's often all that is needed to fix the problems caused by singularities in the input data. – Thomas Klimpel Jan 3 '13 at 9:39
  • The other thing I'm thinking of is modeling part of the input data as boundary conditions. However, the savings from this are often less than a factor two, and boundary conditions are notoriously difficult to get right, at least in my experience. So I would say this is often not worth the effort to do it perfectly (or only worth the effort if the corresponding extension of the problem in that direction is really small, or if you really want high accuracy), and just selecting roughly the right boundary condition and placing the boundary sufficiently far away often works good enough. – Thomas Klimpel Jan 3 '13 at 9:47

Kamil, differential equation solving is global, interpolation is local. In piecewise polynomial interpolation, the accuracy far from the singularity won't be bothered by the singularity. Unfortunately, this isn't at all true for solving an elliptic equation, such as a two-point boundary value problem. The singularity will pollute the approximation globally.

Here's something to try. Solve D(sqrt(x) Du) on [0,1] with homogeneous Dirichlet b.c.s. D is the differentiation operator. Use finite elements or finite differences on an n-point uniform mesh. Compare to a mesh in which the ith point is (1/n)^1.5. Note that the worst error for the uniform mesh is far from the singularity, and much larger than for the graded mesh.

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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