An Unconditionally Monotone C2 Quartic Spline Method with Nonoscillation Derivatives

Show more

1. Introduction

Existing monotone cubic interpolation methods [1] [2] [3] [4] are successful in solving many practical problems. However, they are C^{1} and not fit for certain applications that require a higher degree of smoothness. For a C^{2} continuous monotone interpolation, a quintic polynomial is required [5] [6] [7] , or some subdivision of intervals needs to be performed [8] [9] [10] . An issue with these methods is that the derivative of the interpolation curve can have global oscillations and this limits their usage. There are methods existing for either monotone [11] [12] or with nonoscillation derivatives [13] [14] . In general, a monotone polynomial spline method requires certain constraints on their slope estimate to be satisfied.

In this article, we propose an unconditionally monotone interpolation method that is only quartic, with the C^{2} continuity over the entire domain. The new method has no nonphysical oscillation with its derivative and is 3^{rd} order accurate in space.

2. The Solution Procedure

2.1. Description of the Problem

A set of points $\left\{{x}_{i},{s}_{i}\right\}$ are specified for $i=0,1,2,\cdots ,N$ with $\left\{{x}_{i}\right\}$ ordered increasingly as well as $\left\{{s}_{i}\right\}$ , i.e. ${x}_{i}<{x}_{i+1}$ and ${s}_{i}<{s}_{i+1}$ hold for each integer i in range, as shown in Figure 1. Here N is the number of intervals defined by the $\left(N+1\right)$ given points.

The problem is to construct an explicit polynomial curve $g\left(x\right)$ that passes all the given data points, i.e., $g\left({x}_{i}\right)={s}_{i}$ , for $0\le i\le N$ , and has a continuous second derivative. In addition, the curve must be strictly increasing (being monotonic ${g}^{\prime}\left(x\right)>0$ ). Furthermore, the derivative of the interpolation function, ${g}^{\prime}\left(x\right)$ should have no unnecessary oscillations. Figure 2 demonstrates the oscillation in the derivative space of a monotone interpolation obtained with a simple Hermite-interpolation method applied to a radioactive particle energy distribution problem.

Figure 1. A monotone interpolation of a strictly increasing data set $\left\{{x}_{i},{s}_{i}\right\}$ .

Figure 2. Oscillation in the slope-space of a monotone interpolation. Although the slope of the spline function (red curve) is positive and matches the given areas in each interval defined by the black polylines, thus, exactly passes each data point, such a monotone interpolation is unfit in certain applications like rebinning a data set of radiation energy counts.

2.2. Reduction of the Problem

We consider a reduced problem in the slope space of the original problem. Let the slope function $f\left(x\right)\equiv {g}^{\prime}\left(x\right)$ , then if $f\left(x\right)$ satisfy that

$f\left(x\right)>0,$ (1)

${\int}_{{x}_{i}}^{{x}_{i+1}}f\left(\stackrel{\xaf}{x}\right)\text{d}\stackrel{\xaf}{x}}=\text{\Delta}{s}_{i}\equiv {s}_{i+1}-{s}_{i},$ (2)

and
$f\left(x\right)$ has a continuous first derivative (being C^{1}). Then, we can see that the integral of
$f(\; x\; )$

$g\left(x\right)={s}_{0}+{\displaystyle {\int}_{0}^{x}f\left(\stackrel{\xaf}{x}\right)\text{d}\stackrel{\xaf}{x}},$ (3)

not only passes through all the data points, but also has a well-defined second derivative. Therefore, we pursue the solution of the reduced problem by finding certain $f\left(x\right)$ that satisfy the above given constraints.

The domain $\left[{x}_{0},{x}_{N}\right]$ is divided to N intervals with interval i defined by $\left[{x}_{i},{x}_{i+1}\right]$ . The partial area $\text{\Delta}{s}_{i}$ is bounded by two lines: $y={h}_{i},y=0$ , and the boundaries of interval i, here

${h}_{i}\equiv \text{\Delta}{s}_{i}/\text{\Delta}{x}_{i},$ (4)

where $\text{\Delta}{x}_{i}\equiv {x}_{i+1}-{x}_{i}$ . A solution of the problem must satisfy Equations (1), (2), and (3), and without unnecessary oscillation.

The above description corresponds to certain statistics problems such as re-binning of a radioactive particle energy distribution (see Figure 3 on the next page), with the x-axis being the energy and y-axis being the denisty of particles in each energy bin. A requirement for energy re-binning is that the curve
$f\left(x\right)$ has to be smooth for differentiation, i.e.
$g\left(x\right)$ is C^{2}. Not only that, the solution needs to have no oscillation for a minimal slope variation. Therefore, as the solution is mapped to a different bin-structure, it still makes physical sense.

There can be an infinite number of candidates of the solution. To limit the choices we consider a Hermite cubic-spline between an arbitrarily given pair of data points $\left({x}_{L},{y}_{L}\right)$ and $\left({x}_{R},{y}_{R}\right)$ (here L and R stand for the “left” and “right” boundaries of an interval) such that

$\begin{array}{l}x\left(u\right)={x}_{L}{H}_{00}\left(u\right)+{x}_{R}{H}_{01}\left(u\right)+{p}_{L}{H}_{10}\left(u\right)+{p}_{R}{H}_{11}\left(u\right),\\ y\left(u\right)={y}_{L}{H}_{00}\left(u\right)+{y}_{R}{H}_{01}\left(u\right)+{q}_{L}{H}_{10}\left(u\right)+{q}_{R}{H}_{11}\left(u\right),\end{array}$ (5)

where $0\le u\le 1$ is a non-dimensional parameter. ${p}_{L},{p}_{R}$ are the estimates of $\left(\text{d}x/\text{d}u\right)$ , ${q}_{L},{q}_{R}$ are the estimates of $\left(\text{d}y/\text{d}u\right)$ on the left and the right ends of a given interval. ${H}_{ij}\left(u\right)$ are the Hermit cubic spline base functions that

$\begin{array}{l}{H}_{00}\left(u\right)=1+{u}^{2}\left(2u-3\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{H}_{01}\left(u\right)={u}^{2}\left(3-2u\right),\\ {H}_{10}\left(u\right)=u{\left(u-1\right)}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{H}_{11}\left(u\right)={u}^{2}\left(u-1\right).\end{array}$ (6)

Next, we will show that how $f\left(x\right)$ can be constructed with the above Hermite spline to satisfy Equations (1), (2), and (3), by properly choosing a set of control points.

Figure 3. A good monotone interpolation $g\left(x\right)$ should have no nonphysical oscillations in its slope space $\left(f\left(x\right),x\right)$ besides matching the given areas under the black polylines bounded inside each interval exactly. The green curve is such a positive function $f\left(x\right)$ that satisfies the said constraints, obtained with the proposed interpolation method.

2.3. Area-Matching for Selecting Control Points on Interval-Walls

The slope of $f\left(x\right)$ , ${f}^{\prime}\left(x\right)\equiv {g}^{\u2033}\left(x\right)$ at each inner data point can be estimated numerically. For example, using a quadratic interpolation on the three points $\left({x}_{i-1},{s}_{i-1}\right),\text{\hspace{0.17em}}\left({x}_{i},{s}_{i}\right)$ and $\left({x}_{i+1},{s}_{i+1}\right)$ , one is able to obtain an estimate of ${f}^{\prime}\left({x}_{i}\right)$ , except at the left and right boundaries. We consider the reduced problem as an interface reconstruction problem for volume conservation. The approach is to construct the geometry of the interface contained in interval i and to match the partial volume (area) $\text{\Delta}{s}_{i}$ for each i. The interface piece constructed in interval i in general does not match with the pieces constructed in its neighbor intervals on boundaries. We will apply a Hermit spline later to eliminate the gaps and ensure a global slope continuation of $f\left(x\right)$ .

To start with, we consider a given internal
$\left[{x}_{i},{x}_{i+1}\right]$ as shown in Figure 4. The average of
$f\left(x\right)$ in this interval, h_{i} is defined in Equation (4). We replace the subscripts “i” with “L” (left) and “i + 1” with “R” (right) for generality and let the width of the interval be
$\Delta \equiv {x}_{R}-{x}_{L}$ . We build a local Cartesian coordinate system
$\left(\xi ,\eta \right)$ with its origin at
$\left({x}_{L}+\frac{\text{\Delta}}{2},h\right)$ , and let a quadratic curve represent the interface, such that

$\eta =a{\xi}^{2}+b\xi +c.$

For slope match at the left end $\xi =-\frac{\text{\Delta}}{2}$ and the right end $\xi =+\frac{\text{\Delta}}{2}$ , we have

$-a\text{\Delta}+b={{f}^{\prime}}_{L},\text{\hspace{0.17em}}\text{\hspace{0.17em}}a\text{\Delta}+b={{f}^{\prime}}_{R}.$

The area matching means the integral of $\eta $ over the interval is 0, or

$\frac{a}{12}{\text{\Delta}}^{3}+c\text{\Delta}=0.$

Figure 4. A quadratic fitting in the local coordinate system $\left(\xi ,\eta \right)$ in $\left[{x}_{i},{x}_{i+1}\right]$ . The two end slopes and the area under $y={h}_{i}$ are exactly matched. The blue dots on the ends are to be repositioned by finding the intersection between the interval wall $x={x}_{i}$ and a Hermite spline passing the current yellow control points. The yellow control points in the middle are to be lifted/lowered later with another round of area-matching.

Solving for a, b and c, one arrives at

$a=\frac{{{f}^{\prime}}_{R}-{{f}^{\prime}}_{L}}{2\text{\Delta}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}b=\frac{{{f}^{\prime}}_{R}+{{f}^{\prime}}_{L}}{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}c=\frac{\text{\Delta}}{24}\left({{f}^{\prime}}_{L}-{{f}^{\prime}}_{R}\right).$ (7)

The constant term c carries the position of the interface at the middle of the interval to the 3^{rd} order accuracy (for our fitting here is quadratic).

We will use these mid interval interface positions obtained from the above quadratic fitting as a set of control points to construct our first approximation of the solution. Each interval wall $x={x}_{i}$ intersects the Hermit spline curve passing the set of the mid-interval control points and the intersection is taken as a fixed control point. We have $\left(2N+1\right)$ control points ${C}_{k},\left(k=0,1,2,\cdots ,2N\right)$ . The $\left(N+1\right)$ of them with even subscripts are on the walls of the intervals and are fixed, the rest N of them above the middle points of intervals are to be shifted vertically by matching volumes again.

2.4. Area-Matching for Selecting Mid Interval Control Points

We are to construct a Hermit spline that passes all the control points. For an exact area match, we break each original interval into two subintervals about the control point in the middle of the interval, see Figure 5. With the two neighbor mid interval points, there are three mid interval control points involved in the area-matching. We compute the heights of the mid interval control points by solving a tri-linear linear system.

In Figure 6, the shadowed area ${A}_{H}$ under a Hermite cubic spline in the interval $\left[{x}_{L},{x}_{R}\right]$ can be calculated with

${A}_{H}={\displaystyle {\int}_{0}^{1}y\left(u\right){x}^{\prime}\left(u\right)\text{d}u}.$ (8)

Let the control point immediately left to $x={x}_{L}$ be $LL=\left({x}_{LL},{y}_{LL}\right)$ , and the

Figure 5. The yellow dots at the middle of each interval on their quadratic area fitting curves (green curves) are to be lifted or lowered. The blue control points are the intersection between the interval walls and a cubic spline (red curve) passing the yellow control points. The interval $\left[{x}_{i},{x}_{i+1}\right]$ is broken to two subintervals for another round of area-matching.

Figure 6. The area under a Hermit spline curve defined in a general interior interval $\left[{x}_{L},{x}_{R}\right]$ . To compute the area bounded for the original interval, two areas from the two subintervals are to be added together. In another word, each of the intervals $\left[{x}_{LL},{x}_{L}\right]$ , $\left[{x}_{L},{x}_{R}\right]$ and $\left[{x}_{R},{x}_{RR}\right]$ in this figure should be considered as a subinterval.

one immediately right to $x={x}_{R}$ be $RR=\left({x}_{RR},{y}_{RR}\right)$ . Evaluation of the above area integral provides that

${A}_{H}={\sigma}_{00}{y}_{L}+{\sigma}_{01}{y}_{R}+\frac{1}{2}\left({\sigma}_{10}\left({y}_{R}-{y}_{LL}\right)+{\sigma}_{11}\left({y}_{RR}-{y}_{L}\right)\right),$ (9)

where

$\begin{array}{l}{\sigma}_{00}={x}_{L}{I}_{0000}+{x}_{R}{I}_{0001}+{p}_{L}{I}_{0010}+{p}_{R}{I}_{0011},\\ {\sigma}_{01}={x}_{L}{I}_{0100}+{x}_{R}{I}_{0101}+{p}_{L}{I}_{0110}+{p}_{R}{I}_{0111},\\ {\sigma}_{10}={x}_{L}{I}_{1000}+{x}_{R}{I}_{1001}+{p}_{L}{I}_{1010}+{p}_{R}{I}_{1011},\\ {\sigma}_{11}={x}_{L}{I}_{1100}+{x}_{R}{I}_{1101}+{p}_{L}{I}_{1110}+{p}_{R}{I}_{1111}.\end{array}$ (10)

Our choices of the x-slope terms are

${p}_{L}=\frac{1}{2}\left({x}_{R}-{x}_{LL}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{p}_{R}=\frac{1}{2}\left({x}_{RR}-{x}_{L}\right),$

and the y-slope terms are

${q}_{L}=\frac{1}{2}\left({y}_{R}-{y}_{LL}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{q}_{R}=\frac{1}{2}\left({y}_{RR}-{y}_{L}\right).$ (11)

They have been explicitly substituted in Equation (9) the expression of the physical area. The ${\sigma}_{ij}$ terms depend only on the x-coordinates of the boundaries of an interval. The terms ${I}_{ijkl}$ are area integrals of the Hermit spline defined as

${I}_{ijkl}={\displaystyle {\int}_{0}^{1}{H}_{ij}\left(u\right){{H}^{\prime}}_{kl}\left(u\right)\text{d}u}.$ (12)

Each of $i,j,k,l$ takes the values 0 or 1. These integrals are evaluated as the following

${I}_{0000}=-\frac{1}{2},{I}_{0010}=\frac{1}{10},{I}_{0001}=\frac{1}{2},{I}_{0011}=-\frac{1}{10},$

${I}_{1000}=-\frac{1}{10},{I}_{1010}=0,{I}_{1001}=\frac{1}{10},{I}_{1011}=-\frac{1}{60},$

${I}_{0100}=-\frac{1}{2},{I}_{0110}=-\frac{1}{10},{I}_{0101}=\frac{1}{2},{I}_{0111}=\frac{1}{10},$

${I}_{1100}=\frac{1}{10},{I}_{1110}=\frac{1}{60},{I}_{1101}=-\frac{1}{10},{I}_{1111}=0.$

Now let us consider the i^{th} interval, with which 3 control points are involved, see in Figure 7. They are
${C}_{2i}=\left({x}_{i},{y}_{i}\right),{C}_{2i+1}=\left({x}_{i+\frac{1}{2}},{y}_{i+\frac{1}{2}}\right)$ and
${C}_{2\left(i+1\right)}=\left({x}_{i+1},{y}_{i+1}\right)$ with
${x}_{i+\frac{1}{2}}\equiv \frac{1}{2}\left({x}_{i}+{x}_{i+1}\right)$ , and
${y}_{i+\frac{1}{2}}$ to be determined. For area matching we

Figure 7. The area under the Hermit spline curves on interval i is the sum of two half-interval areas defined with Equation (13). Each area matching involves three mid interval control points (yellow). The solution of a tri-diagonal linear system determines their heights.

need to enforce that the sum of the area contribution from the left subinterval and the right subinterval to equal an known value that

${A}_{H}^{\text{left}}+{A}_{H}^{\text{right}}={h}_{i}\text{\Delta}{x}_{i}\equiv {s}_{i+1}-{s}_{i}.$ (13)

Utilizing Equation (9) one arrives at

${A}_{H}^{\text{left}}={\sigma}_{00}^{\text{left}}{y}_{i}+{\sigma}_{01}^{\text{left}}{y}_{i+\frac{1}{2}}+\frac{1}{2}\left({\sigma}_{10}^{\text{left}}\left({y}_{i+\frac{1}{2}}-{y}_{i-\frac{1}{2}}\right)+{\sigma}_{11}^{\text{left}}\left({y}_{i+1}-{y}_{i}\right)\right),$

${A}_{H}^{\text{right}}={\sigma}_{00}^{\text{right}}{y}_{i+\frac{1}{2}}+{\sigma}_{01}^{\text{right}}{y}_{i+1}+\frac{1}{2}\left({\sigma}_{10}^{\text{right}}\left({y}_{i+1}-{y}_{i}\right)+{\sigma}_{11}^{\text{right}}\left({y}_{i+\frac{3}{2}}-{y}_{i+\frac{1}{2}}\right)\right).$

Because each of the ${\sigma}_{ij}$ terms involves only the x-coordinates of the control points, it is a constant. The sum of areas under each Hermit cubit splines defined on a subinterval is simply a linear combination of ${y}_{i-\frac{1}{2}}$ , ${y}_{i+\frac{1}{2}}$ , and ${y}_{i+\frac{3}{2}}$ . Therefore, we have a tri-diagonal linear system involving all intervals to solve in order to match the area exactly in each interval, with certain boundary conditions provided for the left-most and the right-most intervals.

2.5. Boundary Conditions

For the interval on the left boundary, we must specify the slope terms ${p}_{L},{q}_{L}$ . For the interval on the right boundary, we must provide the slope terms ${p}_{R},{q}_{R}$ as well. Currently we assume two kinds of boundary conditions. The first is a symmetrical boundary condition with which one sets a ghost control points at the reflection point of the nearest inner control point. The other one is a counter-symmetric condition by setting a ghost control point out of the boundary by extending the line-segment defined by the two nearest known control points involved (see Figure 8). These ghost points provide closure of the solution of the tri-diagonal linear system mentioned in the last subsection.

Figure 8. A demonstration of the symmetric (the right side) and the counter-symmetric (the left side) boundary conditions. A ghost control point (pink) is employed in each case.

2.6. Positivity in the Slope-Space

We have obtained $f\left(x\right)$ with a Hermite cubic spline on the control points computed in the previous sections. This solution may not necessarily be positive when ${h}_{i}$ is close or equal to zero. Thus, we need to have a treatment in the possible case this negativity occurs. Fortunately, we have a nearly trivial treatment that is rather easy to perform, at least for isolated cases.

Considering the worst case that $\Delta {s}_{i}=0$ , in the $\left(g\left(x\right),x\right)$ space this means $g\left(x\right)$ must be a constant in the corresponding interval. Equivalently, the slope $f\left(x\right)$ must be zero everywhere in interval i, otherwise any variation would create some negative slopes then the monotone condition is violated. Specifically, the slopes at the two ends must be zeros.

Thus, we choose to enforce the slope terms ${q}_{L},{q}_{R}$ to zeros in the Hermite cubic spline in case an interval contains a point with a negative $f\left(x\right)$ . This means the two control points on interval boundaries are at the same height. Since we also assume a troubled interval is isolated, we lift the neighbor control

points at $\left({x}_{i-\frac{1}{2}},{y}_{i-\frac{1}{2}}\right)$ and $\left({x}_{i+\frac{3}{2}},{y}_{i+\frac{3}{2}}\right)$ to match the areas in the two neighbor

intervals. Because the area match condition is linear for a single variable ${y}_{i-\frac{1}{2}}$ or ${y}_{i+\frac{3}{2}}$ in either neighbor interval so the solution is trivial and does not affect rest of the intervals.

After the above treatment there will be no occurrence of $f\left(x\right)<0$ anymore. Therefore, the monotonicity of $g\left(x\right)$ is satisfied unconditionally, see Figure 9. Not to mention that no unnecessary oscillations are introduced with this local treatment.

Figure 9. With setting the end slopes to zeros (the green curve) and matching the area under the polylines again, the isolated slope-space negativity in a single interval (the red curve) is fixed. The locally modified spline still has a continuous derivative crossing the end points of the interval and does not affect the solution elsewhere.

3. An Unconditional Monotone C^{2} Spline

3.1. Constructing the Proposed Method

We have obtained a Hermite cubic spline in previous sections. The spline is non-negative, differentiable, and the area under it bounded by the walls of intervals and the x-axis exactly matches a specified area for each interval. Therefore, the integration of $f\left(x\right)$ provides a monotone interpolation of the data set $\left\{{x}_{i},{s}_{i}\right\},i=0,1,2,\cdots ,N$ with an excellent quality. We have picked a general parametric form of $x=x\left(u\right)$ and $y=y\left(u\right)$ in Equation (5). In practice one can simply pick that

$u=\frac{x-{x}_{i}}{\text{\Delta}{x}_{i}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}f\left(x\right)=y\left(u\right),$

for interval “i”. All the previous discussions would still hold. This means we have an unconditionally monotone interpolation that is only quartic. It is an integral of the Hermit-cubic splines. Its order is lower than some of the existing monotone interpolation methods. Although we have split each interval to two, thus increased the number of knots. Nevertheless, since our analysis is done in the slope space for a cubic spline, the proposed method may be easier to handle than other monotone interpolation methods. The new monotone interpolation can be explicitly expressed as the follows

$\begin{array}{c}g\left(x\right)={s}_{i}+\left({x}_{i+\frac{1}{2}}-{x}_{i}\right){\displaystyle {\int}_{0}^{u}f\left(u\right)\text{d}u}\\ ={s}_{i}+\left({x}_{i+\frac{1}{2}}-{x}_{i}\right)\left({y}_{i}{G}_{00}\left(u\right)+{y}_{i+\frac{1}{2}}{G}_{01}\left(u\right)+{q}_{i}{G}_{10}\left(u\right)+{q}_{i+\frac{1}{2}}{G}_{11}\left(u\right)\right),\end{array}$ (14)

for the left subinterval of the original interval “i” with $u\equiv \left(x-{x}_{i}\right)/\left({x}_{i+\frac{1}{2}}-{x}_{i}\right)$ .

The quartic functions ${G}_{ij}\left(u\right)$ are integrations of the Hermite base functions ${H}_{ij}\left(u\right)$ and

${G}_{00}\left(u\right)=u\left(1-{u}^{2}-\frac{{u}^{3}}{2}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{G}_{01}\left(u\right)={u}^{3}\left(1-\frac{u}{2}\right),$

${G}_{10}\left(u\right)=\frac{{u}^{2}}{12}\left(3{u}^{2}-8u+6\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{G}_{11}\left(u\right)=\frac{{u}^{3}}{12}\left(3u-4\right).$

Let us define that

${s}_{i+\frac{1}{2}}=g\left({x}_{i+\frac{1}{2}}\right)={s}_{i}+\frac{1}{2}\left({x}_{i+\frac{1}{2}}-{x}_{i}\right)\left({y}_{i}+{y}_{i+1}+\frac{1}{6}\left({q}_{i}+{q}_{i+\frac{1}{2}}\right)\right).$

Then for the right subinterval $\left({x}_{i+\frac{1}{2}},{x}_{i+1}\right)$ , we have

$\begin{array}{l}g\left(x\right)={s}_{i+\frac{1}{2}}+\left({x}_{i+1}-{x}_{i+\frac{1}{2}}\right){\displaystyle {\int}_{0}^{u}f\left(u\right)\text{d}u}\\ ={s}_{i+\frac{1}{2}}+\left({x}_{i+1}-{x}_{i+\frac{1}{2}}\right)\left({y}_{i+\frac{1}{2}}{G}_{00}\left(u\right)+{y}_{i+1}{G}_{01}\left(u\right)+{q}_{i+\frac{1}{2}}{G}_{10}\left(u\right)+{q}_{i+1}{G}_{11}\left(u\right)\right),\end{array}$ (15)

with $u=\left(x-{x}_{i+\frac{1}{2}}\right)/\left({x}_{i+1}-{x}_{i+\frac{1}{2}}\right)$ this time. The slope terms are defined as

${q}_{i}=\frac{1}{2}\left({y}_{i+\frac{1}{2}}-{y}_{i-\frac{1}{2}}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{q}_{i+\frac{1}{2}}=\frac{1}{2}\left({y}_{i+1}-{y}_{i}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}{q}_{i+1}=\frac{1}{2}\left({y}_{i+\frac{3}{2}}-{y}_{i+\frac{1}{2}}\right),$

if interval “i” contains no negativeness. In the case of a possible isolated negativeness the corresponding q terms are taken to zeros. Because the y (thus q) terms in the above equations are obtained with area matching, Equations (1), (2), and (3) are all satisfied. Because
$f\left(x\right)$ is differentiable,
$g\left(x\right)$ is an unconditionally monotone C^{2} quartic spline. Besides all the above, the proposed method does not have nonphysical oscillations on its derivative for we have explicit slope control with a Hermite cubic spline. This is a desirable feature. It is possible to further refine the solution by minimize the curvature (say) of
$f\left(x\right)$ to adjust the level of the control points for minimal variation in slope. However, we are confident that the proposed solution is of third order accuracy with a quadratic interface reconstruction. Therefore, the control points are almost right on the ideal solution assuming which exists. Then, a further refinement is hardly necessary for a practical purpose.

3.2. In Two-Dimensions

Consider a set of data triplets $\left\{{x}_{i},{y}_{j},{s}_{i,j}\right\}$ given for $i=1,2,\cdots ,\left(I-1\right),I$ and $j=1,2,\cdots ,\left(J-1\right),J$ with the properties ${x}_{i+1}>{x}_{i},{y}_{j+1}>{y}_{j}$ , and ${s}_{\left(i+1\right),j}>{s}_{i,j}$ , ${s}_{i,\left(j+1\right)}>{s}_{i,j}$ for an arbitrary pair of integers i, j in the range. Can one construct a polynomial surface $S\left(x,y\right)$ that passes every data point (i.e. $S\left({x}_{i},{y}_{j}\right)={s}_{i,j}$ for all feasible integer pair $\left(i,j\right)$ ), and has positive spatial derivatives

$f\left(x,y\right)\equiv \frac{\partial S\left(x,y\right)}{\partial x}>0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}g\left(x,y\right)\equiv \frac{\partial S\left(x,y\right)}{\partial y}>0,$

with no oscillation for both $f$ and $g$ ?

If the given dataset is also consistent with a sufficient condition for data monotonicity

${v}_{i,j}\equiv {s}_{i+1,j+1}+{s}_{i,j}-{s}_{i,j+1}-{s}_{i+1,j}>0,$ (16)

for all valid i and j, then, a method for interface-reconstruction in three-dimensions using volume conservation may be applied to constructing a two-dimensional monotone spline that is at least twice differentiable, with no oscillation on the first derivatives of $S\left(x,y\right)$ .

Without loss of generality, we assume ${s}_{i,j}>0$ and modify a given data set satisfying Equation (16), by adding a layer of phantom data points by simply choosing ${x}_{0}<{x}_{1}$ , ${y}_{0}<{y}_{1}$ , and assigning ${s}_{i,0}=0,{s}_{0,j}=0$ for all $0\le i\le I$ , $0\le j\le J$ . The monotone feature that ${s}_{\left(i+1\right),j}>{s}_{i,j},{s}_{i,\left(j+1\right)}>{s}_{i,j}$ can be derived from Equation (16) using the boundary values.

The interface reconstruction is done in the $\left({\partial}^{2}s/\partial x\partial y\right)$ space. The “fluid volume” confined in each rectangular cell ${x}_{i}\le x\le {x}_{i+1},{y}_{j}\le y\le {y}_{j+1}$ , is defined by Equation (16), which is a difference expression for the cross-derivative at the cell center. Because slope estimates can be computed with certain difference schemes, there are enough data to support the definition of a local quadratic polynomial (or some other algebraic expression) that bounds ${v}_{i,j}$ exactly in this cell.

The average of the intersections between the straight line $x={x}_{i},y={y}_{i}$ and the reconstructed local interface expressions in the surrounding cells can be used as a fixed control point. The control points above the middle point of each edge of a given cell can be set similarly (see Figure 10). A bicubic spline polynomial can be defined over all the control points and an exact volume match for all the cells can be performed to determine the heights of the control points at the center of cells. Solution of a sparse linear system is required because a volume integral under such a bicubic spline is a linear combination of the heights of neighbor control points above cell centers.

Finally, the monotone spline can be obtained with an integral of the bicubic spline $V\left(x,y\right)$ described above

$S\left(x,y\right)={\displaystyle {\int}_{{x}_{0}}^{x}{\displaystyle {\int}_{{y}_{0}}^{y}V\left(\stackrel{\xaf}{x},\stackrel{\xaf}{y}\right)\text{d}\stackrel{\xaf}{x}\text{d}\stackrel{\xaf}{y}}}.$

Because $V=\left({\partial}^{2}S/\partial x\partial y\right)>0$ holds everywhere by construction and on the boundaries $x={x}_{0}$ and $y={y}_{0}$ the first derivatives of S are nonnegative, as the spatial integrals of $\left({\partial}^{2}S/\partial x\partial y\right)$ , $\partial S/\partial x>0$ and $\partial S/\partial y>0$ must hold everywhere inside the computational domain. Thus the spline is monotone. Because explicit slope control is provided for $V\left(x,y\right)$ with setting the control points by volume conservation, there is no oscillation on the derivatives of $S\left(x,y\right)$ for $\left({\partial}^{2}S/\partial x\partial y\right)$ has no oscillation. In addition, the spatial slopes is continuous for a bicubic polynomial $V\left(x,y\right)\equiv \left({\partial}^{2}S/\partial x\partial y\right)$ , $S\left(x,y\right)$ is then at least twice differentiable.

Figure 10. The top view of an internal cell (solid line) and its neighbors (dashed lines) are shown above. The blue control points are computed from local interface reconstruction and are fixed. The yellow control point at the center of each cell is going to be shifted in the direction perpendicular to the page to match volumes defined by Equation (16) by solving a linear system in order to construct a bicubic spline $V\left(x,y\right)=\left({\partial}^{2}S/\partial x\partial y\right)$ .

4. An Application: Arbitrary Rebinning of Statistical Data

A set of statistical data pair $\left\{{x}_{i},{s}_{i}\right\},\left(i=0,1,2,\cdots ,N-1,N\right)$ is given. Which can be considered as density of radioactive particle counts vs. energy, as shown in Figure 11. It is desired to bin the data in such a way that the number of counts is even in each bin. We are going to demonstrate how to solve this problem with constructing the proposed monotone interpolation.

4.1. Slope Estimate

At each internal knot ${x}_{i}$ , we find its neighbor knots ${x}_{i-1}$ and ${x}_{i+1}$ and fit a quadratic polynomial

$S\left(x\right)={\alpha}_{i}{\left(x-{x}_{i}\right)}^{2}+{\beta}_{i}\left(x-{x}_{i}\right)+{\gamma}_{i}$

that passes the three data points and the solution is explicit that

${\alpha}_{i}=\frac{{h}_{i}-{h}_{i-1}}{{x}_{i+1}-{x}_{i-1}},\text{\hspace{0.17em}}{\beta}_{i}=\left({h}_{i}\left(\text{\Delta}{x}_{i-1}\right)+{h}_{i-1}\left(\text{\Delta}{x}_{i}\right)\right),\text{\hspace{0.17em}}{\gamma}_{i}={s}_{i}.$

For an estimate of slopes we can differentiate $S\left(x\right)$ and take

$f\left({x}_{i}\right)\approx {S}^{\prime}\left({x}_{i}\right)={\beta}_{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{f}^{\prime}\left({x}_{i}\right)\approx {S}^{\u2033}\left({x}_{i}\right)=2{\alpha}_{i}.$

4.2. A Quadratic Polynomial Area-Match to Estimate Mid-Interval Control Points

For each interval $\left[{x}_{i},{x}_{i+1}\right]$ , we fit a quadratic polynomial

$F\left(x\right)=a{\left(x-{x}_{i+\frac{1}{2}}\right)}^{2}+b\left(x-{x}_{i+\frac{1}{2}}\right)+\left(c-{h}_{i}\right)$

with ${F}^{\prime}\left({x}_{i}\right)=2{\alpha}_{i},{F}^{\prime}\left({x}_{i+1}\right)=2{\alpha}_{i+1}$ the slope estimates described above. One

Figure 11. The statistical counts of certain radio-active events are binned. The x-axis is energy density and the y-axis stands for the average event count per energy unit. This plot has the same interpretation as shown in Figure 2.

finds that

$a=\left({S}^{\u2033}\left({x}_{i+1}\right)-{S}^{\u2033}\left({x}_{i}\right)\right)/\left(2\text{\Delta}{x}_{i}\right),\text{\hspace{0.17em}}b=\left({S}^{\u2033}\left({x}_{i}\right)+{S}^{\u2033}\left({x}_{i+1}\right)\right)/2,\text{\hspace{0.17em}}c=-a\frac{{\left(\text{\Delta}{x}_{i}\right)}^{2}}{12}.$

The value of $c=F\left({x}_{i+\frac{1}{2}}\right)$ provides the estimate of the height of the mid interval control point ${C}_{2i+1}$ . The quadratic area fitting function $F\left(x\right)$ described

above has the 3^{rd} order spatial accuracy and provides the initial heights of the mid interval control points.

4.3. A Hermit Spline to Determine Control Points above Knots

A Hermit cubic spline passing the set of mid interval area matching points obtained above intersects the line $x={x}_{i}$ and the intersection is taken as the position of control point ${C}_{2i}$ . In the unlikely case that this control point is below the x-axis, for positivity, we would lift it to $\left({x}_{i},0\right)$ thus move this control-point closer to the true solution which is by definition above the x-axis. These control points are not to be moved.

4.4. Adjusting Mid-Interval Control Points

At this step, each original interval is broken to two subintervals of equal lengths. Each interval $\left[{x}_{i},{x}_{i+1}\right]$ corresponds to a mid-interval control point

${C}_{2i+1}=\left({x}_{i+\frac{1}{2}},{y}_{i+\frac{1}{2}}\right)$ . ${y}_{i+\frac{1}{2}}$ is computed by solving a tri-diagonal linear system

defined by the area-matching Equation (13) for all i’s.

4.5. Constructing the Proposed Monotone Spline

All the control points are now determined. They define a unique Hermit spline curve $f\left(x\right)$ that passes all the control points, has a continuous derivative crossing the bin-walls, and matches the counts in each bin. We are able to integrate $f\left(x\right)$ to construct the monotone interpolation function $g\left(x\right)$ as described with Equation (14), Equation (15).

4.6. Solving for Evenly Dividing Counts in Each New Bins

Because $g\left(x\right)$ is monotone, one is able to evenly divide the vertical axis in the $\left(g\left(x\right),x\right)$ space between $\left(g\left({x}_{0}\right),g\left({x}_{N}\right)\right)$ to M intervals. Let the new accumulated bin counts be

${S}_{j}={s}_{0}+\frac{j\left({s}_{N}-{s}_{0}\right)}{M},\text{\hspace{0.17em}}\left(j=0,1,2,\cdots ,M\right).$

Then one solves for each j that $g\left({x}_{j}\right)={S}_{j}$ for the new bin-system. Because the derivative ${g}^{\prime}\left(x\right)=f\left(x\right)$ is available at any point in the domain, a Newton- Raphson method

${x}_{j}^{k+1}={x}_{j}^{k}-\left(g\left({x}^{k}\right)-{S}_{j}\right)/f\left({x}^{k}\right),$

Figure 12. The statistical counts in Figure 11 are rebinned with the proposed monotone C^{2} non-oscillation quartic interpolation. The black curve is the slope function
$f\left(x\right)$ . The red bin structure is the original and the green bin structure is for even bin counts in each bin.

will generate the solution quickly for
$g\left(x\right)$ is monotone. The first guess can be picked with a search to find the i^{th} bin that contains the solution, i.e.

${s}_{i}<{S}_{j}<{s}_{i+1}$ then, simply taken to ${x}_{j}^{0}={x}_{i+\frac{1}{2}}$ .

A solution in the case of $N=28$ and $M=40$ is shown in Figure 12. One should easily understand the feature of no oscillation with $f\left(x\right)$ is necessary for a sensible solution of rebinning the statistical data.

5. Conclusion

A one dimensional monotone C^{2} quartic interpolation method is proposed. The original problem is reduced to a volume matching problem in the slope space to locate a set of control points on the solution curve. A quadratic fitting in each interval is employed to first approximately locate the control points with area matching. The solution of a tri-diagonal linear system is employed to relocate the control points above the middle of each interval to exactly matching the areas under a Hermit-cubic spline curve while fixing the control points above the original knots. The integration of the solution in the slope-space provides the desired unconditional monotone C^{2} interpolation. The proposed method has the feature of being of lower order (quartic) and with no undesired oscillation in the slope space. An application to the practical problem of rebinning a set of nonnegative statistical data shows the proposed method is effective in one-dimension. The proposed method may be extended to two-dimensions in a similar fashion with a higher smoothness for data sets that satisfy an extra constraint.

Acknowledgements

This work is performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

The authors would like to express appreciation to Simon Labov and Chad Noble of Lawrence Livermore National Laboratory for their encouragement and support.

References

[1] Fritsch, F.N. and Carlson, R.E. (1980) Monotone Piecewise Cubic Interpolation. SIAM Journal on Numerical Analysis, 17, 238-246.

https://doi.org/10.1137/0717021

[2] Fritsch, F.N. and Butland, J. (1984) A Method for Constructing Local Monotone Piecewise Cubic Interpolants. SIAM Journal on Scientific and Statistical Computing, 5, 300-304.

http://epubs.siam.org/doi/abs/10.1137/0905021

[3] Carlson, R.E. and Fritsch, F.N. (1985) Monotone Piecewise Bicubic Interpolation. SIAM Journal on Numerical Analysis, 22, 386-400.

https://doi.org/10.1137/0722023

[4] Steffen, M. (1990) A Simple Method for Monotonic Interpolation in One Dimension. Astronomy and Astrophysics, 239, 443.

[5] Heß, W. and Schmidt, J.W. (1994) Positive Quartic, Monotone Quintic C2-Spline Interpolation in One and Two Dimensions. Journal of Computational and Applied Mathematics, 55, 51-67.

https://doi.org/10.1016/0377-0427(94)90184-8

[6] Ling, N.S., Kong, V.P. and Ong, B.H. (2009) A Class of Bézier-Like Splines in Smooth Monotone Interpolation. Malaysian Journal of Mathematical Sciences, 3, 203-226.

[7] Dougherty, R.L., Edelman, A. and Hyman, J.M. (1989) Nonnegativity-, Monotonicity-, or Convexity-Preserving Cubic and Quintic Hermite Interpolation. Mathematics of Computation (United States), 52, 471-494.

https://doi.org/10.1090/S0025-5718-1989-0962209-1

[8] Jochen, W., Schmidt, J.W. and Heß, W. (1995) An Always Successful Method in Univariate Convex C2 Interpolation. Numerische Mathematik, 71, 237-252.

https://doi.org/10.1007/s002110050143

[9] Sablonniere, P. and Mani, C. (1997) Monotone Interpolation of Order 3 by C2 Cubic Splines. IMA Journal of Numerical Analysis, 17, 305-320.

https://doi.org/10.1093/imanum/17.2.305

[10] Sablonniere, P. and Merrien, J.-L. (2003) Monotone and Convex C1 Hermite Interpolants Generated by a Subdivision Scheme. Constructive Approximation, 19, 279-298.

https://doi.org/10.1007/s00365-002-0512-3

[11] Hyman, J.M. (1983) Accurate Monotonicity Preserving Cubic Interpolation. SIAM Journal on Scientific and Statistical Computing, 4, 645-654.

https://doi.org/10.1137/0904045

[12] Wolberg, G. and Alfy, I. (1999) Monotonic Cubic Spline Interpolation. Proceeding CGI’99 Proceedings of the International Conference on Computer Graphics, Washington DC, 188.

[13] Han, H. and Guo, X. (2018) Cubic Hermite Interpolation with Minimal Derivative Oscillation. Journal of Computational and Applied Mathematics, 331, 82-87.

https://doi.org/10.1016/j.cam.2017.09.049

[14] Wolberg, G. and Alfy, I. (2002) An Energy-Minimization Framework for Monotonic Cubic Spline Interpolation. Journal of Computational and Applied Mathematics, 143, 145-188.

https://doi.org/10.1016/S0377-0427(01)00506-4