Nereis. Interdisciplinary Ibero-American Journal of Methods, Modelling and Simulation.

INTRODUCTION

In Numerical Methods, such as interpolation, numerical integration or numerical solving of differential equations, we have to discretize the domain wherein we are looking for an approximate numerical solution of a given problem. As initial approximation to this discretization, usually uniform meshes are used. However, if some information is known in advance about the function being discretized, such an abrupt change or a violent oscillation in the neighborhood of a point, say a critical point, we can improve the numerical approximation by using non-uniform meshes which concentrate the number of mesh points nearby this critical point.

The scope of this paper is just to provide a general theory of how we can generate this non-uniform meshes for one-dimensional intervals. From this theory, we can derive closed-form expressions for some particular cases, namely, a Loretzian-type and an exponential-type meshes. Finally, we will provide some examples to see how this kind of meshes improves the interpolation of some functions.

Uniform and non-uniform-meshes

When we have a one-dimensional mesh of n points over an interval (x_min, x_max), we generate n – 1 subintervals (x_i, x_i+1), i = 1, 2,… n – 1. If the mesh is uniform, all subintervals are of the same length h,

(1)

Summing up all the subintervals, we obtain the length of the whole interval,

(2)

We can obtain non-uniform meshes with the aid of a positive function f(x) > 0 (hereafter termed as mesh density function), defining the subintervals as,

(3)

For convenience, the mesh density function f(x) is normalized to unity within the whole interval,

(4)

Thereby, the subintervals will be smaller where the mesh density function is greater (see Fig. 1).

Note that we can recover the uniform mesh taking a constant mesh density function f(x) = C. Indeed, in this case, from (3), we have

(5)

Therefore, summing up in (5) all the subintervals, we obtain

(6)

However, note that from (4),

(7)

thus, (6) and (7) yields,

(8)

Finally, combining (5), (7), and (8), taking into account (1), we arrive at (2), which is the formula for a uniform mesh.

PARTICULAR CASES OF NON-UNIFORM MESHES

Lorentzian mesh density function

Let us consider now a Lorentzian-type mesh density function, that for convenience, we will express as follows:

(9)

The above mesh density function (9) is positive in the real domain for á, β > 0, and has got a unique maximum at x = x₀ (see Fig. 2). Therefore, this point is going to be the critical point around which the density of points in the mesh is going to be higher.

Now, in order to calculate the non-uniform mesh with this Lorentzian-type density function, let us define a subinterval centered on x = x₀,and of width h^* = 2δ (see Fig. 2). Hereafter we will term this subinterval as the maximum subinterval. According to (3), the maximum subinterval has to satisfy,

(10)

As aforementioned, if we have n points in the mesh, the length h of each subinterval in a uniform mesh will be given by (2). However, the length of the maximum subinterval must be smaller, thus

(11)

The shrinkage factor ε will tell us how much the points of the mesh gather around the critical point x₀.

The points of the mesh are calculated iteratively as follows. Considering a Lorentzian density (9) in (3), and knowing that [1]

(12)

we have, for i = 1, 2,…, n – 1,

(13)

where for the initial iteration we know that x₁ = x_min.

Let us calculate now the parameters á and β from the condition for the maximum subinterval (10) and the normalization condition (4). For this purpose, apply the following identity [2]

(14)

to (12), setting integration limits, hence

(15)

Therefore, according to (15), the condition for the maximum subinterval (10) is rewritten as,

(16)

Also, the normalization condition (4) is rewritten as,

(17)

where we have set,

(18)

Also, notice that from (2) and (18) we can rewrite (11) as

(19)

Solving for β² in (16) and (17), and equating both results, we can solve for β, obtaining,

(20)

Now, inserting (20) in (17), we arrive at

(21)

The above equation (21) can be solved numerically using Brent’s method [3], taking as initial bracketing (0, π).

Procedure

In summary, knowing the interval I = (x_min, x_max), the critical point x₀ within the interval I, the number of points in the mesh n, and the shrinkage factor ε, we can compute A with (8), Δ and Κ with (18), and δ with (19). Once all these parameters are set, we can solve for á with (21). Next, β is computed using (20). Finally, the points of the mesh are calculated iteratively using (13).

Exponential mesh density function

The Lorentzian-type mesh cannot be used when the critical point x₀ is at one of the extremes of the interval (x_min, x_max). For these cases, consider an exponential mesh density function,

(22)

The exponential function (22) is positive in the real domain for á > 0. Also, in the interval (x_min, x_max), the maximum is at x = x_max when β > 0, and at x = x_min when β < 0. Therefore, the critical point will be found at one of the extremes of the interval (x_min, x_max) depending on the sign of β.

The points of the mesh are calculated as follows. Inserting (22) in (3), we have,

(23)

where the initial iteration is given by x₁ = x_min.

In order to calculate the parameters á and β, we use the normalization condition (4), which according to (18) and (22) reads as

(24)

Now, consider β > 0, thus the maximum subinterval is given by (x_max – δ, x_max). Similar to (11), the length of this maximum subinterval must be shrunk by a factor ε with respect to a uniform mesh, hence

(25)

Therefore, considering this maximum subinterval in (3), we get a second condition for á and β,

(26)

Now, from (24) and (26), we obtain the following equation for β,

(27)

that can be solved using Newton-Raphson’s method. For the initial iteration β₀, we can expand in McLaurin series the exponentials given in (27) up to second order, arriving at the following approximate equation,

(28)

Once β is computed, we can calculate á from (24) as,

(29)

As mentioned before, if we consider the first interval as the maximum subinterval, then we have to change β → – β.

NUMERICAL EXAMPLES

Interpolation with a Lorentzian mesh

In order to test the usefulness of the Lorentzian mesh, we have interpolated the sigmoid-like function:

(30)

taking n = 9 points. Fig. 3 shows the plot of the g(x) function and the interpolated function with second order splines [4] using a uniform mesh within the interval (–2, 1). We can see that the interpolation deviates from the g(x) function where this function changes abruptly.

However, taking the same number of points, n = 9, but using a Lorentzian mesh centered in x₀ = –0.2 with a shrinkage factor ε = 0.4, the second order splines interpolation improves the result, as it is shown in Fig. 4.

We can evaluate numerically the performance of both interpolations shown in Figs. 3 and 4, by using the following distance [5] between two non-negative functions f₁ and f₂, integrable in a certain interval I = (a, b):

(31)

In [5] is proved that 0 ≤ Δ_I(f₁, f₂) ≤ 1, wherein Δ_I(f₁, f₂) = 0 means that f₁ and f₂ are overlapped within the integration interval I = (a, b), and Δ_I(f₁, f₂) = 1 means that both functions are relatively infinitely far one from each other.

Denoting g_unif(x) as the interpolating function with a uniform mesh (Fig. 3), and g_lorentz(x) as the interpolating function with a Lorentzian mesh (Fig. 4), we have obtained,

(32)

which indicates that the Lorentz mesh improves the interpolation quite significantly.

Interpolation with an exponential mesh

Now, let us interpolate the function,

(33)

within the interval (0, π), taking n = 10 points in the mesh. Fig. 5 shows the plot of the g(x) function and the interpolated function with second order splines [4] using a uniform mesh. We can see that the interpolation deviates from the g(x) function nearby the origin since there the function oscillates significantly.

However, taking the same number of points, n = 10, but using an exponential mesh with the critical point at x_min = 0 with a shrinkage factor ε = 0.25, the second order splines interpolation improves the result, as it is shown in Fig. 6.

Denoting g_unif(x) as the interpolating function with a uniform mesh (Fig. 5), and g_exp(x) as the interpolating function with an exponential mesh (Fig. 6), we have obtained,

(34)

which indicates that the exponential mesh improves the interpolation quite significantly.

CONCLUSIONS

We have obtained a general procedure to generate non-uniform meshes in terms of the density of points required (mesh density function) within a given interval, (3) and (4). This is useful when we know in advance that a function changes abruptly or oscillates violently in the neighborhood of a certain point of the interval (i.e. a critical point x₀), thereby more points in the mesh are needed nearby x₀. In particular, we have obtained closed-form expressions for a Lorentzian mesh density function in (13), which can be applied when the central subinterval belongs to the interval. When the critical point is on one extreme of the interval, an exponential mesh is used (23). Finally, numerical examples are given to highlight how this kind of non-uniform meshes improves the interpolation.

ACKNOWLEDGEMENTS

The author wishes to acknowledge the financial support received from Universidad Católica de Valencia under Grant 2017-160-001.

LITERATURE CITED

1. Spiegel, MR. Mathematical handbook of formulas and tables. New York: McGraw-Hill, 1968. Eqn. 14.39.
2. Gradsthteyn IS, Ryzhik IM. Table of integrals, series and products. 7th ed. New York: Academic Press Inc; 2007. Eqn. 1.625.9.
3. Brent, RP. Algorithms for Minimization without Derivatives. Dover, 2002 (Original edition 1973).
4. Hoffman, JD, Frankel, S. Numerical methods for engineers and scientists. 2nd ed. New York: Marcel Dekker Inc, 2001. Sect. 4.9.
5. González-Santander JL, Martín G. Relative distance between two scalar fields. Application to mathematical modeling approximation. Math Method Appl Sci 2014;37:2906-2922.