# A uniform non-linear subdivision scheme
reproducing polynomials at any non-uniform grid

Sergio López-Ureña<EMAIL_ADDRESS>Dept. de Matemàtiques,
Universitat de València, Doctor Moliner Street 50, 46100 Burjassot, Valencia,
Spain.

###### Abstract

In this paper, we introduce a novel non-linear uniform subdivision scheme for
the generation of curves in $\mathbb{R}^{n}$, $n\geq 2$. This scheme is
distinguished by its capacity to reproduce second-degree polynomial data on
non-uniform grids without necessitating prior knowledge of the grid
specificities. Our approach exploits the potential of annihilation operators
to infer the underlying grid, thereby obviating the need for end-users to
specify such information. We define the scheme in a non-stationary manner,
ensuring that it progressively approaches a classical linear scheme as the
iteration number increases, all while preserving its polynomial reproduction
capability.

The convergence is established through two distinct theoretical methods.
Firstly, we propose a new class of schemes, including ours, for which we
establish $\mathcal{C}^{1}$ convergence by combining results from the analysis
of quasilinear schemes and asymptotically equivalent linear non-uniform non-
stationary schemes. Secondly, we adapt conventional analytical tools for non-
linear schemes to the non-stationary case, allowing us to again conclude the
convergence of the proposed class of schemes.

We show its practical usefulness through numerical examples, showing that the
generated curves are curvature continuous.

###### keywords:

Subdivision schemes, polynomial reproduction, non-uniform, non-stationary,
quasilinear, annihilation operators, asymptotically equivalent.

††journal: arXiv

##  1 Introduction

Recursive subdivision is a well-established technique for creating smooth
geometric models from discrete initial data. It has found extensive
application in Computer-Aided Design (CAD) thanks to its exceptional
performance and straightforward implementation (e.g., see Dyn [1992], Ma
[2005]).

In this paper, our focus is on generating curves in $\mathbb{R}^{n}$, $n\geq
2$. To achieve this, we will employ the following type of subdivision schemes.

######  Definition 1 (Vector-valued subdivision scheme).

For a given initial data $f^{0}\in\ell(\mathbb{Z},\mathbb{R}^{n})$, so that
$f^{0}_{i}\in\mathbb{R}^{n}$ $\forall i\in\mathbb{Z}$, more data is generated
by the recursive formula $f^{k+1}=S_{k}f^{k}$,
$k\in\mathbb{N}_{0}:=\mathbb{N}\cup\\{0\\}$, where
$S_{k}:\ell(\mathbb{Z},\mathbb{R}^{n})\longrightarrow\ell(\mathbb{Z},\mathbb{R}^{n})$
is a _subdivision operator_ defined as

$(S_{k}f)_{2i}=\Psi_{2i}^{k}(f_{i-p},\ldots,f_{i+p+1}),\quad(S_{k}f)_{2i+1}=\Psi_{2i+1}^{k}(f_{i-p},\ldots,f_{i+p+1}),\qquad
i\in\mathbb{Z},\quad f\in\ell(\mathbb{Z},\mathbb{R}^{n}),$

for some $p\in\mathbb{N}_{0}$ and some
$\Psi^{k}_{i}:(\mathbb{R}^{n})^{2p+2}\longrightarrow\mathbb{R}^{n}$,
$i\in\mathbb{Z}$. The _subdivision scheme_ is identified as the sequence of
subdivision operators $\mathcal{S}=\\{S_{k}\\}_{k\in\mathbb{N}_{0}}$. The
scheme is _linear_ iff each $\Psi^{k}_{i}$ is a linear function. The scheme is
_uniform_ iff $\Psi^{k}_{2i}=\Psi^{k}_{0}$ and $\Psi^{k}_{2i+1}=\Psi^{k}_{1}$,
$\forall i\in\mathbb{Z}$, $\forall k\in\mathbb{N}_{0}$. The scheme is
_stationary_ iff $\Psi^{k}_{i}=\Psi^{0}_{i}$, $\forall i\in\mathbb{Z}$,
$\forall k\in\mathbb{N}_{0}$. We say that the scheme is _based on a scalar-
valued subdivision scheme_ (or simply _scalar-based_) iff there exists
$\psi^{k}_{i}:\mathbb{R}^{2p+2}\longrightarrow\mathbb{R}$ such that

$\Psi^{k}_{i}(f_{i-p},\ldots,f_{i+p+1})=(\psi^{k}_{i}((f_{i-p})_{j},\ldots,(f_{i+p+1})_{j}))_{j=1}^{n},\qquad\forall
f_{i-p},\ldots,f_{i+p+1}\in\mathbb{R}^{n},\ \forall i\in\mathbb{Z},\ \forall
k\in\mathbb{N}_{0}.$

The properties of linear uniform subdivision schemes have been extensively
studied, including non-stationary schemes, which are valuable for achieving
exact geometric representations of conic sections through exponential
polynomial reproduction (see Dyn et al. [2003], Conti and Romani [2011]).

Non-uniform linear schemes have also been the subject of previous research (as
in Dyn et al. [2014]). However, their practical application in CAD may be
limited, necessitating knowledge of the underlying non-uniform grid or
requiring end-users to adjust certain shape parameters along the grid (e.g.,
see Kuijt and Van Damme [2002]). To address this issue, some scholars have
suggested estimating these parameters from existing data, either at the
initial step (as in Beccari et al. [2011]) or at each step (as in Jeong et al.
[2021], Dyn et al. [2009]), with the latter approach leading to uniform non-
linear schemes. A similar strategy has been employed in this work: we consider
a family of non-uniform linear schemes and replace the grid-dependent
parameters with non-linear functions that solely dependent on the data,
leading to a uniform non-linear scheme.

In this paper, we consider subdivision schemes as general as in Definition 1.
In the absence of some the adjectives “linear”, “uniform”, “stationary” and
“scalar-based”, we do not assume them.

In the following, we introduce some notation and basic definitions. For any
$f\in\ell(\mathbb{Z},\mathbb{R}^{n})$, we denote $\nabla f:=(\nabla
f_{i})_{i\in\mathbb{Z}}\in\ell(\mathbb{Z},\mathbb{R}^{n})$, $\nabla
f_{i}:=f_{i+1}-f_{i}\in\mathbb{R}^{n}$. For bounded sequences, we consider

$\|f\|_{\infty}:=\sup_{i\in\mathbb{Z}}|f_{i}|,\quad
f\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n}),$

for some fixed $|\cdot|$ norm in $\mathbb{R}^{n}$. Abusing the notation, we
write $\langle a,b\rangle:=\sum_{i}a_{i}b_{i}\in\mathbb{R}^{m}$, for $a$ and
$b$ having the same length, $a_{i}\in\mathbb{R}$, $b_{i}\in\mathbb{R}^{m}$,
with $m\in\mathbb{N}$.

######  Definition 2 (Convergence).

We say that a sequence
$\\{f^{k}\\}_{k\in\mathbb{N}_{0}}\subset\ell(\mathbb{Z},\mathbb{R}^{n})$
_converges uniformly_ to a continuous function
$F:\mathbb{R}\longrightarrow\mathbb{R}^{n}$ iff for any compact interval
$I\subset\mathbb{R}$,

$\lim_{k\to+\infty}\sup_{i\in\mathbb{Z}\cap 2^{k}I}|f^{k}_{i}-F(i2^{-k})|=0.$
(1)

A subdivision scheme $\mathcal{S}$ is _uniformly convergent_ iff, for any
$f^{0}\in\ell(\mathbb{Z},\mathbb{R}^{n})$, the sequence
$\\{f^{k}\\}_{k\in\mathbb{N}_{0}}$ defined recursively by $f^{k+1}=S_{k}f^{k}$
converges uniformly to a continuous function, denoted by $S^{\infty}f^{0}$,
and $S^{\infty}\neq 0$. A scheme is $\mathcal{C}^{\mu}$, $\mu\in(0,+\infty)$,
iff it is convergent and
$S^{\infty}f^{0}\in\mathcal{C}^{\mu}(\mathbb{R},\mathbb{R}^{n})$, $\forall
f^{0}\in\ell(\mathbb{Z},\mathbb{R}^{n})$ (Hölder regularity is considered
here). A scheme is $\mathcal{C}^{\mu-}$ iff it is $\mathcal{C}^{\mu-\epsilon}$
for any $\epsilon\in(0,\mu)$.

The concept of reproduction refers to the ability of a subdivision scheme to
exactly generate a specific family of functions, typically polynomials. It has
been thoroughly examined, but to the best of our knowledge, it has been barely
addressed in the context of non-uniform grids (see Levin [2003]).
Consequently, we propose the following reproduction definition, which is
compatible with non-uniform grids.

######  Definition 3 (Reproduction).

For each $k\in\mathbb{N}_{0}$, let
$\xi^{k}=(\xi^{k}_{i})_{i\in\mathbb{Z}}\subset\mathbb{R}$ be a grid, so that
$\xi^{k}_{i}<\xi^{k}_{i+1}$. Let
$\mathcal{F}\subset\mathcal{C}(\mathbb{R},\mathbb{R}^{n})$ be a set of
functions. We say that a subdivision operator $S_{k}$ _reproduces_
$\mathcal{F}$ for $(\xi^{k},\xi^{k+1})$ iff

$S_{k}F|_{\xi^{k}}=F|_{\xi^{k+1}},\qquad\forall F\in\mathcal{F},$

where $F|_{\xi^{k}}=(F(\xi^{k}_{i}))_{i\in\mathbb{Z}}$. If this holds true for
every $k\in\mathbb{N}_{0}$, we say that the subdivision scheme $\mathcal{S}$
_stepwise reproduces_ $\mathcal{F}$ for $\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$.

We say that $\mathcal{S}$ _reproduces_ $\mathcal{F}$ on $\xi^{0}$ iff there
exists some parametrization of $\mathbb{R}$, say
$\sigma:\mathbb{R}\longrightarrow\mathbb{R}$ strictly increasing, such that
the subdivision process converges to $F\circ\sigma$ when the input data is
$f^{0}=F|_{\xi^{0}}$, $\forall F\in\mathcal{F}$.

The typical definition of reproduction does not consider a parametrization
$\sigma$. This is because it is provided for uniform schemes, where $\sigma$
is implicitly selected as the identity function and $\xi^{0}$ is a uniform
(i.e., equispaced) grid.

######  Proposition 4.

Let $\mathcal{F}\subset\mathcal{C}(\mathbb{R},\mathbb{R}^{n})$ be and let
$\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$ be a uniformly convergent sequence. If
$\mathcal{S}$ stepwise reproduces $\mathcal{F}$ for
$\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$, then $\mathcal{S}$ reproduces
$\mathcal{F}$ on $\xi^{0}$.

In addition, $S^{\infty}(F|_{\xi^{0}})=F\circ\sigma$, where
$\sigma:\mathbb{R}\longrightarrow\mathbb{R}$ is the limit function of
$\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$.

###### Proof.

Consider $I\subset\mathbb{R}$ an arbitrary compact interval and
$F\in\mathcal{F}$ an arbitrary function of the space being reproduced. Given
the stepwise reproduction, we know that starting from any $f^{0}=F|_{\xi^{0}}$
the generated iterates are $f^{k}=F|_{\xi^{k}}$. Thus,

$|f^{k}_{i}-(F\circ\sigma)(i2^{-k})|=|F(\xi^{k}_{i})-F(\sigma(i2^{-k}))|.$

Given $\epsilon>0$, since $F$ is uniformly continuous on $I$, there exists
$\delta>0$ such that

$|x-y|<\delta\quad\longrightarrow\quad|F(x)-F(y)|<\epsilon.$

Since $\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$ converges uniformly to $\sigma$,
there exists $k_{0}\in\mathbb{N}_{0}$ such that

$|\xi^{k}_{i}-\sigma(i2^{-k})|<\delta,\qquad\forall i\in\mathbb{Z}\cap
2^{k}I,\quad\forall k\geq k_{0}.$

Combining these results, we get

$|f^{k}_{i}-(F\circ\sigma)(i2^{-k})|<\epsilon,\quad\forall i\in\mathbb{Z}\cap
2^{k}I,\quad\forall k\geq k_{0}.$

Hence, $\\{f^{k}\\}_{k\in\mathbb{N}_{0}}$ converges uniformly to
$F\circ\sigma$ on $I$. Since $I$ and $F$ were arbitrarily chosen, we conclude
that $\mathcal{S}$ reproduces $\mathcal{F}$ on $\xi^{0}$. ∎

Thanks to Dyn et al. [2003], Conti and Romani [2011], the algebraic conditions
for a linear uniform scheme to reproduce exponential polynomials are well-
established. For instance, to reproduce a circle, the number of sampling
points is required to define a reproducing linear non-stationary scheme. In
the linear case, the consideration of non-stationary schemes is mandatory to
reproduce exponential polynomials.

In later studies (see Donat and López-Ureña [2019], López-Ureña and Viscardi
[2021]), we investigated the possibility of defining a single non-linear but
stationary scheme with the capability to simultaneously reproduce multiple
exponential polynomial spaces. In the previous example of a circle, this
approach eliminates the need to know the number of sampling points. This type
of schemes was accomplished by leveraging the existence of annihilation
operators for exponential polynomials (see Dyn et al. [2003], Conti et al.
[2022]).

In this work, we harness annihilation operators to develop a non-linear but
uniform subdivision scheme with the ability to reproduce second-degree
polynomials even if they are sampled on non-uniform grids. This scheme can
handle many non-uniform grids without the need for prior knowledge of those
grids. No preprocessing steps are performed. This property is precisely stated
in Theorem 15.

The capability to represent parabolas in CAD applications, afforded by second-
degree polynomials reproduction, holds little significance. This research
pretends to mark an initial step towards achieving the reproduction of
exponential polynomials (and in particular, conic sections) on non-uniform
grids, even when the grid information is not known in advance. Such
reproduction capabilities can be particularly valuable in practical scenarios
where users provide a set of points, here denoted as $f^{0}$, representing a
common curve like a conic section, without specifying the underlying
parametrization or grid spacings. Consequently, these curves could be
reproduced on a wide range of non-uniform grids without burdening the user
with providing such information.

More in details, we deduce and study a subdivision scheme that is able to
reproduce the set of functions which coordinates are second-degree
polynomials,

$\Pi_{2}^{n}:=\\{F=(F_{1},\ldots,F_{n}):\mathbb{R}\longrightarrow\mathbb{R}^{n},\quad
F_{i}\in\Pi_{2},\quad\forall i=1,\ldots,n\\}.$

######  Remark 5.

From now on, let us consider a sequence of grids
$\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$, so that $\xi^{k}_{i}<\xi^{k}_{i+1}$,
$\forall i\in\mathbb{Z}$, $\forall k\in\mathbb{N}_{0}$. We measure the
uniformity of the grid $\xi^{k}$ by the following quantities:

$\alpha_{i}^{k}:=\frac{\nabla\xi^{k}_{i-1}}{\nabla\xi^{k}_{i}},\quad\beta_{i}^{k}:=\frac{\nabla\xi^{k}_{i+1}}{\nabla\xi^{k}_{i}},\qquad
i\in\mathbb{Z}.$

Observe that $\alpha_{i}^{k},\beta_{i}^{k}\in(0,+\infty)$,
$\beta_{i}^{k}=1/\alpha_{i+1}^{k}$, and that a uniform grid is characterized
by $\alpha_{i}^{k}=1$, $\forall i\in\mathbb{Z}$.

Second-degree polynomials satisfy the following annihilation property.

######  Proposition 6 (Annihilation property).

For any grid $\xi^{k}$ and for all $i\in\mathbb{Z}$, $F\in\Pi_{2}^{n}$, $n\geq
1$,

$A_{i}^{k}(F(\xi^{k}_{i})-F(\xi^{k}_{i-1}))-(F(\xi^{k}_{i+1})-F(\xi^{k}_{i}))+B_{i}^{k}(F(\xi^{k}_{i+2})-F(\xi^{k}_{i+1}))=0\in\mathbb{R}^{n},$
(2)

where

$A_{i}^{k}:=\frac{\beta_{i}^{k}+1}{\alpha_{i}^{k}(\alpha_{i}^{k}+\beta_{i}^{k}+2)},\quad
B_{i}^{k}:=\frac{\alpha_{i}^{k}+1}{\beta_{i}^{k}(\alpha_{i}^{k}+\beta_{i}^{k}+2)}.$
(3)

###### Proof.

By the linearity of (2) respect to $F$, and since $\Pi_{2}^{n}$ is a tensor
product space, it is enough to check this property on a basis of
$\Pi_{2}^{1}$, say $\\{1,x,x^{2}\\}$. It can be confirmed by simplifying the
algebraic expressions that result from replacing (3) in (2) for each
$F\in\\{1,x,x^{2}\\}$. We provide a Mathematica notebook file (see Section
Reproducibility) where this verification, and others, are performed. ∎

The Lagrange interpolation formula can be employed to readily derive a linear
non-uniform subdivision operator that reproduces $\Pi_{3}^{n}$ for
$(\xi^{k},\xi^{k+1})$. It forms the basis of our proposal for a uniform non-
linear scheme, as it did for other schemes (see Dyn et al. [2004, 2009]).

######  Definition 7 (Lagrange subdivision scheme).

Given $\\{(\alpha^{k}_{i},\beta^{k}_{i})\\}_{i\in\mathbb{Z}}$,
$\alpha^{k}_{i},\beta^{k}_{i}\in(0,+\infty)$, $k\in\mathbb{N}_{0}$, we
consider the subdivision operators

$\displaystyle L_{k}$
$\displaystyle:\ell(\mathbb{Z},\mathbb{R}^{n})\longrightarrow\ell(\mathbb{Z},\mathbb{R}^{n})$
$\displaystyle(L_{k}f)_{2i+l}$
$\displaystyle:=\Lambda_{l}[\alpha^{k}_{i},\beta^{k}_{i}](f_{i-1},f_{i},f_{i+1},f_{i+2}),\quad
l=0,1,\ i\in\mathbb{Z},$
$\displaystyle\Lambda_{0}[\alpha,\beta](f_{-1},f_{0},f_{1},f_{2})$
$\displaystyle:=\sum_{j=-1}^{2}a_{j}(\alpha,\beta)\ \cdot\
f_{j},\quad\Lambda_{1}[\alpha,\beta](f_{-1},f_{0},f_{1},f_{2}):=\sum_{i=-1}^{2}a_{1-j}(\beta,\alpha)\
\cdot\ f_{j},$ (4)
$\displaystyle\begin{split}a_{j}:&(0,+\infty)^{2}\longrightarrow\mathbb{R},\quad
j=-1,0,1,2,\\\
a_{-1}(\alpha,\beta)&=-\frac{3(4\beta+3)}{64\left(\alpha^{2}+\alpha\right)(\alpha+\beta+1)},\quad
a_{0}(\alpha,\beta)=\frac{3(4\alpha+1)(4\beta+3)}{64\alpha(\beta+1)},\\\
a_{2}(\alpha,\beta)&=-\frac{3(4\alpha+1)}{64\left(\beta^{2}+\beta\right)(\alpha+\beta+1)},\quad
a_{1}(\alpha,\beta)=\frac{(4\alpha+1)(4\beta+3)}{64\beta(\alpha+1)}.\end{split}$
(5)

We have defined $L_{k}$ using $\Lambda_{0}$ and $\Lambda_{1}$ to emphasize
that the dependence on the grid and the iteration $k$ is exclusively given by
$\alpha^{k}_{i}$ and $\beta^{k}_{i}$. Our proposal of subdivision scheme,
which will be defined in Section 2, consists of replacing
$\alpha^{k}_{i},\beta^{k}_{i}$ by some non-linear functions depending only on
the data $f$.

The reason for considering $L_{k}$ to reproduce $\Pi_{3}^{n}$ instead of
$\Pi_{2}^{n}$, our target, is that $L_{k}$ is uniquely determined by that
4-point stencil and its reproducing condition. The choice of 4 points is due
to it being the minimum number required to locally use Proposition 6 for
computing $\alpha^{k}_{i}$ and $\beta^{k}_{i}$ from
$f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2}$, as we will see in Section 2.

######  Proposition 8.

Let $\xi^{k}$ be a grid and define $\xi^{k+1}$ by applying the Chaikin’s
subdivision operator, i.e.,

$\xi^{k+1}_{2i}=\frac{3}{4}\xi^{k}_{i}+\frac{1}{4}\xi^{k}_{i+1},\qquad\xi^{k+1}_{2i+1}=\frac{1}{4}\xi^{k}_{i}+\frac{3}{4}\xi^{k}_{i+1},\qquad
i\in\mathbb{Z}.$

The subdivision operator in Definition 7 reproduces $\Pi^{n}_{3}$ for
$(\xi^{k},\xi^{k+1})$.

###### Proof.

The Chaikin’s operator is strictly monotonicity preserving, implying that
$\xi^{k+1}$ is strictly increasing because $\xi^{k}$ is. $L_{k}$ was computed
using the Lagrange interpolation formula, ensuring reproduction. The steps
taken to derive it are outlined here, and for a direct check of the
reproduction property, we refer to the Mathematica file (see Section
Reproducibility).

First, given the data $f_{j}=F(\xi_{j})$, $j=-1,0,1,2$, attached to some nodes
$\xi_{-1},\xi_{0},\xi_{1},\xi_{2}$, there is a unique cubic polynomial
interpolating them, that can be expressed as a linear combination of the
Lagrange basis functions:

$F(x)=\sum_{i=-1}^{2}f_{i}\prod_{\begin{subarray}{c}j=-1\\\ j\neq
i\end{subarray}}^{2}\frac{x-\xi_{j}}{\xi_{i}-\xi_{j}}.$

We are only concerned with two values of $x$,
$x=\frac{3}{4}\xi_{0}+\frac{1}{4}\xi_{1}$ and
$x=\frac{1}{4}\xi_{0}+\frac{3}{4}\xi_{1}$, for deriving the subdivision rules.
Here, we present the first choice of $x$, as the outcome for the second one
should mirror the first, with the order of data and nodes reversed.

To simplify the expression, we denote $\alpha:=\nabla\xi_{-1}/\nabla\xi_{0}$,
$\beta:=\nabla\xi_{1}/\nabla\xi_{0}$, and

$F(x)=\sum_{i=-1}^{2}a_{i}f_{i},\qquad
a_{i}:=\prod_{\begin{subarray}{c}j=-1\\\ j\neq
i\end{subarray}}^{2}\frac{x-\xi_{j}}{\xi_{i}-\xi_{j}}.$

To express $a_{i}$ in terms of $\alpha$ and $\beta$, we compute the following:

$\displaystyle\frac{\xi_{2}-\xi_{-1}}{\nabla\xi_{0}}$
$\displaystyle=\frac{\nabla\xi_{1}+\nabla\xi_{0}+\nabla\xi_{-1}}{\nabla\xi_{0}}=\beta+1+\alpha,\quad\frac{\xi_{1}-\xi_{-1}}{\nabla\xi_{0}}=1+\alpha,\quad\frac{\xi_{2}-\xi_{0}}{\nabla\xi_{0}}=\beta+1,$
$\displaystyle\frac{x-\xi_{0}}{\nabla\xi_{0}}$
$\displaystyle=\frac{-\frac{1}{4}\xi_{0}+\frac{1}{4}\xi_{1}}{\nabla\xi_{0}}=\frac{1}{4},\quad\frac{x-\xi_{-1}}{\nabla\xi_{0}}=\frac{x-\xi_{0}+\nabla\xi_{-1}}{\nabla\xi_{0}}=\frac{1}{4}+\alpha,\quad\frac{x-\xi_{1}}{\nabla\xi_{0}}=\frac{\frac{3}{4}\xi_{0}-\frac{3}{4}\xi_{1}}{\nabla\xi_{0}}=-\frac{3}{4},$
$\displaystyle\frac{x-\xi_{2}}{\nabla\xi_{0}}$
$\displaystyle=\frac{x-\xi_{1}-\nabla\xi_{1}}{\nabla\xi_{0}}=-\frac{3}{4}-\beta.$

Then, it is straightforward to see that

$\displaystyle a_{-1}$
$\displaystyle=\prod_{j=0}^{2}\frac{x-\xi_{j}}{\xi_{-1}-\xi_{j}}=\prod_{j=0}^{2}\frac{(x-\xi_{j})/\nabla\xi_{0}}{(\xi_{-1}-\xi_{j})/\nabla\xi_{0}}=\frac{\frac{1}{4}(-\frac{3}{4})(-\frac{3}{4}-\beta)}{-\alpha(-1-\alpha)(-1-\alpha-\beta)}=-\frac{3(3+4\beta)}{64\alpha(1+\alpha)(1+\alpha+\beta)}.$

We can proceed similarly with the other coefficients to obtain (5). ∎

######  Remark 9.

It can be verified that $\sum_{i=-1}^{2}a_{i}(\alpha,\beta)=1$, a consequence
of the reproduction of constant functions, which are included in
$\Pi_{3}^{n}$. Additionally, the scheme is scalar-based. These two properties
imply that $L_{k}$ commutes with affine transformations, i.e.,

$L_{k}(\mathcal{A}_{M,c}f)=\mathcal{A}_{M,c}(L_{k}f),\quad\mathcal{A}_{M,c}f:=(Mf_{i}+c)_{i\in\mathbb{Z}},\quad\forall
M\in\mathbb{R}^{n\times n},\ \forall c\in\mathbb{R}^{n},\ \forall
f\in\ell(\mathbb{Z},\mathbb{R}^{n}).$

This property is useful in CAD applications and it is easily verified by
direct computation.

######  Remark 10.

For uniform (i.e. equispaced) grids, where $\alpha^{k}_{i}=\beta^{k}_{i}=1$,
the Lagrange subdivision operators coincide with the one studied in Dyn et al.
[2004], which rules are

$\displaystyle\Lambda_{0}[1,1](f_{-1},f_{0},f_{1},f_{2})=\frac{1}{128}(-7f_{-1}+105f_{0}+35f_{1}-5f_{2}),$
$\displaystyle\Lambda_{1}[1,1](f_{-1},f_{0},f_{1},f_{2})=\frac{1}{128}(-5f_{-1}+35f_{0}+105f_{1}-7f_{2}).$

This linear scheme is $\mathcal{C}^{2}$.

We have introduced all the necessary components to establish a single non-
linear and uniform subdivision scheme capable of reproducing $\Pi_{2}^{n}$ on
a variety of grids. This scheme will be defined next, in Section 2. It is
important to note that the scheme will be defined in a non-stationary manner,
which is not a requisite for achieving reproduction capability but proves to
be highly advantageous for demonstrating convergence. We adopt two different
approaches to prove convergence, as we find them interesting from a
theoretical development perspective. In Section 3, we combine the results of
_asymptotic equivalence_ from Dyn et al. [2014] for linear non-uniform non-
stationary schemes with the _quasilinear_ concept of Cohen et al. [2003]. This
synthesis gives rise to a new class of subdivision schemes that demonstrate
convergence and, with the addition of an easily met requirement, achieve
$\mathcal{C}^{1}$ continuity. The last is applicable to our scheme. In Section
4, we analyse convergence using standard analytical tools for non-linear
schemes, demonstrating their applicability in the non-stationary case by
easily adapting them. Section 5 illustrates the performance of our proposal in
several numerical examples, highlighting its ability to produce curvature-
continuous curves. Our work concludes with final remarks and considerations
for the future in Section 6.

##  2 Definition of the new subdivision scheme. Reproduction

######  Definition 11.

For a fixed $\rho\geq 0$ and $n\geq 2$, we consider the uniform non-linear
non-stationary subdivision scheme defined by the rules
$\Psi^{k}_{0},\Psi^{k}_{1}:(\mathbb{R}^{n})^{4}\longrightarrow\mathbb{R}^{n}$,

$\displaystyle\begin{split}\Psi^{k}_{l}(f_{-1},f_{0},f_{1},f_{2})&:=\Lambda_{l}[((\alpha^{k},\beta^{k})\circ(A,B))(\nabla
f_{-1},\nabla f_{0},\nabla f_{1})](f_{-1},f_{0},f_{1},f_{2}),\quad
l=0,1,\end{split}$ (6)

where $\Lambda_{0},\Lambda_{1}$ are defined in (4) and

$\displaystyle\alpha^{k},\beta^{k}:(0,$
$\displaystyle+\infty)^{2}\longrightarrow[(1-2^{-k}\rho)^{-1},1+2^{-k}\rho],\qquad
A,B:(\mathbb{R}^{n})^{3}\longrightarrow(0,+\infty)$
$\displaystyle\begin{split}\alpha^{k}(A,B)&:=\max\left\\{(1+2^{-k}\rho)^{-1},\min\left\\{1+2^{-k}\rho,\frac{1}{A+\frac{\sqrt{A(A+1)B(B+1)}}{B+1}}\right\\}\right\\},\\\
\beta^{k}(A,B)&:=\alpha^{k}(B,A),\end{split}$ (7)
$\displaystyle\begin{split}A(\delta_{-1},\delta_{0},\delta_{1})&:=\begin{cases}\frac{1}{2},&\text{if
}\text{num}\cdot\text{den}=0,\\\
\text{abs}(\frac{\text{num}}{\text{den}}),&\text{otherwise}.\end{cases},\quad
B(\delta_{-1},\delta_{0},\delta_{1}):=A(\delta_{1},\delta_{0},\delta_{-1})\\\
\text{num}&=\det\begin{pmatrix}(\delta_{-1})^{T}\delta_{0}&(\delta_{-1})^{T}\delta_{1}\\\
(\delta_{1})^{T}\delta_{0}&|\delta_{1}|_{2}^{2}\end{pmatrix},\quad\text{den}=\det\begin{pmatrix}|\delta_{-1}|_{2}^{2}&(\delta_{-1})^{T}\delta_{1}\\\
(\delta_{-1})^{T}\delta_{1}&|\delta_{1}|_{2}^{2}\end{pmatrix},\end{split}$ (8)

where $\delta_{-1},\delta_{0},\delta_{1}\in\mathbb{R}^{n}$ are written as
columns vectors.

######  Remark 12.

For $n=2$, observe that for any
$\delta_{-1},\delta_{0},\delta_{1}\in\mathbb{R}^{2}$, we have that

num
$\displaystyle=\det\left(\begin{pmatrix}\delta_{-1}&\delta_{1}\end{pmatrix}^{T}\begin{pmatrix}\delta_{0}&\delta_{1}\
\end{pmatrix}\right)=\det\begin{pmatrix}\delta_{-1}&\delta_{1}\end{pmatrix}\det\begin{pmatrix}\delta_{0}&\delta_{1}\end{pmatrix},$
den
$\displaystyle=\det\left(\begin{pmatrix}\delta_{-1}&\delta_{1}\end{pmatrix}^{T}\begin{pmatrix}\delta_{-1}&\delta_{1}\
\end{pmatrix}\right)=\det\begin{pmatrix}\delta_{-1}&\delta_{1}\end{pmatrix}^{2}.$

Hence, $A$ can be rewritten as

$A(\delta_{-1},\delta_{0},\delta_{1}):=\begin{cases}\frac{1}{2},&\text{if
}\prod_{i\neq j}\det\begin{pmatrix}\delta_{i}&\delta_{j}\end{pmatrix}=0,\\\
\text{abs}\left(\frac{\det\begin{pmatrix}\delta_{0}&\delta_{1}\end{pmatrix}}{\det\begin{pmatrix}\delta_{-1}&\delta_{1}\end{pmatrix}}\right),&\text{otherwise}.\end{cases}$
(9)

Note that $\alpha^{k}$ and $\beta^{k}$ involve truncating a number to ensure
it always falls between two values dependent on $2^{-k}\rho$. This approach
has been carefully designed to maintain the reproduction capability intact (as
demonstrated in Theorem 15) while simultaneously ensuring the
$\mathcal{C}^{1}$ convergence of the scheme (as established in Corollary 34).
Observe that $\lim_{k\to+\infty}(\alpha^{k},\beta^{k})=(1,1)$, so that the
proposed subdivision operators converge to the ones of the Lagrange scheme for
uniform grids (see Remark 10). In particular, for $\rho=0$, the proposed
operator is exactly that one, which is linear and possesses good properties.
With this, the asymptotic properties of our scheme, including convergence and
regularity, are linked with the mentioned linear scheme.

Note that, for every $\rho>0$, a distinct scheme is obtained, which is non-
stationary (since $\alpha^{k}$ and $\beta^{k}$ vary among iterations), uniform
(only two subdivision rules per iteration are required, which do not consider
any grid information), and non-linear. As shown in Section 5, $\rho$ is
related with the curve flexibility between points. We have observed that
pleasing shapes emerge for $\rho$ values within the range of $[1,4]$, but this
interval can be even wider, depending on the stylistic preferences of the end-
user.

Observe that if $\alpha^{k},\beta^{k}$ are defined without truncation, their
image set is $(0,+\infty)$, which is unbounded and induces practical problems.
They could be truncated with a constant number, independent of $k$, which
would lead to a stationary scheme. But we did not find any benefit of
considering such stationary version, which introduces several difficulties in
proving convergence and reproduction simultaneously.

Before delving deeper into the analysis, it is essential to consider potential
issues related to the scheme definition.

######  Proposition 13.

For any $\rho\geq 0$, the subdivision scheme in Definition 11 is well-defined.

###### Proof.

For any $\delta_{-1},\delta_{0},\delta_{1}\in\mathbb{R}^{n}$, we have that
$A,B$ are well-defined and positive. Then, $\alpha^{k}$ and $\beta^{k}$, which
are well-defined, can be composed with them and, by definition,

$((\alpha^{k},\beta^{k})\circ(A,B))(\delta_{-1},\delta_{0},\delta_{1})\in[(1+2^{-k}\rho)^{-1},1+2^{-k}\rho]^{2}.$
(10)

The functions $a_{-1},a_{0},a_{1}$ are being evaluated in positive numbers, so
the denominators are positive and we conclude that $\Psi^{k}_{0},\Psi^{k}_{1}$
are well-defined. ∎

The concepts presented in this paper can be applied to create similar schemes
for different purposes. For example, an interpolatory subdivision scheme would
be very practical. However, in our study, we opted to explore a ‘dual’ scheme,
which, in essence, generates new data more akin to the Chaikin’s scheme. The
reason for this choice is revealed in the following lemma, where one can
observe that this particular refinement scheme offers a highly advantageous
property. We leverage this property with the truncated-style definition of
$\alpha^{k}$ and $\beta^{k}$, which becomes key to proving $\mathcal{C}^{1}$
convergence and reproduction simultaneously.

######  Lemma 14.

Let $\xi^{k},\xi^{k+1}$ be two grids fulfilling

$\xi^{k+1}_{2i}=\frac{3}{4}\xi^{k}_{i}+\frac{1}{4}\xi^{k}_{i+1},\qquad\xi^{k+1}_{2i+1}=\frac{1}{4}\xi^{k}_{i}+\frac{3}{4}\xi^{k}_{i+1},\qquad\forall
i\in\mathbb{Z}.$ (11)

Then, for any $\rho\geq 0$,

$(1+\rho)^{-1}\leq\alpha_{i}^{k}\leq 1+\rho,\ \forall
i\in\mathbb{Z}\quad\longrightarrow\quad\left(1+\frac{\rho}{2}\right)^{-1}\leq\alpha_{i}^{k+1}\leq
1+\frac{\rho}{2},\ \forall i\in\mathbb{Z}.$

###### Proof.

We are going to prove that
$\alpha_{2i}^{k+1}=\frac{\nabla\xi^{k+1}_{2i+1}}{\nabla\xi^{k+1}_{2i}}$ and
$\alpha_{2i+1}^{k+1}=\frac{\nabla\xi^{k+1}_{2i+2}}{\nabla\xi^{k+1}_{2i+1}}$
are in between these bounds, $\forall i\in\mathbb{Z}$. First, observe that

$(1+\rho)^{-1}\leq x\leq 1+\rho,\quad\longrightarrow\quad(1+\rho)^{-1}\leq
x^{-1}\leq 1+\rho.$ (12)

Second, we deduce that

$\nabla\xi^{k+1}_{2i+1}=\xi^{k+1}_{2i+2}-\xi^{k+1}_{2i+1}=\frac{1}{4}\xi^{k}_{i+2}-\frac{1}{4}\xi^{k}_{i},\quad\nabla\xi^{k+1}_{2i}=\xi^{k+1}_{2i+1}-\xi^{k+1}_{2i}=\frac{1}{2}\xi^{k}_{i+1}-\frac{1}{2}\xi^{k}_{i}.$

Then,

$\displaystyle\frac{\nabla\xi^{k+1}_{2i+1}}{\nabla\xi^{k+1}_{2i}}-1$
$\displaystyle=\frac{\frac{1}{4}\xi^{k}_{i+2}-\frac{1}{4}\xi^{k}_{i}}{\frac{1}{2}\xi^{k}_{i+1}-\frac{1}{2}\xi^{k}_{i}}-1=\frac{1}{2}\left(\frac{\xi^{k}_{i+2}-\xi^{k}_{i}}{\xi^{k}_{i+1}-\xi^{k}_{i}}-2\right)=\frac{1}{2}\left(\frac{\xi^{k}_{i+2}-\xi^{k}_{i}-\xi^{k}_{i+1}+\xi^{k}_{i}}{\xi^{k}_{i+1}-\xi^{k}_{i}}-1\right)$
$\displaystyle=\frac{1}{2}\left(\frac{\nabla\xi^{k}_{i+1}}{\nabla\xi^{k}_{i}}-1\right),$
$\displaystyle\frac{\nabla\xi^{k+1}_{2i+1}}{\nabla\xi^{k+1}_{2i+2}}-1$
$\displaystyle=\frac{\frac{1}{4}\xi^{k}_{i+2}-\frac{1}{4}\xi^{k}_{i}}{\frac{1}{2}\xi^{k}_{i+2}-\frac{1}{2}\xi^{k}_{i+1}}-1=\frac{1}{2}\left(\frac{\xi^{k}_{i+2}-\xi^{k}_{i}}{\xi^{k}_{i+2}-\xi^{k}_{i+1}}-2\right)=\frac{1}{2}\left(\frac{\xi^{k}_{i+2}-\xi^{k}_{i}-\xi^{k}_{i+2}+\xi^{k}_{i+1}}{\xi^{k}_{i+2}-\xi^{k}_{i+1}}-1\right)$
$\displaystyle=\frac{1}{2}\left(\frac{\nabla\xi^{k}_{i}}{\nabla\xi^{k}_{i+1}}-1\right).$

Applying to the last two formulas the fact that

$x\in[(1+\rho)^{-1},1+\rho]\quad\longrightarrow\quad
1+\frac{1}{2}(x-1)\in\left[\frac{2+\rho}{2+2\rho},1+\frac{\rho}{2}\right]\subset\left[\left(1+\frac{\rho}{2}\right)^{-1},1+\frac{\rho}{2}\right]$

(the latter can be checked in the Mathematica file, see Section
Reproducibility), we arrive to

$\left(1+\frac{\rho}{2}\right)^{-1}\leq\frac{\nabla\xi^{k+1}_{2i+1}}{\nabla\xi^{k+1}_{2i}},\frac{\nabla\xi^{k+1}_{2i+1}}{\nabla\xi^{k+1}_{2i+2}}\leq
1+\frac{\rho}{2}.$

Using (12), we arrive to the claim. ∎

We are now prepared to demonstrate that the scheme reproduces non-degenerate
parabolas (i.e., the image of $F\in\Pi_{2}^{n}$ is not confined to a line).
Additionally, we establish that the scheme generates lines when the data is
collinear.

######  Theorem 15.

Let us consider $\mathcal{L}=\\{F\in\Pi_{2}^{n}:F(\mathbb{R})\text{ is
contained in a line}\\}$. For any $\rho\geq 0$, the subdivision scheme in
Definition 11 reproduces $\mathcal{F}=\Pi_{2}^{n}\setminus\mathcal{L}$ on any
$\xi^{0}$ grid fulfilling

$(1+\rho)^{-1}\leq\alpha^{0}_{i}\leq 1+\rho,\quad i\in\mathbb{Z}.$

Moreover, if $f^{0}$ lies on a line, then $f^{k}$ lies on the same line.

###### Proof.

Using Proposition 4, it is sufficient to prove the stepwise reproduction for
some sequence of grids $\\{\xi^{k}\\}_{k\in\mathbb{N}_{0}}$. We define the
grids using the Chaikin’s subdivision scheme, as in (11).

Denote $f^{k}=F|_{\xi^{k}}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$,
$F\in\mathcal{F}$. Let us focus on four consecutive data
$f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2}$, that indeed fulfil

$F(\xi^{k}_{j})=f^{k}_{j},\quad j=i-1,\ldots,i+2.$ (13)

We have to prove that

$\Psi^{k}_{l}(f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2})=F(\xi^{k+1}_{2i+l}),\quad
l=0,1.$

Considering that the proposed scheme is primarily built upon the Lagrange
scheme of Definition 7, which stepwise reproduces $\Pi_{2}^{n}$ according to
Proposition 8, our focus shifts to confirming that the non-linear scheme
exactly computes $\alpha^{k}_{i}$ and $\beta^{k}_{i}$ from the data. This
requires us to verify that

$(\alpha^{k}_{i},\beta^{k}_{i})=((\alpha^{k},\beta^{k})\circ(A,B))(\nabla
f^{k}_{i-1},\nabla f^{k}_{i},\nabla f^{k}_{i+1}),$

because it would imply

$\Psi^{k}_{l}(f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2})=\Lambda_{l}[\alpha^{k}_{i},\beta^{k}_{i}](f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2})=(L_{k}f^{k})_{2i+l}=F(\xi^{k+1}_{2i+l}),\qquad
l=0,1.$

Recall the annihilation property for $\Pi^{2}_{n}$ in Proposition 6:

$A_{i}^{k}\nabla f^{k}_{i-1}-\nabla f^{k}_{i}+B_{i}^{k}\nabla
f^{k}_{i+1}=0\in\mathbb{R}^{n},$ (14)

where

$A_{i}^{k}=\frac{\beta_{i}^{k}+1}{\alpha_{i}^{k}(\alpha_{i}^{k}+\beta_{i}^{k}+2)},\qquad
B_{i}^{k}=\frac{\alpha_{i}^{k}+1}{\beta_{i}^{k}(\alpha_{i}^{k}+\beta_{i}^{k}+2)}.$
(15)

To ensure that $\alpha^{k}_{i},\beta^{k}_{i}$ can be recovered from
$A^{k}_{i},B^{k}_{i}$, we observe that the function

$(\alpha,\beta)\mapsto\left(\frac{\beta+1}{\alpha(\alpha+\beta+2)},\frac{\alpha+1}{\beta(\alpha+\beta+2)}\right)$

is bijective from $(0,+\infty)^{2}$ to itself, with inverse

$(A,B)\mapsto\left(\frac{1}{A+\frac{\sqrt{A(A+1)B(B+1)}}{B+1}},\frac{1}{B+\frac{\sqrt{A(A+1)B(B+1)}}{A+1}}\right),\quad
A,B>0.$

See the Mathematica file (Section Reproducibility) for a verification. Hence,
$\alpha^{k}_{i}=\alpha^{k}(A^{k}_{i},B^{k}_{i})$ and
$\beta^{k}_{i}=\beta^{k}(A^{k}_{i},B^{k}_{i})$, with $\alpha^{k},\beta^{k}$ as
defined in (7), provided that

$(1+2^{-k}\rho)^{-1}\leq\alpha^{k}_{i},\beta^{k}_{i}\leq 1+2^{-k}\rho.$

We can ensure the last condition for $k=0$, by hypothesis, and for any $k>0$,
by applying Lemma 14 recursively.

The last point of this proof consists in checking that

$A^{k}_{i}=A(\nabla f^{k}_{i-1},\nabla f^{k}_{i},\nabla f^{k}_{i+1}),\quad
B^{k}_{i}=B(\nabla f^{k}_{i-1},\nabla f^{k}_{i},\nabla f^{k}_{i+1}),$ (16)

with $A,B$ as defined in (8). It is done by solving the $n\times 2$ linear
system (14), with unknowns $A_{i}^{k},B_{i}^{k}$. Despite the system size, we
know that there exists a solution, the one in (15) and it is unique provided
that $\text{rank}\begin{pmatrix}\nabla f^{k}_{i-1}&\nabla
f^{k}_{i+1}\end{pmatrix}=2$. In that case, the solution can be written as:

$\begin{pmatrix}A^{k}_{i}\\\
B^{k}_{i}\end{pmatrix}=\left(\begin{pmatrix}\nabla f^{k}_{i-1}&\nabla
f^{k}_{i+1}\end{pmatrix}^{T}\begin{pmatrix}\nabla f^{k}_{i-1}&\nabla
f^{k}_{i+1}\end{pmatrix}\right)^{-1}\begin{pmatrix}\nabla f^{k}_{i-1}&\nabla
f^{k}_{i+1}\end{pmatrix}^{T}\nabla f^{k}_{i},$

where we are writing $\nabla f^{k}_{i}\in\mathbb{R}^{n}$ as column vectors.
Performing the matrices multiplications:

$\begin{pmatrix}A^{k}_{i}\\\ B^{k}_{i}\end{pmatrix}=\begin{pmatrix}|\nabla
f^{k}_{i-1}|_{2}^{2}&(\nabla f^{k}_{i-1})^{T}\nabla f^{k}_{i+1}\\\ (\nabla
f^{k}_{i-1})^{T}\nabla f^{k}_{i+1}&|\nabla
f^{k}_{i+1}|_{2}^{2}\end{pmatrix}^{-1}\begin{pmatrix}(\nabla
f^{k}_{i-1})^{T}\nabla f^{k}_{i}\\\ (\nabla f^{k}_{i+1})^{T}\nabla
f^{k}_{i}\end{pmatrix}.$

Applying the Cramer’s formula, we arrive to the definition of $A,B$ in (8),
provided that $\text{num}\cdot\text{den}\neq 0$. Since
$\text{rank}\begin{pmatrix}\nabla f^{k}_{i-1}&\nabla
f^{k}_{i+1}\end{pmatrix}=2$, then $\text{den}\neq 0$. By (15), we know that
$A^{k}_{i},B^{k}_{i}$ are positive numbers, implying that $\text{num}\neq 0$.

Now we address the case $\text{rank}\begin{pmatrix}\nabla f^{k}_{i-1}&\nabla
f^{k}_{i+1}\end{pmatrix}<2$, that is, when $\nabla f^{k}_{i-1}$ and $\nabla
f^{k}_{i+1}$ are proportional. According to (14),

$\nabla f^{k}_{i}=A_{i}^{k}\nabla f^{k}_{i-1}+B_{i}^{k}\nabla f^{k}_{i+1}.$
(17)

Thus, we have that the three vectors are proportional, meaning that
$f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2}$ are co-linear. This implies
that $F\in\mathcal{L}$. Since $\Psi^{k}_{l}$, $l=0,1$, consists of applying
$\Lambda_{l}[\alpha,\beta]$ for some $\alpha,\beta$ and, by Remark 9, it
always performs convex combinations of the data, then
$f^{k+1}_{2i},f^{k+1}_{2i+1}$ lies on the same line as
$f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2}$. ∎

######  Remark 16.

According to (17), if $f^{k}=F|_{\xi^{k}}$, $F\in\Pi_{2}^{n}$, then the three
vectors $\nabla f^{k}_{i-1},\nabla f^{k}_{i},\nabla f^{k}_{i+1}$ are linearly
dependent and, as a result, $f^{k}_{i-1},f^{k}_{i},f^{k}_{i+1},f^{k}_{i+2}$
lie on the same plane. If $F\in\Pi_{2}^{n}\setminus\mathcal{L}$, then $\nabla
f^{k}_{i-1},\nabla f^{k}_{i+1}$ are linearly independent and they form a basis
of that plane. In addition, since $A_{i}^{k},B_{i}^{k}>0$, then $\nabla
f^{k}_{i}$ falls into the first quadrant defined by this basis. This is
illustrated in Figure 1-(a). We can also see that $\nabla f^{k}_{i}$ cannot
fall into any other quadrant when the data describe a parabola.

$\nabla f^{k}_{i-1}$$\nabla f^{k}_{i}$$\nabla f^{k}_{i+1}$ | $\nabla f^{k}_{i-1}$$\nabla f^{k}_{i}$$\nabla f^{k}_{i+1}$ | $\nabla f^{k}_{i-1}$$\nabla f^{k}_{i}$$\nabla f^{k}_{i+1}$ | $\nabla f^{k}_{i-1}$$\nabla f^{k}_{i}$$\nabla f^{k}_{i+1}$
---|---|---|---
(a) | (b) | (c) | (d)
Figure 1: Four cases where $\nabla f^{k}_{i-1},\nabla f^{k}_{i+1}$ are
linearly independent and $\nabla f^{k}_{i}$ falls into the first, second,
third or fourth quadrant (from the left to the right).

If $\nabla f^{k}_{i-1},\nabla f^{k}_{i+1}$ are linearly dependent, then
$\nabla f^{k}_{i}$ is also proportional to them by (17), and the four data
points are collinear. This happens when $F\in\mathcal{L}$, for instance,
$F(x)=(x^{2},x^{2})$. We would also like to emphasize that if $\nabla
f^{k}_{i}=0$, then $\nabla f^{k}_{i-1},\nabla f^{k}_{i+1}$ are linearly
dependent with opposite directions according to (17).

It becomes impossible to determine $F(\xi^{k+1}_{2i})$ and
$F(\xi^{k+1}_{2i+1})$ in the last case: Since $\nabla f^{k}_{i-1},\nabla
f^{k}_{i+1}$ are linearly dependent, the rank of the system (14) is one, and
only one equation (one coordinate) is sufficient to represent the system
information. This is insufficient for reproducing the data, and neither this
scheme nor any other one can reproduce the data without knowing
$\alpha^{k}_{i},\beta^{k}_{i}$, as we demonstrate in the following example:
Let us consider the family of polynomials

$\left\\{p(x)=\frac{x(x-1)}{\alpha(\alpha+1)}:\alpha>0\right\\}\subset\Pi_{2}^{1}.$

All of them fulfil that

$p(-\alpha)=1,\quad p(0)=0,\quad p(1)=0,\quad p(\alpha)=1,$

for some $\alpha>0$. Consequently, provided the four data points $1,0,0,1,$ it
becomes impossible to determine $\alpha$ and, by extension, $p(\frac{1}{4})$
and $p(\frac{3}{4})$. For an example in $\Pi_{2}^{n}$, $n\geq 2$, just
consider the proposed family of polynomials in each coordinate.

######  Remark 17.

As a consequence of Theorem 15, we can prove that the scheme is not
$\mathcal{C}^{2}$, for any $\rho>0$. For $\rho=0$, it is the linear scheme
proposed in Dyn et al. [2004] which is $\mathcal{C}^{2}$. Let us consider the
initial grid $\xi^{0}_{i}=(1+\rho)^{i}$, $i\in\mathbb{Z}$, so that

$\alpha^{0}_{i}=\frac{\nabla\xi^{0}_{i+1}}{\nabla\xi^{0}_{i}}=1+\rho.$

Let us consider the initial data $f^{0}=F|_{\xi^{0}}$ with $F(x)=(x,x^{2})$.

On the one hand, by Theorem 15 and Proposition 4,
$S^{\infty}(f^{0})=F\circ\sigma$, where $\sigma$ is the limit function of the
Chaikin’s scheme applied to $\xi^{0}$, which is $\mathcal{C}^{2-}$. More
precisely, $\sigma$ is a quadratic spline whose second derivative is piecewise
constant, with jumps at the integers. Thus,
$S^{\infty}(f^{0})=F\circ\sigma=(\sigma,\sigma^{2})$ is also
$\mathcal{C}^{2-}$.

On the other hand, the limit curve is precisely a parabola, which is
infinitely smooth. This example underscores that the perception of
‘smoothness’ depends on whether we treat the limit as a function or a curve.
It emphasizes that a smooth curve does not necessarily imply a smooth function
(its parametrization). Conversely, the cycloidal curve serves as a classical
illustration which demonstrates that a smooth parametrization does not
guarantee a smooth curve.

We show in the numerical experiments of Section 5 that the scheme generates
curves with continuous curvature. We hope that this could be theoretically
established in the future. For the moment, we prove in Section 3 that it is
$\mathcal{C}^{1}$ for any initial data and any choice of $\rho\geq 0$.

This example also demonstrates that the regularity of the scheme is
constrained by the regularity of the implicit mesh refinement. So, replacing
the Chaikin’s scheme with some other smoother and monotonicity preserving
scheme may be advantageous.

In the following proposition, we establish that the scheme commutes with
similarity transformations, including rotations, translations, scalings, and
reflections. This property is desirable for CAD applications, ensuring
consistent outcomes when applying these standard transformations to the data.
It is noteworthy that any similarity transformation is an affine
transformation, so the next result is weaker than the one in Remark 9.

######  Proposition 18.

For any $\rho\geq 0$, the subdivision scheme in Definition 11 commutes with
similarity transformations.

###### Proof.

Similarly to Remark 9, we need to prove that if a similarity transformation is
applied to the data, i.e.

$\bar{f}_{i}:=\mu Mf_{i}+c,\qquad\forall i=-1,0,1,2,$ (18)

being $\mu\in\mathbb{R}$, $c\in\mathbb{R}^{n}$ and $M\in\mathbb{R}^{n\times
n}$ an orthogonal matrix, then

$\Psi_{l}^{k}(\bar{f}_{-1},\bar{f}_{0},\bar{f}_{1},\bar{f}_{2})=\mu
M\Psi^{k}_{l}(f_{-1},f_{0},f_{1},f_{2})+c,\quad l=0,1.$ (19)

The case $\mu=0$ is trivial since $\Lambda_{l}[\alpha,\beta]$ reproduces
constants for any $\alpha,\beta>0$, $l=0,1$. For $\mu\neq 0$, first observe
that

$\nabla\bar{f}_{i}=\mu M\nabla f_{i},\qquad\forall i=-1,0,1,2.$

Then, looking at (8), it is straightforwardly verified that

$(A,B)(\nabla f_{-1},\nabla f_{0},\nabla
f_{1})=(A,B)(\nabla\bar{f}_{-1},\nabla\bar{f}_{0},\nabla\bar{f}_{1}).$

Consequently,

$\Lambda_{l}[((\alpha^{k},\beta^{k})\circ(A,B))(\nabla f_{-1},\nabla
f_{0},\nabla
f_{1})]=\Lambda_{l}[((\alpha^{k},\beta^{k})\circ(A,B))(\nabla\bar{f}_{-1},\nabla\bar{f}_{0},\nabla\bar{f}_{1})],\quad
l=0,1.$

Then, (19) is achieved provided that $\Lambda_{l}[\alpha,\beta]$, $l=0,1$,
also fulfil (19) for any $\alpha,\beta>0$, which certainly holds true by
Remark 9. ∎

##  3 Convergence by asymptotic equivalence and quasilinearity

In this section, we introduce a broad class of non-linear, non-stationary,
uniform subdivision schemes, which includes the proposed scheme. We
demonstrate that these schemes are $\mathcal{C}^{1}$ by utilizing two
concepts: their quasilinearity (refer to Cohen et al. [2003]) and their
asymptotic equivalence (refer to Dyn et al. [2014]) to linear, stationary,
uniform schemes. To apply the findings from these papers, we assume that the
data is bounded, i.e., $f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$,
$S_{k}:\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})\longrightarrow\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$,
$k\in\mathbb{N}_{0}$. This is a common constraint in subdivision theory, which
is readily met in practice due to the finite nature of the data.

Let us introduce each of these concepts.

######  Definition 19 (Quasilinear subdivision scheme. Based on Definition 1
of Cohen et al. [2003]).

A subdivision operator $S$ is _data dependent_ if for each
$f\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$ there exists a linear
subdivision operator $S[f]$ such that $Sf=S[f]f.$ A subdivision scheme is
_quasilinear_ if each subdivision operator is data dependent.

######  Remark 20.

According to the Euler’s theorem, any function
$F:\mathbb{R}^{n}\longrightarrow\mathbb{R}$ that is 1-homogeneous and
continuously differentiable fulfils $F(x)=\langle\nabla F(x),x\rangle$. This
can be applied to the subdivision rules to deduce that any 1-homogeneous
subdivision scheme that admits generalized Jacobian (see Donat et al. [2017]
for further information) is quasilinear. Thus, quasilinearity is a very common
property among the non-linear subdivision schemes in the literature.

######  Remark 21.

Let $\mathcal{S}=\\{S_{k}\\}_{k\in\mathbb{N}_{0}}$ be a quasilinear
subdivision scheme. For a given
$f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$, consider
$f^{k+1}=S_{k}f^{k}$, $k\in\mathbb{N}_{0}$. $S_{k}[f^{k}]$ denotes the linear
subdivision operator such that $f^{k+1}=S_{k}[f^{k}]f^{k}$. We consider the
linear (and possibly non-uniform) subdivision scheme
$\mathcal{S}[f^{0}]:=\\{S_{k}[f^{k}]\\}_{k\in\mathbb{N}_{0}}$. Observe that,
for each $f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$, $\mathcal{S}$ and
$\mathcal{S}[f^{0}]$ generate exactly the same sequences
$\\{f^{k}\\}_{k\in\mathbb{N}_{0}}$. We will exploit this fact to examine the
smoothness of $\mathcal{S}$ through $\mathcal{S}[f^{0}]$. For this analysis,
we can apply the results for linear schemes, particularly those found in Dyn
et al. [2014].

######  Proposition 22.

Let $\mathcal{S}=\\{S_{k}\\}_{k\in\mathbb{N}_{0}}$ be a quasilinear
subdivision scheme and let $m\in\mathbb{N}_{0}$ be. If $\mathcal{S}[f^{0}]$ is
$\mathcal{C}^{m}$, $\forall f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$,
then $\mathcal{S}$ is $\mathcal{C}^{m}$.

###### Proof.

For each $f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$, the subdivision
schemes $\mathcal{S}$ and $\mathcal{S}[f^{0}]$ produce exactly the same
sequence $\\{f^{k}\\}_{k\in\mathbb{N}_{0}}$. Since $\mathcal{S}[f^{0}]$
converges to a $\mathcal{C}^{m}$ function, $\mathcal{S}$ also converges to
that $\mathcal{C}^{m}$ function when starting from $f^{0}$. ∎

We now revisit the property of asymptotic equivalence, which was previously
examined for linear, non-uniform, non-stationary subdivision schemes in Dyn et
al. [2014]. This property has proven to be effective in the analysis of
subdivision schemes that reproduce extended Chebyshev systems. The findings
presented in that paper are exclusively concerned with linear scalar-valued
subdivision schemes, i.e., $n=1$. However, by applying these results to each
component individually, they also hold true for linear vector-valued scalar-
based subdivision schemes.

Similarly to Dyn et al. [2014], given a linear scalar-based subdivision scheme
$\mathcal{S}=\\{S_{k}\\}_{k\in\mathbb{N}_{0}}$, we denote by
$b^{i,k}\in\ell(\mathbb{Z},\mathbb{R})$, $i\in\mathbb{Z}$,
$k\in\mathbb{N}_{0}$, some _masks_ that allow to compute the refinement as
follows:

$(S_{k}f^{k})_{i}=\sum_{j\in\mathbb{Z}}b^{i,k}_{i-2j}f^{k}_{j},\quad
i\in\mathbb{Z},\ f^{k}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n}).$

Observe that only half of the coefficients of each mask are used (the even or
the odd ones, depending on the parity of $i$), so that the unused half of
$b^{i,k}$ can be defined in several ways, for instance, in some convenient
manner for a specific purpose (as in Theorem 30).

The next two properties link the convergence and regularity of two linear
subdivision schemes.

######  Definition 23 (Asymptotically equivalent. Based on Section 3.2 of Dyn
et al. [2014]).

Let
$\mathcal{S}=\\{S_{k}\\}_{k\in\mathbb{N}_{0}},\tilde{\mathcal{S}}=\\{\tilde{S}_{k}\\}_{k\in\mathbb{N}_{0}}$
be two linear scalar-based subdivision schemes. They are _asymptotically
equivalent_ if

$\sum_{k=0}^{\infty}\sup_{i\in\mathbb{Z}}\\{\|b^{i,k}-{\tilde{b}}^{i,k}\|\\}<\infty,$

where $\|\cdot\|$ is any norm for the masks.

######  Definition 24 (Property A. Based on Definition 4.1 of Dyn et al.
[2014]).

Consider the Laurent polynomials

$b^{i,k}(z):=\sum_{j\in\mathbb{Z}}b^{i,k}_{j}z^{j},\quad
d^{i,k}_{m}(z):=\sum_{j=0}^{m}(-1)^{j}\binom{m}{j}z^{j}b^{i-j,k}(z),\quad
m\in\mathbb{N},\ i\in\mathbb{Z},\ k\in\mathbb{N}_{0}.$

A linear scalar-based scheme $\mathcal{S}$ satisfies _property A of order $m$_
if

$\sum_{k=0}^{\infty}2^{k(m-r)}\left|\frac{\partial^{r}d^{i,k}_{m}}{\partial
z^{r}}(\pm 1)\right|<\infty,\qquad 0\leq r<m,\quad i\in\mathbb{Z}.$

Now, we propose a straightforward result that links Property A with the
reproduction of constants.

######  Lemma 25.

If $\sum_{j\in\mathbb{Z}}b^{i,k}_{2j}=\sum_{j\in\mathbb{Z}}b^{i,k}_{1+2j}=1$,
$\forall i\in\mathbb{Z}$, then $d^{i,k}_{m}(\pm 1)=0$, $\forall
i\in\mathbb{Z}$, $\forall m\in\mathbb{N}$.

###### Proof.

The hypothesis implies that the symbols of the masks fulfil that
$b^{i,k}(1)=2$ and $b^{i,k}(-1)=0$. Then

$\displaystyle d^{i,k}_{m}(1)$
$\displaystyle=\sum_{j=0}^{m}(-1)^{j}\binom{m}{j}\cdot 1\cdot 2=0,\quad
d^{i,k}_{m}(-1)=\sum_{j=0}^{m}(-1)^{j}\binom{m}{j}(-1)^{j}\cdot 0=0.$

∎

######  Theorem 26 (Proposition 3.2 of Dyn et al. [2014] and the latter
comment).

Let $\mathcal{S},\tilde{\mathcal{S}}$ be two asymptotically equivalent
subdivision schemes, being $\tilde{\mathcal{S}}$ uniform and stationary. If
$\tilde{\mathcal{S}}$ converges, then $\mathcal{S}$ does.

######  Theorem 27 (Theorem 4.2 of Dyn et al. [2014]).

Let $\mathcal{S},\tilde{\mathcal{S}}$ be two asymptotically equivalent
subdivision schemes, being $\tilde{\mathcal{S}}$ uniform, stationary and
$\mathcal{C}^{m}$, for some $m\in\mathbb{N}$. If $\mathcal{S}$ satisfies
property A of order $m$, then $\mathcal{S}$ is $\mathcal{C}^{m}$.

Now, we show how to establish a connection between the concepts of
quasilinearity and asymptotic equivalence.

######  Corolary 28.

Let $\mathcal{S}$ be a quasilinear subdivision scheme. Let
$\tilde{\mathcal{S}}$ be a convergent, linear, uniform, stationary and scalar-
based subdivision scheme. If $\mathcal{S}[f^{0}]$ is scalar-based and
asymptotically equivalent to $\tilde{\mathcal{S}}$, $\forall
f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$, then $\mathcal{S}$ is
convergent. In addition, if $\mathcal{S}[f^{0}]$ satisfies property $A$ of
order $m$, $\forall f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$, for
some $m\in\mathbb{N}$, and $\tilde{\mathcal{S}}$ is $\mathcal{C}^{m}$, then
$\mathcal{S}$ is $\mathcal{C}^{m}$.

###### Proof.

First, apply Theorem 26 (and Theorem 27) to obtain that $\mathcal{S}[f^{0}]$
is convergent (and $\mathcal{C}^{m}$), $\forall
f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$. Second, apply Proposition
22 to deduce that $\mathcal{S}$ is convergent (and $\mathcal{C}^{m}$) as well.
∎

Next, we propose a broad class of subdivision schemes, ours included, and
apply Corollary 28 to achieve the convergence of any scheme within this class.
We believe this could be useful in the study of future non-linear non-
stationary subdivision schemes.

######  Definition 29 (Class B).

A subdivision scheme $\mathcal{S}$ belongs to the _class B_ iff it is uniform
and there exist $q\in\mathbb{N}$, $D\subset\mathbb{R}^{q}$, $0\in D$,
$\pi_{k}:(\mathbb{R}^{n})^{2p+2}\longrightarrow D$,
$\vec{\Lambda}_{0},\vec{\Lambda}_{1}:D\longrightarrow\mathbb{R}^{2p+2}$,
$\rho,C>0$ such that

$\displaystyle\Psi^{k}_{l}(f_{-p},\ldots,f_{p+1})$
$\displaystyle=\langle\vec{\Lambda}_{l}(\pi_{k}(f_{-p},\ldots,f_{p+1})),(f_{-p},\ldots,f_{p+1})\rangle,$
$\displaystyle l=0,1,\ \forall f_{j}\in\mathbb{R}^{n},\ j=-p,\ldots,i+1,$
$\displaystyle|\pi_{k}(f_{-p},\ldots,f_{p+1})|_{\mathbb{R}^{q}}$
$\displaystyle\leq 2^{-k}\rho,$ $\displaystyle\forall f_{j}\in\mathbb{R}^{n},\
j=-p,\ldots,i+1,$
$\displaystyle|\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0)|_{\mathbb{R}^{2p+2}}$
$\displaystyle\leq C|\pi|_{\mathbb{R}^{q}},$ $\displaystyle\forall\pi\in D,\
l=0,1,$

for some vector norms $|\cdot|_{\mathbb{R}^{q}},|\cdot|_{\mathbb{R}^{2p+2}}$.
We say that $\mathcal{S}$ is _asymptotically equivalent_ to a linear uniform
stationary subdivision scheme, that we denote $\mathcal{S}^{*}$, which rules
are given by the coefficients $\vec{\Lambda}_{l}(0)$, $l=0,1$:

$\Psi^{*}_{l}(f_{-p},\ldots,f_{p+1})=\langle\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle.$

We demand that $\mathcal{S}^{*}$ is convergent.

Next theorem justifies that the term _asymptotically equivalent_ in last
definition is suitable.

######  Theorem 30.

If $\mathcal{S}$ belongs to class B, then

  1. 1.

$\mathcal{S}$ is quasilinear

  2. 2.

The schemes $\mathcal{S}[f^{0}]$ and $\mathcal{S}^{*}$ are scalar-based,
$\forall f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$.

  3. 3.

For a given $f^{k}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$, the masks of
$S_{k}[f^{k}]$, say $\\{b^{i,k}\\}_{i\in\mathbb{Z}}$, can be chosen as

$(b^{2i,k}_{l-2j})_{j=-p}^{p+1}=\vec{\Lambda}_{l}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),\qquad
l=0,1,$ (20)

with $b^{2i,k}_{j}=0$ if $j<-2p-2$ or $j>1+2p$, and $b^{2i+1,k}=b^{2i,k}$,
$\forall i\in\mathbb{Z}$.

  4. 4.

$\mathcal{S}[f^{0}]$ is asymptotically equivalent to $\mathcal{S}^{*}$,
$\forall f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$.

###### Proof.

It is quasilinear since $S_{k}f^{k}=S_{k}[f^{k}]f^{k}$ where $S_{k}[f^{k}]$ is
precisely

$(S_{k}[f^{k}]g)_{2i+l}=\langle\vec{\Lambda}_{l}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),(g_{i-p},\ldots,g_{i+p+1})\rangle,\quad
l=0,1,\ g\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n}).$ (21)

Observe that the subdivision rules of $\mathcal{S}[f^{0}]$ and
$\mathcal{S}^{*}$ consists of linear combinations (with scalar coefficients)
of vector data. Thus, they are scalar-based.

Let us consider an arbitrary but fixed
$f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$. Using the notation in
Remark 21, we need to define $b^{i,k}$, $i\in\mathbb{Z}$, such that

$\langle\vec{\Lambda}_{l}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),(g_{i-p},\ldots,g_{i+p+1})\rangle=(S_{k}[f^{k}]g)_{2i+l}=\sum_{j\in\mathbb{Z}}b^{2i+l,k}_{l+2(i-j)}g_{j}=\sum_{j\in\mathbb{Z}}b^{2i+l,k}_{l-2j}g_{i+j},\quad
l=0,1.$

Hence,

$(b^{2i+l,k}_{l-2j})_{j=-p}^{p+1}:=\vec{\Lambda}_{l}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})).$

We are omitting the dependence of $b^{2i+l,k}$ on $f^{k}$. For $l=0$, this is

$b^{2i,k}_{-2j}=\vec{\Lambda}_{0}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),\quad
j=-p,\ldots,p+1.$

Since the odd positions of $b^{2j,k}$ are not used, we can freely define them
as

$b^{2i,k}_{1-2j}:=\vec{\Lambda}_{1}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),\quad
j=-p,\ldots,p+1.$

We can proceed analogously with $l=1$, obtaining that

$b^{2i+1,k}_{1-2j}=\vec{\Lambda}_{1}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),\quad
j=-p,\ldots,p+1,$

and defining

$b^{2i+1,k}_{-2j}:=\vec{\Lambda}_{0}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),\quad
j=-p,\ldots,p+1.$

By defining the unused coefficients in this manner, we have got (20) and
$b^{2i,k}=b^{2i+1,k},\ \forall i\in\mathbb{Z}.$

Denote by $b^{*}$ the mask of the linear, uniform and stationary scheme
$\mathcal{S}^{*}$, which is

$(b^{*}_{l-2j})_{j=-p}^{p+1}=\vec{\Lambda}_{l}(0,\ldots,0),\qquad l=0,1.$

We are going to prove that $\mathcal{S}[f^{0}]$ is asymptotically equivalent
to $\mathcal{S}^{*}$. Note that

$\|b^{2i,k}-b^{*}\|_{\infty}=\|b^{2i+1,k}-b^{*}\|_{\infty}=\max_{l=0,1}\\{|\vec{\Lambda}_{l}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1}))-\vec{\Lambda}_{l}(0,\ldots,0)|_{\infty}\\}.$

By hypothesis,

$\|b^{2i,k}-b^{*}\|_{\infty}=\|b^{2i+1,k}-b^{*}\|_{\infty}\leq
C|\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})|_{\infty}\leq 2^{-k}C\rho.$

Finally, this implies that

$\sum_{k=0}^{\infty}\sup_{i\in\mathbb{Z}}\\{\|b^{i,k}-b^{*}\|\\}\leq\sum_{k=0}^{\infty}2^{-k}C\rho<\infty.$

∎

The constants $\rho,C$ could be _absorbed_ by the vector norms and, thus,
removed from the definition of class B. We keep them for the sake of clarity.

######  Corolary 31.

Any scheme belonging to class B is convergent.

###### Proof.

Apply Theorem 30 and Corollary 28. ∎

######  Proposition 32.

For any $\rho\geq 0$, the subdivision scheme in Definition 11 belongs to class
B.

###### Proof.

Taking $q=2$, $p=1$,

$\displaystyle\pi_{k}(f_{-1},f_{0},f_{1},f_{2})$
$\displaystyle=((\alpha^{k},\beta^{k})\circ(A,B))(\nabla f_{-1},\nabla
f_{0},\nabla f_{1})-(1,1)$ (22)
$\displaystyle\begin{split}\vec{\Lambda}_{0}(\pi)&=(a_{-1},a_{0},a_{1},a_{2})(\pi_{1}+1,\pi_{2}+1),\\\
\vec{\Lambda}_{1}(\pi)&=(a_{2},a_{1},a_{0},a_{-1})(\pi_{2}+1,\pi_{1}+1).\end{split}$
(23)

Observe that we are using function evaluations to shorten the notation in last
equations.

The values that $\alpha^{k},\beta^{k}$ can reach are contained in
$[(1+2^{-k}\rho)^{-1},1+2^{-k}\rho]$ by definition and, as a consequence,

$|\alpha^{k}(A,B)-1|,|\beta^{k}(A,B)-1|\leq 2^{-k}\rho,\qquad\forall A,B>0.$

Hence, $|\pi_{k}(f_{-1},f_{0},f_{1},f_{2})|_{\infty}\leq 2^{-k}\rho$,
$D=[(1+2^{-k}\rho)^{-1}-1,2^{-k}\rho]^{2}$ and $(0,0)\in D$.

Observe that $a_{i}$ are Lipschitz operators in
$[(1+2^{-k}\rho)^{-1},1+2^{-k}\rho]^{2}$, then $\vec{\Lambda}_{l}$, $l=0,1$,
is Lipschitz in $D$. Thus, taking as $C$ the Lipschitz constant we have that

$|\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0,0)|_{\infty}\leq
C|\pi|_{\infty},\quad\pi\in D,\ l=0,1.$

Finally, as commented in Remark 10, we have that $\mathcal{S}^{*}$ is
convergent. ∎

Note that if $\mathcal{S}$ belongs to class B, it may not be scalar-based, as
is the case with the proposed scheme. However, as demonstrated in Theorem 30,
$\mathcal{S}[f^{0}]$ is scalar-based for any
$f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$.

######  Theorem 33.

Let $\mathcal{S}$ be belonging to class B such that
$\langle\vec{\Lambda}_{l}(\pi),(1,\ldots,1)\rangle=1,\,\forall\pi\in D,\
l=0,1.$ Then $\mathcal{S}$ reproduces constants. In addition, if
$\mathcal{S}^{*}$ is $\mathcal{C}^{1}$, then $\mathcal{S}$ is
$\mathcal{C}^{1}$.

###### Proof.

Let $f^{0}=(\ldots,c,c,c,\ldots)$, $c\in\mathbb{R}^{n}$, be some constant
sequence. We will prove by induction on $k\in\mathbb{N}_{0}$ that
$f^{k}=f^{0}$. Using the definition of class B and the hypothesis, we have
that

$\displaystyle f^{k+1}_{2i+l}$
$\displaystyle=\Psi^{k}_{l}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})=\langle\vec{\Lambda}_{l}(\pi_{k}(c,\ldots,c)),(c,\ldots,c)\rangle$
$\displaystyle=c\
\langle\vec{\Lambda}_{l}(\pi_{k}(c,\ldots,c)),(1,\ldots,1)\rangle=c,$
$\displaystyle l=0,1,\ i\in\mathbb{Z}.$

To prove that $\mathcal{S}$ is $\mathcal{C}^{1}$, by Corollary 28, we have
prove that $\mathcal{S}[f^{0}]$ satisfies property A of order 1, $\forall
f^{0}\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$. Let us consider that
$f^{0}$ is fixed and $f^{k}$ is computed by $\mathcal{S}$. We denote by
$\\{b^{i,k}\\}_{i\in\mathbb{Z}}$ the masks of $S_{k}[f^{k}]$. We have to prove
that

$\sum_{k=0}^{\infty}2^{k}\left|d^{i,k}_{1}(\pm 1)\right|<\infty,\quad
i\in\mathbb{Z}.$

By (20), the odd and even positions of the masks of $b^{i,k}$ sums 1, since

$\sum_{j\in\mathbb{Z}}b^{2i,k}_{l-2j}=\sum_{j=-p}^{p+1}b^{2i,k}_{l-2j}=\langle(b^{2i,k}_{l-2j})_{j=-p}^{p+1},(1,\ldots,1)\rangle=\langle\vec{\Lambda}_{l}(\pi_{k}(f^{k}_{i-p},\ldots,f^{k}_{i+p+1})),(1,\ldots,1)\rangle=1,\qquad
l=0,1,$

and $b^{2i,k}=b^{2i+1,k}$. By Lemma 25, we have that $d^{i,k}_{1}(\pm 1)=0$,
so that the property A is clearly satisfied. ∎

######  Corolary 34.

For any $\rho\geq 0$, the subdivision scheme in Definition 11 is
$\mathcal{C}^{1}$.

###### Proof.

We check that fulfils Theorem 33. It is of class B by Proposition 32 and the
coordinates of $\vec{\Lambda}_{l}(\pi)$ sum 1 by (23) and by Remark 9. ∎

##  4 Convergence by non-linear subdivision theory

The truncation procedure in (7) have been intentionally designed it to make
the scheme non-stationary and ‘asymptotically equivalent’ to a linear scheme.
This property simplifies the process of demonstrating convergence, in contrast
to the potential complexity involved in proving the convergence of the
stationary version of the scheme (consisting of removing the truncation
procedure). The latter has a large non-linear expression that might be
challenging to manage using conventional analytical tools for non-linear
schemes.

In addition, the subdivision operators of the proposed scheme belong to the
class of non-linear operators studied in Amat et al. [2011]. Subsequently, we
adapt [Amat et al., 2011, Theorem 1] to our non-stationary context.

######  Theorem 35.

Let
$\mathcal{S}=\\{S_{k}:\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})\longrightarrow\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})\\}_{k\in\mathbb{N}_{0}}$
be a subdivision scheme such that

$S_{k}=T+\mathcal{F}_{k}\circ\nabla,\qquad k\in\mathbb{N}_{0},$ (24)

where $\mathcal{T}=\\{T\\}_{k\in\mathbb{N}_{0}}$ is a convergent, linear and
stationary subdivision operator and
$\mathcal{F}_{k}:\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})\longrightarrow\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$
is a non-linear operator. If there exists $M>0$, $c\in(0,1)$,
$k_{0},K\in\mathbb{N}_{0}$ such that, $\forall k\geq k_{0}$,

$\displaystyle\|\mathcal{F}_{k}\delta\|_{\infty}$ $\displaystyle\leq
M\|\delta\|_{\infty},$
$\displaystyle\forall\delta\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n}),$ (25)
$\displaystyle\|\nabla S_{k+K}\cdots S_{k}f\|_{\infty}$ $\displaystyle\leq
c\|\nabla f\|_{\infty},$ $\displaystyle\forall
f\in\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n}),$ (26)

then the subdivision scheme converges.

###### Proof.

This result was proven in Amat et al. [2011] for scalar-valued ($n=1$) and
stationary ($S_{k}=S_{0}$, $\forall k\in\mathbb{N}_{0}$, $k_{0}=0$)
subdivision schemes. However, the proof would be fairly the same for $n>1$ and
the non-stationary case. In fact, the theorem in Amat et al. [2011] rely on
the original results in Daubechies et al. [2004], where the non-stationary
case and $n=2$ was considered.

Furthermore, the initial iterations have no impact on the convergence of a
subdivision scheme, as we can regard $f^{k_{0}}$ as the initial data. For the
subsequent iterations, the hypotheses remain valid. ∎

######  Theorem 36.

Let $\mathcal{S}$ be belonging to class B such that
$\langle\vec{\Lambda}_{l}(\pi),(1,\ldots,1)\rangle=1,\,\forall\pi\in D,\
l=0,1.$ Suppose that there exists
$\pi^{[1]}_{k}:(\mathbb{R}^{n})^{2p+1}\longrightarrow D$ such that
$\pi_{k}(f_{-p},\ldots,f_{p+1})=\pi^{[1]}_{k}(\nabla f_{-p},\ldots,\nabla
f_{p})$. Then, $\mathcal{S}$ fulfils the hypothesis of Theorem 35 (thus, it
converges).

###### Proof.

First, we prove that $\mathcal{S}$ can be written as in (24). By hypothesis,
$\langle\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0),(1,\ldots,1)\rangle=0,\,\forall\pi\in
D,\ l=0,1.$ Hence,

$\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0)\in\operatorname*{span}\left\\{(1,1,\ldots,1,1)\right\\}^{\perp}=\operatorname*{span}\left\\{(-1,1,0,0,\ldots),(0,-1,1,0,\ldots),\ldots,(\ldots,0,-1,1)\right\\}.$

Then, there exist some coefficients
$\vec{\Lambda}^{[1]}_{l}(\pi)=((\vec{\Lambda}^{[1]}_{l}(\pi))_{j})_{j=-p}^{p}\in\mathbb{R}^{2p+1}$
such that

$\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0)=\sum_{j=-p}^{p}(\vec{\Lambda}^{[1]}_{l}(\pi))_{j}(0,\ldots,\overset{(j)}{-1},1,0,\ldots,0),\qquad
l=0,1.$

This means that

$\displaystyle\langle\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle$
$\displaystyle=\sum_{j=-p}^{p}(\vec{\Lambda}^{[1]}_{l}(\pi))_{j}\langle(0,\ldots,\overset{(j)}{-1},1,0,\ldots,0),(f_{-p},\ldots,f_{p+1})\rangle$
$\displaystyle=\sum_{j=-p}^{p}(\vec{\Lambda}^{[1]}_{l}(\pi))_{j}\nabla
f_{j}=\langle\vec{\Lambda}^{[1]}_{l}(\pi),(\nabla f_{-p},\ldots,\nabla
f_{p})\rangle.$

Now, for $l=0,1$, by Definition 29 (class B), by the latter equation and by
hypothesis, in this order:

$\displaystyle\Psi^{k}_{l}(f_{-p},\ldots,f_{p+1})$
$\displaystyle=\langle\vec{\Lambda}_{l}(\pi_{k}(f_{-p},\ldots,f_{p+1})),(f_{-p},\ldots,f_{p+1})\rangle$
$\displaystyle=\langle\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle+\langle\vec{\Lambda}_{l}(\pi_{k}(f_{-p},\ldots,f_{p+1}))-\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle$
$\displaystyle=\langle\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle+\langle\vec{\Lambda}^{[1]}_{l}(\pi_{k}(f_{-p},\ldots,f_{p+1})),(\nabla
f_{-p},\ldots,\nabla f_{p})\rangle$
$\displaystyle=\langle\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle+\langle\vec{\Lambda}^{[1]}_{l}(\pi^{[1]}_{k}(\nabla
f_{-p},\ldots,\nabla f_{p})),(\nabla f_{-p},\ldots,\nabla f_{p})\rangle.$

Thus, the subdivision operator can be written as
$S_{k}=T+\mathcal{F}_{k}\circ\nabla$, where $T$ coincides with the subdivision
operator of $\mathcal{S}^{*}$ (its coefficients are given by
$\vec{\Lambda}_{l}(0)$, $l=0,1$) and
$\mathcal{F}_{k}:\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})\longrightarrow\ell_{\infty}(\mathbb{Z},\mathbb{R}^{n})$
is defined as

$\displaystyle(\mathcal{F}_{k}(\delta))_{2i+l}$
$\displaystyle:=\langle\vec{\Lambda}^{[1]}_{l}(\pi^{[1]}_{k}(\delta_{i-p},\ldots,\delta_{i+p})),(\delta_{i-p},\ldots,\delta_{i+p})\rangle,$
$\displaystyle l=0,1.$ (27)

Second, we prove that (25) is satisfied. By definition of class B, taking the
vector norm $|\cdot|_{\mathbb{R}^{q}}:=|\cdot|_{\infty}$ and considering

$(f_{-p},\ldots,f_{p+1}):=(0,\delta_{-p},\delta_{-p}+\delta_{-p+1},\ldots,\sum_{j=-p}^{p}\delta_{j}),$
(28)

so that $(\nabla f_{-p},\ldots,\nabla f_{p})=(\delta_{-p},\ldots,\delta_{p})$,
we have that

$\displaystyle|\pi^{[1]}_{k}(\delta_{-p},\ldots,\delta_{p})|_{\infty}=|\pi^{[1]}_{k}(\nabla
f_{-p},\ldots,\nabla
f_{p})|_{\infty}=|\pi_{k}(f_{-p},\ldots,f_{p+1})|_{\infty}\leq 2^{-k}\rho.$
(29)

That is, $\pi^{[1]}_{k}$ is bounded and converges uniformly to 0 when
$k\to\infty$. We can also deduce that

$\displaystyle\langle\vec{\Lambda}^{[1]}_{l}(\pi),(\delta_{-p},\ldots,\delta_{p})\rangle=\langle\vec{\Lambda}^{[1]}_{l}(\pi),(\nabla
f_{-p},\ldots,\nabla
f_{p})\rangle=\langle\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0),(f_{-p},\ldots,f_{p+1})\rangle.$
(30)

Observe that the last expression is an equality between vectors in
$\mathbb{R}^{n}$. Now we only pay attention to one of the coordinates, the
first one, for example. For each $j\in\\{-p,\ldots,p\\}$, we choose the data
such that $(\delta_{j})_{1}=1$ and $(\delta_{i})_{1}=0$, $i\neq j$. By (28),

$((f_{-p})_{1},\ldots,(f_{p+1})_{1})=(0,0,\ldots,0,\overset{(j)}{1},1,\ldots,1)$

and, by (30),

$\displaystyle|(\vec{\Lambda}^{[1]}_{l}(\pi))_{j}|$
$\displaystyle=|(\langle\vec{\Lambda}^{[1]}_{l}(\pi),(\delta_{-p},\ldots,\delta_{p})\rangle)_{1}|=|(\langle\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0),f_{-p},\ldots,f_{p+1}\rangle)_{1}|$
$\displaystyle=|\langle\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0),((f_{-p})_{1},\ldots,(f_{p+1})_{1})\rangle|=|\langle\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0),(0,0,\ldots,1,1\ldots)\rangle|$
$\displaystyle\leq|\vec{\Lambda}_{l}(\pi)-\vec{\Lambda}_{l}(0)|_{1}\leq
C|\pi|_{\infty}.$

Observe that we have conveniently considered in the last line that
$|\cdot|_{\mathbb{R}^{2p+2}}:=|\cdot|_{1}$. Then
$|\vec{\Lambda}^{[1]}_{l}(\pi)|_{\infty}\leq C|\pi|_{\infty}$, $l=0,1$. With
this, (29) and (27), we conclude that

$\displaystyle|(\mathcal{F}_{k}(\delta))_{2i+l}|$
$\displaystyle\leq|\vec{\Lambda}^{[1]}_{l}(\pi^{[1]}_{k}(\delta_{i-p},\ldots,\delta_{i+p}))|_{1}\|(\delta_{i-p},\ldots,\delta_{i+p})\|_{\infty}$
$\displaystyle\leq(2p+1)|\vec{\Lambda}^{[1]}_{l}(\pi^{[1]}_{k}(\delta_{i-p},\ldots,\delta_{i+p}))|_{\infty}\|(\delta_{i-p},\ldots,\delta_{i+p})\|_{\infty}$
$\displaystyle\leq(2p+1)C|\pi^{[1]}_{k}(\delta_{i-p},\ldots,\delta_{i+p})|_{\infty}\|(\delta_{i-p},\ldots,\delta_{i+p})\|_{\infty}\leq(2p+1)C2^{-k}\rho\|(\delta_{i-p},\ldots,\delta_{i+p})\|_{\infty}.$

That is, $\|\mathcal{F}_{k}\delta\|_{\infty}\leq M_{k}\|\delta\|_{\infty}$
with $M_{k}=(2p+1)C2^{-k}\rho\leq(2p+1)C\rho$. In addition, we conclude that
the sequences of operators $\\{\mathcal{F}_{k}\\}_{k\in\mathbb{N}_{0}}$
converges uniformly to zero.

Finally, we prove (26). Observe that,

$\nabla\circ S_{k}=\nabla\circ
T+\nabla\circ\mathcal{F}_{k}\circ\nabla=T^{[1]}\circ\nabla+\nabla\circ\mathcal{F}_{k}\circ\nabla=(T^{[1]}+\nabla\circ\mathcal{F})\circ\nabla=S_{k}^{[1]}\circ\nabla,$

where $T^{[1]}$ is referred to as the _difference operator_ of $T$, and
$S_{k}^{[1]}:=T^{[1]}+\nabla\circ\mathcal{F}_{k}$ turns out to be the
difference operator of $S_{k}$.

It is well known (see Dyn [1992]) that for a linear, stationary, uniform and
convergent scheme, as the one defined by $T$, there exists $K\in\mathbb{N}$
such that

$\|T^{[1]}\circ\overset{(K)}{\cdots}\circ T^{[1]}\|_{\infty}=c<1.$

Using the uniform convergence of $\mathcal{F}_{k}$ to zero, we deduce that

$\displaystyle\lim_{k\to+\infty}\|\nabla S_{k+K}\cdots S_{k}f\|_{\infty}$
$\displaystyle=\lim_{k\to+\infty}\|S^{[1]}_{k+K}\cdots S^{[1]}_{k}\nabla
f\|_{\infty}=\lim_{k\to+\infty}\|(T^{[1]}+\nabla\circ\mathcal{F}_{k+K})\cdots(T^{[1]}+\nabla\circ\mathcal{F}_{k})\nabla
f\|_{\infty}$ $\displaystyle=\|T^{[1]}\circ\overset{(K)}{\cdots}\circ
T^{[1]}\nabla f\|_{\infty}\leq c\|\nabla f\|_{\infty}.$

Hence, for $c^{\prime}=(1+c)/2<1$, there exists $k_{0}\in\mathbb{N}$ such that

$\|\nabla S_{k+K}\cdots S_{k}f\|_{\infty}\leq c^{\prime}\|\nabla
f\|_{\infty},\qquad\forall k\geq k_{0}.$

∎

##  5 Numerical experiments

In this section, we introduce several numerical experiments. They have been
designed to illustrate the reproduction of parabolas on non-uniform grids and
the generation of 2D and 3D curves using the proposed scheme. We consider
various values of $\rho$ to demonstrate the impact of this parameter on the
final shape. For comparison, we also consider $\rho=0$. In this case, the
proposed scheme coincides with the Lagrange-type linear uniform scheme of
Remark 10. All computations were executed using MATLAB R2022a. The codes
necessary to replicate these experiments are readily available (see section
Reproducibility).

###  5.1 Parabolas reproduction

Figure 2: Left, application of 5 iterations of the scheme with
$\rho\in\\{0,2,6\\}$ (dotted blue, solid green and dashed red, respectively)
to parabolic non-uniform data (black dots). Right, the reproduction error in
log-scale.

We start by confirming the scheme ability to reproduce second degree
polynomials. To show this, we choose a function with a commonly recognized
graph: $F(x)=(x,x^{2})$. It is worth noting that this experiment is applicable
to any $F\in\Pi_{2}^{n}$, where $n\geq 2$, yielding comparable outcomes.

We randomly generated a non-uniform grid $\xi^{0}=(\xi^{0}_{i})_{i=1}^{15}$
within the interval $[-2,2]$,

$\displaystyle\xi^{0}=($
$\displaystyle-1.9381,-1.7893,-1.5751,-1.3313,-1.2075,-0.9235,-0.8321,$
$\displaystyle-0.6420,-0.5104,-0.2734,-0.0412,0.9514,1.6813,1.8065,1.9363),$

for which

$\alpha_{11}^{0}=0.2339\leq\alpha_{i}^{0}=\frac{\nabla\xi_{i-1}}{\nabla\xi_{i}}\leq\alpha_{13}^{0}=5.8299,\quad\beta_{12}^{0}=0.1715\leq\beta_{i}^{0}=\frac{\nabla\xi_{i+1}}{\nabla\xi_{i}}\leq\beta_{10}^{0}=4.2748,\,\forall
i=2,\ldots,13.$

According to Theorem 15, $\rho\geq 4.8299$ must be chosen to achieve an exact
reproduction of the parabola.

In Figure 2-left, we present the results of applying the scheme with
$\rho\in\\{0,2,6\\}$, each for 5 iterations. The boundaries of the interval
$[-2,2]$ are not displayed because the data ‘shrinks’ after each iteration, a
typical behaviour in subdivision. The scheme with $\rho=0$ generates a bump
near the vertex of the parabola. The curve obtained with $\rho=2$ is nearer to
the parabola and no bumps appear. Only for $\rho=6$ is an exact reproduction
achieved up to machine precision.

On the right side of the figure, the errors in the reproduction, computed as
$F|_{\xi^{5}}-S_{5}(F|_{\xi^{0}})$, are displayed on a logarithmic scale. In
some small regions (near $x=1$, for instance), the error is precisely 0, which
cannot be visualized on a log-scale, producing some gaps in the graphic.

###  5.2 2D curve example

|
---|---
|
Figure 3: Example of 2D curves describing a rabbit generated by the presented
scheme with $\rho\in\\{0,2,6\\}$ (dotted blue, solid green and dashed red,
respectively) generated from initial data with varying spacing (black dots).
In the top-left figure, the entire curves can be seen, while the rest of
figures are zooms of interesting regions (ear, tail and paw).

Next, we present the shapes generated by the scheme when a specific drawing is
desired to be created. We have chosen a set of points that outline a rabbit.
We applied the scheme for 5 iterations with $\rho\in\\{0,2,6\\}$. The results
can be examined in Figure 3. Note that the spacing between points varies along
the curve, as there are regions that require more detail, while others do not.
It motivates the use of a non-uniform-based scheme, as ours, which implicitly
examines the distances between the points and suggests the grid spacings for
selecting a Lagrange scheme from Definition 7.

Figure 4: The cumulative arc length versus the curvature of the curve
generated by the scheme with $\rho=2$ in Figure 3 is displayed in the left
side. A zoom over the highest peak is presented in the right side. The
curvature is estimated using the discrete set of points generated by the
scheme after 12 iterations, and the curvature estimation at each generated
point is marked with a black dot. The numerous dots and the continuity of the
curvature along the rabbit curve create the appearance of a continuous
function, as evident in the display.

Note that for $\rho=0$, the generated curve displays noticeable oscillations
in areas such as the ear and paw. Conversely, in other regions like the tail
and back, it appears somewhat ‘stiff’. This highlights the limited
adaptability of the linear scheme. In contrast, the non-linear scheme
($\rho>0$) appears to be more adaptive and flexible. A comparison of the
results for $\rho=2$ and $\rho=6$ clearly demonstrates that this parameter
precisely controls the level of flexibility. It is worth noting the shapes
generated for $\rho=6$ in the ear (which, in our opinion, is also a good
representation of an ear), in the tail (which has two undesirable bumps at the
beginning and end), and in the paw (where the result is quite undesirable due
to an excessive arc length). Good-looking results are obtained for $\rho=2$.

We computed the curvature along the rabbit curve and determined that, for
$\rho\in\\{0,2,6\\}$, it exhibits continuity. Figure 4 specifically depicts
the case for $\rho=2$. It is worth noting that, despite the presence of peaks
where the curvature achieves local maxima, a closer examination reveals that
these peaks are also smooth.

###  5.3 3D curve example

Figure 5: Example of 3D curves describing trefoil knots, taking $v=0$ (top
row) and $v=\pi$ (bottom row) in (31), generated by the presented scheme with
$\rho\in\\{0,1,2\\}$ and with $\rho\in\\{0,2,6\\}$, respectively. The line
aspect for each $\rho$ value is: $\rho=0$, dotted blue; $\rho=1$, dash-dotted
magenta; $\rho=2$, solid green; and $\rho=6$, dashed red. We consider the
front, lateral and top views of these 3D curves (left, central and right
columns, respectively).

Now we proceed with a 3D curve example. Let us consider a trefoil knot surface
described by the parametric equations

$\begin{split}x&=r\frac{\sin(3u)}{2+\cos(v)},\\\
y&=r\frac{\sin(u)+2\sin(2u)}{2+\cos\left(v+\frac{2\pi}{3}\right)},\\\
z&=\frac{r}{8}(\cos(u)-2\cos(2u))(2+\cos(v))\left(2+\cos\left(v+\frac{2\pi}{3}\right)\right),\end{split}$
(31)

with $u,v\in[0,2\pi]$, $r>0$. We look at the two coordinate lines obtained
with $v=0$ and $v=\pi$, $r=1$, and use the points given by $u=2\pi i/9$,
$i\in\mathbb{Z}$, as initial data. We apply the scheme for 5 iterations with
$\rho\in\\{0,1,2\\}$ for $v=0$ (top row) and $\rho\in\\{0,2,6\\}$ for $v=\pi$
(bottom row). We show in Figure 5 the resulting 3D curves from front, lateral,
and top views. We can see that the curves are smooth and the parameter $\rho$
governs the level of flexibility.

##  6 Conclusions

In this study, we introduced a uniform non-linear non-stationary scheme for
generating curves in $\mathbb{R}^{n}$, where $n\geq 2$. It relies on
annihilation operators to automatically estimate grid spacing and select an
appropriate Lagrange-type scheme.

The scheme is dependent on a positive parameter, which allows for control over
both the level of flexibility and the variety of non-uniform grids on which
second-degree polynomials are reproduced. We conducted numerical experiments
to validate these findings.

The scheme has been defined in a non-stationary manner, which facilitated the
proof of its $\mathcal{C}^{1}$ convergence through similarities with a well-
behaved subdivision scheme. We proposed a new class of subdivision schemes, to
which our scheme belongs, and proved their convergence, and even
$\mathcal{C}^{1}$ continuity under a condition that is easily met. We
established convergence in two distinct ways. Firstly, by leveraging the
properties of quasilinear and asymptotically equivalent schemes. Secondly, by
adapting a result from non-linear stationary subdivision theory to our non-
stationary context.

For future research, we propose several avenues. Firstly, the current scheme
exhibited stable behaviour during our numerical testing, meaning that slight
modifications of the data induced small variations on the generated curve.
However, it is worth noting that the piecewise definition of its rules is not
continuous, and thus, stability is not guaranteed. This scheme could be
modified to ensure that its rules are Lipschitz functions, which could result
in a stable scheme. We believe this is achievable, although it might
necessitate sacrificing certain properties, such as limiting its reproduction
capability.

Secondly, the simultaneous reproduction of parabolas on many non-uniform grid
could be extended to higher-degree polynomials or to exponential polynomials.
Another challenging task would be the extension of this work to the
multivariate context.

Thirdly, the regularity was constrained by the regularity of the implicit mesh
refinement, which was based on Chaikin’s scheme. The adoption of a smoother
and strictly monotonicity preserving scheme may be considered.

Fourthly, while the selection of the parameter $\rho$ is straightforward in
practice to achieve good-looking results, finding its optimal value (in a
certain sense) could be interesting. Rather than considering a constant
parameter $\rho$, it could be defined as a function of the data, denoted as
$\rho(f_{i-p},\ldots,f_{i+p+1})$. This would result in a parameter-less scheme
that could potentially offer, for instance, improved adaptability or reduced
distance between the generated curve and the data.

Finally, we demonstrated $\mathcal{C}^{1}$ regularity and also argued that it
is not $\mathcal{C}^{2}$. Nevertheless, the numerical experiments suggest that
it could be curvature continuous. A comprehensive analysis of this topic could
be conducted to reveal the conditions under which the generated curve is
curvature continuous.

## Acknowledgments

This research has been supported by project CIAICO/2021/227 (funded by
Conselleria de Innovación, Universidades, Ciencia y Sociedad digital,
Generalitat Valenciana) and by grant PID2020-117211GB-I00 (funded by
MCIN/AEI/10.13039/501100011033).

## Reproducibility

The MATLAB codes used for the numerical results in this paper, as well as a
Wolfram Mathematica notebook file containing symbolic computations that
complement the proofs presented in the manuscript, are available on Github:
https://github.com/serlou/uniform4nonuniform-subdivision.

## Declaration of Generative AI and AI-assisted technologies in the writing
process

During the preparation of this work the author used ChatGPT in order to
improve readability and language. After using this service, the author
reviewed and edited the content as needed and take full responsibility for the
content of the publication.

## References

  * Amat et al. [2011] Amat, S., Dadourian, K., Liandrat, J., 2011.  Analysis of a class of nonlinear subdivision schemes and associated multiresolution transforms.  Advances in computational Mathematics 34, 253–277.
  * Beccari et al. [2011] Beccari, C.V., Casciola, G., Romani, L., 2011.  Non-uniform interpolatory curve subdivision with edge parameters built upon compactly supported fundamental splines.  BIT Numerical Mathematics 51, 781–808.
  * Cohen et al. [2003] Cohen, A., Dyn, N., Matei, B., 2003.  Quasilinear subdivision schemes with applications to eno interpolation.  Applied and Computational Harmonic Analysis 15, 89–116.
  * Conti et al. [2022] Conti, C., López-Ureña, S., Romani, L., 2022.  Annihilation operators for exponential spaces in subdivision.  Applied Mathematics and Computation 418, 126796.
  * Conti and Romani [2011] Conti, C., Romani, L., 2011\.  Algebraic conditions on non-stationary subdivision symbols for exponential polynomial reproduction.  J. Comput. Appl. Math. 236, 543–556.
  * Daubechies et al. [2004] Daubechies, I., Runborg, O., Sweldens, W., 2004.  Normal multiresolution approximation of curves.  Constructive Approximation 20, 399–463.
  * Donat and López-Ureña [2019] Donat, R., López-Ureña, S., 2019\.  Nonlinear stationary subdivision schemes reproducing hyperbolic and trigonometric functions.  Adv. Comput. Math. 45, 3137–3172.
  * Donat et al. [2017] Donat, R., López-Ureña, S., Santágueda, M., 2017.  A family of non-oscillatory 6-point interpolatory subdivision schemes.  Advances in Computational Mathematics 43, 849–883.
  * Dyn [1992] Dyn, N., 1992.  Subdivision schemes in computer-aided geometric design, in: Advances in numerical analysis, Vol. II (Lancaster, 1990). Oxford Univ. Press, New York. Oxford Sci. Publ., pp. 36–104.
  * Dyn et al. [2004] Dyn, N., Floater, M.S., Hormann, K., 2004.  A $C^{2}$ four-point subdivision scheme with fourth order accuracy and its extensions.  Mathematical Methods for Curves and Surfaces: Tromsø 2005, 145–156.
  * Dyn et al. [2009] Dyn, N., Floater, M.S., Hormann, K., 2009.  Four-point curve subdivision based on iterated chordal and centripetal parameterizations.  Computer Aided Geometric Design 26, 279–286.
  * Dyn et al. [2003] Dyn, N., Levin, D., Luzzatto, A., 2003.  Exponentials reproducing subdivision schemes.  Found. Comput. Math. 3, 187–206.
  * Dyn et al. [2014] Dyn, N., Levin, D., Yoon, J., 2014.  A new method for the analysis of univariate nonuniform subdivision schemes.  Constructive Approximation 40, 173–188.
  * Jeong et al. [2021] Jeong, B., Yang, H., Yoon, J., 2021.  A non-uniform corner-cutting subdivision scheme with an improved accuracy.  Journal of Computational and Applied Mathematics 391, 113446.
  * Kuijt and Van Damme [2002] Kuijt, F., Van Damme, R., 2002\.  Shape preserving interpolatory subdivision schemes for nonuniform data.  Journal of Approximation theory 114, 1–32.
  * Levin [2003] Levin, A., 2003.  Polynomial generation and quasi-interpolation in stationary non-uniform subdivision.  Computer Aided Geometric Design 20, 41–60.
  * López-Ureña and Viscardi [2021] López-Ureña, S., Viscardi, A., 2021\.  An annihilator-based strategy for the automatic detection of exponential polynomial spaces in subdivision.  Computer Aided Geometric Design 85, 101976.
  * Ma [2005] Ma, W., 2005.  Subdivision surfaces for cad—an overview.  Computer-Aided Design 37, 693–709.