An Unconditionally Monotone C 2 Quartic Spline Method with Nonoscillation Derivatives

A one-dimensional monotone interpolation method based on interface reconstruction with partial volumes in the slope-space utilizing the Hermite cubic-spline, is proposed. The new method is only quartic, however is C and unconditionally monotone. A set of control points is employed to constrain the curvature of the interpolation function and to eliminate possible nonphysical oscillations in the slope space. An extension of this method in two-dimensions is also discussed.


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.

Description of the Problem
A set of points { } , i i x s are specified for 0,1, 2, , i N = with { } i x ordered increasingly as well as { } i s , i.e.
x x + < and 1 i i s s + < hold for each integer i in range, as shown in Figure 1.Here N is the number of intervals defined by the ( ) The problem is to construct an explicit polynomial curve ( ) g x that passes all the given data points, i.e., ( ) i i g x s = , for 0 i N ≤ ≤ , and has a continuous second derivative.In addition, the curve must be strictly increasing (being monotonic ( ) 0 g x ′ > ).Furthermore, the derivative of the interpolation function, ( ) 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 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.

Reduction of the Problem
We consider a reduced problem in the slope space of the original problem.Let the slope function ( ) ( ) ( ) and ( ) f x has a continuous first derivative (being C 1 ).Then, we can see that the integral of ( ) 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 x that satisfy the above given constraints. where A solution of the problem must satisfy Equations ( 1), ( , 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 ( ) has to be smooth for differentiation, i.e. ( ) 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 ( ) x y and ( ) x y (here L and R stand for the "left" and "right" boundaries of an interval) such that , where 0 p p are the estimates of ( ) q q are the estimates of ( )

H u u u H u u u H u u u H u u u
Next, we will show that how ( )  that satisfies the said constraints, obtained with the proposed interpolation method.

Area-Matching for Selecting Control Points on Interval-Walls
The slope of ( ) at each inner data point can be estimated numerically.For example, using a quadratic interpolation on the three points , , , , , one is able to obtain an estimate of ( ) i f x ′ , 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) Δ i s 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 ( ) To start with, we consider a given internal [ ] x x + as shown in Figure 4.
The average of ( ) f x 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 R L x x ∆ ≡ − .We build a local Cartesian coordinate system ( ) , ξ η with its origin at Δ , 2 , and let a quadratic curve represent the interface, such that The area matching means the integral of η over the interval is 0, or 3 Δ Δ 0. 12 a c + = 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 i x x = 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 ( ) 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.

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 H A under a Hermite cubic spline in the in- terval [ ] x x can be calculated with ( ) ( ) Let the control point immediately left to where , Our choices of the x-slope terms are and the y-slope terms are ( ) ( ) They have been explicitly substituted in Equation ( 9) the expression of the physical area.The ij σ terms depend only on the x-coordinates of the bounda- ries of an interval.The terms ijkl I are area integrals of the Hermit spline defined as ( ) ( ) Each of , , , i j k l takes the values 0 or 1.These integrals are evaluated as the following 0000 0010 0001 0011 , , , Δ .
Utilizing Equation ( 9) one arrives at ( ) Because each of the 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 , and . 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.

Boundary Conditions
For the interval on the left boundary, we must specify the slope terms , L L p q .
For the interval on the right boundary, we must provide the slope terms ,

Positivity in the Slope-Space
We have obtained ( ) f x with a Hermite cubic spline on the control points computed in the previous sections.This solution may not necessarily be positive when i h 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 0 i s ∆ =, in the ( ) ( ) , g x x space this means ( ) g x must be a constant in the corresponding interval.Equivalently, the slope ( ) f x 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 , L R q q to zeros in the Hermite cubic spline in case an interval contains a point with a negative ( ) 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 to match the areas in the two neighbor intervals.Because the area match condition is linear for a single variable After the above treatment there will be no occurrence of ( ) 0 f x < anymore.
Therefore, the monotonicity of ( ) g x is satisfied unconditionally, see Figure 9.
Not to mention that no unnecessary oscillations are introduced with this local treatment.

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 ( ) with an excellent quality.We have picked a general parametric form of ( ) ( ) y y u = in Equation ( 5).In practice one can simply pick that 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 for the left subinterval of the original interval "i" with ( ) ( ) ( ) Let us define that Then for the right subinterval this time.The slope terms are defined as ( ) 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 x is differentiable, ( ) g x is an un- conditionally 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 ( ) 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.

In Two-Dimensions
Consider a set of data triplets { } ,

S x y S x y f x y g x y x y
with no oscillation for both f and g ?If the given dataset is also consistent with a sufficient condition for data monotonicity , 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 x y .Without loss of generality, we assume , 0 i j s > and modify a given data set satisfying Equation (16), by adding a layer of phantom data points by simply choosing 0 1 x x < , 0 1 y y < , and assigning ,0 0, 0, 0 , i i j j x x x y y y , 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 , i j v exactly in this cell.
The average of the intersections between the straight line , i i x x y y = = 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 ( ) , , d d .

An Application: Arbitrary Rebinning of Statistical Data
A set of statistical data pair { } ( ) 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.

Slope Estimate
At each internal knot i x , we find its neighbor knots that passes the three data points and the solution is explicit that For an estimate of slopes we can differentiate ( ) S x and take

A Quadratic Polynomial Area-Match to Estimate Mid-Interval Control Points
For each interval [ ] 1 , i i x x + , we fit a quadratic polynomial ( ) ( ) 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 The value of provides the estimate of the height of the mid interval control point 2 1 i C + .The quadratic area fitting function ( ) F x described above has the 3 rd order spatial accuracy and provides the initial heights of the mid interval control points.

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 i x x = and the intersection is taken as the position of control point 2i C .In the unlikely case that this control point is below the x-axis, for positivity, we would lift it to ( ) ,0 i x 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.

Adjusting Mid-Interval Control Points
At this step, each original interval is broken to two subintervals of equal lengths.Each interval [ ] 1 , i i x x + corresponds to a mid-interval control point  will generate the solution quickly for ( ) g x is monotone.The first guess can be picked with a search to find the i th bin that contains the solution, i.e.

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.

Figure 1 .
Figure 1.A monotone interpolation of a strictly increasing data set { } , i i x s .
d y u on the left and the right ends of a given interval.( ) ij H u are the Hermit cubic spline base functions that

Figure 3 .
Figure 3.A good monotone interpolation ( ) g x should have no nonphysical oscillations

For slope match at the left end Δ 2 ξ
= − and the right end

Figure 4 .
Figure 4.A quadratic fitting in the local coordinate system ( ) , ξ η in [ ] 1 , i i x x + .The two = , and the Advances in Pure Mathematics

Figure 5 .
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 [ ] 1 , i ix x + is broken to two subintervals for another round of

Figure 6 .
Figure 6.The area under a Hermit spline curve defined in a general interior interval [ ] , L R x x .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 [ ] , LL L x x , Now let us consider the i th interval, with which 3 control points are involved, see in Figure7.They are

+
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.
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 Figure8).These ghost points provide closure of the solution of the tri-diagonal linear system mentioned in the last subsection.

Figure 8 .
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.
interval so the solution is trivial and does not affect rest of the intervals.

Figure 9 .
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.
pair of integers i, j in the range.Can one construct a polynomial surface ( ) , S x y that passes every data point (i.e. ( ) , , integer pair ( ) , i j ), and has positive spatial derivatives ( ) from Equation (16) using the boundary values.The interface reconstruction is done in the ( )2 s x y ∂ ∂ ∂ space.The "fluid volume" confined in each rectangular cell 1 1

2 ,
must hold every- where inside the computational domain.Thus the spline is monotone.Because explicit slope control is provided for ( ) , V x y with setting the control points by volume conservation, there is no oscillation on the derivatives of ( ) has no oscillation.In addition, the spatial slopes is continuous for a bicubic polynomial ( ) () is then at least twice diffe- rentiable.

Figure 10 . 2 ,
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 ( ) ( ) 2 , V x y S x y = ∂ ∂ ∂ .

Figure 12 .
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 x .The red bin structure is the original and the green bin structure is for even bin counts in each bin.
Figure 12.One should easily understand the feature of no oscillation with ( ) f x is necessary for a sensible solution of rebinning the statistical data.