3

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.

Hemmer
  • 265
  • 2
  • 10

1 Answers1

4

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.

Steve
  • 541
  • 2
  • 14
  • 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 https://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: http://i.imgur.com/bYf5LzZ.png – Hemmer Jan 18 '16 at 14:24
  • 1
    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
  • Last point: piecewise linear is convenient if you have pointwise densities to mesh to (think error estimates or the like), but if your smooth mesh density function is integrable, then it's probably best to use that. If your root finding is taking a long time, you could code up potentially quicker Newton solver with $x_{i+1}^{(k+1)} = x_{i+1}^{(k)} - \frac{\int \rho - \text{target}}{\rho}(x_{i+1}^{(k)})$ – Steve Jan 19 '16 at 09:32