# Generating a non-uniform grid

I am interested in generating a 1D non-uniform grid on the interval [0, L] with N points, where a region of width $\sigma$ and centred at $\mu$ is at a higher density and where the transition from low and high densities of grid points occurs over a length $\ell$. This is for a finite-difference code, where particular attention is required to a region with large gradients.

In my current implementation I specify low and high grid spacings and interpolate between these using a pair of tanh functions. However in this scheme the total number of points isn't known a priori. I can then iteratively adjust the low and high densities until the total number of points is N but this is quite convoluted in practice.

My question then is: can anyone help by describing a scheme for generating such a grid in a simpler way, perhaps something analogous to the way you might generate Gauss-Lobato points x(i) = cos(i PI / N) where i is the grid index.

You achieve this with equidistribution to a mesh density function $\rho(x)$.

If you consider $x$ as a continuous map from $\xi \in [0,1]$ into your domain $[0,L]$, then the statement $x$ equidistributes $\rho$ is equivalent to $$\int_0^{x(\xi)} \rho(x') \mathrm{d}x' = \xi \theta\,,$$ where $$\theta = \int_0^L \rho(x) \mathrm{d}x\,.$$ Not sure if this has a name, but I've been calling it MMPDE0, for moving mesh PDE 0 (rather confusing out of context, but from work in Huang, Ren, Russell (1994)) Another way of stating this, is that you're looking for the $x_i$ for $i=0,\dotsc,N$ that satisfy $$\frac1{\theta} \int_{x_{i-1}}^{x_i}\rho(x')\mathrm{d}x' = \frac1N\,.$$

In your case, your mesh density function is piecewise linear, so you can find an explicit expression for $\int \rho(x)\mathrm{d}x$ (which will be piecewise quadratic), which you can solve to find the next $x_i$ from $x_{i-1}$.

(Note that this will give $N+1$ points in $[0,L]$ with $x_0 = 0$, $x_N=L$, so adjust for your application as necessary).

Here is a possible algorithm you could use to find your $x_i$: Here, $y_j$ are the ordered knots of your piecewise linear mesh density function

• $x_0 = 0$.
• For $i$ from 1 to $N-1$
• Find the highest $j$ such that $$\int_0^{y_j}\rho(x) \mathrm{d}x < \frac{i \theta}{N},$$ (trapezium rule)
• Find $x_i$ such that \begin{aligned} \int_0^{x_i}\rho(x')\mathrm{d}x' &= \int_0^{y_j}\rho(x')\mathrm{d}x' + \int_{y_j}^{x_i}\rho(x')\mathrm{d}x'\\ &= \frac{i\theta}{N} \end{aligned} (involves solving a quadratic equation in $x_i$)
• $x_N = L$

It is generally a good idea to smooth your mesh density function before equidistributing to it. It should still be okay to represent it as a piecewise linear mesh density function, just with a lot more knots.

• Thanks this appears to be just what I'm looking for, although I may need a clarification (I'm not familiar with the knots terminology). I've coded a working implementation where $\rho(x) = 1 + \tanh((x-(\mu-\sigma))/\ell) - \tanh((x-(\mu+\sigma))/\ell)$, where the indefinite integral $\int \rho(x)dx$ has an analytic form. Currently this works, although I need to find $x_i$ from $x_{i-1}$ using a root finder (see gist.github.com/hemmer/04ab3689493b86b3bb3d). I feel like the root finding step can be done analytically though? – Hemmer Jan 18 '16 at 14:02
• By knots I mean the non-smooth points of a piecewise linear function. With regards to the analytic form, you might be able to rearrange for $x$, but I don't think so (did a quick check with a symbolic toolbox, but I might have made a mistake). Are you now getting out $N$ points then? Are you just looking for a quicker implementation? – Steve Jan 18 '16 at 14:16
• Yes I'm getting N points out now - your answer has been very useful. If I understand correctly, I could replace the $\tanh$ interpolation between high and low densities with a linear interpolation - this would allow me to have a simpler analytic expression for finding $x_i$ from $x_{i-1}$. See: i.imgur.com/bYf5LzZ.png – Hemmer Jan 18 '16 at 14:24
• Piecewise linear $\rho$ is a convenient monitor function since you can just use the quadratic formula, but for a more general $\rho$ it might not be possible to get an explicit $x$. Another useful one is a Witch of Agnesi $\rho_{\epsilon}(x) = \frac{\epsilon}{\epsilon^2 + x^2}$, which gives $x(\xi) = \epsilon \tan\left(\theta \xi + \text{arctan}\left(\frac{a}{\epsilon}\right)\right)$ (interval $[a,b]$) – Steve Jan 18 '16 at 14:25
• Hmm, it might be simpler but if you're going to use it for finite difference it's advisable to have a smooth $\rho$. You could approximate your smooth function by a high(ish) resolution piecewise linear function – Steve Jan 18 '16 at 14:29