2013 January 56 1 1 XX Dehbi L et al.

Corresponding author

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Algebraic Curve Interpolation for Intervals via Symbolic-Numeric Computation

Lydia Dehbi  Zhengfeng Yang  Chao Peng  Yaochen Xu  Zhenbing Zeng  School of
Software Engineering, East China Normal University, Shanghai 200062, China;
Department of Mathematics, Shanghai University, Shanghai 200444, China;

###### Abstract

Algebraic curve interpolation is described by specifying the location of $N$
points in the plane and constructing an algebraic curve of a function $f$ that
should pass through them. In this paper, we propose a novel approach to
construct the algebraic curve that interpolates a set of data (points or
neighborhoods). This approach aims to search the polynomial with the smallest
degree interpolating the given data. Moreover, the paper also presents an
efficient method to reconstruct the algebraic curve of integer coefficients
with the smallest degree and the least monomials that interpolates the
provided data. The problems are converted into optimization problems and are
solved via Lagrange multipliers methods and symbolic computation. Various
examples are presented to illustrate the proposed approaches.

###### keywords:

Algebraic Curve, Mathematics Mechanization, Symbolic-Numeric Computation,
Sparse interpolation.

03C40, 65D05, 65D10

##  1 Introduction

Polynomial interpolation is a rather old mathematical problem [16] and has
significant applications in various practical engineering problems [2, 14, 24,
7, 6, 8, 20, 9]. It has become a fundamental research tool in approximation
theory and numerical analysis [15, 13, 4]. Compared to univariate polynomial
interpolation, the study of multivariate polynomial interpolation and implicit
function interpolation is relatively new. In recent years, many scholars have
focused on these topics and proposed new algorithms [5, 11, 12, 17, 22, 23,
30, 29].

In the study of general interpolation problems, researchers seek to construct
functions of the form $y=f(x_{1},x_{2},\ldots,x_{n})$ that minimize the
overall deviation from the given data (such as [30, 29]). In contrast to
existing research, the problem discussed in this paper is as follows: the
given data consists of small neighborhoods in space (such as square or
spherical neighborhoods), and the goal of interpolation is to construct
polynomial implicit functions that pass through these small neighborhoods. The
measure of interpolation quality is not the deviation between the interpolated
function and the known data, but rather the properties of the degree and
coefficients of the constructed polynomial. For example, the goal is to
minimize the degree of the polynomial and ensure that the coefficients are
integers and have the smallest possible absolute values. Implicit function
interpolation problems, where $f(x_{i},y_{i})=0,1\leq i\leq N$, seek the
function $f(x,y)$ whose curve passes through these points or satisfies
$f(x_{i},y_{i})\approx 0$. Such problems often arise in engineering or
experimental data processing, such as constructing empirical formulas based on
experimental data (as seen in [32, 18]) or problems in computer graphics that
involve reconstructing implicit surfaces from discrete points (see [28, 25,
21, 19, 1, 26, 8, 10]).

In previous research, if there were errors in the data, it was assumed that
the given data represented small neighborhoods, and the goal was for the
interpolated implicit function to pass through the centers of these
neighborhoods as closely as possible. In our problem, the size of the data
neighborhoods is fixed, and the interpolated function only needs to pass
through any point within the neighborhood. The implicit functions obtained
directly using existing methods generally do not satisfy our optimization
objectives regarding the degree and coefficients of the implicit functions. In
practical applications, we often need to fit an empirical formula based on a
limited set of data obtained from experiments or observations. For these
problems, if the degree of the interpolated polynomial implicit function is
too high or the range of coefficient values is too large, the predictive
capability of the formula may be reduced. Therefore, it is desirable for the
interpolated formula to avoid overfitting and be as simple as possible in
terms of its form (degree and coefficients). This increases the likelihood of
capturing the essential connections between the relevant variables.

The rest of this paper is organized as follows. In the next section, we will
first study the computation problem of finding the lowest degree of an
implicit polynomial $f(x,y)$ passing through a given set of accurate points
$(x_{1},y_{1}),(x_{2},y_{2}),\ldots,(x_{N},y_{N})$ on the plane. We construct
a discriminant matrix $A^{[d]}$ of order $(d+1)(d+2)/2$ and prove the
relationship between the lowest degree of the implicit polynomial and the
determinant of the discriminant matrix being zero. Then, we investigate the
computation problem of finding the lowest degree of a polynomial passing
through a given set of rectangular regions $B_{1},B_{2},\ldots,B_{N}$ on the
plane, extending the discriminant method from accurate points to the case of
rectangular neighborhoods. In section 3, we study related sparse interpolation
problems. By analyzing the motion data of 38 celestial bodies in the solar
system, including planets, asteroids, and dwarf planets, we search for the
minimal number of terms in the lowest degree polynomial curve. In section 4,
we investigate the interval interpolation problem of rational-coefficient or
integer-coefficient implicit functions with given monomials. We provide an
algorithm to compute a positive number for the interpolation problem of
integer-coefficient polynomials, such that among the interpolations with
coefficients whose absolute values do not exceed this number, there exists a
unique one passing through all the given neighborhoods. In addition to
presenting rigorous proofs for several theorems, the paper also provides
concrete computational processes for several meaningful theoretical and
practical problems, aiming to assist readers in solving similar problems or
extending the methods presented in this paper. Finally, Section 5 concludes
the paper with a brief summary.

##  2 The lowest degree of approximate implicit interpolation

Let’s consider a set of points $P_{1},P_{2},\cdots,P_{N}$ in the two-
dimensional plane $\mathbb{R}^{2}$. Our objective is to find the simplest
polynomial equation, denoted as $f(x,y)$, that has the lowest possible degree
and passes through all the given points. In the case where the coefficients of
the polynomial are in $\mathbb{Q}$, we also aim to determine the equivalent
polynomial equation with integer coefficients.

First of all, we show the existence of such a polynomial. We define the set
$\mathcal{P}$ as

$\mathcal{P}:=V(P_{1},P_{2},\cdots,P_{N})=\\{f\in\mathbb{R}[x,y]\,|\,f(x_{i},y_{i})=0,~{}~{}i=1,2,\cdots,N\\},$

and we assume a polynomial template of degree $d$

$f(x,y)=c_{0}+c_{1}\cdot x+c_{2}\cdot y+\cdots+c_{{p}-2}\cdot
x^{d-1}y+c_{{p}-1}\cdot y^{d},$ (2.1)

where $f$ has $p={d+2\choose 2}$ coefficients. Let ${\mathbf{c}}_{d}$ and
${\mathbf{v}}_{d}$ be the coefficient vector and the monomial vector of $f$,
i.e., ${\mathbf{c}}_{d}=[c_{0},c_{1},c_{2},\ldots,c_{p-1}]^{T}$ and
${\mathbf{v}}_{d}=[1,x,y,x^{2},xy,y^{2},\ldots,xy^{d-1},y^{d}]^{T}$, $f$ can
be written as $f={\mathbf{c}}_{d}^{T}\,{\mathbf{v}}_{d}$. Let
${\mathbf{v}}_{d}(x_{i},y_{i})$ be the vector by evaluating the monomial
vector ${\mathbf{v}}_{d}$ with the given point $(x_{i},y_{i})$. Then a system
of homogeneous linear equations is constructed as follows:

$\displaystyle{\mathbf{c}}_{d}^{T}\,{\mathbf{v}}_{d}(x_{j},y_{j})=c_{0}+c_{1}\cdot
x_{j}+c_{2}\cdot y_{j}+\cdots+c_{p-2}\cdot x_{j}y_{j}^{d-1}+c_{p-1}\cdot
y_{j}^{d}=0,\quad j=1,2,\ldots,N.$

When $d$ is large enough, e. g., verifying, $d\geq
u(N):=\lfloor\sqrt{2N+1/4}-1/2\rfloor\approx\sqrt{2N}-1$, we have the
following inequality

$2p=(d+1)(d+2)\geq\left(\lfloor\sqrt{2N+\frac{1}{4}}-\frac{1}{2}\rfloor+1\right)\left(\lfloor\sqrt{2N+\frac{1}{4}}-\frac{1}{2}\rfloor+2\right)>2N.$

Therefore, the system of equations (2.1) must have a non-zero solution
${\mathbf{c}}_{d}^{\ast}=[c_{0}^{\ast},c_{1}^{\ast},c_{2}^{\ast},\ldots,c_{p-1}^{\ast}]^{T}$,
and when $(x_{i},y_{i})\,(1\leq i\leq N)$ is a rational number, each component
of the solution vector can also be taken as an integer. This proves the
existence of a polynomial of degree $d$:

$f^{\ast}(x,y):=C\cdot({{\mathbf{c}}_{d}^{\ast}}^{T}\,{\mathbf{v}}_{d})=C\cdot(c_{0}^{\ast}+c_{1}^{\ast}\cdot
x+c_{2}^{\ast}\cdot y+c_{p-2}^{ast}\cdot xy^{d-1}++c_{p-1}^{\ast}\cdot
y^{d})\in\mathcal{P}.$

Obviously, when $N$ points are in some special positions, and $d$ smaller than
$\lfloor\sqrt{2N}\rfloor-1$ the system of homogeneous equations (2.1) may also
have non-zero solutions, so $f(x,y)$ can be constructed according to
(LABEL:poly-k), and the implicit function with smallest degree and passing
through the given points can be obtained.

The number of coefficients of an algebraic curve of degree $d$ grows
dramatically as $d$ increases. Hence, keeping the degree of a curve in a
reasonable range is very important for fast computation and fewer numerical
errors. In what follows, we present a way to determine the smallest degree of
polynomials that pass through a given dataset.

Next, we present a method to determine whether $d$ is the smallest degree
based on the coordinates of $P_{i}$ without solving the system of linear
equations. Given $N$ points and a positive integer $d$, the interpolation of
the implicit function matrix ($d-th$ implicit interpolation aggregate matrix)
is defined as follows.

######  Definition 2.1.

Given $N$ points ($P_{i}=(x_{i},y_{i})$), an integer $d\geq 1$ , and
$p={d+2\choose 2}$, we denote $A^{[d]}$ as a $p\times p$ symmetric matrix
defined as $A^{[d]}=\sum_{i=1}^{N}A_{i}$, where
$A_{i}={\mathbf{v}}_{d}(x_{i},y_{i}){\mathbf{v}}_{d}(x_{i},y_{i})^{T}$, that
is,

$A_{i}=[1,x_{i},y_{i},\cdots,x_{i}^{d},x_{i}^{d-1}y_{i},\cdots,y_{i}^{d}]^{T}\,[1,x_{i},y_{i},\cdots,x_{i}^{d},x_{i}^{d-1}y_{i},\cdots,y_{i}^{d}].$

We have the following theorem.

######  Theorem 2.2.

Given $N$ points and an integer $d$, if the curve $f(x,y)=0$ passes through
the given points, then

$\min_{f\in\mathcal{P}}~{}~{}\mbox{\rm degree}~{}(f)\leq
d~{}~{}\Leftrightarrow~{}~{}\det(A^{[d]})=0.$ (2.2)

######  Proof 2.3.

Let us first prove $\Leftarrow$.

Assuming $(x_{1},y_{1}),\cdots,$ $(x_{N},y_{N})$ satisfy $\det(A^{[d]})=0$, we
have the homogeneous linear equation system

$A^{[d]}\,{\mathbf{c}}_{d}=[0,0,\cdots,0]^{T}.$

This equation has a non-zero solution
${\mathbf{c}}_{d}^{\ast}=[c_{0}^{\ast},c_{1}^{\ast},\cdots,c_{p-1}^{\ast}]^{T}\in\mathbb{R}^{p}\setminus\\{[0,0,\cdots,0]^{T}\\}$,
leading to:

$\sum_{i=1}^{N}A_{i}{\mathbf{c}}_{d}^{\ast}=[0,0,\cdots,0]^{T}.$

Thus, we can derive the following expression:

$\displaystyle\sum_{i=1}^{N}(c_{0}^{\ast}+c_{1}^{\ast}x_{i}+c_{2}^{\ast}y_{i}+\cdots+c_{p-1-d}^{\ast}x_{i}^{d}+\cdots+c_{p-1}^{\ast}y_{i}^{d})^{2}$
(2.3) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}({{\mathbf{c}}_{d}^{\ast}}^{T}\,{\mathbf{v}}_{d}(x_{i},y_{i}))({\mathbf{v}}_{d}(x_{i},y_{i})^{T}\,{\mathbf{c}}_{d}^{\ast})=\sum_{i=1}^{N}{{\mathbf{c}}_{d}^{\ast}}^{T}\,({\mathbf{v}}_{d}(x_{i},y_{i})\,{\mathbf{v}}_{d}(x_{i},y_{i})^{T}\,{\mathbf{c}}_{d}^{\ast})=\sum_{i=1}^{N}{{\mathbf{c}}_{d}^{\ast}}^{T}\,(A_{i}{\mathbf{c}}_{d}^{\ast})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}{{\mathbf{c}}_{d}^{\ast}}^{T}\,[0,0,\cdots,0]^{T}~{}~{}=0,$

which implies that

$f(x_{i},y_{i})=c_{0}^{\ast}+c_{1}^{\ast}x_{i}+c_{2}^{\ast}y_{i}+\cdots+c_{p-1-d}^{\ast}x_{i}^{d}+\cdots+c_{p-1}^{\ast}y_{i}^{d}=0,~{}~{}i=1,2,\cdots,N.$

Therefore, we can conclude that $f(x,y)\in\mathcal{P}$, satisfying the degree
constraint (2.2).

To prove the reverse implication, we Suppose that $f(x,y)$ is a polynomial of
degree $d$ passing through the points $(x_{i},y_{i})$ for $i=1,2,\cdots,N$,
with unknown coefficients $c_{0},\ldots,c_{p-1}$. Consider the quadratic
optimization problem:

$\displaystyle\left\\{\begin{array}[]{ccl}&{\displaystyle\min_{{\mathbf{c}}_{d}\in{\mathbb{R}}^{p}}}&{\displaystyle\sum_{i=1}^{N}}({\mathbf{v}}_{d}(x_{i},y_{i})^{T}\,{\mathbf{c}}_{d})^{2}\\\
&s.t.&\|{\mathbf{c}}_{d}\|_{2}^{2}=1.\end{array}\right.$

we can use the Lagrange multiplier method and the necessary optimality
conditions to derive the following equality:

$\sum_{i=1}^{N}\left(\hat{c}_{0}+\hat{c}_{1}x_{i}+\cdots+\hat{c}_{p-1}y_{i}^{d}\right)^{2}+\lambda(\hat{c}_{0}^{2}+\hat{c}_{1}^{2}+\cdots+\hat{c}_{p-1}^{2})=0,$

From this equation, together with the fact that $(x_{i},y_{i})\in Z(f)$ and
$[\hat{c}_{0},\hat{c}_{1},\cdots,\hat{c}_{p-1}]^{T}\in\mathbb{R}^{p}\setminus\\{[0,0,\cdots,0]^{T}\\},$
we can conclude that $\lambda=0$. Thus,
$[\hat{c}_{0},\hat{c}_{1},\cdots,\hat{c}_{p-1}]^{T}$ is a non-zero solution of
the following homogeneous linear equations

$\displaystyle\sum_{i=1}^{N}(c_{0}+x_{i}c_{1}+y_{i}c_{2}+\cdots+x_{i}^{d}c_{p-1-d}+\cdots+y_{i}^{d}c_{p-1})=0,$
$\displaystyle\sum_{i=1}^{N}x_{i}(c_{0}+x_{i}c_{1}+y_{i}c_{2}+\cdots+x_{i}^{d}c_{p-1-d}+\cdots+y_{i}^{d}c_{p-1})=0,$
$\displaystyle\vdots$
$\displaystyle\sum_{i=1}^{N}y_{i}^{d}(c_{0}+x_{i}c_{1}+y_{i}c_{2}+\cdots+x_{i}^{d}c_{p-1-d}+\cdots+y_{i}^{d}c_{p-1})=0.$

In matrix form, this can be written as:

$\left[\begin{bmatrix}1\\\ x_{1}\\\ \vdots\\\
y_{1}^{d}\end{bmatrix}[1,x_{1},\cdots,y_{1}^{d}]+\cdots+\begin{bmatrix}1\\\
x_{N}\\\ \vdots\\\
y_{N}^{d}\end{bmatrix}[1,x_{N},\cdots,y_{N}^{d}]\right]\;\begin{bmatrix}c_{0}\\\
c_{1}\\\ \vdots\\\
c_{p-1}\end{bmatrix}=(\sum_{i=1}^{N}A_{i})\cdot\begin{bmatrix}c_{0}\\\
c_{1}\\\ \vdots\\\ c_{p-1}\end{bmatrix}=\begin{bmatrix}0\\\ 0\\\ \vdots\\\
0\end{bmatrix}$

which can be further simplified as:

$A^{[d]}\,[c_{0},c_{1},\cdots,c_{p-1}]^{T}=[0,0,\cdots,0]^{T}.$

Thus, we have $\det(A^{[d]})=0.$

######  Corollary 2.4.

For any integer $d$ and any points
$(x_{1},y_{1}),\cdots,(x_{q},y_{q})\in\mathbb{R}^{2}$ with $q<{{d+2}\choose
2}$, we have

$\det\left(\sum_{i=1}^{q}A_{i}\right)=0.$

######  Proof 2.5.

From Definition 2.1, $\mbox{ \rm rank}(A_{i})=1$ for each $i$, which implies

$\mbox{\rm rank}(\sum_{i=1}^{q}A_{i})\leq\sum_{i=1}^{q}\mbox{\rm
rank}(A_{i})=q.$

Since $q<(d+1)(d+2)/2$, it leads to $\det\left(\sum_{i=1}^{q}A_{i}\right)=0$.

The following result provides a lower bound on the degree of polynomials in
$V(P_{1},P_{2},\cdots,P_{N})$, based on the well-known Bézout theorem
(cf.[31]).

######  Corollary 2.6.

For any points
$P_{1}=(x_{1},y_{1}),\cdots,P_{N}=(x_{N},y_{N})\in\mathbb{R}^{2}$. If there
exists a polynomial $f\in\mathbb{R}[x,y]$ of degree $d$ that passes through
$(x_{1},y_{1}),\cdots,(x_{q},y_{q})$ but not through
$(x_{q+1},y_{q+1}),\cdots,$ $(x_{N},y_{N}),$ where ${d+2\choose 2}-1\leq q<N,$
then

$\min_{f\in\mathcal{P}}~{}~{}\mbox{\rm
degree}~{}(f)\geq~{}\max~{}(d+1,~{}\lceil\frac{q}{d}\rceil).$

In other words, if we find a subset of $q$ points in
$P_{1},P_{2},\cdots,P_{N}$ that lie on a straight line, then it implies that
either all the points $P_{1},P_{2},\cdots,P_{N}$ lie on that same line, or any
polynomial passing through $P_{1},P_{2},\cdots,P_{N}$ must have a degree of at
least $q$.

The following example explain the usage of Theorem 2.2 and Corollary 2.4, 2.6

######  Example 2.7.

Consider the following example to illustrate our approach. We have 22 points
in the plane, as shown in Figure 1. The coordinates of these points are given
in Table 1.

Table 1: The coordinates of the $22$ points.
$\begin{array}[]{|c|c|c||c|c|c||c|c|c|}\hline\cr\phantom{PPP}&\phantom{xxxx}x\phantom{xxxx}&\phantom{yy}y\phantom{yy}&\phantom{PPP}&\phantom{xxxx}x\phantom{xxxx}&\phantom{yy}y\phantom{yy}&\phantom{PPP}&\phantom{xxxx}x\phantom{xxxx}&\phantom{yyyi}y\phantom{yyyi}\\\
\hline\cr
P_{1}&0&1&P_{9}&1&2&P_{17}&\frac{\sqrt{18+2\,\sqrt{73}}}{6}&\frac{2}{3}\\\
\hline\cr
P_{2}&-\frac{2}{3}&\frac{1}{3}&P_{10}&\frac{2\,\sqrt{5}}{3}&\frac{5}{3}&P_{18}&-\frac{\sqrt{9-\sqrt{73}}}{3}&\frac{4}{3}\\\
\hline\cr
P_{3}&\frac{2}{3}&\frac{1}{3}&P_{11}&\frac{\sqrt{10}}{6}&\frac{5}{3}&P_{19}&\frac{\sqrt{9+\sqrt{73}}}{3}&\frac{4}{3}\\\
\hline\cr
P_{4}&-\frac{\sqrt{2}}{6}&\frac{1}{3}&P_{12}&\frac{\sqrt{12-6\,\sqrt{2}}}{8}&\frac{1}{4}&P_{20}&-\frac{\sqrt{9+\sqrt{73}}}{3}&\frac{4}{3}\\\
\hline\cr
P_{5}&-\frac{2\,\sqrt{5}}{3}&\frac{5}{3}&P_{13}&\frac{\sqrt{12+6\,\sqrt{2}}}{8}&\frac{1}{4}&P_{21}&-\frac{\sqrt{12+9\,\sqrt{2}}}{4}&1+\frac{3\sqrt{2}}{4}\\\
\hline\cr
P_{6}&-\sqrt{2}&2&P_{14}&-\frac{\sqrt{18-6\,\sqrt{7}}}{4}&\frac{3}{2}&P_{22}&\frac{1}{4}&1+\frac{3\sqrt{2}}{4}\\\
\hline\cr
P_{7}&-\frac{\sqrt{6}}{2}&1&P_{15}&\frac{\sqrt{18-6\,\sqrt{7}}}{4}&\frac{3}{2}&&&\\\
\hline\cr
P_{8}&\frac{\sqrt{6}}{2}&1&P_{16}&-\frac{\sqrt{18-2\,\sqrt{73}}}{6}&\frac{2}{3}&&&\\\
\hline\cr\end{array}$

We aim to compute the minimal polynomial ${f}(x,y)=0$ which defines the
implicit equation of the curve $f(x,y)=0$ passing through these 22 points. We
have $N=22$ and

${{5+2}\choose 2}=21<22,\quad{{6+2}\choose 2}=28>22.$

Thus, the smallest degree of a polynomial that passes through the 22 points at
most $6$,

$\mbox{\rm degree}~{}{(f)}\leq 6.$

Therefore, it is evident that the line $y=1$ passes through the points
$P_{1},~{}P_{7}$ and $P_{8}$ but the other $19$ points do not lie on this
line, so according to Corollary 2.6, this implies that:

$\mbox{\rm degree}~{}({f})\geq\max\\{1+1,\lceil\frac{3}{1}\rceil\\}=3.$

So, we first search for a polynomial of degree $d=3$ in the following form

$f(x,y)=c_{0}+c_{1}x+c_{2}y+c_{3}x^{2}+c_{4}xy+c_{5}y^{2}+c_{6}x^{3}+c_{7}x^{2}y+c_{8}xy^{2}+c_{9}y^{3}.$

Figure 1: The 22 points (small black circles) on the plane selected from the
quartic curve ${y}^{4}-2\,{y}^{3}+{y}^{2}-3\,{x}^{2}y+2\,{x}^{4}=0$ for
recovering the polynomial.

Since $(d+1)(d+2)/2=10$, we first choose 10 points from the 22 points, for
instance, $P_{1},\cdots,P_{10}$. . We then check whether there are cubic
curves passing through these 10 points. For this purpose, we substitute the
coordinates $x=x_{i},y=y_{i}$ of each point into a $10\times 10$ symmetric
matrix generated by

$[1,x,y,x^{2},xy,y^{2},x^{3},x^{2}y,xy^{2},y^{3}]^{T}[1,x,y,x^{2},xy,y^{2},x^{3},x^{2}y,xy^{2},y^{3}],$

as follows:

$A_{i}=\left(\begin{array}[]{cccccccccc}1&x_{i}&y_{i}&{x_{i}^{2}}&x_{i}y_{i}&{y_{i}^{2}}&{x_{i}^{3}}&{x_{i}^{2}}y_{i}&x_{i}{y_{i}^{2}}&{y_{i}^{3}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
x_{i}&{x_{i}^{2}}&x_{i}y_{i}&{x_{i}^{3}}&{x_{i}^{2}}y_{i}&x_{i}{y_{i}^{2}}&{x_{i}^{4}}&{x_{i}^{3}}y_{i}&{x_{i}^{2}}{y_{i}^{2}}&{y_{i}^{3}}x_{i}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
y_{i}&x_{i}y_{i}&{y_{i}^{2}}&{x_{i}^{2}}y_{i}&x_{i}{y_{i}^{2}}&{y_{i}^{3}}&{x_{i}^{3}}y_{i}&{x_{i}^{2}}{y_{i}^{2}}&{y_{i}^{3}}x_{i}&{y_{i}^{4}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{x_{i}^{2}}&{x_{i}^{3}}&{x_{i}}^{2}y_{i}&{x^{4}}&{x_{i}^{3}}y_{i}&{x_{i}^{2}}{y_{i}^{2}}&{x_{i}^{5}}&{x_{i}^{4}}y_{i}&{x_{i}^{3}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{2}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
x_{i}y_{i}&{x_{i}^{2}}y_{i}&x_{i}{y_{i}^{2}}&{x_{i}^{3}}y_{i}&{x_{i}^{2}}{y_{i}^{2}}&{y_{i}^{3}}x_{i}&{x_{i}^{4}}y_{i}&{x_{i}^{3}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{2}}&{y_{i}^{4}}x_{i}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{y_{i}^{2}}&x_{i}{y_{i}^{2}}&{y_{i}^{3}}&{x_{i}^{2}}{y_{i}^{2}}&{y_{i}^{3}}x_{i}&{y_{i}^{4}}&{x_{i}^{3}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{2}}&{y_{i}^{4}}x_{i}&{y_{i}^{5}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{x_{i}^{3}}&{x_{i}^{4}}&{x_{i}^{3}}y_{i}&{x_{i}^{5}}&{x_{i}^{4}}y_{i}&{x_{i}^{3}}{y_{i}^{2}}&{x_{i}^{6}}&{x_{i}^{5}}y_{i}&{x_{i}^{4}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{3}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{x_{i}^{2}}y_{i}&{x_{i}^{3}}y_{i}&{x_{i}^{2}}{y_{i}^{2}}&{x_{i}^{4}}y_{i}&{x_{i}^{3}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{2}}&{x_{i}^{5}}y_{i}&{x_{i}^{4}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{3}}&{y_{i}^{4}}{x_{i}^{2}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
x_{i}{y_{i}^{2}}&{x_{i}^{2}}{y_{i}^{2}}&{y_{i}^{3}}x_{i}&{x_{i}^{3}}{y_{i}}^{2}&{y_{i}^{3}}{x_{i}^{2}}&{y_{i}^{4}}x_{i}&{x_{i}^{4}}{y_{i}^{2}}&{y_{i}^{3}}{x_{i}^{3}}&{y_{i}^{4}}{x_{i}}^{2}&{y_{i}^{5}}x_{i}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{y_{i}^{3}}&{y_{i}^{3}}x_{i}&{y_{i}^{4}}&{y_{i}^{3}}{x_{i}^{2}}&{y_{i}^{4}}x_{i}&{y_{i}^{5}}&{y_{i}^{3}}{x_{i}^{3}}&{y_{i}^{4}}{x_{i}^{2}}&{y_{i}^{5}}x_{i}&{y_{i}^{6}}\end{array}\right).$

By summing up all the matrices $A_{i}$, we obtain the matrix
$A^{[3]}=\sum_{i=1}^{10}A_{i}$, which is a $10\times 10$ symmetric matrix
given by:

$A^{[3]}=\left(\begin{array}[]{cc}U&V\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr V^{T}&W\end{array}\right),$

where $U,V,W$ are $5\times 5$ matrices and $U,W$ are symmetric, as follows:

$\begin{array}[]{l}U=\left(\begin{array}[]{ccccc}\scalebox{0.7071}{$10$}&\scalebox{0.7071}{$-$}\frac{7}{6}\scalebox{0.7071}{$\sqrt{2}+1$}&{\frac{34}{3}}&{\frac{205}{18}}&\scalebox{0.7071}{$-$}{\frac{37}{18}}\scalebox{0.7071}{$\sqrt{2}+2$}\\\\[7.0pt]
\scalebox{0.7071}{$-$}\frac{7}{6}\scalebox{0.7071}{$\sqrt{2}+1$}&{\frac{205}{18}}&\scalebox{0.7071}{$-$}{\frac{37}{18}}\scalebox{0.7071}{$\sqrt{2}+2$}&\scalebox{0.7071}{$-$}{\frac{217}{108}}\scalebox{0.7071}{$\sqrt{2}+1$}&{\frac{301}{18}}\\\\[7.0pt]
{\frac{34}{3}}&\scalebox{0.7071}{$-$}{\frac{37}{18}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{152}{9}}&{\frac{301}{18}}&\scalebox{0.7071}{$-$}{\frac{217}{54}}\scalebox{0.7071}{$\sqrt{2}+4$}\\\\[7.0pt]
{\frac{205}{18}}&\scalebox{0.7071}{$-$}{\frac{217}{108}}\scalebox{0.7071}{$\sqrt{2}+1$}&{\frac{301}{18}}&{\frac{6407}{324}}&\scalebox{0.7071}{$-$}{\frac{1297}{324}}\scalebox{0.7071}{$\sqrt{2}+2$}\\\\[7.0pt]
\scalebox{0.7071}{$-$}{\frac{37}{18}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{301}{18}}&\scalebox{0.7071}{$-$}{\frac{217}{54}}\scalebox{0.7071}{$\sqrt{2}+4$}&\scalebox{0.7071}{$-$}{\frac{1297}{324}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{4447}{162}}\end{array}\right),\end{array}$
$\begin{array}[]{l}V=\left(\begin{array}[]{ccccc}{\frac{152}{9}}&\scalebox{0.7071}{$-$}{\frac{217}{108}}\scalebox{0.7071}{$\sqrt{2}+1$}&{\frac{301}{18}}&\scalebox{0.7071}{$-$}{\frac{217}{54}}\scalebox{0.7071}{$\sqrt{2}+4$}&{\frac{766}{27}}\\\\[7.0pt]
\scalebox{0.7071}{$-$}{\frac{217}{54}}\scalebox{0.7071}{$\sqrt{2}+4$}&{\frac{6407}{324}}&\scalebox{0.7071}{$-$}{\frac{1297}{324}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{4447}{162}}&\scalebox{0.7071}{$-$}{\frac{1297}{162}}\scalebox{0.7071}{$\sqrt{2}+8$}\\\\[7.0pt]
{\frac{766}{27}}&\scalebox{0.7071}{$-$}{\frac{1297}{324}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{4447}{162}}&\scalebox{0.7071}{$-$}{\frac{1297}{162}}\scalebox{0.7071}{$\sqrt{2}+8$}&{\frac{4088}{81}}\\\\[7.0pt]
{\frac{4447}{162}}&\scalebox{0.7071}{$-$}{\frac{7777}{1944}}\scalebox{0.7071}{$\sqrt{2}+1$}&{\frac{30223}{972}}&\scalebox{0.7071}{$-$}{\frac{7777}{972}}\scalebox{0.7071}{$\sqrt{2}+4$}&{\frac{857}{18}}\\\\[7.0pt]
\scalebox{0.7071}{$-$}{\frac{1297}{162}}\scalebox{0.7071}{$\sqrt{2}+8$}&{\frac{30223}{972}}&\scalebox{0.7071}{$-$}{\frac{7777}{972}}\scalebox{0.7071}{$\sqrt{2}+4$}&{\frac{857}{18}}&\scalebox{0.7071}{$-$}{\frac{7777}{486}}\scalebox{0.7071}{$\sqrt{2}+16$}\end{array}\right),\end{array}$
$\begin{array}[]{l}W=\left(\begin{array}[]{ccccc}{\frac{4088}{81}}&\scalebox{0.7071}{$-$}{\frac{7777}{972}}\scalebox{0.7071}{$\sqrt{2}+4$}&{\frac{857}{18}}&\scalebox{0.7071}{$-$}{\frac{7777}{486}}\scalebox{0.7071}{$\sqrt{2}+16$}&{\frac{22534}{243}}\\\\[7.0pt]
\scalebox{0.7071}{$-$}{\frac{7777}{972}}\scalebox{0.7071}{$\sqrt{2}+4$}&{\frac{220879}{5832}}&\scalebox{0.7071}{$-$}{\frac{46657}{5832}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{151571}{2916}}&\scalebox{0.7071}{$-$}{\frac{46657}{2916}}\scalebox{0.7071}{$\sqrt{2}+8$}\\\\[7.0pt]
{\frac{857}{18}}&\scalebox{0.7071}{$-$}{\frac{46657}{5832}}\scalebox{0.7071}{$\sqrt{2}+2$}&{\frac{151571}{2916}}&\scalebox{0.7071}{$-$}{\frac{46657}{2916}}\scalebox{0.7071}{$\sqrt{2}+8$}&{\frac{124375}{1458}}\\\\[7.0pt]
\scalebox{0.7071}{$-$}{\frac{7777}{486}}\scalebox{0.7071}{$\sqrt{2}+16$}&{\frac{151571}{2916}}&\scalebox{0.7071}{$-$}{\frac{46657}{2916}}\scalebox{0.7071}{$\sqrt{2}+8$}&{\frac{124375}{1458}}&\scalebox{0.7071}{$-$}{\frac{46657}{1458}}\scalebox{0.7071}{$\sqrt{2}+32$}\\\\[7.0pt]
{\frac{22534}{243}}&\scalebox{0.7071}{$-$}{\frac{46657}{2916}}\scalebox{0.7071}{$\sqrt{2}+8$}&{\frac{124375}{1458}}&\scalebox{0.7071}{$-$}{\frac{46657}{1458}}\scalebox{0.7071}{$\sqrt{2}+32$}&{\frac{126752}{729}}\end{array}\right).\end{array}$

By calling the computer algebra software Maple, we can verify that:

$\det(A^{[3]})={\frac{341196800000}{282429536481}}+{\frac{280985600000}{847288609443}}\,\sqrt{2}\approx
1.67707205.$

According to Theorem 2.2, we can conclude that no cubic curve $f(x,y)$ passes
through the selected 10 points. and thus, it cannot pass through the given 22
points.

Next, we consider the case $d=4$. Since ${4+2\choose 2}=15$ we select 5
additional interpolating points from the remaining points, for example,
$P_{11},\cdots,P_{15}$. Similarly to before, we calculate the matrix
$A_{i}\,(1\leq i\leq 15)$ for each of the 15 points and then sum them to
obtain the $15\times 15$ matrix $A^{[4]}:=\sum_{i=1}^{15}A_{i}$. Then, we
calculate the determinant of $A^{[4]}$ using Maple, and the result is
$\det(A^{[4]})=0$. According to Theorem 2.2, this shows that there is a
quartic curve passing through $P_{1},P_{2},\cdots,P_{15}$

To determine the coefficients of the quartic curve, we solve the system of
equations:

$A^{[4]}\,[c_{0},c_{1},\cdots,c_{14}]^{T}=[0,0,\cdots,0]^{T}$

where $[c_{0},c_{1},\ldots,c_{14}]^{T}$ represents the coefficient vector. The
non-zero coefficients are found to be:

$c_{5}=1,\quad c_{7}=-3,\quad c_{9}=-2,\quad c_{10}=2,\quad c_{14}=1.$

Hence, $P_{1},P_{2},\cdots,P_{15}$ lie on the following quartic curve:

$f(x,y)=y^{2}-3x^{2}y-2y^{3}+2x^{4}+y^{4}.$

By substituting the coordinates of the other 7 points into this polynomial, we
observe that all the 22 points lie on the quadratic curve. Finally, it can be
asserted that the lowest degree of polynomial curves passing through 22 points
simultaneously is 4.

If the data contains errors, we can consider that the position of each point
$(x_{i},y_{i})$ lies within a hemisphere defined by
$(x-x_{i})^{2}+(y-y_{i})^{2}\leq\delta^{2}$, or within a rectangle defined by
$[x_{i}-\delta,x_{i}+\delta]\times[y_{i}-\delta,y_{i}+\delta]$. In this case,
implicit function interpolation can be formulated as the following
optimization problem:

$\displaystyle\left\\{\begin{array}[]{ccl}&{\displaystyle\min_{{\mathbf{c}}_{d},\,x_{i},\,y_{i}}}&{\displaystyle\sum_{i=1}^{N}}(f(x_{i},y_{i},{\mathbf{c}}_{d}))^{2}\\\
&s.t.&(x_{i},y_{i})\in B_{i},\quad i=1,2,\cdots,N,\\\
&&\|{\mathbf{c}}_{d}\|_{2}^{2}=1.\end{array}\right.$ (2.8)

In this case, if we find a real-coefficient polynomial of the desired degree,
we can perturb it to obtain a rational-coefficient polynomial that passes
through the given neighborhood $B_{i}$. In some practical problems, the given
data $(x_{i},y_{i})$ may have random errors that follow a certain distribution
(such as a normal distribution). The goal is to construct a polynomial curve
$f(x,y)=0$ that approximately satisfies the equation.

$|f(x_{i},y_{i})|<\epsilon,\quad i=1,2,\cdots,N,$

where $\epsilon$ is a small positive number. The optimization problem
corresponding to this type of problem is equivalent to adding $N$ linear
inequalities of the following form to the constraints of problem (2.8):

$-\epsilon<{{\mathbf{c}}_{d}}^{T}\,{\mathbf{v}}_{d}=c_{0}+c_{1}\cdot
x+c_{2}\cdot y+c_{p-2}\cdot xy^{d-1}++c_{p-1}\cdot y^{d}<\epsilon.$

Now consider the case in problem (2.8) where $B_{i}$ represents a rectangular
neighborhood. By applying the Lagrange multiplier method, the condition for
the optimum can be transformed into the following semi-algebraic decision
problem:

$\displaystyle\left(\exists(x_{1},y_{1})\in B_{1},\cdots,(x_{N},y_{N})\in
B_{N}\right)\wedge\,\exists{\mathbf{c}}_{d}\in{\mathbb{R}}^{p}\,\left(\|{\mathbf{c}}_{d}\|_{2}=1\right),$
$\displaystyle\left(A(x_{1},y_{1},\cdots,x_{N},y_{N})\right)_{p\times
p}{\mathbf{c}}_{d}=[0,0,\cdots,0]^{T}.$ (2.9)

For this problem, applying Theorem 2.2 yields the result of the lowest degree
of polynomials that pass through all neighborhoods.

######  Corollary 2.8.

For $N$ rectangles on the plane
$B_{i}=(\underline{a_{i}},\overline{a_{i}})\times(\underline{b_{i}},\overline{b_{i}})\,(i=1,2,\cdots,N)$
and a positive integer $d$, let

$V(B_{1},B_{2},\cdots,B_{N})=\\{f\in\mathbb{Z}[x,y]\,|\,\mbox{Zero}(f)\cap
B_{i}\not=\emptyset\\},$

where

$\mbox{Zero}(f)=\\{(x,y)\in\mathbb{R}^{2}|f(x,y)=0\\},$

and $D_{d}:\mathbb{R}^{2N}\rightarrow R$ a polynomial function defined by

$D_{d}(x_{1},y_{1},\cdots,x_{N},y_{N})=\det(A^{[d]}(x_{1},y_{1},\cdots,x_{N},y_{N})).$

If $D_{d}$ is positive definite on the box $B_{1}\times
B_{2}\times\cdots\times B_{N}\subset\mathbb{R}^{2N}$ then

$\min~{}\\{\mbox{\rm degree}~{}(f)\,|\,~{}f\in
V(B_{1},B_{2},\cdots,B_{N})\\}~{}\geq~{}d+1.$ (2.10)

######  Proof 2.9.

Since for any ${\mathbf{c}}_{d}\in{\mathbb{R}}^{p}$, we have

$\displaystyle{\mathbf{c}}_{d}\cdot
A_{d}(x_{1},y_{1},\cdots,x_{N},y_{N})\cdot{{\mathbf{c}}_{d}}^{T}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left({\mathbf{c}}_{d}\cdot[1,x_{i},y_{i},\cdots,x_{i}^{d},\cdots,y_{i}^{d}]^{T}(1,x_{i},y_{i},\cdots,x_{i}^{d},\cdots,y_{i}^{d})\cdot{{\mathbf{c}}_{d}}^{T}\right),$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}([1,x_{i},y_{i},\cdots,x_{i}^{d},\cdots,y_{i}^{d}]\cdot{\mathbf{c}}_{d})^{2}~{}~{}\geq~{}~{}0,$

the matrix $A^{[d]}$ is a $p\times p$ semi-definite matrix and $D_{d}\geq 0$
for all $(x_{1},y_{1},\cdots,x_{N},y_{N})\in\mathbb{R}$. Thus, if (2.10) is
not valid, then there exists
$[c_{0},c_{1},\cdots,c_{p-1}]\in\mathbb{R}^{p}\setminus\\{[0,0,\cdots,0]^{T}\\}$
such that

$f=c_{0}+c_{1}x+c_{2}y+\cdots+c_{p-1-d}x^{d}+\cdots+c_{p-1}y^{d}\in
V(B_{1},B_{2},\cdots,B_{N}),$

that is,

$f(x_{i},y_{i})=c_{0}+c_{1}x_{i}+c_{2}y_{i}+\cdots+c_{p-1-d}x_{i}^{d}+\cdots+c_{p-1}y_{i}^{d}=0,\quad
i=1,2,\cdots,N,$

and therefore, for some $(x_{i},y_{i})\in B_{i}\ (i=1,2,\cdots,N)$, we have

$D_{d}(x_{1},x_{2},\cdots,x_{N},y_{N})=0.$

According to Corollary 2.4, when $N<{{d+2}\choose 2}$, for any
$(x_{1},y_{1},\cdots,x_{N},y_{N})\in B_{1}\times B_{2}\times\cdots\times
B_{N}$, we always have:

$D_{d}(x_{1},y_{1},\cdots,x_{N},y_{N})=\det(A^{[d]}(x_{1},y_{1},\cdots,x_{N},y_{N}))=0,$

which means that $D_{d}$ cannot be positive definite over $B_{1}\times
B_{2}\times\cdots\times B_{N}$. When $N\geq{{d+2}\choose 2}$, the
determination of whether $D_{d}$ is positive definite over $B_{1}\times
B_{2}\times\cdots\times B_{N}$, as stated in Corollary 2.10, can be
transformed into the following optimization problem:

$\min\\{D_{d}(x_{1},y_{1},\cdots,x_{N},y_{N})\,|\,\underline{a}_{i}\leq
x_{i}\leq\overline{a}_{i},\underline{b}_{i}\leq
y_{i}\leq\overline{b}_{i},i=1,2,\cdots,N\\}.$ (2.11)

Alternatively, it can be expressed as the following quantifier elimination
problem:

$\begin{array}[]{ll}\exists(x_{1},y_{1},\cdots,x_{N},y_{N})\\\
\phantom{\exists}(\underline{a}_{i}\leq
x_{i}\leq\overline{a}_{i},\underline{b}_{i}\leq
y_{i}\leq\overline{b}_{i},i=1,\cdots,N)\wedge
D_{d}(x_{1},y_{1},\cdots,x_{N},y_{N})=0.\end{array}$

Note that the number of terms in $D_{2}$ increases from $125865$ for $N=6$ to
$782775$ for $N=7$. To simplify $D_{d}$, we introduce the matrix $M_{d}$ of
order $p\times N$ defined as follows.

######  Definition 2.10.

Let $d$ be a positive integer, $p={d+2\choose 2}$, and
$(x_{1},y_{1}),(x_{2},y_{2}),\cdots,(x_{N},y_{N})$ be $N$ points in the plane.
The $p\times N$ matrix $M_{d}$ defined below is referred to as the $d$th order
implicit function interpolation basis matrix for these $N$ points:

$M_{d}=M(x_{1},y_{1},x_{2},y_{2},\cdots,x_{N},y_{N}):=\left(\begin{array}[]{llcr}1&1&\cdots&1\\\
x_{1}&x_{2}&\cdots&x_{N}\\\ y_{1}&y_{2}&\cdots&y_{N}\\\
\vdots&\vdots&\ddots&\vdots\\\ x_{1}^{d}&x_{2}^{d}&\cdots&x_{N}^{d}\\\
x_{1}^{d-1}y_{1}&x_{2}^{d-1}y_{2}&\cdots&x_{N}^{d-1}y_{N}\\\
\vdots&\vdots&\ddots&\vdots\\\
y_{1}^{d}&y_{2}^{d}&\cdots&y_{N}^{d}\end{array}\right).$

Given $N$ points, there is a connection between the $d$-th order implicit
function interpolation basis matrix $M_{d}$ and the matrix $A^{[d]}$.

######  Lemma 2.11.

For any $N$ points $\\{(x_{i},y_{i}),i=1\cdots N\\}\in\mathbb{R}^{2N}$, we
have:

$A^{[d]}=M_{d}\cdot M_{d}^{T}.$

n particular, when $N=p$, we also have:

$\det(A^{[d]})=\det(M_{d}^{2}).$

######  Proof 2.12.

$M_{d}\cdot M_{d}^{T}=\left(\begin{array}[]{llcr}1&1&\cdots&1\\\
x_{1}&x_{2}&\cdots&x_{N}\\\ y_{1}&y_{2}&\cdots&y_{N}\\\
\vdots&\vdots&\ddots&\vdots\\\ x_{1}^{d}&x_{2}^{d}&\cdots&x_{N}^{d}\\\
x_{1}^{d-1}y_{1}&x_{2}^{d-1}y_{2}&\cdots&x_{N}^{d-1}y_{N}\\\
\vdots&\vdots&\ddots&\vdots\\\
y_{1}^{d}&y_{2}^{d}&\cdots&y_{N}^{d}\end{array}\right).\left(\begin{array}[]{cccccccc}1&x_{1}&y_{1}&\cdots&x_{1}^{d}&x_{1}^{d-1}y_{1}&\cdots&y_{1}^{d}\\\
1&x_{2}&y_{2}&\cdots&x_{2}^{d}&x_{2}^{d-1}y_{2}&\cdots&y_{2}^{d}\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
1&x_{N}&y_{N}&\cdots&x_{N}^{d}&x_{N}^{d-1}y_{N}&\cdots&y_{N}^{d}\end{array}\right),$

$=\scalebox{0.9}{$\left(\begin{array}[]{cccccccc}N&\left[x_{i}\right]&\left[y_{i}\right]&\cdots&\left[x_{i}^{d}\right]&\left[x_{i}^{d-1}y_{i}\right]&\cdots&\left[y_{i}^{d}\right]\\\\[5.0pt]
\left[x_{i}\right]&\left[x_{i}\,x_{i}\right]&\left[x_{i}\,y_{i}\right]&\cdots&\left[x_{i}\,x_{i}^{d}\right]&\left[x_{i}\,x_{i}^{d-1}y_{i}\right]&\cdots&\left[x_{i}\,y_{i}^{d}\right]\\\\[5.0pt]
\left[y_{i}\right]&\left[y_{i}\,x_{i}\right]&\left[y_{i}\,y_{i}\right]&\cdots&\left[y_{i}\,x_{i}^{d}\right]&\left[y_{i}\,x_{i}^{d-1}y_{i}\right]&\cdots&\left[y_{i}\,y_{i}^{d}\right]\\\\[5.0pt]
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\\[5.0pt]
\left[x_{i}^{d}\right]&\left[x_{i}^{d}\,x_{i}\right]&\left[x_{i}^{d}\,y_{i}\right]&\cdots&\left[x_{i}^{d}\,x_{i}^{d}\right]&\left[x_{i}^{d}\,x_{i}^{d-1}y_{i}\right]&\cdots&\left[x_{i}^{d}\,y_{i}^{d}\right]\\\\[5.0pt]
\left[x_{i}^{d-1}y_{i}\right]&\left[x_{i}^{d-1}y_{i}\,x_{i}\right]&\left[x_{i}^{d-1}y_{i}\,y_{i}\right]&\cdots&\left[x_{i}^{d-1}y_{i}\,x_{i}^{d}\right]&\left[x_{i}^{d-1}y_{i}\,x_{i}^{d-1}y_{i}\right]&\cdots&\left[x_{i}^{d-1}y_{i}\,y_{i}^{d}\right]\\\\[5.0pt]
\vdots&\vdots&\vdots&\ddots&\cdots&\vdots&\ddots&\vdots\\\\[5.0pt]
\left[y_{i}^{d}\right]&\left[y_{i}^{d}\,x_{i}\right]&\left[y_{i}^{d}\,y_{i}\right]&\cdots&\left[y_{i}^{d}\,x_{i}^{d}\right]&\left[y_{i}^{d}\,x_{1}^{d-1}y_{i}\right]&\cdots&\left[y_{1}^{d}\,y_{i}^{d}\right]\end{array}\right).$}\\\
$ (2.12)

Here, we use the notation
$[x_{i}],[y_{i}],[x_{i}y_{i}],\ldots,[x_{i}^{d}],[x_{i}^{d-1}y_{i}],\ldots,[y_{i}^{d}]$
to record

$\displaystyle[x_{i}]$ $\displaystyle=x_{1}+x_{2}+\cdots+x_{N},$
$\displaystyle[y_{i}]$ $\displaystyle=y_{1}+y_{2}+\cdots+y_{N},$
$\displaystyle[x_{i}y_{i}]$
$\displaystyle=x_{1}y_{1}+x_{2}y_{2}+\cdots+x_{N}y_{N},$ $\displaystyle\cdots$
$\displaystyle[x_{i}^{d}]$
$\displaystyle=x_{1}^{d}+x_{2}^{d}+\cdots+x_{N}^{d},$
$\displaystyle[x_{i}^{d-1}y_{i}]$
$\displaystyle=x_{1}^{d-1}y_{1}+x_{2}^{d-1}y_{2}+\cdots+x_{N}^{d-1}y_{N},$
$\displaystyle\cdots$ $\displaystyle[y_{i}^{d}]$
$\displaystyle=y_{1}^{d}+y_{2}^{d}+\cdots+y_{N}^{d}.$

Comparing these definitions with the definition of the matrix $A^{[d]}$ in
Definition 2.1, we can observe that the matrix (2.12) is exactly equal to the
expanded form of $A^{[d]}$.

The corollary below follows directly from Corollary 2.10.

######  Corollary 2.13.

Consider rectangles
$B_{i}=(\underline{a_{i}},\overline{a_{i}})\times(\underline{b_{i}},\overline{b_{i}})$
for $i=1,2,\cdots,N)$, and $\Delta_{d}$ a polynomial function defined on
$\mathbb{R}^{2p}$ by

$\Delta_{d}=\det(M_{d}(x_{1},y_{1},\cdots,x_{p},y_{p})).$

If $\Delta_{d}$ is positive definite or negative definite on the set
$B_{j_{1}}\times B_{j_{2}}\times\cdots\times B_{j_{p}}\subset\mathbb{R}^{2p}$
for at least one $(j_{1},j_{2},\cdots,j_{p})$, with $1\leq
j_{1}<j_{2}<\cdots<j_{p}\leq 2N$, then

$\min~{}\\{\mbox{\rm degree}(f)\,|\,f\in V(B_{1},B_{2},\cdots,B_{N})\\}\geq
d+1.$

The degree of the polynomial $\Delta_{d}$ is given by

$\mbox{\rm
degree}(\Delta_{d})=\sum_{j=1}^{d}j(j+1)=\frac{1}{3}\,d\left(d+1\right)\left(d+2\right).$

The degrees with respect to $x_{i}$ and $y_{i}$ are

$\quad\mbox{\rm degree}(\Delta_{d},x_{i})=d,\quad\mbox{\rm
degree}(\Delta_{d},y_{i})=d,\text{for all}i=1,2,\cdots,p.$

The expanded form of $\Delta_{d}$ contains $p!$ monomials. For example, when
$d=2$, the minimum value of $p$ is ${4+2\choose 2}=6$. Therefore, the total
degree of the polynomial $\Delta_{2}$ is 8, and the number of monomials it
contains is at least $p!$, which is 720.

######  Example 2.14.

The following data represents the approximate values of the average distance
from the Sun (semi-major axis, which is the mean of perihelion and aphelion
distances) and the orbital period of the six closest planets to the Sun in the
solar system: Mercury, Venus, Earth, Mars, Jupiter, and Saturn. Each data
point $(r,p)$ is represented by a rectangle on a plane, with the side length
of the square being $0.0006$. Here is an example for the data of Mercury:
$\mbox{Mer}:=[0.3868,0.3874]\times[0.2405,0,2411]$. The unit of $r$ is in
astronomical units ($A.U.$), where $1A.U.=149,597,870$ kilometers. The unit of
$p$ is the Earth’s orbital period, which is used as a reference for time.

$\left.\begin{array}[]{cll}&\mbox{Mercury):}&\mbox{\phantom{xx}
Mer}=[0.3868,0.3874]\times[0.2405,0,2411],\\\
&\mbox{Venus):}&\mbox{\phantom{xx}
Ven}=[0.7230,0.7236]\times[0.6149,0.6155],\\\
&\mbox{Earth:}&\mbox{\phantom{xx} Erd}=[1,1]\times[1,1],\\\
&\mbox{Mars:}&\mbox{\phantom{xx} Mar}=[1.5234,1.5240]\times[1.8805,1.8811],\\\
&\mbox{Jupiter:}&\mbox{\phantom{xx}
Jup}=[5.2041,5.2047]\times[11.8617,11.8623],\\\
&\mbox{Saturn:}&\mbox{\phantom{xx}
Sat}=[9.5823,9.5829]\times[29.4568,29.4574].\end{array}\right\\}$ (2.13)

Determine whether there exists a quadratic curve on the plane passing through
the six given rectangles.

First, we take the centers of the six rectangles as
$(x_{i}^{\ast},y_{i}^{\ast})$, and represent each rectangle using
$x_{i}=x_{i}^{\ast}+\Delta x_{i}$ and $y_{i}=y_{i}^{\ast}+\Delta y_{i}$.

$\displaystyle\mbox{Mer}:$ $\displaystyle x_{1}=\frac{3871}{10000}+\Delta
x_{1},\quad y_{1}=\frac{2408}{10000}+\Delta y_{1},$ $\displaystyle\mbox{Ven}:$
$\displaystyle x_{2}=\frac{7233}{10000}+\Delta x_{2},\quad
y_{2}=\frac{6152}{10000}+\Delta y_{2},$ $\displaystyle\mbox{Erd}:$
$\displaystyle x_{3}=1,\quad y_{3}=1,$ $\displaystyle\mbox{Mar}:$
$\displaystyle x_{4}=\frac{15237}{10000}+\Delta x_{4},\quad
y_{4}=\frac{18808}{10000}+\Delta y_{4},$ $\displaystyle\mbox{Jup}:$
$\displaystyle x_{5}=\frac{52044}{10000}+\Delta x_{5},\quad
y_{5}=\frac{11862}{10000}+\Delta y_{5},$ $\displaystyle\mbox{Sat}:$
$\displaystyle x_{6}=\frac{95826}{10000}+\Delta x_{6},\quad
y_{6}=\frac{294571}{10000}+\Delta y_{6}.$

Secondly, substitute $x_{i},y_{i}$ to calculate the determinant below

$M_{2}=\left|\begin{array}[]{cccccc}1&{x_{1}}&{y_{1}}&{{x_{1}}}^{2}&{x_{1}}\,{y_{1}}&{{y_{1}}}^{2}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
1&{x_{2}}&{y_{2}}&{{x_{2}}}^{2}&{x_{2}}\,{y_{2}}&{{y_{2}}}^{2}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr
1&{x_{3}}&{y_{3}}&{{x_{3}}}^{2}&{x_{3}}\,{y_{3}}&{{y_{3}}}^{2}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr
1&{x_{4}}&{y_{4}}&{{x_{4}}}^{2}&{x_{4}}\,{y_{4}}&{{y_{4}}}^{2}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr
1&{x_{5}}&{y_{5}}&{{x_{5}}}^{2}&{x_{5}}\,{y_{5}}&{{y_{5}}}^{2}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr
1&{x_{6}}&{y_{6}}&{{x_{6}}}^{2}&{x_{6}}\,{y_{6}}&{{y_{6}}}^{2}\end{array}\right|.$

Then $M_{2}(\Delta x_{1},\Delta y_{1},\cdots,\Delta x_{6},\Delta y_{6})$ is a
polynomial of degree 8 with 720 monomials:

$M_{2}={\frac{728327184304772538780278933931}{781250000000000000000000000000}}+\sum_{j=1}^{8}H_{k}(\Delta
x_{1},\Delta y_{1},\ldots,\Delta x_{6},\Delta y_{6}),$

Here, $H_{k}(\Delta x_{1},\Delta y_{1},\ldots,\Delta x_{6},\Delta
y_{6}),(k=1,2,\ldots,8)$ represents a $k$th-degree homogeneous polynomial. For
example:

$\displaystyle H_{1}$ $\displaystyle=c_{11}\Delta x_{1}+c_{12}\Delta
y_{1}+\ldots+c_{112}\Delta y_{6},\;$ ($12$ monomials) $\displaystyle H_{2}$
$\displaystyle=c_{21}\Delta x_{1}^{2}+\ldots+c_{26}\Delta
x_{6}^{2}+c_{27}\Delta y_{1}^{2}+\ldots+c_{212}\Delta y_{6}^{2}$
$\displaystyle\phantom{xxxxxx}+c_{213}\Delta x_{1}y_{1}+c_{214}\Delta
x_{1}y_{2}+\ldots+c_{248}\Delta x_{6}y_{6},\;$ ($48$ monomials)

where $c_{11},\ldots,c_{112},c_{21},\ldots,c_{248}$ are specific rational
numbers. For brevity, their exact values are not shown here. It is easy to see
that for any positive number $\delta>0$ and variables $\Delta x_{1},\Delta
y_{1},\cdots,\Delta x_{5},\Delta y_{5}\in[-\delta,\delta]$, the following
inequalities hold:

$\displaystyle H_{1}$
$\displaystyle\geq-(|c_{11}|+|c_{12}|+\ldots+|c_{112}|)\delta,\;$
$\displaystyle H_{2}$
$\displaystyle=-(|c_{21}|+\ldots+|c_{27}|+\ldots+|c_{213}|+\ldots+|c_{248}|)\delta^{2}.$

Similarly, for $H_{3},\ldots,H_{8}$, when $\Delta x_{1},\Delta
y_{1},\cdots,\Delta x_{5},\Delta y_{5}\in[-\delta,\delta]$, the following
inequalities hold:

$H_{k}\geq-C_{k}\times\delta^{k},\quad 3\leq k\leq 8,$

where $C_{k}$ is a positive number determined by the coefficients of the
polynomial $H_{k}$. By substituting specific numerical values, we obtain the
inequality:

$\displaystyle M_{2}$ $\displaystyle\geq$
$\displaystyle{\frac{728327184304772538780278933931}{781250000000000000000000000000}}-{\frac{15486384648337660843877367483}{6250000000000000000000000}}\,{\delta}$
$\displaystyle-{\frac{4025428097433897230769192269}{125000000000000000000000}}\,{\delta}^{2}-{\frac{2430088906604087412321393}{12500000000000000000}}\,{\delta}^{3}$
$\displaystyle-{\frac{264230612082325102439}{625000000000000}}\,{\delta}^{4}-{\frac{155328292156034629}{500000000000}}\,{\delta}^{5}$
$\displaystyle-{\frac{2943856054691}{50000000}}\,\delta^{6}-{\frac{11662287}{2500}}\,{\delta}^{7}-120\,{\delta}^{8}:=k_{2}(\delta)$

It is easy to see that $k_{2}(0)=0.932258\cdots>0$, the smallest positive real
root of $k_{2}(\delta)=0$ is $0.000374\cdots$, and the largest negative real
root is $-1.024959\cdots$. This implies that when $-1.024960<\delta<0.000374$

$M_{2}>k_{2}(\delta)>0$

which means that for $(x_{i},y_{i})$, $1\leq i\leq 6$, within the following
rectangle:

$[x_{i}^{\ast}-\frac{36}{100000},x_{i}^{\ast}+\frac{36}{100000}]\times[y_{i}^{\ast}-\frac{36}{100000},y_{i}^{\ast}+\frac{36}{100000}]$
(2.14)

where $(x_{i}^{\ast},y_{i}^{\ast})$ represents the center of each of the six
small rectangles, the polynomial $M_{2}(x_{1},y_{1},\ldots,x_{6},y_{6})$ is
positive definite.

Finally, it is verified that the square regions $(r,p)$ of the six planets
given in (2.13) (with side length $0.0003$) are all contained within the six
rectangular regions described in (2.14) (with side length $0.00036$ and
identical centers). In other words, no quadratic curves in terms of $r$ and
$p$ pass through the six square regions in (2.13).

##  3 The sparsity aspect of approximate implicit interpolation

Consider a set of $N$ points $P_{i}=(x_{i},y_{i})$ on the two-dimensional
plane $\mathbb{R}^{2}$ and a given positive integer $d$. In the previous
section, we learned that if $d>\lfloor\sqrt{2N}\rfloor$, there always exists
an algebraic curve of degree $d$ passing through the given points. However, if
$d<\lfloor\sqrt{2N}\rfloor$, there may not be such a curve passing through all
of these points simultaneously. In this section, we discuss how to find a
polynomial curve $f(x,y)=0$ of degree at most $d$ such that the difference
between these points and the curve, represented by
$\sum_{i=1}^{N}f(x_{i},y_{i})^{2}$, is minimized.

To approach this problem, we assume that $f(x,y)$ IS expressed as:

$f(x,y)=c_{0}+c_{1}\,x+c_{2}\,y+\cdots+c_{p-d}\,x^{k}+\cdots+c_{p-1}\,y^{d},$

where $c_{0},c_{1},\cdots,c_{p-1}$ are real numbers satisfying the condition:

$c_{0}^{2}+c_{1}^{2}+\cdots+c_{p-1}^{2}=1.$

To solve this optimization problem, we employ the Lagrange multiplier method.
The problem can be formulated as follows:

$\displaystyle\left\\{\begin{array}[]{ccl}&{\displaystyle\min}&{\displaystyle
g(c_{0},c_{1},\ldots,c_{p-1}):=\sum_{i=1}^{N}}(c_{0}+c_{1}x_{i}+c_{2}y_{i}+\cdots+c_{p-1}y_{i}^{d})^{2}\\\
&s.t.&c_{0}^{2}+c_{1}^{2}+\cdots+c_{p-1}^{2}=1.\end{array}\right.$ (3.3)

The Lagrangian function for this problem is defined as:

$\mathcal{L}:=\sum_{i=1}^{N}\left(c_{0}+c_{1}x_{i}+c_{2}y_{i}+\cdots+c_{p-1}y_{i}^{d}\right)^{2}+\lambda\cdot(c_{0}^{2}+c_{1}^{2}+\cdots+c_{p-1}).$

By calculating the partial derivatives of $\mathcal{L}$ with respect to
$c_{0},c_{1},\cdots,c_{p-1}$, we obtain the first-order conditions for the
optimization problem (3.3):

$\frac{\partial\mathcal{L}}{\partial
c_{0}}=0,\quad\frac{\partial\mathcal{L}}{\partial
c_{1}}=0,\quad\frac{\partial\mathcal{L}}{\partial
c_{2}}=0,\quad\cdots,\quad\frac{\partial\mathcal{L}}{\partial c_{p-1}}=0.$

These equations can be expressed in matrix form as:

$\left(\frac{\partial^{2}\mathcal{L}}{\partial c_{i}\partial
c_{j}}\right)_{i,j=0,1,\cdots,p-1}\left(\begin{array}[]{l}c_{0}\\\ c_{1}\\\
\vdots\\\ c_{p-1}\end{array}\right)=\left(\begin{array}[]{c}0\\\ 0\\\
\vdots\\\ 0\end{array}\right),$ (3.4)

where

$\frac{\partial^{2}\mathcal{L}}{\partial c_{i}\partial
c_{j}}\,=\left\\{\begin{array}[]{ll}\dfrac{\partial^{2}g}{\partial
c_{i}^{2}}+2\lambda,&\mbox{if $j=i$},\\\\[7.0pt]
\dfrac{\partial^{2}g}{\partial c_{i}\partial c_{j}},&\mbox{if
$j\not=i$},\end{array}\right.\quad(i,j=0,1,\cdots,p-1).$

In this way, by substituting the right-hand sides of the two equations into
the equation

$\det\left(\frac{\partial^{2}\mathcal{L}}{\partial c_{i}\partial
c_{j}}\right)=0$

we obtain a polynomial equation in $\lambda$:

$a_{0}+a_{1}\lambda+a_{2}\lambda^{2}+\cdots+a_{p-1}\lambda^{p-1}+2^{p}\lambda^{p}=0,$
(3.5)

which is actually the characteristic equation of the matrix

$M(x_{1},y_{1},\ldots,x_{N},y_{N})=\left(-\frac{1}{2}\frac{\partial^{2}g}{\partial
c_{i}\partial c_{j}}\right)$

Since the matrix $M$ is a real symmetric matrix, the equation (3.5) has at
least one real root. After solving for $\lambda$, we substitute it into the
homogeneous linear equation (3.4) to solve for the desired coefficients
$c_{0},c_{1},\ldots,c_{p-1}$.

when $N<p-1$, the coefficients $a_{0},\cdots,a_{p-1-N}$in (3.5) are zero,
resulting in the solution to (3.3) forming a linear subspace of dimension
$p-N$ in $\mathbb{R}^{p}$. The intersection of this linear space with the unit
sphere, defined as:

$S^{p-1}:=\\{(c_{0},c_{1},\cdots,c_{p-1})\,|\,c_{0}^{2}+c_{1}^{2}+\cdots+c_{p-1}^{2}=1\\}$

is an infinite set. Therefore, further exploration is required to find the
sparse polynomial $f(x,y)$ that passes through the given points. For instance,
when considering the 10 points $P_{1},P_{2},\cdots,P_{10}$ as provided in
Example 2.7, any quartic curve

$f(x,y)=c_{0}+c_{1}x+c_{2}y+\cdots+c_{10}x^{4}+c_{11}x^{3}y+\cdots+c_{14}y^{4}$

passing through these 10 points satisfies the equation:

$(c_{0},c_{1},\cdots,c_{9})^{T}=M_{10\times
5}(c_{10},c_{11},\cdots,c_{14})^{T},$ (3.6)

where

$M_{10\times
5}=\left(\begin{array}[]{clcrc}\scalebox{0.8}{$0$}&-{\frac{95}{47}}-{\frac{20\,\sqrt{2}}{47}}&-{\frac{56}{47}}+{\frac{8\,\sqrt{2}}{47}}&{\frac{40}{141}}+{\frac{16\,\sqrt{2}}{47}}&\scalebox{0.8}{$0$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\scalebox{0.8}{$0$}&-{\frac{25\,\sqrt{2}}{282}}-{\frac{10}{47}}&-{\frac{35}{141}}+{\frac{5\,\sqrt{2}}{141}}&-{\frac{70}{141}}+{\frac{10\,\sqrt{2}}{141}}&\scalebox{0.8}{$0$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\scalebox{0.8}{$0$}&-{\frac{15\,\sqrt{2}}{94}}-{\frac{242}{141}}&-{\frac{21}{47}}+{\frac{3\,\sqrt{2}}{47}}&-{\frac{136}{47}}+{\frac{6\,\sqrt{2}}{47}}&\scalebox{0.8}{$0$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\scalebox{0.8}{$0$}&-{\frac{119\,\sqrt{2}}{282}}-{\frac{10}{47}}&{\frac{59}{141}}-{\frac{14\,\sqrt{2}}{47}}&{\frac{10\,\sqrt{2}}{141}}+{\frac{8}{47}}&\scalebox{0.8}{$0$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\scalebox{0.8}{$0$}&{\frac{125\,\sqrt{2}}{141}}+{\frac{488}{141}}&{\frac{350}{141}}-{\frac{50\,\sqrt{2}}{141}}&{\frac{277}{141}}-{\frac{100\,\sqrt{2}}{141}}&\scalebox{0.8}{$0$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr-\frac{3}{2}&{\frac{119\,\sqrt{2}}{282}}+{\frac{10}{47}}&-{\frac{200}{141}}+{\frac{14\,\sqrt{2}}{47}}&-{\frac{10\,\sqrt{2}}{141}}-{\frac{8}{47}}&\scalebox{0.8}{$0$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr-{\frac{34}{9}}&-{\frac{427\,\sqrt{2}}{423}}-{\frac{445}{141}}&-{\frac{623\,\sqrt{2}}{141}}-{\frac{1252}{423}}&-{\frac{824\,\sqrt{2}}{423}}+{\frac{356}{141}}&{\frac{77}{9}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr-\frac{5}{9}&-{\frac{35\,\sqrt{2}}{423}}-{\frac{25}{47}}&-{\frac{380}{423}}-{\frac{35\,\sqrt{2}}{47}}&-{\frac{160\,\sqrt{2}}{423}}+{\frac{20}{47}}&{\frac{10}{9}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\frac{3}{2}&{\frac{7\,\sqrt{2}}{47}}+{\frac{45}{47}}&{\frac{40}{141}}+{\frac{63\,\sqrt{2}}{47}}&{\frac{32\,\sqrt{2}}{47}}-{\frac{36}{47}}&-\scalebox{0.8}{$5$}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{\frac{17}{6}}&{\frac{133\,\sqrt{2}}{141}}+{\frac{385}{141}}&{\frac{168}{47}}+{\frac{539\,\sqrt{2}}{141}}&-{\frac{308}{141}}+{\frac{232\,\sqrt{2}}{141}}&-{\frac{17}{3}}\end{array}\right).$

And therefore, the search of the sparse quartic polynomial $f(x,y)=0$ that
passes the 10 points leads to an optimization problem as follows:

$\displaystyle\min$ $\displaystyle||(c_{0},c_{1},c_{2},\cdots,c_{14})||_{0},$
s.t. $\displaystyle(c_{0},c_{1},\cdots,c_{9})^{T}=M_{10\times
5}(c_{10},c_{11},\cdots,c_{14})^{T},$ $\displaystyle
c_{0}^{2}+c_{1}^{2}+\cdots+c_{14}^{2}=1,$

where the zero-norm $||\cdot||_{0}$ represents the zero-norm of a vector
(i.e., the number of its non-zero elements). This optimization problem aims to
find a sparse quartic polynomial by minimizing the number of non-zero
coefficients.

To solve this problem, we can employ symbolic computation software such as
Maple. By using Maple, we can verify the following facts:

  * •

For any selection $j_{1},j_{2},\cdots,j_{11}$ with $0\leq
j_{1}<j_{2}<\cdots<j_{11}\leq 14$ the equation system formed by (3.6) has a
non-zero solution satisfying

$c_{j_{1}}=c_{j_{2}}=\cdots=c_{j_{11}}=0$

There are ${15\choose 11}=1365$ such selections.

  * •

Among all ${15\choose 10}=3003$ subsets $S\subset\\{0,1,2,\cdots,14\\}$ with
$\\#(S)=10$, there exist exactly two subsets $S_{1},S_{2}$ that satisfy (3.6)
and $\\{c_{j}=0\,|\,j\in S\\}$ with non-zero solutions $N\\!Z$. They are:

$S_{1}=\\{0,1,2,3,4,6,8,11,12,13\\},$
$N\\!Z_{1}=\\{c_{5}=\frac{1}{\sqrt{19}},c_{7}=-\frac{3}{\sqrt{19}},c_{9}=-\frac{2}{\sqrt{19}},c_{10}=\frac{2}{\sqrt{19}},c_{14}=\frac{1}{\sqrt{19}}\\}$

$S_{2}=\\{1,3,4,6,7,8,10,11,12,13\\},$
$N\\!Z_{2}=\\{c_{0}=\frac{5/2}{\sqrt{671}},c_{2}=-\frac{51/4}{\sqrt{671}},c_{5}=\frac{77/4}{\sqrt{671}},c_{9}=-\frac{45/4}{\sqrt{671}},c_{14}=\frac{9/4}{\sqrt{671}}\\}.$

The corresponding curves are given by the equations

${y}^{2}-3\,{x}^{2}y-2\,{y}^{3}+2\,{x}^{4}+{y}^{4}=0,$

and

$\frac{5}{2}-{\frac{51}{4}}\,y+{\frac{77}{4}}\,{y}^{2}-{\frac{45}{4}}\,{y}^{3}+\frac{9}{4}\,{y}^{4}=0.$

The previous quartic curve is consistent with the curve obtained by implicit
function interpolation using 15 points in Example 2.7. The latter curve is
actually composed of four lines:

$(y-1)(y-2)(3y-1)(3y-5)=0.$

When the degree $d$ is large, solving the optimization problem (3.6) directly
through exhaustive search becomes computationally challenging. In such cases,
we can approximate the solution by solving the following convex optimization
problem:

$\displaystyle\min$ $\displaystyle\sum_{0\leq i<j\leq
14}c_{i}^{2}c_{j}^{2}=\frac{1}{2}~{}(-c_{0}^{4}-c_{1}^{4}-\cdots-c_{14}^{4}),$
s.t. $\displaystyle(c_{0},\cdots,c_{9})^{T}=M_{10\times
5}(c_{10},\cdots,c_{14})^{T},$ $\displaystyle
c_{0}^{2}+c_{1}^{2}+\cdots+c_{14}^{2}=1.$

By formulating the problem in this way, we transform it into a convex
optimization problem that can be efficiently solved using available
optimization algorithms. The resulting solution provides an approximate sparse
quartic polynomial that satisfies the given conditions. For some practical
applications, obtaining an approximate solution with a certain level of
accuracy within the allowed time is sufficient. For theoretical problems, an
approximate solution can also serve as a foundation for further research.

To illustrate our approach and demonstrate its practicality and effectiveness
in real-world scenarios, we further investigate the analysis of planetary
motion data in Example 2.7.

######  Example 3.1.

Analysis of the actual motion data of planets and asteroids in the solar
system.

Consider the data of the eight major planets, 23 minor planets, and 7 dwarf
planets orbiting the solar system, including the semi-major axis lengths of
their orbits and the orbital periods (data source:
https://en.wikipedia.org/wiki/Planet). In the previous example (Example 2.14)
of the previous section, data for the first 6 planets, which are closer to
Earth, were provided. The data for Uranus, Neptune, minor planets, and dwarf
planets are given in Table 2. The observational errors in the data, denoted by
$r$, are within $\pm 0.001AU$, and the errors in the orbital periods, denoted
by $p$, are within $\pm 0.01$ years.

Table 2: Orbit data coordinates for Uranus, Neptune, 23 asteroids, and 7 dwarf
planets.

asteroid |  $r$ (in Au) |  $p$ (in yr) | asteroid |  $r$ (in Au) |  $p$ (in yr)
---|---|---|---|---|---
Uranus($\ast\ast$) | 19.2184 | 84.0205 |  Neptune($\ast\ast$) | 30.1104 | 164.7693
Ceres($\ast$) | 2.768 | 4.60 | Pallas | 2.772 | 4.61
Juno | 2.671 | 4.36 | Vesta | 2.362 | 3.63
Astraea | 1.574 | 1.97 | Hebe | 2.426 | 3.78
Iris | 2.385 | 3.68 | Flora | 2.202 | 3.27
Metis | 2.387 | 3.69 | Hygiea | 3.1421 | 5.57
Parthenope | 2.453 | 3.84 | Victoria | 2.334 | 3.57
Egeria | 2.576 | 4.13 | Irene | 2.150 | 3.16
Eunomia | 2.643 | 4.30 | Psyche | 2.921 | 4.99
Thetis | 2.470 | 3.88 | Melpomene | 2.296 | 3.48
Fortune | 2.441 | 3.81 | Proserpina | 2.656 | 4.33
Bellona | 2.777 | 4.63 | Amphitrite | 2.554 | 4.08
Leukothea | 2.990 | 5.17 | Fides | 2.641 | 4.29
Chiron($\ast$) | 13.708 | 50.76 |  Pluto($\ast$) | 39.48 | 247.94
Haumea($\ast$) | 43.218 | 284.12 |  Makemake($\ast$) | 45.715 | 309.09
Eris($\ast$) | 67.781 | 558.04 |  Sedna ($\ast$) | 506.8 | 11411

In the table, the symbol $*$ represents a dwarf planet, and the symbol $**$
represents a planet.

Figure 2 presents the orbital data for planets in the solar system and some
selected minor and dwarf planets shown in Table 2. The theoretical dimensions
of each rectangle in the figure are $0.02\times 0.002$. For clarity, the
actual display uses the $\blacksquare$ symbol from the tiny size package in
LaTeX, which is horizontally and vertically magnified by a factor of $2$ and
$20$ respectively.

From the previous section’s Example 2.14, we know that no quadratic curve
passes through these rectangular neighborhoods simultaneously. Our goal is to
determine whether there exist cubic curves that pass through the $6$ small
rectangular neighborhoods mentioned in (2.13), as well as $32$ small
rectangular neighborhoods centered around the data in Table 2 with dimensions
$0.002\times 0.01$.

If such cubic curves exist, we aim to find a curve with the fewest possible
terms. The ideal scenario would involve re-discovering Kepler’s Third Law
using rigorous analytical methods based on this actual data.

$p$$r$$1$$2$$3$$4$$5$$1$$2$$3$$\blacksquare$Mercury$\blacksquare$Venus$\blacksquare$Earth$\blacksquare$Mars$\blacksquare$Astraea$\blacksquare$Irene/Flora$\blacksquare$$\blacksquare$Melpomene/Victoria/Vesta/Metis$\blacksquare$$\blacksquare$$\blacksquare$$\blacksquare$Hebe/Fortune/Parthenope/Thetis$\blacksquare$$\blacksquare$$\blacksquare$$\blacksquare$Amphitrite/Egeria$\blacksquare$$\blacksquare$Fides/Proserpina/Juno$\blacksquare$$\blacksquare$$\blacksquare$Ceres/Bellona$\blacksquare$$\blacksquare$Psyche$\blacksquare$Leukothea$\blacksquare$Hygiea

Figure 2: The orbital periods of some planets and asteroids in the solar
system and the semi-major axes of their orbits. The size of $\blacksquare$ is
adjusted to reflect the real error of the data.

To achieve this goal, our focus is on constructing cubic curves of the form:

$k(r,p)=c_{0}+c_{1}r+c_{2}p+c_{3}r^{2}+c_{4}rp+c_{5}p^{2}+c_{6}r^{3}+c_{7}r^{2}p+c_{8}rp^{2}+c_{9}p^{3}=0,$
(3.9)

such that they pass through the point $(1,1)$ and the 37 other small square
regions. We aim to maximize the number of coefficients
$c_{0},c_{1},\ldots,c_{9}\in\mathbb{R}$ that are zero, while ensuring that at
least one of $c_{6},c_{7},c_{8},c_{9}$ is non-zero.

Without loss of generality, we assume that the first non-zero coefficient
$c_{i}$ is positive, and the sum of squares of all coefficients is equal to 1,
i.e.,

$c_{0}^{2}+c_{1}^{2}+\cdots+c_{9}^{2}=1.$

To avoid accumulation of errors caused by rounding during calculations, we
represent the 38 data points as rational points $P_{i}=(x_{i},y_{i})$, where
$i=1,2,\ldots,38$, as shown in Table 3.

Table 3: Rational representation of $38$ data points.
$\begin{array}[]{|c|c|c|c|c|c|c|c|c|c|}\hline\cr\phantom{jj}P_{1}\phantom{hh}&\phantom{jj}P_{2}\phantom{hh}&\phantom{jj}P_{3}\phantom{hh}&\phantom{jj}P_{4}\phantom{hh}&\phantom{jj}P_{5}\phantom{hh}&\phantom{jj}P_{6}\phantom{hh}&\phantom{jj}P_{7}\phantom{hh}&\phantom{jj}P_{8}\phantom{hh}&\phantom{jj}P_{9}\phantom{hh}&\phantom{j}P_{10}\phantom{hh}\\\
\hline\cr\frac{3871}{10000}&\frac{7233}{10000}&{\small
1}&\frac{15237}{10000}&\frac{13011}{2500}&\frac{47913}{5000}&\frac{24023}{1250}&\frac{3011}{100}&\frac{987}{25}&\frac{346}{125}\\\
\hline\cr\frac{301}{1250}&\frac{769}{1250}&{\small
1}&\frac{2351}{1250}&\frac{5931}{500}&\frac{294571}{10000}&\frac{168041}{2000}&\frac{1647693}{10000}&\frac{12397}{50}&\frac{23}{5}\\\
\hline\cr
P_{11}&P_{12}&P_{13}&P_{14}&P_{15}&P_{16}&P_{17}&P_{18}&P_{19}&P_{20}\\\
\hline\cr\frac{693}{250}&\frac{2671}{1000}&\frac{1181}{500}&\frac{787}{500}&\frac{1213}{500}&\frac{477}{200}&\frac{1101}{500}&\frac{2387}{1000}&\frac{31421}{10000}&\frac{2453}{1000}\\\
\hline\cr\frac{461}{100}&\frac{109}{25}&\frac{363}{100}&\frac{197}{100}&\frac{189}{50}&\frac{92}{25}&\frac{327}{100}&\frac{369}{100}&\frac{557}{100}&\frac{96}{25}\\\
\hline\cr
P_{21}&P_{22}&P_{23}&P_{24}&P_{25}&P_{26}&P_{27}&P_{28}&P_{29}&P_{30}\\\
\hline\cr\frac{1167}{500}&\frac{322}{125}&\frac{43}{20}&\frac{2643}{1000}&\frac{2921}{1000}&\frac{247}{100}&\frac{287}{125}&\frac{2441}{1000}&\frac{332}{125}&\frac{2777}{1000}\\\
\hline\cr\frac{357}{100}&\frac{413}{100}&\frac{79}{25}&\frac{43}{10}&\frac{499}{100}&\frac{97}{25}&\frac{87}{25}&\frac{381}{100}&\frac{433}{100}&\frac{463}{100}\\\
\hline\cr P_{31}&P_{32}&P_{33}&P_{34}&P_{35}&P_{36}&P_{37}&P_{38}&&\\\
\hline\cr\frac{1277}{500}&\frac{299}{100}&\frac{2641}{1000}&\frac{3427}{250}&\frac{67781}{1000}&\frac{2534}{5}&\frac{21609}{500}&\frac{9143}{200}&&\\\
\hline\cr\frac{102}{25}&\frac{517}{100}&\frac{429}{100}&\frac{1269}{25}&\frac{13951}{25}&{\small
11411}&\frac{7103}{25}&\frac{30909}{100}&&\\\ \hline\cr\end{array}$

It can be easily seen that if any four points
$[x_{1},y_{1}],\ldots,[x_{4},y_{4}]$ from the set of 38 data points lie on a
cubic curve given by
$f(x,y)=c_{0}+c_{1}x+c_{2}y+\cdots+c_{6}x^{3}+\cdots+c_{9}y^{3}$, then the
coefficients $c_{0},c_{1},\ldots,c_{9}$ satisfy the following system of linear
equations:

$\left(\begin{array}[]{cccccccccc}1&x_{1}&y_{1}&x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}&x_{1}^{3}&x_{1}^{2}y_{1}&x_{1}y_{1}^{2}&y_{1}^{3}\\\
1&x_{2}&y_{2}&x_{2}^{2}&x_{2}y_{2}&y_{2}^{2}&x_{2}^{3}&x_{2}^{2}y_{2}&x_{2}y_{2}^{2}&y_{2}^{3}\\\
1&x_{3}&y_{3}&x_{3}^{2}&x_{3}y_{3}&y_{3}^{2}&x_{3}^{3}&x_{3}^{2}y_{3}&x_{3}y_{3}^{2}&y_{3}^{3}\\\
1&x_{4}&y_{4}&x_{4}^{2}&x_{4}y_{4}&y_{4}^{2}&x_{4}^{3}&x_{4}^{2}y_{4}&x_{4}y_{4}^{2}&y_{4}^{3}\end{array}\right)\left(\begin{array}[]{c}c_{0}\\\
c_{1}\\\ \vdots\\\ c_{9}\end{array}\right)=\left(\begin{array}[]{c}0\\\ 0\\\
\vdots\\\ 0\end{array}\right).$ (3.10)

First, we construct nine quadruples $G_{1},G_{2},\ldots,G_{9}$ using the 38
data points. Among them, $G_{1}$ and $G_{2}$ contain the eight major planets
of the solar system, $G_{2}$ to $G_{8}$ include 23 minor planets and the dwarf
planet Ceres (the largest known minor planet in the asteroid belt). $G_{9}$
consists of four dwarf planets.

In Table 3, two dwarf planets (Eris and Sedna) are not included in the
grouping and are reserved for validation purposes. The specific construction
of the nine quadruples is as follows.

Giant planets $\displaystyle G_{1}:=\\{\mbox{\rm Mercury, Venus, Earth,
Mars}\\};$ $\displaystyle G_{2}:=\\{\mbox{\rm Jupiter, Saturn, Uranus,
Neptune}\\};$ Asteroids $\displaystyle G_{3}:=\\{\mbox{\rm Pallas, Juno.
Vesta, Astraea}\\};$ $\displaystyle G_{4}:=\\{\mbox{\rm Hebe, Iris, Flora,
Metis}\\};$ $\displaystyle G_{5}:=\\{\mbox{\rm Hygiea, Parthenope, Victoria,
Egeria}\\};$ $\displaystyle G_{6}:=\\{\mbox{\rm Irene, Eunomia, Psyche,
Thetis}\\};$ $\displaystyle G_{7}:=\\{\mbox{\rm Melpomene, Fortune,
Proserpina, Bellona}\\};$ $\displaystyle G_{8}:=\\{\mbox{\rm Amphitrite,
Leukothea, Fides, Ceres}\\};$ Dwarf planets $\displaystyle G_{9}:=\\{\mbox{\rm
Chiron, Pluto, Haumea, Makemake}\\}.$

Next, for each group of four stars, we perform the following calculations:
Construct a set $S_{(10,5)}$ of five integers ${j_{1},j_{2},\cdots,j_{5}}$
satisfying $0\leq j_{1}<j_{2}<\cdots<j_{5}\leq 5$.

Consider the set $J={j_{1},j_{2},\cdots,j_{5}}$ consisting of five distinct
integers such that $0\leq j_{1}<j_{2}<\cdots<j_{5}\leq 9$, and the set
${6,7,8,9}$ is not a subset of $J$. We construct the set $S_{(10,5)}$ of all
such combinations. The number of elements in $S_{(10,5)}$ is 246. For each
such set of five integers (Integer Quintet) $J$, we define:

$\\{j_{6},j_{7},\cdots,j_{10}\\}=\\{0,1,2,\cdots,9\\}\setminus
J,~{}~{}j_{6}<j_{7}<\cdots<j_{10}.$

Then by substituting $c_{j_{1}}=c_{j_{2}}=c_{j_{3}}=c_{j_{4}}=c_{j_{5}}=0$
into the equation system (3.10) we get four homogeneous equations with $5$
unknowns $j_{6},j_{7},\cdots,j_{10}$. For each non-zero solution
$c_{j_{6}},c_{j_{7}},\ldots,c_{j_{10}}$ of this derived system of equations,
we obtain a non-zero solution $\phi(J):=(c_{0},c_{1},\cdots,c_{9})$ for
(3.10), where at least five of the coefficients are zero. For example, for the
four-star group $G_{1}={\text{Mercury, Venus, Earth, Mars}}$ and the next two
sets of five integers from $S_{(10,5)}$:

$J_{1}=\\{0,1,2,3,7\\},~{}J_{2}=\\{[0,3,5,6,9]\\}$

the solutions to the system of equations (3.10) are respectively:

$\phi(J_{1})=\scalebox{0.866}{
$(0,0,0,0,-0.005115\cdots,-.701381\cdots,.712633\cdots,0,-0.012353\cdots,0.006218\cdots),$
}$ $\phi(J_{2})=\scalebox{0.866}{
$(0,0.154246\cdots,-0.407426\cdots,0,0.582026\cdots,0,0,-0.621253\cdots,0.292406\cdots,0)$.
}$

It took only $0.514$ seconds to generate $\phi(J)$ for all $J\in S_{(10,5)}$.

For each $\phi(J)=(c_{0},c_{1},\cdots,c_{9}),$ we define

$\psi(J)=\\{i~{}|~{}0\leq i\leq
9~{}\text{and}~{}|c_{i}|<\frac{1}{100\sqrt{10}}\\}.$

Here, we take $1/100\sqrt{10}$ as the threshold to set $c_{i}$ to zero if its
absolute value is smaller than this value, based on the following
consideration: When the system of equations (3.10) and $J$ are generated using
a random number generator, $(c_{0},c_{1},\ldots,c_{9})$ is uniformly
distributed on the 10-dimensional unit sphere, and the mean value of each
$|c_{i}|$ is $1/\sqrt{10}$. The probability that the absolute value is smaller
than $1/(100\sqrt{10})$ is approximately given by

$p=\frac{S_{10}(-\frac{1}{100\sqrt{10}},-\frac{1}{100\sqrt{10}})}{S_{10}}\approx\frac{\frac{2}{100\sqrt{10}}\times
V_{9}}{S_{10}}=\frac{64}{7875\sqrt{10}\cdot\pi}\approx 0.000818\cdots,$

where $S_{10}$ is the surface area of the 10-dimensional unit sphere and
$V_{9}$ is the volume of the 9-dimensional unit sphere, given by

$S_{10}=\frac{2\pi^{10/2}}{\Gamma(\frac{10}{2})}=\frac{\pi^{5}}{12},\quad
V_{9}=\frac{\pi^{9/2}}{\Gamma(\frac{9}{2}+1)}=\frac{32\pi^{4}}{945},$

and $S_{10}(-\delta,\delta)$ is the domain consisting of points on the
10-dimensional unit sphere within a distance of $\delta$ from the equator.
When $\delta\approx 0$, we have

$S_{10}(-\delta,\delta)\approx 2\delta\times V_{9}.$

It is evident that $J\subset\psi(J)$. When $J=\psi(J)$, the solutions of the
linear system (3.10) are all non-zero elements. When $J\subsetneq\psi(J)$, the
solutions of the linear system (3.10) contain new zero elements, i.e.,
$\\#\psi(J)\geq 6$. Since we are looking for cubic curves with coefficients
$c_{0},c_{1},\ldots,c_{9}$ that have as many zero elements as possible, for
the set $S_{(10,5)}$ containing 246 elements, we first identify the subset of
index sets $J$ in $\mathcal{S}^{\prime}$ that satisfy the condition
$J\subsetneq\psi(J)$, and then compute
$\mathcal{S}^{\prime\prime}={\psi(J)|J\in\mathcal{S}^{\prime}}$ and discuss
them accordingly.

For $G_{1}$, computation shows that $\mathcal{S}^{\prime}$ contains 57
elements, while $\mathcal{S}^{\prime\prime}$ contains 18 elements. Based on
the number of zero coefficients, they can be further divided into three
categories:

$\displaystyle S_{(10,6)}=$
$\displaystyle\left\\{\scalebox{0.9}{$[0,1,2,3,4,6],[0,1,2,3,5,6],[0,1,2,3,5,9],[0,1,2,3,6,8],[0,1,2,3,6,9]$}\right.;$
(3.11)
$\displaystyle\phantom{xxx}\left.\scalebox{0.9}{$[0,1,2,3,4,5],[0,1,2,3,4,7],[0,1,2,3,4,9],[0,1,2,3,5,8]$};\right.$
$\displaystyle\phantom{xxx}\left.\scalebox{0.9}{$[0,1,2,3,8,9],[0,1,2,4,8,9],[0,1,2,7,8,9],[0,1,3,7,8,9]$}\right\\};$
$\displaystyle S_{(10,7)}=$
$\displaystyle\\{\scalebox{0.9}{$[0,1,2,3,4,7,9],[0,1,2,3,4,8,9],[0,1,2,3,7,8,9],[0,1,2,4,7,8,9]$}\\};$
(3.12) $\displaystyle S_{(10,8)}=$
$\displaystyle\\{\scalebox{0.9}{$[0,1,2,3,4,7,8,9]$}\\}.$ (3.13)

Note that each element $J=[j_{1},j_{2},\ldots,j_{k}]$ in
$S_{(10,k)},(k=6,7,8)$ contains $k$ numbers from $0,1,\ldots,9$, arranged in
ascending order. $S_{(10,6)}$ shows that among all cubic equations (3.9)
passing through the four rectangular regions of
$G_{1}={\text{Mercury},\text{Venus},\text{Earth},\text{Mars}}$, there are
precisely 6 zero elements in the coefficients. Disregarding the specific
numerical values of the non-zero coefficients, these equations must belong to
one of the following 13 forms:

$\displaystyle p^{2}+r^{2}p+rp^{2}+p^{3},\quad rp+r^{2}p+rp^{2}+p^{3},\quad
rp+r^{3}+r^{2}p+rp^{2},\phantom{xxx}$ $\displaystyle
rp+p^{2}+r^{2}p+rp^{2},\quad rp+p^{2}+r^{2}p+p^{3},\phantom{xxx}$
$\displaystyle r^{3}+r^{2}p+rp^{2}+p^{3},\quad p^{2}+r^{3}+rp^{2}+p^{3},\quad
p^{2}+r^{3}+r^{2}p+rp^{2},\quad rp+r^{3}+r^{2}p+p^{3},\phantom{xxx}$ (3.14)
$\displaystyle rp+p^{2}+r^{3}+r^{2}p,\quad r^{2}+p^{2}+r^{3}+r^{2}p,\quad
r^{2}+rp+p^{2}+r^{3},\quad p+rp+p^{2}+r^{3}.\phantom{xxx}$ (3.15)

Among these 13 forms, the first 5 forms are reducible and can be factored out
with either $r$ or $p$. According to the conclusion from Example 2.7 in the
previous section, these forms are actually impossible to occur. $S_{(10,7)}$
shows that among all cubic equations (3.9) passing through $G_{1}$, there are
exactly 7 zero elements in the coefficients, with 4 possible forms:

$c_{5}p^{2}+c_{6}r^{3}+c_{8}rp^{2},\quad
c_{5}p^{2}+c_{6}r^{3}+c_{7}r^{2}p,\quad c_{4}rp+c_{5}p^{2}+c_{6}r^{3},\quad
c_{3}r^{2}+c_{5}p^{2}+c_{6}r^{3}.$ (3.16)

Finally, $S_{(10,8)}$ shows that among all cubic equations (3.9) passing
through $G_{1}$, there are exactly 8 zero elements in the coefficients, and
they must have the form:

$c_{5}p^{2}+c_{6}r^{3}=0.$ (3.17)

Let $\chi({G_{1}})$ be the set of all $13$ forms of cubic polynomials included
in equations (3.14), (3.15), (3.16), and (3.17), namely:

$\displaystyle\chi({G_{1}})=\left\\{r^{3}+r^{2}p+rp^{2}+p^{3},\quad
p^{2}+r^{3}+rp^{2}+p^{3},\quad p^{2}+r^{3}+r^{2}p+rp^{2},\right.$
$\displaystyle\phantom{\chi({G_{1}})=xxx}rp+r^{3}+r^{2}p+p^{3},\quad
rp+p^{2}+r^{3}+r^{2}p,\quad r^{2}+p^{2}+r^{3}+r^{2}p,$
$\displaystyle\phantom{\chi({G_{1}})=xxx}r^{2}+rp+p^{2}+r^{3},\quad
p+rp+p^{2}+r^{3},\quad p^{2}+r^{3}+rp^{2},$
$\displaystyle\phantom{\chi({G_{1}})=xxx}\left.p^{2}+r^{3}+r^{2}p,\quad
rp+p^{2}+r^{3},\quad r^{2}+p^{2}+r^{3},\quad p^{2}+r^{3}\right\\}.$

We can now assert that among the cubic polynomials determined by the four
small rectangular regions of $G_{1}$, if at least $6$ of their coefficients
are zero, then these cubic polynomials, disregarding the specific values of
the non-zero coefficients, must belong to one of the $13$ forms included in
$\chi({G_{1}})$. Hence, we refer to $\chi({G_{1}})$ as the set of sparse cubic
curve types supported by the four-star configuration $G_{1}$. Currently, we
hope that the results obtained from the data through implicit function
interval interpolation, which corresponds formally to the polynomial
$c_{5}p^{2}+c_{6}r^{3}=0$, are also included in the set $\chi({G_{1}})$.

Apply the same method to $G_{2},G_{3},\ldots,G_{9}$. We aim to extract the
common supported sparse cubic curve types/ from $\chi(G_{i})$ for $1\leq i\leq
9$. The results reveal that:

$\chi(G_{1})=\chi(G_{2})=\cdots=\chi(G_{9}).$

Since the construction of $\chi(G_{i})$ for a selected four-star configuration
$G_{i}$ can be completed on average in $0.5$ seconds, we employ the Monte
Carlo method to randomly generate four-star configurations $G_{X}$ among the
$38$ data points listed in Table 3. Our expectation is to obtain a distinct
sparse type $\chi(G_{X})$ from $\chi(G_{i})$ for $i=1,2,\ldots,9$. After
conducting $1000$ experiments, we consistently obtained
$\chi(G_{X})=\chi(G_{1})$. It is worth noting that there are a total of
${38\choose 4}=73,815$ possible four-star configurations that can be
constructed from the $38$ data points. While it is feasible to perform a
thorough computation by traversing through all possibilities, our objective is
to find the simplest form of a cubic curve. Therefore, in this case, we
directly compute the best results obtained by using the binomial expression:

$k_{j}(r,p)=c_{5}p^{2}+c_{6}r^{3}$

to construct the implicit functions passing through all $38$ data points
$P_{i}=[x_{i},y_{i}]$ listed on page 3.

Substituting $x_{i}$ and $y_{i}$ into the equation, we obtain:

$\sum_{i=1}^{38}(y_{i}^{2}\cdot c_{5}+x_{i}^{3}c_{6})^{2}=\alpha
c_{5}^{2}+\beta c_{5}c_{6}+\gamma c_{6}^{2},$ (3.18)

where

$\displaystyle\alpha=33910002280102731124343238312183/(2\times 10^{15}),$
$\displaystyle\beta=42374075982683010624551482072423759/(125\times 10^{16}),$
$\displaystyle\gamma=4236065336850550592183712032020535827017/(25\times
10^{22}).$

The values of $c_{5},c_{6}$ that minimize $\alpha c_{5}^{2}+\beta
c_{5}c_{6}+\gamma c_{6}^{2}$ satisfy

$-\frac{c_{5}}{c_{6}}=\scalebox{0.9}{
${\displaystyle\frac{42374075982683010624551482072423759}{42387502850128413905429047890228750}}$
}\\\ =\scalebox{0.9}{
$\cfrac{1}{1+\cfrac{1}{3155+\cfrac{1}{1+\cfrac{1}{11+\cfrac{1}{75+\cdots}}}}}$}\approx\scalebox{0.9}{$-0.999683\cdots$},$

and the approximate fractions are

$1,\quad\frac{3155}{3156},\quad\frac{3156}{3157},\quad\frac{37871}{37883},\quad\frac{2843481}{2844382},\cdots.$

Using these values, we obtain the following cubic curves:

$k_{1}(r,p)=p^{2}-r^{3},\quad k_{2}(r,p)=p^{2}-\frac{3155}{3156}r^{3},\quad
k_{3}(r,p)=p^{2}-\frac{3156}{3157}r^{3},\quad\cdots$

Given that the errors associated with the data are smaller than $0.001$ AU for
$r$ and $0.01$ years for $p$, a box $B_{i}$ can be assigned to each point
$[x_{i},y_{i}]$ as follows:

$B_{i}:=\\{(x,y)\,|\,x_{i}-\frac{1}{1000}\leq x\leq
x_{i}+\frac{1}{1000},\,y_{i}-\frac{1}{100}\leq y\leq
y_{i}+\frac{1}{100}\\},\quad i=1,2,\cdots,38.$

The following conclusions can be verified:

  * •

The curve $k_{1}(r,p)=p^{2}-r^{3}$ intersects with each $B_{i}$, except for
$i\neq 5,6,7,8,38$.

  * •

The curves $k_{2}(r,p)=0$ and $k_{3}(r,p)=0$ intersect with each $B_{i}$,
except for $i=5,6,7,8,35,36,37,38$. Note that $P_{5},P_{6},P_{7},P_{8}$
correspond to the planets Jupiter, Saturn, Uranus, and Neptune in the solar
system, while $P_{35},P_{36},P_{37},P_{38}$ represent the dwarf planets
Haumea, Makemake, Eris, and Sedna, respectively.

It should be noted that the analysis of orbital data shows that the formulas
$k_{2}(r,p)=0$ and $k_{3}(r,p)=0$ are more accurate than the Kepler’s Third
Law $k_{1}(r,p)=p^{2}-r^{3}$. This is because in reality, the precise formula
relating $r$ and $p$ is given by:

$\frac{p^{2}}{r^{3}}=\frac{4\pi^{2}}{G(M_{\odot}+m)},$

where $M_{\odot}$ is the mass of the Sun, $m$ is the mass of the planet
(asteroid or dwarf planet), and $G$ is the gravitational constant. The
equation $p^{2}-r^{3}=0$ is only correct when $m<<M_{\odot}$.

##  4 The uniqueness problem of algebraic curve interval interpolation with
specified monomials and integer coefficients.

In some proof problems, symbolic computation methods face difficulties in
controlling the space complexity during intermediate calculations. In such
cases, it is possible to instantiate some variables and use symbolic-numeric
hybrid computations to obtain a large number of simplified solutions to the
problem. Then, utilizing implicit function interpolation techniques such as
algebraic curves or surfaces, the solutions of the instantiated problem can be
integrated into solutions of the original problem. Often, through theoretical
analysis, we can anticipate that the integrated solutions form a polynomial or
surface with integer coefficients.

In this context, we not only aim to obtain sparse solutions with fewer terms,
but also desire the coefficients of the resulting polynomials to be reasonably
small. Additionally, we would like to determine the maximum range of
coefficient values for which the solutions remain unique as integer-
coefficient polynomials.

In this section, we discuss the following problem of implicit polynomial
interval interpolation: constructing a polynomial curve that passes through a
given set of vertical short lines on the plane. We focus on studying the
uniqueness of integer-coefficient polynomial interpolation.

We have the following theorem.

######  Theorem 4.1.

Let $m_{1}(x,y)=x^{k_{1}}y^{l_{1}},m_{2}(x,y)=x^{k_{2}}y^{l_{2}},$ and
$m_{q}(x,y)=x^{k_{q}}y^{l_{q}}$ be monomials of $x,y$ and
$f(x,y)=c_{1}m_{1}(x,y)+c_{2}m_{2}(x,y)+\cdots+c_{q}m_{q}(x,y)$ a polynomial
that satisfies $|c_{j}|\leq 1,j=1,2,\cdots,q$ and passes through the line
segments

$S_{i}=\\{(x,y)|x=x_{i}>0,0<y_{i}-\delta\leq y\leq
y_{i}+\delta\\},i=1,2,\cdots,N.$

Then

$\sum_{i=1}^{N}(f(x_{i},y_{i}))^{2}\leq\delta^{2}\,q\cdot\sum_{i=1}^{N}h(x_{i},y_{i}+\delta),$

where

$h(x,y)=\sum_{j=1}^{q}\left(\frac{\partial m_{j}(x,y)}{\partial
y}\right)^{2}.$

######  Proof 4.2.

Assume that $f(x,y)$ passes the following points on each of the given line
segments at point $P_{i}:(x_{i},y_{i}+\Delta y_{i})$, where $|\Delta
y_{i}|\leq\delta$ for $i=1,2,\cdots,N$. Then, for $i=1,2,\cdots,N$ we have

$\displaystyle(f(x_{i},y_{i}))^{2}$ $\displaystyle=$
$\displaystyle(f(x_{i},y_{i})-f(x_{i},y_{i}+\Delta y_{i})^{2}$
$\displaystyle=$
$\displaystyle{\displaystyle\sum_{j=1}^{q}\left(c_{j}(x_{i}^{k_{j}}\,y_{i}^{l_{j}}-x_{i}^{k_{j}}\,(y_{i}+\Delta
y_{i})^{l_{j}})\right)^{2}}$ $\displaystyle\leq$
$\displaystyle{\displaystyle(c_{1}^{2}+c_{2}^{2}+\cdots+c_{q}^{2})\cdot\sum_{j=1}^{q}x_{i}^{2k_{j}}\left(y_{i}^{l_{j}}-(y_{i}+\Delta
y_{i})^{l_{j}}\right)^{2}}$ $\displaystyle=$ $\displaystyle{\displaystyle
q\cdot\sum_{j=1}^{q}x_{i}^{2k_{j}}\left(\Delta y_{i}\cdot
l_{j}\cdot(y_{i}+\theta_{l_{j}}\Delta
y_{i})^{l_{j}-1}\right)^{2}\,\quad(0<\theta_{l_{j}}<1)}$ $\displaystyle\leq$
$\displaystyle{\displaystyle
q\cdot\sum_{j=1}^{q}x_{i}^{2k_{j}}((l_{j})\delta\cdot(y_{i}+\delta)^{l_{j}-1})^{2}}$
$\displaystyle=$
$\displaystyle{\displaystyle\delta^{2}\,q\cdot\sum_{j=1}^{q}l_{j}^{2}\,x_{i}^{2k_{j}}\left(y_{i}+\delta\right)^{2l_{j}-2}=\delta^{2}\,q\cdot\sum_{j=1}^{q}\left(\frac{\partial
m_{j}}{\partial y}(x_{i},y_{i}+\delta)\right)^{2}},$

######  Corollary 4.3.

Let $p=(d+1)(d+2)/2$. If a polynomial

$f(x,y)=c_{0}+c_{1}x^{2}+c_{2}y^{2}+\cdots+c_{p-d}x^{2d}+c_{p-d+1}x^{2d-2}y^{2}+\cdots+c_{p-1}y^{2d},$

with $|c_{j}|\leq 1$ for $j=0,1,\cdots,p$ passes through the line segments

$S_{i}=\\{(x,y)|x=x_{i}>0,0<y_{i}-\delta\leq y\leq
y_{i}+\delta\\},i=1,2,\cdots,N.$

Then

$\sum_{i=1}^{N}(f(x_{i},y_{i}))^{2}\leq
2(d+1)(d+2)\delta^{2}\cdot\sum_{i=1}^{N}h_{d}(x_{i},y_{i}+\delta),$

where

$h_{d}(x,y)=y^{2}\sum_{l,j\geq 0,l+j\leq d-1}(l+1)^{2}x^{4j}y^{4l}.$

In particular,

$\displaystyle
h_{2}(x,y)={y}^{2}\left(1+{x}^{4}+4\,{y}^{4}\right),~{}~{}h_{3}(x,y)={y}^{2}\left(1+{x}^{4}+4\,{y}^{4}+{x}^{8}+4\,{x}^{4}{y}^{4}+9\,{y}^{8}\right),$
$\displaystyle
h_{4}(x,y)={y}^{2}\left(1+{x}^{4}+4\,{y}^{4}+{x}^{8}+4\,{x}^{4}{y}^{4}+9\,{y}^{8}+{x}^{12}+4\,{x}^{8}{y}^{4}+9\,{x}^{4}{y}^{8}+16\,{y}^{12}\right).$

We investigate the problem of rational solutions for the following
optimization problem: given points $P_{i}=(x_{i},y_{i})$, $i=1,2,\cdots,N$, a
small positive number $\delta>0$, and a set of monomials
$m_{j}=x^{k_{j}}y^{l_{j}}$, $j=1,2,\cdots,q$, the objective is to find
rational numbers
$(c_{1},c_{2},\cdots,c_{q})\in[-1,1]^{q}\subset\mathbb{Q}^{q}$ that satisfy
the following conditions:

  * •

The polynomial curve $f(x,y)=c_{1}m_{1}+c_{2}m_{2}+\cdots+c_{q}m_{q}=0$ passes
through the line segment $[x_{i},x_{i}]\times[y_{i}-\delta,y_{i}+\delta]$, for
$i=1,2,\ldots,N$.

  * •

$\max{|c_{j}|,j=1,2,\cdots,q}=1$.

  * •

The sum of squares of $f(x_{i},y_{i})$, for $i=1,2,\ldots,N$, is minimized.

It is clear that if $(c_{1},c_{2},\cdots,c_{q})$ is a local minimum, then
$(-c_{1},-c_{2},\cdots,-c_{q})$ is also a local minimum, and
$\min\\{|c_{j}|,j=1,2,\cdots,q\\}=1$. This optimization problem can be
decomposed into $q$ individual optimization problems $(Q_{j})$, for
$j=1,2,\cdots,q$.

$(Q_{j}):\left\\{\begin{array}[]{rcl}\min&&{\displaystyle\sum_{i=1}^{N}\left(c_{1}x_{i}^{k_{1}}y_{i}^{l_{1}}+c_{2}x_{i}^{k_{2}}y_{i}^{l_{2}}+\cdots+c_{q}x_{i}^{k_{q}}y_{i}^{l_{q}}\right)^{2},}\\\
\mbox{s.t.}&&(c_{1}x_{i}^{k_{1}}\underline{y_{i}}^{l_{1}}+c_{2}x_{i}^{k_{2}}\underline{y_{i}}^{l_{2}}+\cdots+c_{q}x_{i}^{k_{q}}\underline{y_{i}}^{l_{q}})\\\
&&\phantom{x}\times(c_{1}x_{i}^{k_{1}}\overline{y_{i}}^{l_{1}}+c_{2}x_{i}^{k_{2}}\overline{y_{i}}^{l_{2}}+\cdots+c_{q}x_{i}^{k_{q}}\overline{y_{i}}^{l_{q}})<0,\\\
&&\phantom{here}(\underline{y_{i}}=y_{i}-\delta,\;\overline{y_{i}}=y_{i}+\delta),\\\
&&-1<c_{1},c_{2},\cdots,c_{q}<+1,\quad c_{j}=1.\end{array}\right.$ (4.1)

Each sub-problem $(Q_{j})$ can be transformed into a linear inequality problem
with constraints $c_{1},c_{2},\ldots,c_{q}$, where the objective function is a
quadratic function. Therefore, it can be converted into a problem of solving a
system of linear inequalities

If the coordinates of points $P_{1},P_{2},\ldots,P_{N}$ are all rational
numbers, it is possible to obtain a rational solution to the original problem
using symbolic computation methods. Subsequently, asymptotic analysis can be
applied to find the ”optimal” integer coefficient solution to the original
problem.

######  Example 4.4.

Given any real number $b\,(0<b\leq 1)$, let
$E(b):=\\{(x,y)|x^{2}+y^{2}/b^{2}=1\\}$ be an ellipse with $b$ as minor-axis
and the major-axis equal to 1, and $L(b)$ be the maximal perimeter of
triangles inscribed in the ellipse. The aim is to find an explicit formula
that connects $L$ and $b$.

For a given value of $b$, the maximum inscribed triangle can be approximated
using numerical methods to obtain an approximate value of $L(b)$. Let
$\underline{L}(b)$ and $\overline{L}(b)$ be the lower and upper bounds of
$L(b)$, respectively. By employing numerical methods, we can determine the
upper and lower bounds to make the interval
$[\underline{L}(b),\overline{L}(b)]$ arbitrarily small.

For rational numbers $b_{i}=i/16,i=1,2,\cdots,16$, we calculate the upper and
lower bound intervals such that $\overline{L}(b_{i})-\underline{L}(b_{i})\leq
1/10^{14}$. This results in a dataset containing $16$ very short vertical line
segments (see Figure 1). The dataset
$S_{i}:[b_{i},b_{i}]\times[\underline{L_{i}},\overline{L_{i}}]$ is as follows:

$\displaystyle S_{1}:[1/16,4.00391006468634],$ $\displaystyle
S_{2}:[1/8,4.01568603234052],$ $\displaystyle S_{3}:[3/16,4.03546516762171],$
$\displaystyle S_{4}:[1/4,4.06347582517416],$ $\displaystyle
S_{5}:[5/16,4.10003596888417],$ $\displaystyle S_{6}:[3/8,4.14554927249083],$
$\displaystyle S_{7}:[7/16,4.20049672260437],$ $\displaystyle
S_{8}:[1/2,4.26542082282451],$ $\displaystyle S_{9}:[9/16,4.34089883607274],$
$\displaystyle S_{10}:[5/8,4.42750157636573],$ $\displaystyle
S_{11}:[11/16,4.52573604652852],$ $\displaystyle
S_{12}:[3/4,4.63597477802543],$ $\displaystyle
S_{13}:[13/16,4.75838211258069],$ $\displaystyle
S_{14}:[7/8,4.89285523313257],$ $\displaystyle
S_{15}:[15/16,5.03899987143207],$ $\displaystyle S_{16}:[1,5.19615242270664].$

Here, for simplicity, we use $[b_{i},L_{i}]$ to represent the line segment
$[b_{i},L_{i}-\delta]\times[b_{i},L_{i}-\delta]$, where $\delta=0.5\times
10^{-14}$. For example, $S_{1}:=[1/16,4.00391006468634]$ represents the
interval

$[1/16,4.00391006468634-0.5\times
10^{-14}]\times[1/16,4.00391006468634+0.5\times 10^{-14}].$

The current objective is to reconstruct $f(L,b)=0$ using the $16$ vertical
short lines and implicit function interpolation. From the references [3, 27],
it is known that the solution to this problem is given by:

$f(b,L)=\left(1-{b^{2}}\right)^{2}{L^{4}}-8\,\left(2-{b^{2}}\right)\left(1-2\,{b^{2}}\right)\left(1+{b^{2}}\right){L^{2}}-432\,{b^{4}}=0.$
(4.2)

$b$$L$0.514.55.0 Figure 3: The 16 vertical segments determining a minimal
polynomial $f(L,b)=0$. The segment’s length in this picture is enlarged
$10^{9}$ times for a better visualization effect.

We will reconstruct this relationship through implicit function interval
interpolation. The challenge lies in finding the most likely integer
coefficient polynomial from the floating-point results.

Before interpolation, we first analyze the function between the perimeter $L$
of the maximum inscribed triangle and the semi-axis $b$ of the ellipse. This
function can be represented by a polynomial equation $f(L,b)=0$.

To do so, let us assume that $ABC$ is the triangle inscribed in the ellipse
$E(1,b):x^{2}+y^{2}/b^{2}=1$ with the maximum perimeter $L$. The coordinates
of the vertices $A=(x_{1},y_{1})$, $B=(x_{2},y_{2})$, and $C=(x_{3},y_{3})$ of
this maximum inscribed triangle are known. The side lengths of the triangle
are denoted as $BC=u$, $CA=v$, and $AB=w$. We can establish the following
equations:

$\displaystyle f_{0}:=L-(u+v+w)=0,$ $\displaystyle
f_{4}:=x_{1}^{2}+y_{1}^{2}/b^{2}-1=0,$ $\displaystyle
f_{1}:=(x_{2}-x_{3})^{2}+(y_{2}-y_{3})^{2}-u^{2}=0,$ $\displaystyle
f_{5}:=x_{2}^{2}+y_{2}^{2}/b^{2}-1=0,$ $\displaystyle
f_{2}:=(x_{3}-x_{1})^{2}+(y_{3}-y_{1})^{2}-v^{2}=0,$ $\displaystyle
f_{6}:=x_{3}^{2}+y_{3}^{2}/b^{2}-1=0.$ $\displaystyle
f_{3}:=(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}-w^{2}=0,$

Then applying the Lagrange multiplier method we have

$\mathcal{L}=u+v+w-\sum_{j=1}^{6}k_{j}f_{j},$

and

$\displaystyle{\displaystyle f_{7}:=\frac{\partial\mathcal{L}}{\partial
u}=0,\quad f_{8}:=\frac{\partial\mathcal{L}}{\partial v}=0,\quad
f_{9}:=\frac{\partial\mathcal{L}}{\partial w}=0,}$ $\displaystyle
f_{10}:={\displaystyle\frac{\partial\mathcal{L}}{\partial x_{1}}=0},\quad
f_{11}:={\displaystyle\frac{\partial\mathcal{L}}{\partial x_{2}}=0},\quad
f_{12}:={\displaystyle\frac{\partial\mathcal{L}}{\partial x_{3}}=0},$ (4.4)
$\displaystyle f_{13}:={\displaystyle\frac{\partial\mathcal{L}}{\partial
y_{1}}=0},\quad f_{14}:={\displaystyle\frac{\partial\mathcal{L}}{\partial
y_{2}}=0},\quad f_{15}:={\displaystyle\frac{\partial\mathcal{L}}{\partial
y_{3}}=0},$

The system of equations (4.4) and (4.4) consists of $16$ polynomial equations
$f_{i}=0$ $(i=0,1,\ldots,15)$ and $17$ unknowns, namely:

$L,\;b,\;u,\;v,\;w,\;x_{1},\;y_{1},\;x_{2},\;y_{2};\;x_{3},\;y_{3},\;k_{1},\;\cdots,\;k_{6}.$

In theory, it is possible to eliminate $15$ variables from the system
${f_{i}=0,1\leq i\leq q}$. For example, by eliminating all variables except
$L$ and $b$, we obtain a polynomial in $L$ and $b$ denoted as $\mbox{\it
Res/}(L,b)=0$. Therefore, equation (4.2) must be a factor of $\mbox{\it
Res/}(L,b)$.

Previous work has shown that direct symbolic computation can eliminate $12$
variables from ${f_{0},f_{1},\ldots,f_{15}=0}$, specifically
${x_{1},\ldots,y_{6},k_{1},\ldots,k_{6}}$, resulting in a system of $4$
equations only involving $L$, $b$, $u$, $v$, and $w$. In this reduced system,
both $L$ and $b$ appear only with even powers. To save space, we denote these
$4$ intermediate results as:

$g_{1}(L^{2},b^{2},u,v,w)=0,\quad g_{2}(L^{2},b^{2},u,v,w)=0,\quad
g_{3}(L^{2},b^{2},u,v,w)=0,\quad g_{4}(L^{2},b^{2},u,v,w)=0.$ (4.5)

Further elimination of the variables $u$, $v$, and $w$ encounters storage
space issues. To obtain the final result (4.2), a combination of symbolic
computation and sparse interpolation is required. In this paper, we attempt to
use implicit function interval interpolation to construct a polynomial curve
$F(x,y)=0$ passing through the $16$ small line segments in Figure 3. We then
prove that substituting $x=L$ and $y=b$ into $F(x,y)$ yields the solution to
the original problem.

According to equation (4.5), we are searching for curves of the following form
that pass through the given 16 short lines:

$F(x,y)=\,c_{0}+c_{1}x^{2}+c_{2}y^{2}+c_{3}x^{4}+c_{4}x^{2}y^{2}+c_{5}y^{4}+\cdots+c_{p-d-1}x^{2d}+\cdots+c_{p-1}y^{2d},$
(4.6)

Although the upper bound of $d$ can be estimated based on the degrees of
$L,b,u,v,w$ in equation (4.5) and the form of the Sylvester resultant, the
obtained value of $d$ is generally larger than the actual final result.
Additionally, if the interpolation target polynomial has $d>5$, $F(x,y)$
contains 21 monomials:

$1;\quad x^{2},\,y^{2};\quad x^{4},\,x^{2}y^{2},\,y^{4};\quad
x^{6},\,x^{4}y^{2},\,x^{2}y^{4},\,y^{6};$
$x^{8},\,x^{6}y^{2},\,x^{4}y^{4},\,x^{2}y^{6},\,y^{8};\quad
x^{10},\,x^{8}y^{2},\,x^{6}y^{4},\,x^{4}y^{6},\,x^{2}y^{8},\,y^{10},$

which exceeds the current number of data points ($16$). Therefore, we will try
$d=2,3,4$ in ascending order. If $d=2$ is successful, we will then test the
case of $d=1$. If $d=4$ is not successful, we will use numerical optimization
to supplement more data and test the cases of $d=5,\ldots$. In the actual
computations, we found that $d=2,3$ do not have polynomials of the form (4.6)
that simultaneously pass through the 16 line segments in Figure 3, but there
is such a polynomial curve for $d=4$. The specific calculation process is
provided below. According to Theorem 4.1 and Corollary 4.3, we have the
following result:

  1. 1.

If a polynomial
$f(x,y)=c_{0}+c_{1}x^{2}+c_{2}y^{2}+c_{3}x^{4}+c_{4}x^{2}y^{2}+c_{5}y^{4}$
passes the line segments $S_{1},S_{2},\cdots,S_{16}$ and
$c_{0},c_{1},\cdots,c_{5}\in[-1,1]$, then

$\sum_{i=1}^{16}f(x_{i},y_{i})^{2}\leq r_{2}:=1.2740\times 10^{-21};$ (4.7)

  2. 2.

If
$f(x,y)=c_{0}+c_{1}x^{2}+c_{2}y^{2}+c_{3}x^{4}+c_{4}x^{2}y^{2}+c_{5}^{4}+c_{5}y^{4}+c_{6}x^{6}+c_{7}x^{4}y^{2}+c_{8}x^{2}y^{4}+c_{9}y^{6}$
passes the line segments $S_{1},S_{2},\cdots,S_{16}$ and
$c_{0},c_{1},\cdots,c_{9}\in[-1,1]$ , then

$\sum_{i=1}^{16}f(x_{i},y_{i})^{2}\leq r_{3}:=2.2799\times 10^{-18};$ (4.8)

  3. 3.

If
$f(x,y)=c_{0}+c_{1}x^{2}+c_{2}y^{2}+\cdots+c_{10}x^{8}+c_{11}x^{6}y^{2}+\cdots+c_{14}y^{8}$
passes the line segments $S_{1},S_{2},\cdots,S_{16}$ and
$c_{0},c_{1},\cdots,c_{14}\in[-1,1]$, then

$\sum_{i=1}^{16}f(x_{i},y_{i})^{2}\leq r_{4}:=3.2326\times 10^{-15}.$ (4.9)

Let us denote the set of monomials as

$\displaystyle{\cal P}_{2}:=\\{1,x^{2},y^{2},x^{4},x^{2}y^{2},y^{4}\\},$
$\displaystyle{\cal
P}_{3}:=\\{1,x^{2},y^{2},x^{4},x^{2}y^{2},y^{4},x^{6},x^{4}y^{2},x^{2}y^{4},y^{6}\\},$
$\displaystyle{\cal
P}_{4}:=\\{1,x^{2},y^{2},x^{4},x^{2}y^{2},y^{4},x^{6},x^{4}y^{2},x^{2}y^{4},y^{6},x^{8},x^{6}y^{2},x^{4}y^{4},x^{2}y^{6},y^{8}\\},$

and $P_{i}=(x_{i},y_{i})\,(i=1,2,\cdots,16)$ the barycenters of the segments
$S_{i}$ as presented in Table 4.

Table 4:

The barycenters of the 16 segments

$\begin{array}[]{|c|c|c||c|c|c|}\hline\cr\phantom{PP}P_{1}\phantom{PP}&\phantom{x}1/16\phantom{x}&\phantom{x}4.003910064687\phantom{x}&\phantom{PP}P_{9}\phantom{PP}&\phantom{x}9/16\phantom{x}&\phantom{x}4.340898836073\phantom{x}\\\
\hline\cr P_{2}&1/8&4.015686032341&P_{10}&5/8&4.427501576366\\\ \hline\cr
P_{3}&3/16&4.035465167622&P_{11}&11/16&4.525736046529\\\ \hline\cr
P_{4}&1/4&4.063475825175&P_{12}&3/4&4.635974778026\\\ \hline\cr
P_{5}&5/16&4.100035968885&P_{13}&13/16&4.758382112581\\\ \hline\cr
P_{6}&3/8&4.145549272491&P_{14}&7/8&4.892855233133\\\ \hline\cr
P_{7}&7/16&4.200496722605&P_{15}&15/16&5.038999871433\\\ \hline\cr
P_{8}&1/2&4.265420822825&P_{16}&{\small 1}&5.196152422707\\\
\hline\cr\end{array}$

For $d=2,3,4$ and $2\leq q\leq(d+1)(d+2)/2$, we define

$r(d,q):=\min_{|c_{j}|\leq 1,j=1,\cdots,q\atop m_{1},\cdots,m_{q}\in{\cal
P}_{d}}{\sum_{i=1}^{16}\left(\sum_{j=1}^{q}c_{j}m_{j}(x_{i},y_{i})\right)^{2}}.$

Then

$q_{1}<q_{2}\Longrightarrow r(d,q_{1})\leq r(d,q_{2}).$

Let $r_{2},r_{3},r_{4}$ be the real numbers defined in (4.7), (4.8), and
(4.9). It can be inferred from

$r(d,q)>\frac{q}{(d+1)(d+2)/2}\cdot r_{d}$

that there does not exist $m_{1}(x,y),\cdots,m_{q}(x,y)\in{\cal P}d$ and real
numbers $c_{1},\cdots,c_{q}$ such that
$f(x,y)=c_{1}m_{1}(x,y)+\cdots+c_{q}m_{q}(x,y)(\not\equiv 0)$ simultaneously
passes through $S_{1},S_{2},\ldots,S{16}$. By using symbolic computation, the
following results can be obtained

$\begin{array}[]{lclcl}r(2,2)=0.4695\cdots,&&r(3,2)=0.2023,\cdots&&r(4,2)=0.1263\cdots,\\\
r(2,3)=0.001129\cdots,&&r(3,3)=0.001129,\cdots&&r(4,3)=6.4239\times
10^{-4},\\\ r(2,4)=1.0123\times 10^{-6},&&r(3,4)=1.2093\times
10^{-7},&&r(4,4)=1.2093\times 10^{-7},\\\ r(2,5)=3.1611\times
10^{-8},&&r(3,5)=1.4273\times 10^{-8},&&r(4,5)=2.0731\times 10^{-9},\\\
r(2,6)=3.9317\times 10^{-8},&&r(3,6)=4.2241\times
10^{-11},&&r(4,6)=6.9501\times 10^{-12},\\\
\cline{1-1}\cr&&r(3,7)=1.5547\times 10^{-13}&&r(4,7)=1.5547\times 10^{-13},\\\
\cline{5-5}\cr&&r(3,8)=6.6283\times 10^{-15},&&r(4,8)=4.939\times 10^{-30};\\\
&&r(3,9)=5.0264\times 10^{-16},&&r(4,9)\leq 4.939\times 10^{-30},\\\
&&r(3,10)=1.2133\times 10^{-16},&&r(4,10)\leq 4.939\times 10^{-30},\\\
\cline{3-3}\cr&&&&r(4,11)\leq 4.939\times 10^{-30},\\\ &&&&r(4,12)\leq
4.939\times 10^{-30},\\\ &&&&r(4,13)\leq 4.939\times 10^{-30},\\\
&&&&r(4,14)\leq 4.939\times 10^{-30},\\\ &&&&r(4,15)\leq 4.939\times
10^{-30}.\end{array}$

The above results show that when $d=2,3$ and $d=4,q=8$, we have
$r(d,q)>r_{d}$. This indicates that in order for the curve $F(x,y)=0$ in the
form of equation (4.6) to pass through the 16 line segments, the degree of $F$
must be at least 8, and it must have at least 8 monomials.

Furthermore, we have:

$\min\sum_{i=1}^{16}\sum_{j=1}^{8}\left(c_{j}m_{j}(x_{i},y_{i})\right)^{2}=\frac{\overbrace{249227\cdots
358763}^{492\text{digits}}}{\underbrace{754846\cdots\cdots
047616}_{517\text{digits}}}\approx 3.301701\cdots\times 10^{-26},$

and the coefficients of the monomials
${x}^{4},{y}^{2},{y}^{4},{x}^{2}{y}^{2},{x}^{2}{y}^{4},{x}^{4}{y}^{2},{x}^{4}{y}^{4},{x}^{6}{y}^{2},$
can be expressed using continued fractions as follows:

$c_{1}=1,\phantom{alargespacehere}c_{2}=\cfrac{1}{26+\cfrac{1}{1+\cfrac{1}{3329579+\cdots}}},$

$c_{3}=-\cfrac{1}{431+\cfrac{1}{1+\cfrac{1}{208102+\cdots}}},\phantom{space}c_{4}=-\cfrac{1}{18+\cfrac{1}{14773286+\cfrac{1}{36+\cdots}}},$

$c_{5}=\cfrac{1}{216+\cfrac{1}{121148115+\cfrac{1}{10+\cdots}}},\phantom{space}c_{6}=-\cfrac{1}{18+\cfrac{1}{14773286+\cfrac{1}{1+\cdots}}},$

$c_{7}=-\cfrac{1}{431+\cfrac{1}{1+\cfrac{1}{937254+\cdots}}},\phantom{space}c_{8}=\cfrac{1}{26+\cfrac{1}{1+\cfrac{1}{21212373+\cdots}}}.$

The corresponding polynomial curve is given by:

$\displaystyle
f(x,y)=431.9999989\cdots{x}^{4}+16.00000013\cdots{y}^{2}-1.000000008\cdots{y}^{4}-23.99999985\cdots{x}^{2}{y}^{2}$
$\displaystyle+1.999999994\cdots{x}^{2}{y}^{4}-23.99999986\cdots{x}^{4}{y}^{2}-{x}^{4}{y}^{4}+15.99999998\cdots{x}^{6}{y}^{2},$
(4.10)

In this polynomial, the coefficients are multiplied by a certain factor so
that ${\displaystyle\min_{1\leq j\leq q}|c_{j}|=1}$, meaning the coefficient
with the smallest absolute value is set to 1. Naturally, we conjecture that
the corresponding integer coefficient polynomial for this polynomial is given
by:

$F(x,y)=432\,{x}^{4}+16\,{y}^{2}-\,{y}^{4}-24\,{x}^{2}{y}^{2}+2\,{x}^{2}{y}^{4}-24\,{x}^{4}{y}^{2}-{x}^{4}{y}^{4}+16\,{x}^{6}{y}^{2}.$
(4.11)

There is a high possibility that this polynomial can also be expressed as a
linear combination of the line segments $S_{1},S_{2},\ldots,S_{16}$. To
confirm this conjecture, we only need to substitute the values
$x_{i}=i/16,i=1,2,\ldots,16$ into (4.11) and verify, using numerical methods,
that the signs of the resulting univariate polynomial equation $f(x_{i},y)=0$
are opposite at the two endpoints of $S_{i}$. Computational calculations have
indeed confirmed this.

Finally, based on the continued fraction representations of
$c_{1},c_{2},\cdots,c_{8}$ (page 4.4), we find a number $R$ such that, under
the assumption that $c_{1},c_{2},\cdots,c_{8}\in\mathbb{Z}$ and $|z_{i}|\leq
R$, equation (4.11) is the unique eighth-degree curve that passes through
$S_{1},S_{2},\cdots,S_{16}$. To do this, assume that
$c_{1},c_{2},\cdots,c_{8}\in\mathbb{Z}$ and $\max{|c_{j}|,j=1,2,\cdots,8}=R$.
Then we have:

$\frac{|c_{2}|}{R}\in(\frac{1}{27},\frac{1}{27}+\frac{1}{2427263793}),$

which implies,

$\frac{|c_{2}|}{R}-\frac{1}{27}\leq\frac{1}{2427263793},$

and further,

$\frac{1}{27R}\leq\frac{27|c_{2}|-R}{27R}\leq\frac{1}{2427263793},~{}~{}\mbox{\it
if\,}~{}|c_{2}|/R\not=1/27.$

This means

$R\geq 89898659.$

After applying the same procedure to several other continued fractions, we
have

$\displaystyle
R\geq\max\left\\{89900495,531838299/2,130839964313/5,\right.\phantom{x}$
$\displaystyle\phantom{xxx}\left.265919167,404894159,572734097\right\\}$
$\displaystyle\,>1.3083\times 10^{11},$

This indicates that under the condition $c_{j}<1.3088\times
10^{11},(j=1,2,\cdots,8)$, the polynomial (4.11) is the unique algebraic curve
with integer coefficients of degree not exceeding $8$ that simultaneously
passes through $S_{1},S_{2},\ldots,S_{16}$.

Note that the solutions we obtain using interval implicit function
interpolation may only be a subset of the solutions to the original system of
equations. Adding these solutions to the original system
${f_{0},f_{1},\ldots,f_{15}}$ can help us find additional solutions to the
original system or prove that $F(L,b)=0$ is the unique solution to the
original system.

##  5 Conclusion

This article discusses the problem of interval interpolation for planar
algebraic curves. The input consists of several rectangular regions and
vertical short lines on the plane. We discuss the solution methods for the
following three problems and describe the solving process using several
meaningful examples: (1) By simultaneously determining the lowest degree of
the curves in the given small neighborhoods and constructing algorithms for
curve construction, we transform the problem of determining the lowest degree
into a positive-definiteness determination problem on a rectangular region in
a high-dimensional space of multivariate polynomials. (2) By simultaneously
constructing sparse solutions for algebraic curves in the given small
neighborhoods, we study the analysis of observational data for the orbits of
38 celestial bodies in the solar system, including planets, asteroids, and
dwarf planets. From this analysis, we rediscover Kepler’s Third Law. (3) By
simultaneously constructing algorithms for integral coefficient algebraic
curves in the given small neighborhoods, we explain the uniqueness analysis
method for integral coefficient polynomial curves using an interesting set
problem as an example.

The methods constructed in this article, as symbolic-numeric hybrid methods,
help overcome the issue of intermediate process inflation that may occur in
pure symbolic computation methods. They also assist in finding possible
invariants and useful empirical formulas in data with errors. These methods
are suitable for analyzing datasets obtained from practical observations with
large sample sizes but limited accuracy (e.g., Example 2.14 in this article),
as well as handling datasets with high precision and small sample sizes
generated for theoretical research problems (e.g., Example LABEL:ex3) or data
consisting of accurate but non-rational numbers proposed in pure theoretical
research (e.g., Example 2.7). The methods presented in this article can also
be applied to algebraic surface interpolation problems in small neighborhoods
in high-dimensional spaces.

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 2171159 and 12071282).

## References

  *   * [1] Anjyo K, Lewis J P, Pighin F. Scattered data interpolation for computer graphics. In: ACM SIGGRAPH 2014 Courses, New York: Association for Computing Machinery, 2014, 1-69.
  * [2] Bagheri S, Konen W, Emmerich M, et al. Self-adjusting parameter control for surrogate-assisted constrained optimization under limited budgets. Appl. Soft Comput, 2017, 61: 377-393.
  * [3] Chen L, Leng T, Shen L, et al. Finding the symbolic solution of a geometric problem through numerical computations. In: Computer Mathematics, Springer, 2014, 223-241.
  * [4] Datta B. Numerical Methods for Linear Control Systems, volume 1. New York: Academic Press, 2004.
  * [5] De Boor C, Ron A. On multivariate polynomial interpolation. Constr. Approx,1990, 6(3): 287-302.
  * [6] Dell’Accio F, Di Tommaso F, Nouisser O, et al. Solving Poisson equation with Dirichlet conditions through multinode shepard operators. Comput. Math. Appl, 2021, 98: 254-260.
  * [7] Dell’Accio F, Di Tommaso F, Siar N. On the numerical computation of bivariate Lagrange polynomials. Applied Mathematics Letters, 2021, 112:106845(9pp).
  * [8] Dell’Accio F, Di Tommaso F, Siar N, et al. Numerical differentiation on scattered data through multivariate polynomial interpolation. BIT Numer. Math, 2022, 62(3): 773-801.
  * [9] Dini P, Begni A, Ciavarella S, et al. Design and testing novel one-class classifier based on polynomial interpolation with application to networking security. IEEE Access, 2022, 10: 67910-67924.
  * [10] Drake K P, Fuselier E J, Wright G B. Implicit surface reconstruction with a curl-free radial basis function partition of unity method. SIAM J. Sci. Comput, 2022, 44(5): A3018-A3040.
  * [11] Errachid M, Essanhaji A, Messaoudi A. Rmvpia: a new algorithm for computing the Lagrange multivariate polynomial interpolation. Numer. Algorithms, 2020, 84(4): 1507-1534.
  * [12] Essanhaji A, Errachid M. Lagrange multivariate polynomial interpolation: A random algorithmic approach. J. Appl. Math, 2022, 2022: 8227086(8pp).
  * [13] Fausett L V. Numerical Methods: Algorithms and Applications. United States: Prentice Hall, 2003.
  * [14] Feng R, Peng S. Quasi-interpolation scheme for arbitrary dimensional scattered data approximation based on natural neighbors and rbf interpolation. J. Comput. Appl. Math, 2018, 329: 95-105.
  * [15] Forsythe G E. Malcolm M A, Moler C B. Computer methods for mathematical computations. Prentice-Hall series in automatic computation, 1977, 259.
  * [16] Gasca M, Sauer T. On the history of multivariate polynomial interpolation. In: Numerical Analysis: Historical Developments in the 20th Century, Elsevier, 2001, 135-147.
  * [17] Gasca M, Sauer T. Polynomial interpolation in several variables. Adv Comput Math, 2000, 12(4): 377-410.
  * [18] Giménez N. An algorithm for implicit interpolation. Appl. Algebra Eng. Commun. Comput, 2014, 25(1): 119-157.
  * [19] Lewis J P, Pighin F, Anjyo K. Scattered data interpolation and approximation for computer graphics. In: ACM SIGGRAPH ASIA 2010 Courses, New York: Association for Computing Machinery, 2010, 1-73.
  * [20] Meletiou G C, Tasoulis D K, Vrahatis M N. Cryptography through interpolation approximation and computational intelligence methods. Bull. Greek Math. Soc, 2003, 48: 61-75.
  * [21] Morse B S, Yoo T S, Rheingans P, et al. Interpolating implicit surfaces from scattered surface data using compactly supported radial basis functions. In: ACM SIGGRAPH 2005 Courses, New York: Association for Computing Machinery, 2005, 78–es.
  * [22] Neidinger R D. Multivariable interpolating polynomials in newton forms. In: Joint Mathematics Meetings, Washington, DC, 2009, 5-8.
  * [23] Neidinger R D. Multivariate polynomial interpolation in newton forms. SIAM Review, 2019, 61(2): 361-381.
  * [24] Nuñez L, Regis R G, Varela K. Accelerated random search for constrained global optimization assisted by radial basis function surrogates. J. Comput. Appl. Math, 2018, 340: 276-295.
  * [25] Shen C, O’Brien J F, Shewchuk J R. Interpolating and approximating implicit surfaces from polygon soup. In: ACM SIGGRAPH 2004 Papers, 2004, 896-904.
  * [26] Smolik M, Skala V. Large scattered data interpolation with radial basis functions and space subdivision. Integr. Comput.-Aided Eng, 2018, 25(1): 49-62.
  * [27] Tang M, Yang Z, Zeng Z. Resultant elimination via implicit equation interpolation. J. Syst. Sci. Complex, 2016, 29(5): 1411-1435.
  * [28] Thévenaz P, Blu T, Unser M. Image interpolation and resampling. Handbook of medical imaging, processing and analysis, 2000, 1(1): 393-420.
  * [29] Varsamis D N, Karampetakis N P. On the newton bivariate polynomial interpolation with applications. Multidimens. Syst. Signal Process, 2014, 25(1): 179-209.
  * [30] Varsamis D N, Karampetakis N P, Mastorocostas P. An optimal bivariate polynomial interpolation basis for the application of the evaluation-interpolation technique. Appl. Math. Inf. Sci, 2014, 8(1): 117-125.
  * [31] Walker R J. Algebraic Curves. Princeton: Princeton University Press, 1950.
  * [32] Zhang S, Ren H. Hyperfinite interpolation, wu’s method and blending of implicit algebraic surfaces. J. Comput. Sci. Technol, 1999, 14(5): 518-529.