Shape preserving piecewise polynomial interpolation

The basis values are $H_{00}(0)=1$, $H_{01}(0)=H_{10}(0)=H_{11}(0)=0$, and $H_{01}(1)=1$, $H_{00}(1)=H_{10}(1)=H_{11}(1)=0$. Substitution into [F1] gives $P_i(x_i)=y_i$ and $P_i(x_{i+1})=y_{i+1}$.

Differentiating the displayed basis using [F2] gives $(H_{00}^{\mathrm{\prime }},H_{01}^{\mathrm{\prime }},H_{10}^{\mathrm{\prime }},H_{11}^{\mathrm{\prime }})(0)=(0,0,1,0)$ and $(H_{00}^{\mathrm{\prime }},H_{01}^{\mathrm{\prime }},H_{10}^{\mathrm{\prime }},H_{11}^{\mathrm{\prime }})(1)=(0,0,0,1)$. Since $dt\mathrm{/}dx=1\mathrm{/}{h}_i$, differentiation of [F1] therefore gives $P_i^{\mathrm{\prime }}(x_i)=m_i$ and $P_i^{\mathrm{\prime }}(x_{i+1})=m_{i+1}$.

At a shared node $x_i$, the left and right pieces both have value $y_i$ by step 1.1 and both have derivative $m_i$ by step 1.2. Thus the piecewise function and its first derivative match there; the one-sided endpoint statements follow from the same step 1.2 calculation.

Write a degree-at-most-three polynomial in the normalized coordinate as $q(t)=a_0+a_{1}t+a_2{t}^2+a_3{t}^3$. For $q=1$, its endpoint values are $(1,1)$ and scaled slopes are $(0,0)$; for $q=t$, they are $(0,1)$ and $(1,1)$; for $q=t^2$, they are $(0,1)$ and $(0,2)$; and for $q=t^3$, they are $(0,1)$ and $(0,3)$, by [F2] and the chain rule under $t=(x-x_i)\mathrm{/}{h}_i$. Substitution in [F1] gives $H_{00}+H_{01}=1$, $H_{01}+H_{10}+H_{11}=t$, $H_{01}+2H_{11}=t^2$, and $H_{01}+3H_{11}=t^3$. The formula is linear in endpoint values and slopes, so it reproduces their linear combination $q$.

For $d\mathrm{\ne }0$, differentiating [F1] by [F3], using $y_1-y_0=hd$ and $dt\mathrm{/}dx=1\mathrm{/}h$, gives $p^{\mathrm{\prime }}(x)\mathrm{/}d=\alpha +(6-4\alpha -2\beta )t+(3\alpha +3\beta -6)t^2$. Indeed, $H_{01}^{\mathrm{\prime }}=6t-6t^2$, $H_{10}^{\mathrm{\prime }}=1-4t+3t^2$, and $H_{11}^{\mathrm{\prime }}=-2t+3t^2$; substitution produces those coefficients. Direct substitution in [F1] also gives $p(x_0)=y_0$ and $p(x_1)=y_1$.

For $0\le t\le 1$ put ${\lambda }_{00}=(1-\alpha \mathrm{/}3)(1-\beta \mathrm{/}3)$, ${\lambda }_{30}=(\alpha \mathrm{/}3)(1-\beta \mathrm{/}3)$, ${\lambda }_{03}=(1-\alpha \mathrm{/}3)(\beta \mathrm{/}3)$, and ${\lambda }_{33}=(\alpha \mathrm{/}3)(\beta \mathrm{/}3)$. These are nonnegative and sum to one. The polynomial in step 1.1 is affine in $(\alpha ,\beta )$, so it equals the convex combination with these weights of its values at $(0,0),(3,0),(0,3),(3,3)$, namely $6t(1-t)$, $3(1-t{)}^2$, $3t^2$, and $3(1-2t{)}^2$. Each of these four polynomials is nonnegative on $[0,1]$.

step 2.1 gives $p^{\mathrm{\prime }}(x)\mathrm{/}d\ge 0$ throughout the cell. If $d>0$, $p^{\mathrm{\prime }}\ge 0$; if $d<0$, $p^{\mathrm{\prime }}\le 0$. For any two points $u<v$ in the cell, [F2] transfers that derivative sign to $p(v)-p(u)$, so $p$ is monotone; step 1.1 gives its endpoint values $y_0,y_1$, so every value lies between them. If $d=0$, then $y_1=y_0$ and [F1] reduces to $p=y_0(H_{00}+H_{01})=y_0$, since $m_0=m_1=0$ and $H_{00}+H_{01}=1$.

Fix a cell and $t\in [0,1]$. Both $t$ and $1-t$ are nonnegative, so their natural powers are nonnegative by induction using the recursive definition [F10] and closure of nonnegative reals under multiplication. By [F3]–[F4], $\left(\genfrac{}{}{0px}{}{n}{k}\right)$ is a natural that is either zero or positive; [F5] sends its zero case to zero, and [F6] sends its positive case to a positive real. Thus every $B_{k,n}(t)\ge 0$.

Apply [F2] with $x=t$ and $y=1-t$. Since $t+(1-t)=1$, [F10] and induction give $1^n=1$, so the displayed Bernstein basis terms sum to one.

By [F8]–[F9], let $c_{\min }$ and $c_{\max }$ be the minimum and maximum of the finite control list. For each $k$, [F1], step 1.1, and [F7] give $c_{\min }{B}_{k,n}(t)\le c_{i,k}{B}_{k,n}(t)\le c_{\max }{B}_{k,n}(t)$: the control differences and basis weights are nonnegative, so their products are nonnegative. Summing with [F7] and using step 1.2 gives $c_{\min }\le p_i(x)\le c_{\max }$. Therefore nonnegative controls imply $p_i(x)\ge 0$. Each $p_i$ is a polynomial in $x$ after substituting the affine local coordinate, so it is continuous on its cell by [F11]. For $n_i\ge 1$, [F3], [F5], [F10], step 1.2, and [F12] give $B_{0,n_i}(0)=B_{n_i,n_i}(1)=1$ and all other endpoint basis values zero; for $n_i=0$ the sole basis term is constantly one. Thus at a shared node $p_i(x_{i+1})=c_{i,n_i}=c_{i+1,0}=p_{i+1}(x_{i+1})$. Given $\epsilon >0$, the minimum of the two positive continuity radii for the adjacent cell restrictions then gives continuity of the piecewise function at that node; at cell interiors continuity follows from [F11], and at the outer endpoints it follows from the single adjacent cell restriction. Applying the bound cellwise gives nonnegativity throughout the partition.

Each $h_i=x_{i+1}-x_i$ is positive, so every secant $d_i$ is defined and finite. In the two-node branch both slopes equal the already defined secant, and neither endpoint formula is evaluated. When $n\ge 3$, the endpoint denominators $h_0+h_1$ and $h_{n-2}+h_{n-3}$ exist and are positive.

The interior harmonic-mean branch is entered only when $d_{k-1}{d}_k>0$. Thus both secants are nonzero and have the same sign; the weights $w_1=2h_k+h_{k-1}$ and $w_2=h_k+2h_{k-1}$ are positive. The two terms $w_1\mathrm{/}{d}_{k-1}$ and $w_2\mathrm{/}{d}_k$ therefore have that same nonzero sign, so their sum is nonzero and finite. The endpoint clip and cap branches only compare finite real quantities and introduce no division. Hence every output slope is defined and finite.

At an interior node $k$ whose adjacent secants have nonpositive product, [F1] sets $m_k=0$. If their product is positive, write $r=d_k\mathrm{/}{d}_{k-1}>0$ and $W=w_1+w_2=3(h_{k-1}+h_k)$. Dividing the harmonic formula in [F1] by $d_k$ gives $m_k\mathrm{/}{d}_k=W\mathrm{/}(w_{1}r+w_2)$, which is positive and at most $W\mathrm{/}{w}_2\le 3$, because $w_2=h_k+2h_{k-1}\ge h_k+h_{k-1}$. Dividing instead by $d_{k-1}$ gives $m_k\mathrm{/}{d}_{k-1}=W\mathrm{/}(w_1+w_2\mathrm{/}r)\le W\mathrm{/}{w}_1\le 3$, since $w_1=2h_k+h_{k-1}\ge h_k+h_{k-1}$.

For a nonzero first secant, put $q=d_1\mathrm{/}{d}_0$. The raw endpoint ratio is ${\widehat{m}}_0\mathrm{/}{d}_0=1+\frac{h_0}{h_0+h_1}(1-q)$. If $q\ge 0$, this ratio is at most $1+h_0\mathrm{/}(h_0+h_1)<2$; clipping makes it zero if it is nonpositive. If $q<0$, the ratio is greater than $1$, and the cap branch reduces it to $3$ exactly when it exceeds $3$. Thus the final ratio is in $[0,3]$. If $d_0=0$, the raw endpoint value is either zero or has sign different from $d_0$, so the clip sets it to zero. The reflected argument proves the same bound at the right endpoint.

For two nodes both selected slopes equal $d_0$, so they are zero when $d_0=0$ and their ratios to $d_0$ are $1$ otherwise. For $n\ge 3$, if $d_i=0$, its interior endpoint slopes are set to zero by the nonpositive-product branch, and a first or last endpoint slope is zero by the clip rule. If $d_i\mathrm{\ne }0$, each of its endpoint slopes is either an interior slope covered by step 1.1 or an endpoint slope covered by step 1.2. These cases prove both alternatives in the statement.

By [F3], $f^{\mathrm{\prime }}$ is continuous. It cannot take both a positive and a negative value on $[a,b]$, since [F5] would then give a zero of $f^{\mathrm{\prime }}$, contradicting $\mathrm{\mid }{f}^{\mathrm{\prime }}\mathrm{\mid }\ge c$. Thus $f^{\mathrm{\prime }}$ has one sign throughout. On each grid cell [F4] gives $d_j=f^{\mathrm{\prime }}({\xi }_j)$ for some interior ${\xi }_j$, so all secants have that sign and $\mathrm{\mid }{d}_j\mathrm{\mid }\ge c$. Since $a<b$ and $h=(b-a)\mathrm{/}N>0$, the grid is admissible for [F1], and [F8] makes $[a,b]$ compact. Since $f^{\mathrm{\prime }}$, $f^{\mathrm{\prime }\mathrm{\prime }}$, and $f^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}$ are continuous by [F3], [F7] supplies $M_1,M_2,M_3\ge 0$ with $\mathrm{\mid }{f}^{\mathrm{\prime }}\mathrm{\mid }\le M_1$, $\mathrm{\mid }{f}^{\mathrm{\prime }\mathrm{\prime }}\mathrm{\mid }\le M_2$, and $\mathrm{\mid }{f}^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}\mathrm{\mid }\le M_3$ on $[a,b]$.

Fix an interior node $x_k$ and set $s=f^{\mathrm{\prime }}(x_k)$ and $q=f^{\mathrm{\prime }\mathrm{\prime }}(x_k)$. Applying [F6] with degree two at $x_k$ to $x_k-h$ and $x_k+h$ and dividing the remainders by $h$ gives $d_{k-1}=s-hq\mathrm{/}2+e_{-}$ and $d_k=s+hq\mathrm{/}2+e_{+}$, where $\mathrm{\mid }{e}_{-}\mathrm{\mid },\mathrm{\mid }{e}_{+}\mathrm{\mid }\le M_3{h}^2\mathrm{/}6$. Since $\mathrm{\mid }s\mathrm{\mid }\ge c$, for sufficiently small $h$ both secants stay bounded away from zero. The identity $1\mathrm{/}(s+u)=1\mathrm{/}s-u\mathrm{/}{s}^2+u^2\mathrm{/}(s^2(s+u))$ then gives, with a uniform constant, $1\mathrm{/}{d}_{k-1}=1\mathrm{/}s+hq\mathrm{/}(2s^2)+O(h^2)$ and $1\mathrm{/}{d}_k=1\mathrm{/}s-hq\mathrm{/}(2s^2)+O(h^2)$.

On a uniform grid the two weights in [F1] are both $3h$, so the interior slope equals $2\mathrm{/}(1\mathrm{/}{d}_{k-1}+1\mathrm{/}{d}_k)$. The first-order terms in step 1.2 cancel, giving $1\mathrm{/}{d}_{k-1}+1\mathrm{/}{d}_k=2\mathrm{/}s+O(h^2)$. Since $\mathrm{\mid }s\mathrm{\mid }\le M_1$, the leading term has magnitude at least $2\mathrm{/}{M}_1$; thus the denominator stays bounded away from zero for sufficiently small $h$, and its reciprocal yields $m_k=s+O(h^2)$ with a constant uniform in $k$.

At the left endpoint, [F6] at $a$ gives $d_0=f^{\mathrm{\prime }}(a)+h{f}^{\mathrm{\prime }\mathrm{\prime }}(a)\mathrm{/}2+O(M_3{h}^2)$ and $d_1=f^{\mathrm{\prime }}(a)+3h{f}^{\mathrm{\prime }\mathrm{\prime }}(a)\mathrm{/}2+O(M_3{h}^2)$, by expanding $f(a+h)$ and $f(a+2h)$. Therefore the raw endpoint estimate in [F1] is ${\widehat{m}}_0=(3d_0-d_1)\mathrm{/}2=f^{\mathrm{\prime }}(a)+O(M_3{h}^2)$. For small $h$ it has the sign of $f^{\mathrm{\prime }}(a)$; step 1.1 gives that same sign for $d_0,d_1$, so neither the sign clip nor the opposite-sign cap changes it. Reflecting the expansions at $b$ gives the identical $O(h^2)$ error for $m_N$.

The interior and endpoint estimates above use the bounds $M_1,M_2,M_3$ and $c^{-1}$ on the fixed interval; their constants and smallness thresholds are therefore uniform in the node index. Taking a common maximum for the error constants and a common positive minimum for the thresholds gives $C,h_0$ valid at every interior node and both endpoints. The finite list $\mathrm{\mid }{m}_j-f^{\mathrm{\prime }}(x_j)\mathrm{\mid }$, $0\le j\le N$, has a maximum by [F9], and the uniform pointwise estimate bounds it by $C{h}^2$.

By [F2], the algorithm defines one finite slope at each input node. Its output consists of Hermite pieces sharing the same value and selected slope at each common node, so [F3] gives interpolation and $C^1$ matching.

Fix a cell with $d_i\mathrm{\ne }0$. By [F5], both slopes normalized by $d_i$ lie in $[0,3]$; [F4] therefore makes this piece monotone from $y_i$ to $y_{i+1}$ and keeps it between those values.

If $d_i=0$, [F5] gives zero slopes at both endpoints, and [F4] makes that Hermite piece constant at its common endpoint value. For nondecreasing input, every $d_i\ge 0$, so each nonflat cell is nondecreasing by step 1.2, the flat cells are constant by the calculation above, and the pieces agree at nodes by step 1.1; concatenating the cells gives a nondecreasing function. The nonincreasing case is identical with all inequalities reversed.

Since $a<b$ and $h=(b-a)\mathrm{/}N>0$, the grid is admissible and $[a,b]$ is compact by [F9]. By [F7] and [F8], choose $M\ge 0$ with $\mathrm{\mid }{f}^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}\mathrm{\mid }\le M$ on $[a,b]$. Fix a cell $[x_i,x_i+h]$ and let $T_2$ be the degree-two Taylor polynomial of $f$ at $x_i$. By [F6], the residual $r=f-T_2$ satisfies $r(x_i)=r^{\mathrm{\prime }}(x_i)=0$ and, for every $x\in [x_i,x_i+h]$, $\mathrm{\mid }r(x)\mathrm{\mid }\le M{h}^3\mathrm{/}6$; at the right endpoint it also satisfies $\mathrm{\mid }r(x_i+h)\mathrm{\mid }\le M{h}^3\mathrm{/}6$ and $\mathrm{\mid }{r}^{\mathrm{\prime }}(x_i+h)\mathrm{\mid }\le M{h}^2\mathrm{/}2$. The derivative bound is [F6] applied to $f^{\mathrm{\prime }}$ with degree one.

Let $H$ be the exact-derivative Hermite piece for $f$ on this cell. The PCHIP piece $P_h$ and $H$ have identical endpoint values. By [F3], their difference is $h(m_i-f^{\mathrm{\prime }}(x_i))H_{10}(t)+h(m_{i+1}-f^{\mathrm{\prime }}(x_i+h))H_{11}(t)$. Since [F2] guarantees these PCHIP slopes are defined and $\mathrm{\mid }{H}_{10}(t)\mathrm{\mid }=t(1-t{)}^2\le 1$, $\mathrm{\mid }{H}_{11}(t)\mathrm{\mid }\le 1$, [F5] gives $\mathrm{\mid }{P}_h(x)-H(x)\mathrm{\mid }\le 2C_s{h}^3$ for every $x$ in the cell.

For this exact-derivative Hermite piece $H$, linearity of the formula [F3] and reproduction [F4] give $H=T_2+H_r$, where $H_r(t)=r(x_i+h)H_{01}(t)+h{r}^{\mathrm{\prime }}(x_i+h)H_{11}(t)$. On $0\le t\le 1$, $\mathrm{\mid }{H}_{01}(t)\mathrm{\mid }=t^2(3-2t)\le 1$ and $\mathrm{\mid }{H}_{11}(t)\mathrm{\mid }=t^2(1-t)\le 1$, so $\mathrm{\mid }{H}_r(t)\mathrm{\mid }\le (1\mathrm{/}6+1\mathrm{/}2)M{h}^3=2M{h}^3\mathrm{/}3$. Together with $\mathrm{\mid }r(x)\mathrm{\mid }\le M{h}^3\mathrm{/}6$ from step 1.1, this yields $\mathrm{\mid }H(x)-f(x)\mathrm{\mid }\le 5M{h}^3\mathrm{/}6$ for every $x$ in the cell.

Combining steps 2.1 and 1.2 gives $\mathrm{\mid }{P}_h(x)-f(x)\mathrm{\mid }\le (5M\mathrm{/}6+2C_s)h^3$ for every $x\in [a,b]$, since the cells cover the interval. The constants do not depend on $i$ or $N$; put $C=5M\mathrm{/}6+2C_s$. To justify the maximum in the statement, admissibility of the grid for [F1] gives $x_0<\cdots <x_N$, so $a<b$. By [F4], $P_h$ is $C^1$ and hence continuous on $[a,b]$; by [F7], $f$ is continuous there. Set $E(x)=\mathrm{\mid }{P}_h(x)-f(x)\mathrm{\mid }$. For $x,z\in [a,b]$, applying the triangle inequality twice gives $\mid \mathrm{\mid }{P}_h(x)-f(x)\mathrm{\mid }-\mathrm{\mid }{P}_h(z)-f(z)\mathrm{\mid }\mid \le \mathrm{\mid }(P_h(x)-f(x))-(P_h(z)-f(z))\mathrm{\mid }\le \mathrm{\mid }{P}_h(x)-P_h(z)\mathrm{\mid }+\mathrm{\mid }f(x)-f(z)\mathrm{\mid }$, so $E$ is continuous on $[a,b]$. By [F9], $[a,b]$ is compact, and it is nonempty; [F10] therefore gives a point where $E$ attains its maximum. The pointwise bound just proved bounds that maximum by $C{h}^3$.