mylargesymbols”08 mylargesymbols”09 mylargesymbols”00 mylargesymbols”01
mylargesymbols”02 mylargesymbols”03

# Analytic Implicit Functions

Kyung Soo Rim This work was supported by a grant from the National Research
Foundation of Korea, Grant No. NRF-2017R1E1A1A03070307.

###### Abstract

In this paper, we introduce a method of converting implicit equations to the
usual forms of functions locally without differentiability. For a system of
implicit equations which are equipped with continuous functions, if there are
unique analytic implicit functions, that satisfies the system in some
rectangle, then each analytic function is represented as a power series which
is the weak-star limit of partial sums in the space of essentially bounded
functions. We also provide numerical examples in order to demonstrate how the
theoretical results in this article can be applied in practice and to show the
effectiveness of the suggested approaches.

###### keywords:

implicit function, implicit function theorem, inverse function theorem,
analytic function, continuous function, power series, weak-star limit

##  1 Introduction

Setting up implicit equations and solving them has long been so important that
it has virtually come to describe what mathematical analysis and its
applications are all about. Many useful mathematical models have expressions
of implicit functions, even for problems of minimizing or maximizing functions
subject to constraints. From information of a single point solution, the
implicit function theorem allows for understanding the relation between
variables but, in spite of that, the implicit function theorems do not provide
a usual function form defined explicitly. A central issue in this subject is
how to derive a power series expansion from an equation involving parameters
because the representation as a function form enables us to estimate the
mathematical models more easily, in economics, physics, engineering, etc., as
well as mathematical development.

In 1669, Issac Newton introduced one of the first instances of analyzing the
behavior of an implicitly defined function originated from an equation jointed
by mutually correlated variables ([9]). In 1684, Gottfried Leibniz applied
implicit differentiation to calculate partial derivatives, which is a way to
take the derivative of a term with respect to another variable without having
to isolate either variable ([11]). Joseph Lagrange, in 1770, derived an
inversion formula which is one of the fundamental formulas of combinatorics.
The formula is closely related to an inverse function theorem, by formulating
a formal power series expansion for an implicit function [4, Theorem 2.3.1].
In the 19th century, Augustin-Louis Cauchy was the first to state and solve an
implicit equation using a rigorous mathematical form. He gave proof of the
implicit function theorem having a form of power series by using Hardamard’s
estimate for the Taylor coefficient of a given function, which is induced
formally ([4, Theorems 2.4.6 and 6.1.2]). Since then, there have been many
improvements on the existence of usual function expressions from implicit
equations under suitable assumptions.

Most importantly, in 1877, Ulisse Dini was the first to prove the real
variable result of an implicit function theorem ([1]). Conceived over two
hundred years ago as a tool for studying mechanics and physics, the implicit
function theorem has many formulations and is used in many aspects of
mathematics. In addition, there are many particular types of implicit function
theorems which are extended to Banach spaces, even under degenerate or non-
smooth situations (e.g., [8, 7]). Implicit function theorems, even those that
are quite sophisticated, are fundamental and powerful parts of the foundation
of modern mathematics. Almost all studies of an implicit function are related
to its existence rather than showing how they behave.

This research is devoted to determining a form of function of an implicit
function which does not adopt the Taylor series of a function given, as used
in Cauchy’s method, but instead of differentiability, we need the
integrability of the given function. This article highlights the formulations
of a power series representation from an implicit equation such that its
partial sum converges weak-star in the space of essentially bounded functions.
This method relies essentially on the continuity of a given function. In
addition, we present several examples for the computational validity of this
study, along with a numerical calculation for each. While there have been too
many contributions to the theory of implicit functions to mention them all
here, we would like to refer to books [4] and [2] as good overviews.

Throughout the article we use multi-index notations. Let
$\alpha=(\alpha_{1},\ldots,\alpha_{n})$, $\beta=(\beta_{1},\ldots,\beta_{n})$
in $\mathbb{Z}^{n}$ and $c$ a scalar. We denote $\alpha\leq\beta$, $\alpha\leq
c$, and $\alpha\pm c$ if $\alpha_{k}\leq\beta_{k}$, $\alpha_{k}\leq c$, and
$(\alpha_{1}\pm c,\ldots,\alpha_{n}\pm c)$ for every $k$, respectively, where
the plus–minus sign is replaced by either the plus or minus sign in the same
order. Moreover, for a vector-valued function $f$, $f_{k}$ denotes the $k$th
component of $f$. For indexed quantities $d_{\alpha}$,
$\left(d_{\alpha}\right)_{\alpha}$ is called a tensor (or a matrix only when
$\alpha$ has two components) and the determinant of a matrix
$\left(d_{\alpha}\right)_{\alpha}$ is defined by
$|\left(d_{\alpha}\right)_{\alpha}|$ or
$\det\left(d_{\alpha}\right)_{\alpha}$. The subsets $R$ and $I$ of
$\mathbb{R}^{n}$ indicate rectangles which are the Cartesian products of $n$
compact intervals and $|R|$ defines the Lebesgue measure of $R$. Especially,
$I$ represents the codomain of an implicit function. Sometimes the integer $N$
is considered a multi-index, in which case all components are defined as $N$.
With these the main results are formulated in Theorems 3.2, 4.2, and 5.2.

Before understanding the proposed methods, we must first state the classical
implicit function theorem as follows.

######  Implicit function theorem.

Let $f(x,y):\mathbb{R}^{n+m}\to\mathbb{R}^{m}$ be continuously differentiable
with $f(a,b)=\mathbf{0}$. Suppose that
$|J_{f,y}(a,b)|=\det\big{(}\partial_{y_{j}}f_{i}(a,b)\big{)}\neq 0$. Then
there is an $R\times I$ of $(a,b)$ and a unique continuously differentiable
function $g:R\to I$ such that $f(x,g(x))=\mathbf{0}$ in $R$. Moreover, the
partial derivative of $g$ is given by
$\partial_{x_{j}}g(x)=-J_{f,y}(\partial_{y_{j}}f_{i}(x,g(x)))_{m\times
m}^{-1}(\partial_{x_{j}}f(x,g(x)))_{m\times 1}$.

This implicit function theorem has been extended to that of analytic
functions. If every $f_{k}$ is analytic in the implicit function theorem, then
every $g_{k}$ is also analytic (refer to [4, Theorem 6.1.2]).

##  2 Implicit functions

In this section, we derive the integral representation of an implicit function
to isolate a dependent variable even though it is neither transcendental nor
algebraic. Let $\operatorname{sgn}$ be the signum function on $\mathbb{R}$ and
$f(x,y)$ a function on $\mathbb{R}^{n+1}$. For each $x$, define a function
$\operatorname{sgn}_{y}f(x,y)$ by $1$ if $f(x,y)\geq 0$ and $-1$ otherwise. We
begin with a fundamental condition for the local existence of an implicit
function.

######  Assumption 1.

Let $f:\mathbb{R}^{n+1}\to\mathbb{R}$ be continuous. Then there is an $R\times
I\subset\mathbb{R}^{n+1}$ such that for each $x\in R$,
$\operatorname{sgn}_{y}f(x,y)$ has only one jump discontinuity on $I$.

From Assumption 1, by the intermediate value theorem, by the continuity of
$f$, and by the uniqueness of jump discontinuity, there is a continuous
function $g:R\to I$ such that $f(x,g(x))=0$ in $R$. Moreover, if $f$ is
continuously differentiable such that $\partial_{y}f(a,b)\neq 0$ with
$f(a,b)=0$, then by the implicit function theorem there is an $R\times
I\ni(a,b)$, which has a unique continuously differentiable function $y=g(x)$
such that $f(x,g(x))=0$ in $R$. Also, the Jacobian matrix is invertible in
$R$. This implies that for each $x$, $\operatorname{sgn}_{y}f(x,y)$ must have
only one jump discontinuity on $I$. Thus, the sufficient condition of the
implicit function theorem guarantees Assumption 1.

######  Example 2.1.

  * $(i)$

Let $f(x,y)=x^{2}+y^{2}-1$ be a function with $f(0,1)=0$. We want to find a
function for $y$ such that $f=0$ in some closed interval of $x=0$. Let us
consider $R\times I=[-1,1]\times[0,2]\ni(0,1)$. Then
$\operatorname{sgn}_{y}f(x,y)$ has only one jump discontinuity on $I$. (On the
other hand, since $\partial_{y}f(0,1)=2\neq 0$, by application of the implicit
function theorem there is an $R\times I$ of $(0,1)$ such that for each $x\in
R$, $\operatorname{sgn}_{y}f(x,y)$ has the only one jump discontinuity on
$I$.)

  * $(ii)$

Let $f(x,y)=y^{2}-x^{4}$ be a function with $f(0,0)=0$. Although
$\partial_{y}f(0,0)=0$, by taking $R\times I=[-1,1]\times[0,1]$ or
$[-1,1]\times[-1,0]\ni(0,0)$, we can justify that $f$ satisfies Assumption 1.

For each $x\in R$, if we define $n_{y}(f)(x)=1$ if
$\operatorname{sgn}_{y}f(x,y)$ is increasing and $-1$ if
$\operatorname{sgn}_{y}f(x,y)$ is decreasing, then, $n_{y}(f)$ is constant.
Indeed, by the intermediate value theorem we may let $y(x)\in I$ be a unique
solution such that $f(x,y(x))=0$ for each $x\in R$. Suppose that there is
$\xi_{1}$ and $\eta_{1}$ such that $n_{y}(f)(\xi_{1})n_{y}(f)(\eta_{1})=-1$.
We may assume that $n_{y}(f)(\xi_{1})=1$ and $n_{y}(f)(\eta_{1})=-1$. We
choose $r>0$ such that $y(x)+r\in I$ or $y(x)-r\in I$ for every $x\in R$. Let
$x_{0}$ be the midpoint of the line segment $\overline{\xi_{1}\eta_{1}}$
between $\xi_{1}$ and $\eta_{1}$. If $f(x_{0},y(x_{0})+r)>0$ or
$f(x_{0},y(x_{0})-r)<0$, then select $\overline{x_{0}\eta_{1}}$, whereas if
$f(x_{0},y(x_{0})+r)<0$ or $f(x_{0},y(x_{0})-r)>0$, then select
$\overline{\xi_{1}x_{0}}$. Write the selected one as
$\overline{\xi_{2}\eta_{2}}$. This procedure reaches the same situation of
$n_{y}(f)(\xi_{2})=1$ and $n_{y}(f)(\eta_{2})=-1$ as the beginning. Repeating
the above process inductively by taking the midpoint of
$\overline{\xi_{2}\eta_{2}}$, we have the monotone sequences $\xi_{k}$ and
$\eta_{k}$, so that $\xi_{k}\nearrow\zeta$ and $\eta_{k}\searrow\zeta$ for
some $\zeta\in R$. The convergence of $\xi_{k}$ and $\eta_{k}$ and the
continuity of $f$ yield the convergence of $y(\xi_{k})$ and $y(\eta_{k})$,
i.e., $y(\xi_{k})$ and $y(\eta_{k})$ go to $y(\zeta)$ as $k\to\infty$.

It follows that $f(\xi_{k},y(\xi_{k})+r)>0$ (or $f(\xi_{k},y(\xi_{k})-r)<0$)
and $f(\eta_{k},y(\eta_{k})+r)<0$ (or $f(\eta_{k},y(\eta_{k})-r)>0$) for every
$k$. Taking the limit $k\to\infty$, we have $f(\zeta,y(\zeta)+r)>0$ (or
$f(\zeta,y(\zeta)-r)<0$) and $f(\zeta,y(\zeta)+r)<0$ (or
$f(\zeta,y(\zeta)-r)>0$). This leads contradiction. Therefore, we conclude the
following lemma.

######  Lemma 2.1.

In $R$, $n_{y}(f)$ is constant (either $1$ or $-1$).

Let $\Theta$ be the Heaviside step function. For each $x$, define a function
$\Theta_{y}f(x,y)$ that assigns $1$ if $f(x,y)\geq 0$ and $0$ otherwise. For a
rectangle $R^{\prime}\subset R$, we define a quantity as

$\iint_{R^{\prime}\times I}\Theta_{y}f(x,y)\,dx\,dy=v_{R^{\prime}\times I}.$
(1)

The next result provides a necessary condition for the existence of a function
form of an implicit equation.

######  Lemma 2.2.

If $R\times I$ contains a unique function $y=g(x)$ such that $f(x,y)=0$, then

$v_{R^{\prime}\times I}=\frac{n_{y}(f)+1}{2}|R^{\prime}|\max
I+\frac{n_{y}(f)-1}{2}|R^{\prime}|\min I-n_{y}(f)\int_{R^{\prime}}g\,dx$ (2)

for every rectangle $R^{\prime}\subset R$.

By Lemma 2.1, the right-hand side of (2) is well defined.

###### Proof of Lemma 2.2.

Let $R^{\prime}\subset R$ be a rectangle. Suppose that $n_{y}(f)=1$. The
integration (1) is calculated as

$\iint_{R^{\prime}\times I}\Theta_{y}f(x,y)\,dx\,dy=|R^{\prime}|\max
I-\int_{R^{\prime}}g\,dx.$ (3)

If $n_{y}(f)=-1$, then by putting $F(x,y)=f(x,-y)$, we have $n_{y}(F)=1$ on
$-I$ with $y=-g(x)$ which satisfies $F(x,y)=0$. By (3) with $F$ and $y=-g(x)$
on $R\times(-I)$ instead of $f$ and $y=g(x)$,

$\iint_{R^{\prime}\times(-I)}\Theta_{y}F(x,y)\,dx\,dy=|R^{\prime}|\max(-I)+\int_{R^{\prime}}g\,dx.$
(4)

By the change of variables ($-y\mapsto y$), (4) equals

$\iint_{R^{\prime}\times
I}\Theta_{y}f(x,y)\,dx\,dy=-|R^{\prime}|\min(I)+\int_{R^{\prime}}g\,dx.$ (5)

According to (3) and (5), we conclude (2). Therefore, the proof is complete. ∎

Now, we have the following integral representation of an implicit function.

######  Theorem 2.3.

The continuous implicit function $y=g(x)$ such that $f=0$ in $R$, is given by

$g(x)=\frac{1+n_{y}(f)}{2}\max I+\frac{1-n_{y}(f)}{2}\min
I-n_{y}(f)\int_{I}\Theta_{y}f(x,y)\,dy\qquad(x\in R).$ (6)

Equation (6) does not depend on the choice of a rectangle $R\times I$ whenever
it contains $y=g(x)$.

###### Proof of Theorem 2.3.

Let $x\in R$ and take cubes $Q\subset R$ that shrink to a singleton $\\{x\\}$.
Divide $|Q|$ into both sides of (1) and (2) after replacing $R^{\prime}$ with
$Q$. By Fubini’s theorem and by the Lebesgue differentiation theorem,

$\int_{I}\Theta_{y}f(x,y)\,dy=\frac{n_{y}(f)+1}{2}\max
I+\frac{n_{y}(f)-1}{2}\min I-n_{y}(f)g(x)$

for every point of continuity of $f$. Therefore, the proof is complete. ∎

Now, we set up an algebraic operator that will appear in the main theorems.
Let $C=(c_{i_{1}\cdots i_{N}})$ and $A=(a_{j_{1}j_{2}})$ be an
$m_{1}\times\cdots\times m_{N}$ tensor and $n_{1}\times n_{2}$ matrix,
respectively. If $n_{1}=m_{k}$, then we define a tensor contraction between
$A$ and $C$ by

$A\stackrel{{\scriptstyle 1\to
k}}{{\cdot}}C=\left(\sum_{i_{k}=1}^{m_{k}}a_{i_{k}j_{2}}c_{\cdots
i_{k}\cdots}\right)$

which produces an $m_{1}\times\cdots\times m_{k-1}\times n_{2}\times
m_{k+1}\times\cdots\times m_{N}$ tensor. If $n_{2}=m_{k}$, then
$A\stackrel{{\scriptstyle 2\to k}}{{\cdot}}C$ is also defined, and it becomes
an $m_{1}\times\cdots\times m_{k-1}\times n_{1}\times
m_{k+1}\times\cdots\times m_{N}$ tensor. For another $n_{1}^{\prime}\times
n_{2}^{\prime}$ matrix $B$, if $n_{1}=m_{i}$ and $n_{2}^{\prime}=m_{j}$
$(i\neq j)$, then we have the commutative property of

$A\stackrel{{\scriptstyle 1\to i}}{{\cdot}}(B\stackrel{{\scriptstyle 2\to
j}}{{\cdot}}C)=B\stackrel{{\scriptstyle 2\to
j}}{{\cdot}}(A\stackrel{{\scriptstyle 1\to i}}{{\cdot}}C).$

Here, the parentheses are omitted providing the commutative property.

##  3 Polynomial implicit functions

In this section, we derive the implicit function of $y=g(x)$ such that
$f(x,y)=0$, provided $y=g(x)$ is a unique multivariate polynomial on $R\times
I$. Put $R=\prod_{k=1}^{n}[\xi_{k},\eta_{k}]$. Let $N_{k}$ be a partition
number and
$\xi_{k}<\xi_{k}+\Delta_{k}<\xi_{k}+2\Delta_{k}<\cdots<\xi_{k}+N_{k}\Delta_{k}=\eta_{k}$
a partition on $[\xi_{k},\eta_{k}]$, where
$\Delta_{k}=(\eta_{k}-\xi_{k})/N_{k}$ is a partition size. We define a grid
block $R_{\alpha}$ by

$R_{\alpha}=\prod_{k=1}^{n}\left[\xi_{k}+(\alpha_{k}-1)\Delta_{k},\xi_{k}+\alpha_{k}\Delta_{k}\right]\quad(1\leq\alpha_{k}\leq
N_{k}),$ (7)

where $1\leq\alpha_{k}\leq N_{k}$ $(1\leq k\leq n)$. The union of all
$R_{\alpha}$ is $R$ and
$|R_{\alpha}|=\prod_{k=1}^{n}(\eta_{k}-\xi_{k})/N_{k}$. For
$a\in\mathbb{R}^{n}$, we define a matrix by

$V_{N_{k}}=\begin{pmatrix}\left(\xi_{k}+\alpha_{k}\Delta_{k}-a_{k}\right)^{j}-\left(\xi_{k}+(\alpha_{k}-1)\Delta_{k}-a_{k}\right)^{j}\end{pmatrix}_{N_{k}\times
N_{k}}\quad(1\leq\alpha_{k},\,j\leq N_{k}),$ (8)

where $\alpha_{k}$ and $j$ denote a row and column number, respectively.

######  Lemma 3.1.

For every $N_{k}$ $(1\leq k\leq n)$, $|V_{N_{k}}|$ is positive and independent
of $a$. Consequently, $V_{N_{k}}$ is invertible.

###### Proof.

Fix $k$ $(1\leq k\leq n)$. Writing $\gamma_{k}=\xi_{k}-a_{k}$, we get

$|V_{N_{k}}|=\left|\begin{matrix}\left(\left(\gamma_{k}+\alpha_{k}\Delta_{k}\right)^{j}-\left(\gamma_{k}+(\alpha_{k}-1)\Delta_{k}\right)^{j}\right)_{(k,j)}\end{matrix}\right|.$
(9)

From the invariant properties of determinants, add the first row to the
second, the second row to the third, and so on until the end. Then (9) equals

$\left|\;\begin{matrix}\Delta_{k}&\left(\gamma_{k}+\Delta_{k}\right)^{2}-\gamma_{k}^{2}&\cdots&\left(\gamma_{k}+\Delta_{k}\right)^{N_{k}}-\gamma_{k}^{N_{k}}\\\
2\Delta_{k}&\left(\gamma_{k}+2\Delta_{k}\right)^{2}-\gamma_{k}^{2}&\cdots&\left(\gamma_{k}+2\Delta_{k}\right)^{N_{k}}-\gamma_{k}^{N_{k}}\\\
\vdots&\vdots&\ddots&\vdots\\\
N_{k}\Delta_{k}&\left(\gamma_{k}+N_{k}\Delta_{k}\right)^{2}-\gamma_{k}^{2}&\cdots&\left(\gamma_{k}+N_{k}\Delta_{k}\right)^{N_{k}}-\gamma_{k}^{N_{k}}\\\
\end{matrix}\right|$ (10)

which is also equal to

$\left|\;\begin{matrix}1&\gamma_{k}&\gamma_{k}^{2}&\cdots&\gamma_{k}^{N_{k}}\\\
0&\Delta_{k}&\left(\gamma_{k}+\Delta_{k}\right)^{2}-\gamma_{k}^{2}&\cdots&\left(\gamma_{k}+\Delta_{k}\right)^{N_{k}}-\gamma_{k}^{N_{k}}\\\
0&2\Delta_{k}&\left(\gamma_{k}+2\Delta_{k}\right)^{2}-\gamma_{k}^{2}&\cdots&\left(\gamma_{k}+2\Delta_{k}\right)^{N_{k}}-\gamma_{k}^{N_{k}}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&N_{k}\Delta_{k}&\left(\gamma_{k}+N_{k}\Delta_{k}\right)^{2}-\gamma_{k}^{2}&\cdots&\left(\gamma_{k}+N_{k}\Delta_{k}\right)^{N_{k}}-\gamma_{k}^{N_{k}}\\\
\end{matrix}\right|$ (11)

by appending the first row and column of (11) to (10). By adding the first row
of (11) to all other rows, we have

$\mbox{(\ref{determinant2})}=\left|\;\begin{matrix}1&\gamma_{k}&\gamma_{k}^{2}&\cdots&\gamma_{k}^{N_{k}}\\\
1&\gamma_{k}+\Delta_{k}&\left(\gamma_{k}+\Delta_{k}\right)^{2}&\cdots&\left(\gamma_{k}+\Delta_{k}\right)^{N_{k}}\\\
1&\gamma_{k}+2\Delta_{k}&\left(\gamma_{k}+2\Delta_{k}\right)^{2}&\cdots&\left(\gamma_{k}+2\Delta_{k}\right)^{N_{k}}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
1&\gamma_{k}+N_{k}\Delta_{k}&\left(\gamma_{k}+N_{k}\Delta_{k}\right)^{2}&\cdots&\left(\gamma_{k}+N_{k}\Delta_{k}\right)^{N_{k}}\\\
\end{matrix}\right|$ (12)

which is

$\Delta_{k}\prod_{0\leq i<j\leq N_{k}}(j-i)>0,$

where the quantity is strictly positive for any $N_{k}$ and independent of
$a$. Therefore, the proof is complete. ∎

For notational simplicity, write $v_{R_{\alpha}\times I}=v_{\alpha}$ in (1)
and put

$d_{\alpha}=\frac{1+n_{y}(f)}{2}|R_{\alpha}|\max
I+\frac{1-n_{y}(f)}{2}|R_{\alpha}|\min I-n_{y}(f)v_{\alpha}$ (13)

which comes only from $f$. By Lemma 2.2, $d_{\alpha}$ can also be calculated
as

$d_{\alpha}=\int_{R_{\alpha}}g\,dx.$ (14)

Let $W$ be an $N_{1}\times\cdots\times N_{n}$ tensor whose $\alpha$th
component is $\alpha_{1}\cdots\alpha_{n}$. From now on, set
$N=(N_{1},\ldots,N_{n})$.

######  Theorem 3.2.

If $y=g(x)$ is a multivariate polynomial from $R$ to $I$ such that $f(x,y)=0$,
then for $a\in R$, $g(x)=\sum_{0\leq\beta<N}c_{\beta}(x-a)^{\beta}$, where
$c_{\beta}$ is the $(\beta+1)$th component of

$W\circ\Big{[}V_{N_{n}}^{-1}\stackrel{{\scriptstyle 2\to n}}{{\cdot}}\cdots
V_{N_{2}}^{-1}\stackrel{{\scriptstyle 2\to
2}}{{\cdot}}V_{N_{1}}^{-1}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}\begin{pmatrix}d_{\alpha}\end{pmatrix}_{\alpha}\Big{]}\qquad(1\leq\alpha\leq
N),$ (15)

$N_{k}>$ the largest exponent of $x_{k}$, and $\circ$ denotes the Hadamard
product between two tensors.

In Theorem 3.2, every component of (15) vanishes if its index contains a
number larger than the largest exponent of $x_{k}$ for some $k$. Now we will
call (15) the coefficient tensor for $g$.

###### Proof of Theorem 3.2.

For each $k$, take a sufficiently large $N_{k}$ so that $N_{k}>$ the largest
exponent of $x_{k}$ in $g$. By (7),

$\begin{split}\int_{R_{\alpha}}g\,dx&=\sum_{0\leq\beta<N}c_{\beta}\int_{R_{\alpha}}(x-a)^{\beta}\,dx\\\
&=\sum_{0\leq\beta<N}c_{\beta}\prod_{k=1}^{n}\frac{1}{\beta_{k}+1}\left[\left(\xi_{k}+\alpha_{k}\Delta_{k}-a_{k}\right)^{\beta_{k}+1}-\left(\xi_{k}+(\alpha_{k}-1)\Delta_{k}-a_{k}\right)^{\beta_{k}+1}\right].\end{split}$
(16)

By (14), the sum of (16) is equal to $d_{\alpha}$. Moreover, (16) is the
$\alpha$th component of

$V_{N_{1}}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}V_{N_{2}}\stackrel{{\scriptstyle 2\to 2}}{{\cdot}}\cdots
V_{N_{n}}\stackrel{{\scriptstyle 2\to
n}}{{\cdot}}\begin{pmatrix}c_{\beta}/(\beta_{1}+1)\cdots(\beta_{n}+1)\end{pmatrix}_{\beta},$
(17)

where $0\leq\beta\leq N-1$. Since $V_{N_{k}}$ is invertible by Lemma 3.1, we
have

$\begin{pmatrix}c_{\beta}/(\beta_{1}+1)\cdots(\beta_{n}+1)\end{pmatrix}_{\beta}=V_{N_{n}}^{-1}\stackrel{{\scriptstyle
2\to n}}{{\cdot}}\cdots V_{N_{2}}^{-1}\stackrel{{\scriptstyle 2\to
2}}{{\cdot}}V_{N_{1}}^{-1}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}\begin{pmatrix}d_{\alpha}\end{pmatrix}_{\alpha}.$

So,

$\begin{pmatrix}c_{\beta}\end{pmatrix}_{\beta}=W\circ\Big{[}V_{N_{n}}^{-1}\stackrel{{\scriptstyle
2\to n}}{{\cdot}}\cdots V_{N_{2}}^{-1}\stackrel{{\scriptstyle 2\to
2}}{{\cdot}}V_{N_{1}}^{-1}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}\begin{pmatrix}d_{\alpha}\end{pmatrix}_{\alpha}\Big{]}.$ (18)

For the $N^{\prime}$-partition of $R$ $(N^{\prime}\geq N)$, according to (7),
(8), (13), (14), and (18), let $\big{(}c_{\beta}^{\prime}\big{)}_{\beta}$ be
the coefficient tensor for
$g(x)=\sum_{0\leq\beta<N^{\prime}}c_{\beta}^{\prime}(x-a)^{\beta}$ which
satisfies $f=0$. From the uniqueness of the implicit function, the nontrivial
components of two tensors should be identical, but they are different only in
their sizes by adding zero components. This implies that $c_{\beta}$ does not
depend on the choice of $N$ whenever $N_{k}$ is greater than the largest
exponent of $x_{k}$ of $g(x)$. Therefore, (18) results in the desired
conclusion (15). ∎

Since the necessary condition of the implicit function theorem with $m=1$
implies Assumption 1, as mentioned before, we have a corollary.

######  Corollary 3.2.1.

Suppose that $f$ is continuously differentiable such that
$\partial_{y}f(a,b)\neq 0$ with $f(a,b)=0$. If $y=g(x)$ is a multivariate
polynomial such that $f(x,y)=0$ near $(a,b)$, then there is an $R\times I$ so
that for $a\in R$, $g(x)=\sum_{0\leq\alpha<N}c_{\alpha}(x-a)^{\alpha}$ with
the coefficient tensor of (15).

The following example is about a polynomial implicit function with two
independent variables.

######  Example 3.1.

Let $f:\mathbb{R}^{3}\to\mathbb{R}$ be given by

$\begin{split}f(x,y,z)&=0.5x^{4}+0.5x^{3}y+0.5x^{3}+2x^{2}y+0.5x^{2}z+0.5xy^{2}\\\
&\qquad-0.5xyz+1.5xy-0.5xz+x+1.5y^{2}-0.5yz+2y-z^{2}+3z-2\end{split}$

with $(2,0,-2)$ and $(1,1,4)$ at which $f$ vanishes. We want to solve $f=0$
for $z$ as a function of $x$ and $y$.

First, since $\partial_{z}f(2,0,-2)=8\neq 0$, by application of the implicit
function theorem there is a rectangle of $(2,0,-2)$ on which
$\operatorname{sgn}_{y}f(x,y)$ has the only one jump discontinuity. Choose
$R=[1.5,2.5]\times[-0.5,0.5]\ni(2,0)$ and $I=[-5,1]\ni-2$, so that
$n_{z}(f)=1$ in $R\times I$. The surface of $f=0$ and $R\times I$ are depicted
in Figure 1. With $N_{1}=N_{2}=3$, $V_{N_{1}}$, $V_{N_{2}}$ and $W$ are
calculated as

$\begin{pmatrix}\scalebox{0.8}{$\Delta_{x}$}&\scalebox{0.8}{$(-0.5+\Delta_{x})^{2}-(-0.5)^{2}$}&\scalebox{0.8}{$(-0.5+\Delta_{x})^{3}-(-0.5)^{3}$}\\\
\scalebox{0.8}{$\Delta_{x}$}&\scalebox{0.8}{$(-0.5+2\Delta_{x})^{2}-(-0.5+\Delta_{x})^{2}$}&\scalebox{0.8}{$(-0.5+2\Delta_{x})^{3}-(-0.5+\Delta_{x})^{3}$}\\\
\scalebox{0.8}{$\Delta_{x}$}&\scalebox{0.8}{$(-0.5+3\Delta_{x})^{2}-(-0.5+2\Delta_{x})^{2}$}&\scalebox{0.8}{$(-0.5+3\Delta_{x})^{3}-(-0.5+2\Delta_{x})^{3}$}\end{pmatrix},$
$\begin{pmatrix}\scalebox{0.8}{$\Delta_{y}$}&\scalebox{0.8}{$(-0.5+\Delta_{y})^{2}-(-0.5)^{2}$}&\scalebox{0.8}{$(-0.5+\Delta_{y})^{3}-(-0.5)^{3}$}\\\
\scalebox{0.8}{$\Delta_{y}$}&\scalebox{0.8}{$(-0.5+2\Delta_{y})^{2}-(-0.5+\Delta_{y})^{2}$}&\scalebox{0.8}{$(-0.5+2\Delta_{y})^{3}-(-0.5+\Delta_{y})^{3}$}\\\
\scalebox{0.8}{$\Delta_{y}$}&\scalebox{0.8}{$(-0.5+3\Delta_{y})^{2}-(-0.5+2\Delta_{y})^{2}$}&\scalebox{0.8}{$(-0.5+3\Delta_{y})^{3}-(-0.5+2\Delta_{y})^{3}$}\end{pmatrix},$

and

$\begin{pmatrix}\scalebox{0.8}{$1\cdot 1$}&\scalebox{0.8}{$1\cdot
2$}&\scalebox{0.8}{$1\cdot 3$}\\\ \scalebox{0.8}{$2\cdot
1$}&\scalebox{0.8}{$2\cdot 2$}&\scalebox{0.8}{$2\cdot 3$}\\\
\scalebox{0.8}{$3\cdot 1$}&\scalebox{0.8}{$3\cdot 2$}&\scalebox{0.8}{$3\cdot
3$}\\\ \end{pmatrix},$

respectively. By (13), the matrix of $d_{(i,j)}$ is given by

$\begin{pmatrix}\scalebox{0.8}{$-0.049896967860677$}&\scalebox{0.8}{$-0.136316792736493$}&\scalebox{0.8}{$-0.222736578394150$}\\\
\scalebox{0.8}{$-0.130143979013675$}&\scalebox{0.8}{$-0.222736592404210$}&\scalebox{0.8}{$-0.315329196303353$}\\\
\scalebox{0.8}{$-0.222736603139157$}&\scalebox{0.8}{$-0.321502042403124$}&\scalebox{0.8}{$-0.420267478899573$}\end{pmatrix}.$

By Corollary 3.2.1, the desired function
$z=g(x,y)=\sum_{0\leq\alpha<3}c_{\alpha}(x-2)^{\alpha_{1}}y^{\alpha_{2}}$ is
determined by the coefficient matrix of

$W\circ\Big{[}V_{N_{2}}^{-1}\stackrel{{\scriptstyle 2\to
2}}{{\cdot}}V_{N_{1}}^{-1}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}\begin{pmatrix}d_{\alpha}\end{pmatrix}_{\alpha}\Big{]}\\\
=\begin{matrix}&\\\ \left(\hskip 4.5pt\vphantom{\begin{matrix}12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y$}&\scalebox{0.8}{$y^{2}$}\\\
\hline\cr\\\\[-11.38109pt]
\scalebox{0.8}{$-2$}&\scalebox{0.8}{$-2.5$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-2.500001$}&\boxed{\scalebox{0.8}{$-0.499998$}}&\scalebox{0.8}{$-2\times
10^{-6}$}\\\ \scalebox{0.8}{$-0.499999$}&\scalebox{0.8}{$-3\times
10^{-6}$}&\scalebox{0.8}{$4\times 10^{-6}$}\end{matrix}\begin{matrix}&\\\
\left.\vphantom{\begin{matrix}12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x-2$}\\\
\scalebox{0.8}{$(x-2)^{2}$}\end{matrix}\end{matrix}\\\ ,$ (19)

where $-0.499998$, for example, denotes the coefficient of $(x-2)y$ in the
summation of $g$.

On the other hand, since $\partial_{z}f(1,1,4)=-6\neq 0$, by application of
the implicit function theorem there is a rectangle, for example,
$R=[0.5,1.5]\times[0.5,1.5]\ni(1,1)$ and $I=[2,7]\ni 4$, on which
$n_{z}(f)=-1$. For $N_{1}=N_{2}=4$ (as shown in Figure 1), we derive the
coefficient matrix of

$\begin{matrix}&\\\ \left(\hskip 4.5pt\vphantom{\begin{matrix}12\\\ 12\\\
12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y-1$}&\scalebox{0.8}{$(y-1)^{2}$}&\scalebox{0.8}{$(y-1)^{3}$}\\\
\hline\cr\\\\[-11.38109pt]
\scalebox{0.8}{$4$}&\scalebox{0.8}{$1.000001$}&\scalebox{0.8}{$-1\cdot
10^{-6}$}&\scalebox{0.8}{$1\cdot 10^{-6}$}\\\
\scalebox{0.8}{$2.000001$}&\scalebox{0.8}{$-1\cdot
10^{-6}$}&\scalebox{0.8}{$2\cdot 10^{-6}$}&\scalebox{0.8}{$-2\cdot
10^{-6}$}\\\ \scalebox{0.8}{$0.999999$}&\scalebox{0.8}{$2\cdot
10^{-6}$}&\scalebox{0.8}{$-3\cdot 10^{-6}$}&\scalebox{0.8}{$4\cdot
10^{-6}$}\\\ \scalebox{0.8}{$0$}&\scalebox{0.8}{$-1\cdot
10^{-6}$}&\scalebox{0.8}{$3\cdot 10^{-6}$}&\scalebox{0.8}{$-4\cdot
10^{-6}$}\end{matrix}\begin{matrix}&\\\ \left.\vphantom{\begin{matrix}12\\\
12\\\ 12\\\ 12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x-1$}\\\ \scalebox{0.8}{$(x-1)^{2}$}\\\
\scalebox{0.8}{$(x-1)^{3}$}\end{matrix}\end{matrix}$ (20)

of
$z=g(x,y)=\sum_{0\leq\alpha<4}c_{\alpha}(x-1)^{\alpha_{1}}(y-1)^{\alpha_{2}}$
that satisfies $f=0$ in $R$ by the same method shown above.

In fact, $f$ is factorized by

$f_{1}=\tfrac{5}{2}(x-2)+\tfrac{5}{2}y+z+2+\tfrac{1}{2}(x-2)^{2}+\tfrac{1}{2}(x-2)y$

and

$f_{2}=2(x-1)+(y-1)-z+4+(x-1)^{2}.$

The coefficient matrices of functions for $z$ such that $f_{1}=0$ and
$f_{2}=0$, respectively, are shown as

$\begin{matrix}&\\\ \left(\vphantom{\begin{matrix}12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y$}&\scalebox{0.8}{$y^{2}$}\\\
\hline\cr\\\\[-11.38109pt]
\scalebox{0.8}{$-2$}&\scalebox{0.8}{$-2.5$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-2.5$}&\scalebox{0.8}{$-0.5$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-0.5$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\end{matrix}\hskip
1.79997pt\begin{matrix}&\\\ \left.\vphantom{\begin{matrix}12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x-2$}\\\ \scalebox{0.8}{$(x-2)^{2}$}\end{matrix}\end{matrix}$

and

$\begin{matrix}&\\\ \left(\hskip 4.5pt\vphantom{\begin{matrix}12\\\ 12\\\
12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y-1$}&\scalebox{0.8}{$(y-1)^{2}$}&\scalebox{0.8}{$(y-1)^{3}$}\\\
\hline\cr\\\\[-11.38109pt]
\scalebox{0.8}{$4$}&\scalebox{0.8}{$1$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$2$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$1$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\end{matrix}\begin{matrix}&\\\
\left.\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x-1$}\\\ \scalebox{0.8}{$(x-1)^{2}$}\\\
\scalebox{0.8}{$(x-1)^{3}$}\end{matrix}\end{matrix}$ (21)

which are comparable to (19) and (20), respectively.

Figure 1: Illustration showing that $f(x,y,z)=0$. Two rectangles of
$[1.5,2.5]\times[-0.5,0.5]\times[-5,1]$ and
$[0.5,1.5]\times[0.5,1.5]\times[2,7]$ denote closed neighborhoods of
$(2,0,-1)$ and $(1,1,4)$, respectively.

##  4 Analytic implicit functions

In this section, we derive a power series expansion of an implicit function
which has one dependent variable, and any number of independent variables as
before. The proposed method does not require the size estimate of Taylor
coefficients and even differentiability of $f$ as we have seen in Section 3.
We start with an assumption that the implicit function $y=g(x)$ for
$f(x,y)=0$, is analytic on a rectangle, which means that $g$ is analytic on
some open neighborhood containing the rectangle. By translation and dilation
on the independent variables of $f$, we may assume that the implicit function
is analytic on $[-1-\delta,1+\delta]^{n}$ for some $\delta>0$ and write
$g(x)=\sum_{\beta}c_{\beta}x^{\beta}$ which converges in
$[-1-\delta,1+\delta]^{n}$, where the sum is the power series expansion of $g$
at the origin. Let $R=[-1,1]^{n}$ and put

$-1<-1+\Delta_{x}<-1+2\Delta_{x}<\cdots-1+(N-1)\Delta_{x}<-1+N\Delta_{x}=1$

be a partition on $[-1,1]$ $(\Delta_{x}=2/N)$. According to (13) and (1),
prepare $d_{\alpha}$ from $f$ over $R_{\alpha}\times I$, where

$R_{\alpha}=\prod_{k=1}^{n}[-1+(\alpha_{k}-1)\Delta_{x},-1+\alpha_{k}\Delta_{x}]$

for $\alpha_{k}=1,2,\ldots,N$.

We split $g$ into two parts: one is a partial sum in which the exponent of a
monomial is dominated by an multi-index $N$, another is the remainder of the
series of $g$, e.g.,

$\begin{split}g(x)&=\sum_{0\leq\beta<N}c_{\beta}x^{\beta}+\sum_{\text{remainder
}\beta}c_{\beta}x^{\beta}\\\ &=g_{N}(x)+r_{N}(x),\quad\mbox{say}.\end{split}$
(22)

From analyticity, the series of $g$ converges uniformly and absolutely. There
is a constant $C$ such that $|c_{\beta}x^{\beta}|\leq C$ uniformly in
$x\in[-1-\delta,1+\delta]^{n}$ for every $\beta$. So,
$|c_{\beta}|(1+\delta)^{|\beta|}\leq C$ for every $\beta$ (for the properties
of real analytic functions of several variables, refer to [5]).

In $r_{N}$ of (22), there is an index $\beta$ such that its component is
greater than or equal to $N$. Assume that $\beta_{1}\geq N$. Then, for $x\in
R$,

$\begin{split}\sum_{\text{remainder }\beta}|c_{\beta}||x^{\beta}|&\leq
C\sum_{\text{remainder }\beta}\frac{|x^{\beta}|}{(1+\delta)^{|\beta|}}\\\
&\leq\sum_{\beta_{1}=N}^{\infty}\sum_{\begin{subarray}{c}\beta_{k}=0\\\ k\neq
1\end{subarray}}^{\infty}\frac{1}{(1+\delta)^{|\beta|}}\\\
&=\frac{1}{\delta^{n}(1+\delta)^{N-n}}\end{split}$ (23)

which vanishes as $N\to\infty$. Let $\epsilon>0$. By the triangle inequality
and by (23), we have $\|r_{N}\|_{L^{\infty}(R)}<\epsilon$ and

$\int_{R_{\alpha}}|r_{N}|\,dx\leq|R_{\alpha}|\epsilon\leq 2^{n}\epsilon$ (24)

for every sufficiently large $N$.

For a fixed $N$, we put
$\tilde{g}_{N}(x)=\sum_{0\leq\alpha<N}\tilde{c}_{\alpha}x^{\alpha}$, where
$\tilde{c}_{\alpha}$ is calculated by solving

$\int_{R_{\alpha}}\tilde{g}_{N}\,dx=d_{\alpha}$ (25)

as the derivation of (15). By (14) and by (25),

$\begin{split}\left|\int_{R_{\alpha}}g_{N}-\tilde{g}_{N}\,dx\right|&=\left|\int_{R_{\alpha}}g_{N}\,dx-
d_{\alpha}\right|\\\ &=\left|\int_{R_{\alpha}}g_{N}-g\,dx\right|\\\
&=\left|\int_{R_{\alpha}}r_{N}\,dx\right|,\end{split}$ (26)

where the last equality follows from (22). By the Riesz representation theorem
for the Lebesgue spaces (the duality argument), we realize the next lemma.

######  Lemma 4.1.

The function $\tilde{g}_{N}$ converges to $g$ weak-star in $L^{\infty}(R)$ as
$N\to\infty$.

###### Proof.

We recall (26). By the triangle inequality,

$\int_{R_{\alpha}}|\tilde{g}_{N}-g|dx\leq\int_{R_{\alpha}}|\tilde{g}_{N}-g_{N}|dx+\int_{R_{\alpha}}|r_{N}|dx\to
0$ (27)

uniformly on $\alpha$ as $N\to\infty$. Here, the convergence of (27) follows
from (24). Since the collection of finite linear combinations of
characteristic functions $\chi_{R_{\alpha}}$ which are supported on grid
blocks $R_{\alpha}$ is dense in $L^{1}(R)$, the duality argument ([3, 12])
from (26) and (27) for Lebesgue spaces yield that $\tilde{g}_{N}$ goes to $g$
weak-star in $L^{\infty}(R)$ as $N\to\infty$. Therefore, the proof is
complete. ∎

Let $V_{N}$ and $W$ be a matrix and tensor as in Theorem 3.2 with
$V_{N}=V_{N_{k}}$ and $N=N_{k}$. By Lemma 4.1, we conclude the main theorem of
this section.

######  Theorem 4.2.

If $y=g(x)$ is analytic on $R$ to $I$ such that $f(x,y)=0$, then
$\tilde{g}_{N}(x)=\sum_{0\leq\alpha<N}\tilde{c}_{\alpha}x^{\alpha}$ with the
coefficient tensor of

$W_{N}\circ\Big{[}V_{N}^{-1}\stackrel{{\scriptstyle 2\to n}}{{\cdot}}\cdots
V_{N}^{-1}\stackrel{{\scriptstyle 2\to
2}}{{\cdot}}V_{N}^{-1}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}\begin{pmatrix}d_{\alpha}\end{pmatrix}_{\alpha}\Big{]},$ (28)

converges to $g$ weak-star in $L^{\infty}(R)$ as $N\to\infty$.

As we have seen in Corollary 3.2.1, the next corollary follows.

######  Corollary 4.2.1.

Suppose that $f$ is continuously differentiable such that
$\partial_{y}f(0,0)\neq 0$ with $f(0,0)=0$. If $y=g(x)$ is analytic on $R$ to
$I$ such that $f(x,y)=0$, then then $\tilde{g}_{N}$ converges to $g$ weak-star
in $L^{\infty}(R)$ as $N\to\infty$.

If $f$ is analytic such that $\partial_{y}f(0,0)\neq 0$ with $f(0,0)=0$, then
by the analytic version of the implicit function theorem, Corollary 4.2.1
holds. We now present a numerical approximation of an analytic implicit
function.

######  Example 4.1.

Let $f:\mathbb{R}^{3}\to\mathbb{R}$ be given by

$\begin{split}f(x,y,z)=x^{2}+y^{2}+z^{2}-1\end{split}$

with $f(0,0,1)=0$. We want to find a function of $x$ and $y$ that satisfies
$f=0$ in some rectangle of $(0,0,1)$. Since $f$ is analytic on a neighborhood
of $(0,0,1)$ and $\partial_{z}f(0,0,1)=2>0$, by the analytic version of the
implicit function theorem there is a rectangle of $(0,0,1)$ and a unique
analytic implicit function for $z$ such that $f=0$ in the rectangle.

Choose a rectangle $R\times I$, for example,
$R=[-1/2,1/2]\times[-1/2,1/2]\ni(0,0)$ and $I=[0,3/2]\ni 1$, which are
depicted in Figure 2. In $R\times I$, readily $n_{z}(f)=1$ and with $N=6$, it
follows that

$V_{N}=\begin{pmatrix}\scalebox{0.8}{$\Delta$}&\scalebox{0.8}{$(-1/2+\Delta)^{2}-(-1/2)^{2}$}&\scalebox{0.8}{$\cdots$}&\scalebox{0.8}{$(-1/2+\Delta)^{6}-(-1/2)^{6}$}\\\
\scalebox{0.8}{$\Delta$}&\scalebox{0.8}{$(-1/2+2\Delta)^{2}-(-1/2+\Delta)^{2}$}&\scalebox{0.8}{$\cdots$}&\scalebox{0.8}{$(-1/2+2\Delta)^{6}-(-1/2+\Delta)^{6}$}\\\
\scalebox{0.8}{$\vdots$}&\scalebox{0.8}{$\vdots$}&\scalebox{0.8}{$\ddots$}&\scalebox{0.8}{$\vdots$}\\\
\scalebox{0.8}{$\Delta$}&\scalebox{0.8}{$(-1/2+6\Delta)^{2}-(-1/2+5\Delta)^{2}$}&\scalebox{0.8}{$\cdots$}&\scalebox{0.8}{$(-1/2+6\Delta)^{6}-(-1/2+5\Delta)^{6}$}\\\
\end{pmatrix},$

where $\Delta=1/6$. From (13), the $6\times 6$ matrix of $d_{i,j}$ is given by

$\left(\begin{array}[]{cccccc}\scalebox{0.68}{$\rotatebox{0.0}{$0.022341010165457$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.022341010165457$}$}\\\
\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025909540264518$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025909540264518$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}\\\
\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.027518589197521$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.027518589197521$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}\\\
\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.027518589197521$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.027518589197521$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}\\\
\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025909540264518$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.026726310029805$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025909540264518$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}\\\
\scalebox{0.68}{$\rotatebox{0.0}{$0.022341010165457$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.025065830482660$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.024192273341742$}$}&\scalebox{0.68}{$\rotatebox{0.0}{$0.022341010165457$}$}\end{array}\right).$

By Corollary 4.2.1, the coefficient matrix of
$z=\tilde{g}_{N}(x,y)=\sum_{0\leq\alpha<6}c_{\alpha}x^{\alpha_{1}}y^{\alpha_{2}}$
which approximates $z=z(x,y)$ such that $f(x,y,z(x,y))=0$, is calculated as

$W\circ\Big{[}V_{N}^{-1}\stackrel{{\scriptstyle 2\to
2}}{{\cdot}}V_{N}^{-1}\stackrel{{\scriptstyle 2\to
1}}{{\cdot}}\begin{pmatrix}d_{\alpha}\end{pmatrix}_{\alpha}\Big{]}\\\
=\begin{matrix}&\\\ \left(\hskip 3.59995pt\vphantom{\begin{matrix}12\\\ 12\\\
12\\\ 12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y$}&\scalebox{0.8}{$y^{2}$}&\scalebox{0.8}{$y^{3}$}&\scalebox{0.8}{$y^{4}$}&\scalebox{0.8}{$y^{5}$}\\\
\hline\cr\\\\[-8.53581pt]
\scalebox{0.8}{$0.99997$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.49880$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.14537$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-0.49880$}&\scalebox{0.8}{$0$}&\boxed{\scalebox{0.8}{$-0.24363$}}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.23108$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-0.14537$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.23108$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.50484$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\end{matrix}\begin{matrix}&\\\
\left.\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\ 12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x$}\\\ \scalebox{0.8}{$x^{2}$}\\\ \scalebox{0.8}{$x^{3}$}\\\
\scalebox{0.8}{$x^{4}$}\\\ \scalebox{0.8}{$x^{5}$}\end{matrix}\end{matrix},$
(29)

where $W=\big{(}ij\big{)}_{1\leq i,j\leq 6}$ and, for example, $-0.24363$ in
(29) is the coefficient of $x^{2}y^{2}$ a term of $g$. Both $\tilde{g}_{N}$
and $z(x,y)=(1-x^{2}-y^{2})^{1/2}$ are plotted in Figure 3.

On the other hand, the coefficient matrix of the partial sum $g_{N}$ of the
Taylor series of $z=(1-x^{2}-y^{2})^{1/2}$ is as shown as

$\begin{matrix}&\\\ \left(\hskip 3.59995pt\vphantom{\begin{matrix}12\\\ 12\\\
12\\\ 12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y$}&\scalebox{0.8}{$y^{2}$}&\scalebox{0.8}{$y^{3}$}&\scalebox{0.8}{$y^{4}$}&\scalebox{0.8}{$y^{5}$}\\\
\hline\cr\\\\[-8.53581pt]
\scalebox{0.8}{$1$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.5$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.125$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-0.5$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.25$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.1875$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$-0.125$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-1875$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$-0.234375$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\end{matrix}\begin{matrix}&\\\ \left.\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\\\ 12\\\ 12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x$}\\\ \scalebox{0.8}{$x^{2}$}\\\ \scalebox{0.8}{$x^{3}$}\\\
\scalebox{0.8}{$x^{4}$}\\\ \scalebox{0.8}{$x^{5}$}\end{matrix}\end{matrix}.$

As it is indicated on the graph of Figure 4, in the accuracy comparison
between $\tilde{g}_{N}$ and $g_{N}$, the former is more accurate.

Figure 2: Surface of $f(x,y,z)=x^{2}+y^{2}+z^{2}-1=0$. The rectangle $R\times
I=[-0.5,0.5]\times[-0.5,0.5]\times[0,1.5]$ contains $(0,0,1)$.

(a) $z=\tilde{g}_{N}(x,y)$ with $N=6$

(b) $z(x,y)=(1-x^{2}-y^{2})^{1/2}$

Figure 3: The comparison of $z=\tilde{g}_{N}(x,y)$ and $z=z(x,y)$ in
$[-0.5,0.5]\times[-0.5,0.5]$.

Figure 4:  The comparison of $|z(x,y)-g_{N}(x,y)|$ and
$|z(x,y)-\tilde{g}_{N}(x,y)|$. The former is drawn on the top side and the
latter on the bottom side.

In general, the implicit function and inverse function theorems can be thought
of as equivalent formulations of a similar basic idea. The next example leads
to the polynomial approximation of an inverse function for Kepler’s equation.
When an initial value is given, the methods using iterative calculation of a
trajectory that satisfies Kepler’s equation are widely used in practice.
Although the methods using the formula by a function are very useful, which
are formal infinite sums or depend on the expansion method of the Taylor
series, there is no proper approach to achieve high-performance computing. The
following example provides more accurate numerical values than the other
results on Kepler’s equation.

######  Example 4.2.

Let $M=E-\epsilon\sin(E)$ be the standard Kepler equation which fixes
$(\pi,\pi)$ for every $0\leq\epsilon\leq 1$, where $M$, $E$, and $\epsilon$
are the mean anomaly, the eccentric anomaly, and the eccentricity,
respectively. Kepler’s equation has a unique inverse function for every
eccentricity, since $M$ is increasing monotonically.

Let $K(E,M)=E-\epsilon\sin(E)-M:\mathbb{R}^{2}\to\mathbb{R}$. The Taylor
series does not give a good approximation of the inverse function (the
implicit function for $E$ of $K=0$) in the open neighborhood of $M=0$ and
$\pi$, especially when $\epsilon=1$ because $\partial_{E}K$ vanishes. (To
approximate the inverse function, for example, [4] adopts the Lagrange inverse
theorem to obtain a formal series expansion and [6] provides the usability of
a solution of Kepler’s equation for $E$ as the most recent result.)

Choose a rectangle of $(\pi,\pi)$, for example, $[0,2\pi]\times[-\pi,3\pi]$ on
which $n_{E}(K)=1$. By Theorem 4.2, for $N=28$ we have an approximation of $E$
by $\tilde{E}(M)=\sum_{n=0}^{N-1}c_{n}(M-\pi)^{n}$ calculated as

$\tilde{E}(M)=\left(\begin{array}[]{cccccccccccccccccccccccccccc}\scalebox{0.8}{$\rotatebox{90.0}{$3.14159265358979$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$0.500000286048401$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-6.39953880001189\times
10^{-13}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$0.0104167098436714$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$3.52445662156362\times
10^{-12}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$0.000520826826664063$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-8.33322713405475\times
10^{-12}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$3.35293978405063\times
10^{-5}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$8.99171081830078\times
10^{-12}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$1.71069043293255\times
10^{-6}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-6.35421688275359\times
10^{-12}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$1.35045779641222\times
10^{-6}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$9.47892775420489\times
10^{-13}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-1.13558515793038\times
10^{-6}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$3.97758303345101\times
10^{-13}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$7.16518651568157\times
10^{-7}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-3.64509053298631\times
10^{-13}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-2.89790511992230\times
10^{-7}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$1.17074864690306\times
10^{-13}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$7.74501977110289\times
10^{-8}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-1.92321275158235\times
10^{-14}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-1.34780939066348\times
10^{-8}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$2.14518516845533\times
10^{-15}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$1.46671954613640\times
10^{-9}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-1.16947164649083\times
10^{-16}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-9.03766696253492\times
10^{-11}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$3.13882342334049\times
10^{-18}$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$2.40613536296425\times
10^{-12}$}$}\end{array}\right)_{\kern-4.0pt\scalebox{0.6}{$1\\!\times\\!28$}}\kern-4.0pt\begin{pmatrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$M-\pi$}\\\ \scalebox{0.8}{$(M-\pi)^{2}$}\\\
\scalebox{0.8}{$\vdots$}\\\ \scalebox{0.8}{$(M-\pi)^{26}$}\\\
\scalebox{0.8}{$(M-\pi)^{27}$}\end{pmatrix}_{\kern-4.0pt\scalebox{0.6}{$28\\!\times\\!1$}}.$

In Figures 5 and 6, $\tilde{E}$ and $K(M,\tilde{E}(M))$ are shown for
$\epsilon=0.8$, $0.9$, and $1$ (for comparison with related work, e.g., refer
to [10, 6]).

Figure 5: Polynomial approximation $\tilde{E}=\sum_{n=0}^{27}c_{n}(M-\pi)^{n}$
for the implicit function for $E$ of $K=0$ in $[0,2\pi]$ with $\epsilon=0.8$,
$0.9$, and $1$.

Figure 6: The non-vanishing values of $K(M,\tilde{E}(M))$ denote errors.

From Example 4.2, we have seen that the inverse function for $y=f(x)$ is
obtained by finding an implicit function for $x$ such that $F(x,y)=0$ after
setting $F(x,y)=y-f(x)$. For a system of implicit equations, Theorem 5.2 will
be applicable for constructing a vector-valued inverse function.

##  5 System of analytic implicit functions

In this section, we extend the results for a real-valued function in Section 4
to a vector-valued function. We will prove it by a variable-reduction
technique which is the process of eliminating the dependent variables one by
one. Let $\phi$ be a permutation on $\\{1,2,\ldots,m\\}$ that consists of all
axis numbers for dependent variables and put
$I=\prod_{i=1}^{m}I_{i}\subset\mathbb{R}^{m}$, where $I_{i}$ is an interval.
For a vector $b$ and rectangle $I$, the notations $b^{\\{j\\}}$ and
$I^{\\{j\\}}$ denote the $j$th component and edge deletions of $b$ and $I$,
respectively.

######  Assumption 2.

For a continuous vector-valued function $f:\mathbb{R}^{n+m}\to\mathbb{R}^{m}$,
the following conditions are satisfied in descending induction over the number
of dependent variables: there is a permutation $\phi$ and an $R\times
I\subset\mathbb{R}^{n+m}$ such that

  * $(i)$

for $(x,y^{\\{m\\}})\in R\times I^{\\{m\\}}$,

$\operatorname{sgn}_{y_{m}}f_{\phi(m)}(x,y^{\\{m,\\}},y_{m})$

has only one jump discontinuity as a function of $y_{m}$ on $I_{m}$,

  * $(ii)$

for $(x,y^{\\{m,m-1\\}})\in R\times I^{\\{m,m-1\\}}$,

$\operatorname{sgn}_{y_{m-1}}f_{\phi(m-1)}(x,y^{\\{m,m-1\\}},y_{m-1},h_{\phi(m)}(x,y^{\\{m\\}}))$

has only one jump discontinuity as a function of $y_{m-1}$ on $I_{m-1}$, where
$y_{m}=h_{\phi(m)}$ solves $f_{\phi(m)}(x,y^{\\{m\\}},y_{m})=0$ in $R\times
I^{\\{m,m-1\\}}$.

Using $(i)$ and $(ii)$, by the intermediate value theorem, by the continuity
of $f$, and by the uniqueness of jump discontinuity, we know that
$y_{m}=h_{\phi(m)}$ and $y_{m-1}=h_{\phi(m-1)}$ are continuous, uniquely
determined. Hence, Assumption 2 is well explained by descending induction.

######  Example 5.1.

Let $f=(f_{1},f_{2}):\mathbb{R}^{3}\to\mathbb{R}^{2}$ be given by
$f_{1}(x,y,z)=y-1$, $f_{2}(x,y,z)=x$ with $f(0,1,1)=\mathbf{0}$. By
observation of $f_{1}$ and $f_{2}$, we easily know that there is no rectangle
of $(0,1,1)$ on which $\operatorname{sgn}_{x}f_{1}(x,y,z)$,
$\operatorname{sgn}_{z}f_{1}(x,y,z)$, $\operatorname{sgn}_{y}f_{2}(x,y,z)$,
and $\operatorname{sgn}_{z}f_{2}(x,y,z)$ have only one jump discontinuity.
This means that there is no $\phi$ on a set containing the $z$-axis number,
which satisfies Assumption 2. On the other hand, since
$\partial_{y}f_{1}(0,1,1)=1\neq 0$, by the analytic version of the implicit
function theorem, there is a rectangle of $(0,1,1)$ and a unique analytic
function $y=y(x,z)$ such that $f_{1}=0$ in the rectangle on which
$\operatorname{sgn}_{y}f_{1}(x,y,z)$ has only one jump discontinuity.
Moreover,

$\partial_{x}f_{2}(x,y(x,z),z)\big{|}_{(x,z)=(0,1)}=\partial_{x}f_{2}(0,1,1)+\partial_{y}f_{2}(0,1,1)\partial_{x}y(0,1)=1\neq
0.$

Thus, $\operatorname{sgn}_{x}f_{2}(x,y(x,z),z)$ has only one jump
discontinuity in the rectangle. Now the permutation $\phi$ is defined by
$\phi(1)=2$ and $\phi(2)=1$, where $1$ and $2$ denote the axis numbers of $x$
and $y$, respectively, that satisfies Assumption 2 in some rectangle of
$(0,1,1)$.

######  Example 5.2.

Let $f=(f_{1},f_{2}):\mathbb{R}^{3}\to\mathbb{R}^{2}$ be defined by
$f_{1}(x,y,z)=y-1$ and $f_{2}(x,y,z)=(z-1)^{2}$ with $f(0,1,1)=\mathbf{0}$.
Since $f_{2}$ is a squared function, $\operatorname{sgn}_{z}f_{2}(x,y,z)$ has
no jump discontinuity on any rectangle of $(0,1,1)$. However, if we consider
$f_{1}$ and $\partial_{z}f_{2}$ instead of $f_{2}$, then
$\operatorname{sgn}_{y}f_{1}$ and $\operatorname{sgn}_{z}\partial_{z}f_{2}$
have only one jump discontinuity, for example, on
$[-1,1]\times[-1,2]\times[-1,2]\ni(0,1,1)$. We can find $\phi$ by defining
$\phi(1)=1$ and $\phi(2)=2$, which satisfies Assumption 2, where $1$ and $2$
denote the axis numbers of $y$ and $z$, respectively. Note that
$(f_{1},f_{2})$ and $(f_{1},\partial_{z}f_{2})$ have the same implicit
function for $y$ and $z$ in the rectangle.

We need a lemma to prove the main theorem of this section, which will be shown
by using a method of variable elimination.

######  Lemma 5.1.

If $y_{i}=h_{\phi(i)}$ satisfies

$f_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}},h_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}}),h_{\phi(i+1)}(x,y^{\\{m,m-1,\ldots,i+1\\}}),\ldots,h_{\phi(m)}(x,y^{\\{m\\}}))=0$

in $R\times I^{\\{m,m-1,\ldots,i\\}}$ $(i=m,m-1,\ldots,1)$, then $y=g(x)$ such
that $f=\mathbf{0}$, is continuous and determined uniquely by

$\begin{split}g_{1}(x)&=h_{\phi(1)}(x),\\\
g_{2}(x)&=h_{\phi(2)}(x,g_{1}(x)),\\\ &\;\;\vdots\\\
g_{m}(x)&=h_{\phi(m)}(x,g_{1}(x),g_{2}(x),\ldots,g_{m-1}(x))\end{split}$ (30)

in $R$

###### Proof.

For convenience, we assume that $\phi=id$. By regarding $(x,y^{\\{m\\}})$ as
the independent variables of $f$, from $(i)$ there is a unique continuous
function $y_{m}=h_{m}(x,y^{\\{m\\}}):R\times I^{\\{m\\}}\to I_{m}$ such that
$f_{m}(x,y^{\\{m\\}},y_{m})=0$ in $R\times I^{\\{m\\}}$.

Substitute $y_{m}=h_{m}(x,y^{\\{m\\}})$ into $f$ and put
$f(x,y^{\\{m\\}},h_{m}(x,y^{\\{m\\}}))=f^{[1]}(x,y^{\\{m\\}})$. Precisely,

$\begin{split}f_{1}^{[1]}(x,y^{\\{m\\}})&=f_{1}(x,y^{\\{m\\}},h_{m}(x,y^{\\{m\\}})),\\\
&\;\;\vdots\\\
f_{m-1}^{[1]}(x,y^{\\{m\\}})&=f_{m-1}(x,y^{\\{m\\}},h_{m}(x,y^{\\{m\\}})),\\\
f_{m}^{[1]}(x,y^{\\{m\\}})&=0\end{split}$

on $R\times I^{\\{m\\}}$, where the $m-1$ nontrivial components of $f^{[1]}$
do not depend on the $y_{m}$-variable.

For every $(x,y^{\\{m.m-1\\}})\in R\times I^{\\{m,m-1\\}}$, the property of
$(ii)$ shows that

$\operatorname{sgn}_{y_{m-1}}f_{m-1}^{[1]}(x,y^{\\{m\\}})=\operatorname{sgn}_{y_{m-1}}f_{m-1}(x,y^{\\{m\\}},h_{m}(x,y^{\\{m\\}}))$

has one jump discontinuity on $I_{m-1}$. This also produces a unique
continuous function $y_{m-1}=h_{m-1}(x,y^{\\{m,m-1\\}}):R\times
I^{\\{m,m-1\\}}\to I_{m-1}$ such that
$f^{[1]}_{m-1}(x,y^{\\{m,m-1\\}},y_{m-1})=0$.

Substitute $y_{m-1}=h_{m-1}(x,y^{\\{m,m-1\\}})$ into $f^{[1]}$ and similar to
the argument above, put
$f^{[1]}(x,y^{\\{m,m-1\\}},h_{m-1}(x,y^{\\{m,m-1\\}}))=f^{[2]}(x,y^{\\{m,m-1\\}})$,
i.e.,

$\begin{split}f_{1}^{[2]}(x,y^{\\{m,m-1\\}})&=f_{1}^{[1]}(x,y^{\\{m,m-1\\}},h_{m-1}(x,y^{\\{m,m-1\\}})),\\\
&\;\;\vdots\\\
f_{m-2}^{[2]}(x,y^{\\{m,m-1\\}})&=f_{m-2}^{[1]}(x,y^{\\{m,m-1\\}},h_{m-1}(x,y^{\\{m,m-1\\}})),\\\
f_{m-1}^{[2]}(x,y^{\\{m,m-1\\}})&=0,\\\
f_{m}^{[2]}(x,y^{\\{m,m-1\\}})&=0\end{split}$

on $R\times I^{\\{m,m-1\\}}$. Then the $m-2$ nontrivial components of
$f^{[2]}$ clearly do not depend on the variables of $y_{m}$ and $y_{m-1}$.

Similarly, for $(x,y^{\\{m,m-1,m-2\\}})\in R\times I^{\\{m,m-1,m-2\\}}$, the
inductively assumed property of $(ii)$ shows that

$\begin{split}\operatorname{sgn}_{m-2}f_{m-2}^{[2]}(x,y^{\\{m.m-1\\}})&=\operatorname{sgn}_{m-2}f_{m-2}^{[1]}(x,y^{\\{m,m-1\\}},h_{m-1}(x,y^{\\{m,m-1\\}}))\\\
&=\operatorname{sgn}_{m-2}f_{m-2}(x,y^{\\{m,m-1\\}},h_{m-1}(x,y^{\\{m,m-1\\}}),h_{m}(x,y^{\\{m\\}}))\end{split}$

has one jump discontinuity on $I_{m-2}$. Again, we gain a unique continuous
function $y_{m-2}=h_{m-2}(x,y^{\\{m,m-1,m-2\\}}):R\times
I^{\\{m,m-1,m-2\\}}\to I_{m-2}$ which satisfies the implicit equation of
$f^{[2]}_{m-2}(x,y^{\\{m,m-1,m-2\\}},y_{m-2})=0$.

Repeating the above variable-reduction process inductively until reaching
$y_{1}$, we get a unique continuous function $y_{1}=h_{1}(x):R_{n}\to I_{1}$
such that $f^{[m-1]}_{1}(x,h_{1}(x))=0$. Eventually, the unique continuous
function $y=g(x)$ for $f(x,y)=\mathbf{0}$ is calculated as

$\begin{split}y_{1}&=h_{1}(x),\\\ y_{2}&=h_{2}(x,y_{1}),\\\ &\;\;\vdots\\\
y_{m}&=h_{m}(x,y_{1},\ldots,y_{m-1})\end{split}$ (31)

in $R$. Therefore, we obtain (30) and the proof is complete. ∎

Although $y_{i}=h_{\phi(i)}$ in Lemma 5.1 has an integral representation as
shown in (6), which does not reveal itself algebraically. From (30) and (31),
if every $h_{\phi(i)}$ is analytic, then every $g_{i}$ is also analytic, and
vice versa. If these are analytic, then every $g_{i}$ is calculated as the
limit of a sequence of multivariate polynomials in the weak-star topology on
$L^{\infty}$. The necessity of the implicit function theorem and the
analyticity of implicit functions yield the following main theorem and
corollary.

######  Theorem 5.2.

If $f$ is continuously differentiable such that $|J_{f,y}(a,b)|\neq 0$ with
$f(a,b)=\mathbf{0}$ and every component of an implicit function for
$f=\mathbf{0}$ is analytic on a neighborhood of $(a,b)$, then the implicit
function is calculated uniquely in some rectangle of $(a,b)$ as (30).

Since the analytic version of the implicit function theorem guarantees a
unique existence of a system of analytic implicit functions in some rectangle,
we have a corollary.

######  Corollary 5.2.1.

Suppose that every component of $f$ is analytic such that $|J_{f,y}(a,b)|\neq
0$ with $f(a,b)=\mathbf{0}$. Then the implicit function for $f=\mathbf{0}$ is
calculated uniquely in some rectangle of $(a,b)$ as (30) .

###### Proof of Theorem 5.2.

It suffices to show that there is an $R\times I$ and a permutation $\phi$ with
which $y_{i}=h_{\phi(i)}$ is analytic on $R\times I^{\\{m,m-1,\ldots,i\\}}$
$(i=m,m-1,\ldots,1)$, which satisfies Assumption 2. Then by Lemma 5.1, the
desired conclusion follows.

Since that $|J_{f,y}(a,b)|\neq 0$, the column vectors of $J_{f,y}(a,b)$ are
not all zero. There is $i$ such that $\partial_{y_{m}}f_{i}(a,b)\neq 0$ and
define $\phi(m)=i$. Thus there is a unique continuous function
$y_{m}=h_{i}(x,y^{\\{m\\}})$ such that
$f_{i}(x,y^{\\{m\\}},h_{i}(x,y^{\\{m\\}}))=0$ for some rectangle of $(a,b)$.
In addition, by Corollary 4.2.1 with analyticity, $y_{m}=h_{i}$ is analytic on
$R^{(1)}\times\prod_{k=1}^{m-1}I_{k}^{(1)}$ to $I_{m}^{(1)}$ for some
rectangle $R^{(1)}\times\prod_{k=1}^{m}I_{k}^{(1)}\ni(a,b)$.

Put $f^{[1]}(x,y^{\\{m\\}})=f^{\\{i\\}}(x,y^{\\{m\\}},h_{i}(x,y^{\\{m\\}}))$,
where $f^{\\{i\\}}$ is the $i$th component removal from $f$. First, we prove
that $|J_{f^{[1]},y^{\\{m\\}}}(a,b^{\\{m\\}})|\neq 0$. Since
$\partial_{y_{m}}f_{i}(a,b)\neq 0$, at $(a,b^{\\{m\\}})$ the determinant of

$|J_{f^{[1]},y^{\\{m\\}}}|=\begin{vmatrix}\partial_{y_{1}}f_{1}^{[1]}&\cdots&\partial_{y_{m-1}}f_{1}^{[1]}\\\
\vdots&\ddots&\vdots\\\
\partial_{y_{1}}f_{i-1}^{[1]}&\cdots&\partial_{y_{m-1}}f_{i-1}^{[1]}\\\
\partial_{y_{1}}f_{i+1}^{[1]}&\cdots&\partial_{y_{m-1}}f_{i+1}^{[1]}\\\
\vdots&\ddots&\vdots\\\
\partial_{y_{1}}f_{m}^{[1]}&\cdots&\partial_{y_{m-1}}f_{m}^{[1]}\end{vmatrix}$

is equal to

$\frac{(-1)^{i+m}}{\partial_{y_{m}}f_{i}(a,b)}\begin{vmatrix}\partial_{y_{1}}f_{1}^{[1]}&\cdots&\partial_{y_{m-1}}f_{1}^{[1]}&\partial_{y_{m}}f_{1}(a,b)\\\
\vdots&\ddots&\cdots&\vdots&\\\
\partial_{y_{1}}f_{i-1}^{[1]}&\cdots&\partial_{y_{m-1}}f_{i-1}^{[1]}&\partial_{y_{m}}f_{i-1}(a,b)\\\
0&\cdots&0&\partial_{y_{m}}f_{i}(a,b)\\\
\partial_{y_{1}}f_{i+1}^{[1]}&\cdots&\partial_{y_{m-1}}f_{i+1}^{[1]}&\partial_{y_{m}}f_{i+1}(a,b)\\\
\vdots&\ddots&\vdots&\vdots\\\
\partial_{y_{1}}f_{m}^{[1]}&\cdots&\partial_{y_{m-1}}f_{m}^{[1]}&\partial_{y_{m}}f_{m}(a,b)\end{vmatrix}.$
(32)

Let $1\leq k\leq m-1$. Subtract $\partial_{y_{k}}h_{i}$ times the last column
of (32) from the other columns. Then the $k$th column of (32) is calculated as

$\begin{pmatrix}\partial_{y_{k}}f_{1}^{[1]}-\partial_{y_{m}}f_{1}(a,b)\partial_{y_{k}}h_{i}\\\
\vdots\\\
\partial_{y_{k}}f_{i-1}^{[1]}-\partial_{y_{m}}f_{i-1}(a,b)\partial_{y_{k}}h_{i}\\\
-\partial_{y_{m}}f_{i}(a,b)\partial_{y_{k}}h_{i}\\\
\partial_{y_{k}}f_{i+1}^{[1]}-\partial_{y_{m}}f_{i+1}(a,b)\partial_{y_{k}}h_{i}\\\
\vdots\\\
\partial_{y_{k}}f_{m}^{[1]}-\partial_{y_{m}}f_{m}(a,b)\partial_{y_{k}}h_{i}\end{pmatrix}=\begin{pmatrix}\partial_{y_{k}}f_{1}(a,b)\\\
\vdots\\\ \partial_{y_{k}}f_{i-1}(a,b)\\\
-\partial_{y_{m}}f_{i}(a,b)\partial_{y_{k}}h_{i}\\\
\partial_{y_{k}}f_{i+1}(a,b)\\\ \vdots\\\
\partial_{y_{k}}f_{m}(a,b)\end{pmatrix}$ (33)

by the chain rule,
$\partial_{y_{k}}f_{j}^{[1]}=\partial_{y_{k}}f_{j}+\partial_{y_{m}}f_{j}\partial_{y_{k}}h_{i}$.
Since $f_{i}(x,y^{\\{m\\}},h_{i}(x,y^{\\{m\\}}))=0$ for $(x,y^{\\{m\\}})\in
R^{1}\times\prod_{k=1}^{m-1}I_{k}^{1}$, we get

$0=\partial_{y_{k}}f_{i}(x,y^{\\{m\\}},h_{i}(x,y^{\\{m\\}}))\Big{|}_{(x,y^{\\{m\\}})=(a,b^{\\{m\\}})}=\big{[}\partial_{y_{k}}f^{[1]}_{i}-\partial_{y_{m}}f_{i}\partial_{y_{k}}h_{i}\big{]}_{(x,y)=(a,b)}.$

This identity and the change of variables of
$\partial_{y_{k}}f^{[1]}_{i}=\partial_{y_{k}}f_{i}+\partial_{y_{m}}f_{i}\partial_{y_{k}}h_{i}$
give

$\begin{split}-\partial_{y_{m}}f_{i}(a,b)\partial_{y_{k}}h_{i}(a,b^{\\{m\\}})&=\partial_{y_{k}}f_{i}(a,b)-\partial_{y_{k}}f_{i}(x,y^{\\{m\\}},h_{i}(x,y^{\\{m\\}}))\Big{|}_{(x,y^{\\{m\\}})=(a,b^{\\{m\\}})}\\\
&=\partial_{y_{k}}f_{i}(a,b).\end{split}$

So, (33) equals

$\begin{pmatrix}\partial_{y_{k}}f_{1}(a,b)\\\ \vdots\\\
\partial_{y_{k}}f_{i-1}(a,b)\\\ \partial_{y_{k}}f_{i}(a,b)\\\
\partial_{y_{k}}f_{i+1}(a,b)\\\ \vdots\\\
\partial_{y_{k}}f_{m}(a,b)\end{pmatrix}$

for $1\leq k\leq m$. Hence,

$|J_{f^{[1]},y^{\\{m\\}}}(a,b^{\\{m\\}})|=\frac{(-1)^{i+m}}{\partial_{y_{m}}f_{i}(a,b)}|J_{f,y}(a,b)|\neq
0.$

Next, since $f^{[1]}$ satisfies $|J_{f^{[1]},y^{\\{m\\}}}(a,b^{\\{m\\}})|\neq
0$ with $f^{[1]}(a,b^{\\{m\\}})=\mathbf{0}$, we can take $i^{\prime}$ such
that $\partial_{y_{m-1}}f^{[1]}_{i^{\prime}}(a,b^{\\{m\\}})\neq 0$. Here,
$i^{\prime}\neq i$; therefore, put $\phi(m-1)=i^{\prime}$. As before, by
Corollary 4.2.1 there is a unique analytic function
$y_{m-1}=h_{i^{\prime}}(x,y^{\\{m,m-1\\}}):R^{(2)}\times\prod_{k=1}^{m-2}I_{k}^{(2)}\to
I_{m-1}^{(2)}$ for some rectangle
$R^{(2)}\times\prod_{k=1}^{m-1}I_{k}^{(2)}\ni(a,b^{\\{m\\}})$ in which
$f^{[1]}_{i^{\prime}}(x,y^{\\{m,m-1\\}},h_{i^{\prime}}(x,y^{\\{m,m-1\\}}))=0$.

Put
$f^{[2]}(x,y^{\\{m,m-1\\}})=f^{[1]\\{i^{\prime}\\}}(x,y^{\\{m,m-1\\}},h_{i^{\prime}}(x,y^{\\{m,m-1\\}}))$,
where $f^{[1]\\{i^{\prime}\\}}$ indicates the $i^{\prime}$th component removal
from $f^{[1]}$. By the exact same method above,

$|J_{f^{[2]},y^{\\{m,m-1\\}}}(a,b^{\\{m,m-1\\}})|\neq 0.$

In descending induction, $\phi$ is well defined as a permutation on
$\\{1,2,\ldots,m\\}$ and take

$R=\bigcap_{k=1}^{m}R^{(k)}\quad\mbox{and}\quad
I=\bigcap_{k=1}^{m}I_{1}^{(k)}\times\bigcap_{k=1}^{m-1}I_{2}^{(k)}\times\cdots\times\bigcap_{k=1}^{2}I_{m-1}^{(k)}\times
I_{m}^{(1)}$

which contain $a$ and $b$, respectively. Finally, we will prove that $f$
satisfies Assumption 2. Since $y_{i}=h_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}})$
is a unique analytic function such that

$\begin{split}f_{\phi(i)}(x,&y^{\\{m,m-1,\ldots,i\\}},h_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}}),\ldots,h_{\phi(m)}(x,y^{\\{m\\}}))\\\
&=f^{[m-i]}_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}},h_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}}))=0\end{split}$

in $R\times I^{\\{m,m-1,\ldots,i\\}}$, for every
$(x,y^{\\{m,m-1\ldots,i\\}})\in R\times I^{\\{m,m-1,\ldots,i\\}}$ there is

$(y_{i+1},\ldots,y_{m})=(h_{\phi(i+1)}(x,y^{\\{m,m-1,\ldots,i+1\\}}),\ldots,h_{\phi(m)}(x,y^{\\{m\\}}))\in
I^{\\{i,i-1,\ldots,1\\}}$

so that

$\operatorname{sgn}_{y_{i}}f_{\phi(i)}(x,y^{\\{m,m-1,\ldots,i\\}},y_{i},y_{i+1},\ldots,y_{m})$

has only one jump discontinuity on $I_{i}^{(m-i+1)}$ for $i=m,m-1,\ldots,1$.
Therefore, the proof is complete. ∎

The following two examples of implicit functions having three variables and
two equations will serve to illustrate Theorem 5.2 and Corollary 5.2.1. In the
first example, an implicit function will be approximated when the Jacobian
matrix is non-degenerate, while the second is an example which has a
degenerate Jacobian matrix.

######  Example 5.3.

Let $f=(f_{1},f_{2}):\mathbb{R}^{3}\to\mathbb{R}^{2}$ be defined by
$f_{1}(x,y,z)=x+y^{2}+z^{3}-6$, $f_{2}(x,y,z)=x^{3}y-z-1$ with
$f(1,2,1)=\mathbf{0}$. We focus on calculating $g(x)=(g_{1}(x),g_{2}(x))$ with
$y=g_{1}(x)$ and $z=g_{2}(x)$ such that $f=\mathbf{0}$ in an interval of
$x=1$. Since $f_{1}$ and $f_{2}$ are analytic near $(x,y)=(1,2)$ and
$|J_{f,(y,z)}(1,2,1)|=-7\neq 0$, by Corollary 5.2.1 there is a rectangle, for
example, $R=[0.5,1.5]\ni x=1$ and $I=[1.5,2.5]\times[-2,8]\ni(y,z)=(2,1)$ on
which both $\operatorname{sgn}_{z}f_{1}(x,y,z)$ and
$\operatorname{sgn}_{z}f_{2}(x,y,z)$ have only one jump discontinuity. (The
surfaces $f_{1}=0$ and $f_{2}=0$ are illustrated in Figure 7, in which
$R\times I$ is also depicted.)

First, by observation, $n_{z}(f_{2})=-1$. For $N=4$,
$z=h_{2}(x,y)=\sum_{0\leq\alpha<4}c_{\alpha}(x-1)^{\alpha_{1}}(y-2)^{\alpha_{2}}$
approximates $z=z(x,y)$ such that $f_{2}(x,y,z(x,y))=0$, where $c_{\alpha}$ is
the $(\alpha+1)$th component of

$\begin{matrix}&\\\ \left(\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y-2$}&\scalebox{0.8}{$(y-2)^{2}$}&\scalebox{0.8}{$(y-2)^{3}$}&\\\
\hline\cr\\\\[-11.38109pt]
\scalebox{0.8}{$0.999999999999997$}&\scalebox{0.8}{$0.999999999999976$}&\scalebox{0.8}{$-0.000000000000017$}&\scalebox{0.8}{$0.000000000000259$}\\\
\scalebox{0.8}{$6.000000000000016$}&\scalebox{0.8}{$3.000000000000230$}&\scalebox{0.8}{$0.000000000000490$}&\scalebox{0.8}{$-0.000000000002899$}\\\
\scalebox{0.8}{$6.000000000000000$}&\scalebox{0.8}{$3.000000000000279$}&\scalebox{0.8}{$0.000000000000160$}&\scalebox{0.8}{$-0.000000000002686$}\\\
\scalebox{0.8}{$1.999999999999848$}&\scalebox{0.8}{$0.999999999998846$}&\scalebox{0.8}{$-0.000000000003098$}&\scalebox{0.8}{$0.000000000015994$}\end{matrix}\begin{matrix}&\\\
\left.\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x-1$}\\\ \scalebox{0.8}{$(x-1)^{2}$}\\\
\scalebox{0.8}{$(x-1)^{3}$}\end{matrix}\end{matrix}$ (34)

($z=h_{2}(x,y)$ is drawn in Figure 8). Furthermore, the coefficients of
$z=z(x,y)$ (which is a polynomial in two variables) are given by

$\begin{matrix}&\\\ \left(\kern-4.0pt\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.8}{$1$}&\scalebox{0.8}{$y-2$}&\scalebox{0.8}{$(y-2)^{2}$}&\scalebox{0.8}{$(y-2)^{3}$}&\\\
\hline\cr\\\\[-11.38109pt]
\scalebox{0.8}{$1$}&\scalebox{0.8}{$1$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$6$}&\scalebox{0.8}{$3$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$6$}&\scalebox{0.8}{$3$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\\\
\scalebox{0.8}{$2$}&\scalebox{0.8}{$1$}&\scalebox{0.8}{$0$}&\scalebox{0.8}{$0$}\end{matrix}\begin{matrix}&\\\
\left.\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.8}{$1$}\\\
\scalebox{0.8}{$x-1$}\\\ \scalebox{0.8}{$(x-1)^{2}$}\\\
\scalebox{0.8}{$(x-1)^{3}$}\end{matrix}\end{matrix}\kern-6.0pt$

which is comparable to (34).

According to (30), put $f_{1}^{[1]}(x,y)=f_{1}(x,y,h_{2}(x,y))=0$ (as shown in
Figure 9). Since $\operatorname{sgn}_{y}f_{1}^{[1]}(x,y)$ has only one jump
discontinuity and $n_{y}(f_{1}^{[1]})=1$, with $N=25$ we gain
$y=h_{1}(x)=\sum_{k=0}^{24}c_{k}(x-1)^{k}$ in $[-0.5,2]\times[0,5]$, which is
an approximation of $y=y(x)$ such that $f_{1}^{[1]}(x,y(x))=0$, where

$h_{1}(x)=\left(\begin{array}[]{ccccccccccccccccccccccccc}\scalebox{0.8}{$\rotatebox{90.0}{$2.0004078854$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-2.7149208367$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-4.9633925567$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$19.6874353356$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$25.2217132975$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-235.2544997108$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-40.2428746167$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$1987.0933621889$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-94.7436991654$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-11429.1583398359$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-1087.1586128315$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$44047.6659489883$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$15624.2539582495$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-110502.7914796270$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-71424.8415543300$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$170165.9282511103$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$167091.0376411914$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-138031.7528541524$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-213411.3308704084$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$22951.9362353907$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$137736.7774057497$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$43043.3433440683$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-30478.2871965461$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-22061.1876771915$}$}&\scalebox{0.8}{$\rotatebox{90.0}{$-3982.2743839354$}$}\end{array}\right)_{\kern-4.0pt\scalebox{0.6}{$1\\!\times\\!25$}}\kern-4.0pt\begin{pmatrix}1\\\
x-1\\\ (x-1)^{2}\\\ \vdots\\\
(x-1)^{24}\end{pmatrix}_{\kern-3.0pt\scalebox{0.6}{$25\\!\times\\!1$}}$

(shown in Figure 10). Finally,

$\begin{split}y=g_{1}(x)&=h_{1}(x)\\\
z=g_{2}(x)&=h_{2}(x,h_{1}(x))\end{split}$

are the approximations of $y=y(x)$ and $z=z(x)$ such that
$f(x,y(x),z(x))=\mathbf{0}$, respectively, in $R$. In addition, the quantities
$f_{1}(x,g_{1}(x),g_{2}(x))$ and $f_{2}(x,g_{1}(x),g_{2}(x))$ are illustrated
in Figure 11 and the approximated curve $x\mapsto(g_{1}(x),g_{2}(x))$ for
$f=\mathbf{0}$ is shown in Figure 12.

Figure 7: The surfaces of $f_{1}(x,y,z)=x+y^{2}+z^{3}-6=0$ and
$f_{2}(x,y,z)=x^{3}y-z-1=0$ and $R\times I$.

Figure 8: The surface
$z=h_{2}(x,y)=\sum_{0\leq\alpha<4}c_{\alpha}(x-1)^{\alpha_{1}}(y-2)^{\alpha_{2}}$.

Figure 9: A graphical plot of the implicit function of
$f_{1}(x,y,h_{2}(x,y))=0$ in $[-0.5,2]\times[0,5]$.

Figure 10: The function $y=h_{1}(x)$ is denoted in blue in $R\times[1.5,2.5]$
containing $(1,2)$, which comes from $f_{1}(x,y,h_{2}(x,y))=0$. The red
indicates a graphical plot of $f_{1}(x,y,h_{2}(x,y))=0$.

Figure 11: The blue and red curves are drawn by $f_{1}(x,g_{1}(x),g_{2}(x))$
and $f_{2}(x,g_{1}(x),g_{2}(x))$, respectively.

Figure 12: The black curve denotes $x\mapsto(y,z)$, where $y=g_{1}(x)$ and
$z=g_{2}(x)$, which approximate $y=y(x)$ and $z=z(x)$ such that
$f_{1}(x,y,z)=0$ and $f_{2}(x,y,z)=0$ in $R\times I$, respectively.

The rectangle $R$ in Assumption 2 may contain a degenerate critical point.
Regardless of the case, Theorem 5.2 still works. We present such an example to
show that a system of implicit functions is approximated by multivariate
polynomials.

######  Example 5.4.

Consider the pair of equations

$\begin{split}f_{1}(x,y,z)&=x^{2}+y^{2}+(x-1)z^{2},\\\
f_{2}(x,y,z)&=x^{2}+2y^{2}+(y-1)z^{2}\end{split}$ (35)

with $(0,0,0)$ at which $f_{1}$ and $f_{2}$ vanish. Since there is no linear
term in (35), at $(0,0,0)$ the rank of the Jacobian matrix is zero (as
depicted in Figure 13). However, confining a range, for example, to $0\leq
x\leq 1$, $0\leq y\leq 1$, and $0\leq z<\infty$ and taking a rectangle
$[0,0.5]\times[0,0.5]\times[0,3]$ which has a corner at $(0,0,0)$. On the
rectangle, by observation, $\operatorname{sgn}_{z}f_{k}(x,y,z)$ has only one
jump discontinuity and $n_{z}(f_{k})=-1$, $n_{x}(f_{k})=n_{y}(f_{k})=1$
$(k=1,2)$.

By Theorem 5.2 with $N=8$, $z=h_{1}(x,y)$ and $z=h_{2}(x,y)$ approximate
$z=z_{1}(x,y)$ and $z=z_{2}(x,y)$ such that $f_{1}(x,y,z_{1})=0$ and
$f_{2}(x,y,z_{2})=0$ whose coefficient matrices are

$\begin{matrix}&\\\ \left(\kern-4.0pt\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\\\ 12\\\ 12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.75}{$1$}&\scalebox{0.75}{$y$}&\scalebox{0.75}{$y^{2}$}&\scalebox{0.75}{$y^{3}$}&\scalebox{0.75}{$y^{4}$}&\scalebox{0.75}{$y^{5}$}&\scalebox{0.75}{$y^{6}$}&\scalebox{0.75}{$y^{7}$}&\\\
\hline\cr\\\\[-5.69054pt]
\scalebox{0.55}{$\rotatebox{0.0}{0.408241400667676}$}&\scalebox{0.55}{$\rotatebox{0.0}{0.816531359019761}$}&\scalebox{0.55}{$\rotatebox{0.0}{0.816117756871955}$}&\scalebox{0.55}{$\rotatebox{0.0}{-1.63446077805188}$}&\scalebox{0.55}{$\rotatebox{0.0}{2.50792658208951}$}&\scalebox{0.55}{$\rotatebox{0.0}{-1.48883805443688}$}&\scalebox{0.55}{$\rotatebox{0.0}{-8.02452447454164}$}&\scalebox{0.55}{$\rotatebox{0.0}{21.2048659921929}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{1.08868964456673}$}&\scalebox{0.55}{$\rotatebox{0.0}{-1.08903120663835}$}&\scalebox{0.55}{$\rotatebox{0.0}{2.19322111746780}$}&\scalebox{0.55}{$\rotatebox{0.0}{2.15725284820275}$}&\scalebox{0.55}{$\rotatebox{0.0}{-21.7192534727769}$}&\scalebox{0.55}{$\rotatebox{0.0}{57.4809959627986}$}&\scalebox{0.55}{$\rotatebox{0.0}{22.3087274847608}$}&\scalebox{0.55}{$\rotatebox{0.0}{-292.817460147257}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{1.63286148893869}$}&\scalebox{0.55}{$\rotatebox{0.0}{1.10162416860967}$}&\scalebox{0.55}{$\rotatebox{0.0}{-10.2243383275022}$}&\scalebox{0.55}{$\rotatebox{0.0}{25.8424726647171}$}&\scalebox{0.55}{$\rotatebox{0.0}{55.6146230887942}$}&\scalebox{0.55}{$\rotatebox{0.0}{-522.134552371370}$}&\scalebox{0.55}{$\rotatebox{0.0}{345.242287373544}$}&\scalebox{0.55}{$\rotatebox{0.0}{2375.12070190542}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{-0.243427885794297}$}&\scalebox{0.55}{$\rotatebox{0.0}{3.86133661190640}$}&\scalebox{0.55}{$\rotatebox{0.0}{18.5593250904849}$}&\scalebox{0.55}{$\rotatebox{0.0}{-161.347274796296}$}&\scalebox{0.55}{$\rotatebox{0.0}{318.132884269998}$}&\scalebox{0.55}{$\rotatebox{0.0}{1280.58538056068}$}&\scalebox{0.55}{$\rotatebox{0.0}{-5803.03015776689}$}&\scalebox{0.55}{$\rotatebox{0.0}{4194.37168512031}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{2.88860850137300}$}&\scalebox{0.55}{$\rotatebox{0.0}{-20.2239288354349}$}&\scalebox{0.55}{$\rotatebox{0.0}{78.8833088283334}$}&\scalebox{0.55}{$\rotatebox{0.0}{185.824867341472}$}&\scalebox{0.55}{$\rotatebox{0.0}{-5658.62904533275}$}&\scalebox{0.55}{$\rotatebox{0.0}{21857.4809412704}$}&\scalebox{0.55}{$\rotatebox{0.0}{61487.2615726907}$}&\scalebox{0.55}{$\rotatebox{0.0}{-321312.022583376}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{0.810126065016658}$}&\scalebox{0.55}{$\rotatebox{0.0}{45.9761921136848}$}&\scalebox{0.55}{$\rotatebox{0.0}{-519.510533394060}$}&\scalebox{0.55}{$\rotatebox{0.0}{1708.18678563788}$}&\scalebox{0.55}{$\rotatebox{0.0}{19624.6879190755}$}&\scalebox{0.55}{$\rotatebox{0.0}{-127143.126595965}$}&\scalebox{0.55}{$\rotatebox{0.0}{-153229.075397079}$}&\scalebox{0.55}{$\rotatebox{0.0}{1382061.71528484}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{-5.34208712455283}$}&\scalebox{0.55}{$\rotatebox{0.0}{38.6432668335212}$}&\scalebox{0.55}{$\rotatebox{0.0}{134.504396172128}$}&\scalebox{0.55}{$\rotatebox{0.0}{-4977.91637626210}$}&\scalebox{0.55}{$\rotatebox{0.0}{65268.6651956214}$}&\scalebox{0.55}{$\rotatebox{0.0}{-186901.326165857}$}&\scalebox{0.55}{$\rotatebox{0.0}{-1100349.48234761}$}&\scalebox{0.55}{$\rotatebox{0.0}{4454278.48681176}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{18.7678308525790}$}&\scalebox{0.55}{$\rotatebox{0.0}{-276.726744898086}$}&\scalebox{0.55}{$\rotatebox{0.0}{2817.28742085553}$}&\scalebox{0.55}{$\rotatebox{0.0}{-246.968670877888}$}&\scalebox{0.55}{$\rotatebox{0.0}{-314854.684257636}$}&\scalebox{0.55}{$\rotatebox{0.0}{1459642.19080275}$}&\scalebox{0.55}{$\rotatebox{0.0}{4208677.39377466}$}&\scalebox{0.55}{$\rotatebox{0.0}{-22904429.3420260}$}\\\
\end{matrix}\begin{matrix}&\\\\[4.04997pt] \left.\vphantom{\begin{matrix}12\\\
12\\\ 12\\\ 12\\\ 12\\\ 12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.65}{$1$}\\\
\scalebox{0.65}{$x$}\\\ \scalebox{0.65}{$x^{2}$}\\\
\scalebox{0.65}{$x^{3}$}\\\ \scalebox{0.65}{$x^{4}$}\\\
\scalebox{0.65}{$x^{5}$}\\\ \scalebox{0.65}{$x^{6}$}\\\
\scalebox{0.65}{$x^{7}$}\end{matrix}\end{matrix}$

and

$\begin{matrix}&\\\ \left(\kern-4.0pt\vphantom{\begin{matrix}12\\\ 12\\\ 12\\\
12\\\ 12\\\ 12\\\ 12\\\
12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.75}{$1$}&\scalebox{0.75}{$y$}&\scalebox{0.75}{$y^{2}$}&\scalebox{0.75}{$y^{3}$}&\scalebox{0.75}{$y^{4}$}&\scalebox{0.75}{$y^{5}$}&\scalebox{0.75}{$y^{6}$}&\scalebox{0.75}{$y^{7}$}&\\\
\hline\cr\\\\[-5.69054pt]
\scalebox{0.55}{$\rotatebox{0.0}{0.499998460565297}$}&\scalebox{0.55}{$\rotatebox{0.0}{1.66668593709208}$}&\scalebox{0.55}{$\rotatebox{0.0}{2.10817219588402}$}&\scalebox{0.55}{$\rotatebox{0.0}{-0.511194801196446}$}&\scalebox{0.55}{$\rotatebox{0.0}{6.36464990241272}$}&\scalebox{0.55}{$\rotatebox{0.0}{-7.09380895560241}$}&\scalebox{0.55}{$\rotatebox{0.0}{-4.67859770787413}$}&\scalebox{0.55}{$\rotatebox{0.0}{57.9648809959836}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{0.666653645507011}$}&\scalebox{0.55}{$\rotatebox{0.0}{-1.33330925771400}$}&\scalebox{0.55}{$\rotatebox{0.0}{2.84808687276566}$}&\scalebox{0.55}{$\rotatebox{0.0}{-1.66216527457875}$}&\scalebox{0.55}{$\rotatebox{0.0}{-18.2420025212197}$}&\scalebox{0.55}{$\rotatebox{0.0}{104.884228084066}$}&\scalebox{0.55}{$\rotatebox{0.0}{-165.532336874666}$}&\scalebox{0.55}{$\rotatebox{0.0}{-91.7951419714518}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{0.888956582839726}$}&\scalebox{0.55}{$\rotatebox{0.0}{0.591804280528688}$}&\scalebox{0.55}{$\rotatebox{0.0}{-8.91692976651670}$}&\scalebox{0.55}{$\rotatebox{0.0}{42.5967553060441}$}&\scalebox{0.55}{$\rotatebox{0.0}{-110.442206417017}$}&\scalebox{0.55}{$\rotatebox{0.0}{-112.581487186787}$}&\scalebox{0.55}{$\rotatebox{0.0}{2237.66959943601}$}&\scalebox{0.55}{$\rotatebox{0.0}{-5202.74439326453}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{-1.18333322661317}$}&\scalebox{0.55}{$\rotatebox{0.0}{5.53432073339570}$}&\scalebox{0.55}{$\rotatebox{0.0}{-7.65765640294483}$}&\scalebox{0.55}{$\rotatebox{0.0}{-83.2590332723340}$}&\scalebox{0.55}{$\rotatebox{0.0}{1050.98717515900}$}&\scalebox{0.55}{$\rotatebox{0.0}{-3369.24411323074}$}&\scalebox{0.55}{$\rotatebox{0.0}{-10968.6904511511}$}&\scalebox{0.55}{$\rotatebox{0.0}{52113.7334662549}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{0.755515127011205}$}&\scalebox{0.55}{$\rotatebox{0.0}{-14.1281347165804}$}&\scalebox{0.55}{$\rotatebox{0.0}{117.141832651004}$}&\scalebox{0.55}{$\rotatebox{0.0}{-446.923762231644}$}&\scalebox{0.55}{$\rotatebox{0.0}{-3051.84753534045}$}&\scalebox{0.55}{$\rotatebox{0.0}{25541.2822154715}$}&\scalebox{0.55}{$\rotatebox{0.0}{5367.23753030004}$}&\scalebox{0.55}{$\rotatebox{0.0}{-208929.001410320}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{1.15562806583154}$}&\scalebox{0.55}{$\rotatebox{0.0}{7.31554371402922}$}&\scalebox{0.55}{$\rotatebox{0.0}{-250.790503985844}$}&\scalebox{0.55}{$\rotatebox{0.0}{2273.89071435435}$}&\scalebox{0.55}{$\rotatebox{0.0}{-5099.53557921255}$}&\scalebox{0.55}{$\rotatebox{0.0}{-37753.2238770695}$}&\scalebox{0.55}{$\rotatebox{0.0}{199854.775802096}$}&\scalebox{0.55}{$\rotatebox{0.0}{-229837.302557044}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{-3.42193556701504}$}&\scalebox{0.55}{$\rotatebox{0.0}{45.9408234577807}$}&\scalebox{0.55}{$\rotatebox{0.0}{-642.740625181014}$}&\scalebox{0.55}{$\rotatebox{0.0}{1883.35634160718}$}&\scalebox{0.55}{$\rotatebox{0.0}{64000.0674666749}$}&\scalebox{0.55}{$\rotatebox{0.0}{-360910.140011759}$}&\scalebox{0.55}{$\rotatebox{0.0}{-774407.965648020}$}&\scalebox{0.55}{$\rotatebox{0.0}{5036249.52507314}$}\\\
\scalebox{0.55}{$\rotatebox{0.0}{2.04331521739714}$}&\scalebox{0.55}{$\rotatebox{0.0}{-95.5669077741437}$}&\scalebox{0.55}{$\rotatebox{0.0}{2664.22271348121}$}&\scalebox{0.55}{$\rotatebox{0.0}{-17979.3545608622}$}&\scalebox{0.55}{$\rotatebox{0.0}{-131490.726833923}$}&\scalebox{0.55}{$\rotatebox{0.0}{1166833.18379635}$}&\scalebox{0.55}{$\rotatebox{0.0}{782649.427964863}$}&\scalebox{0.55}{$\rotatebox{0.0}{-11779204.4616064}$}\\\
\end{matrix}\begin{matrix}&\\\\[4.04997pt] \left.\vphantom{\begin{matrix}12\\\
12\\\ 12\\\ 12\\\ 12\\\ 12\\\ 12\\\
12\end{matrix}}\right)\begin{matrix}\scalebox{0.65}{$1$}\\\
\scalebox{0.65}{$x$}\\\ \scalebox{0.65}{$x^{2}$}\\\
\scalebox{0.65}{$x^{3}$}\\\ \scalebox{0.65}{$x^{4}$}\\\
\scalebox{0.65}{$x^{5}$}\\\ \scalebox{0.65}{$x^{6}$}\\\
\scalebox{0.65}{$x^{7}$}\end{matrix}\end{matrix},$

respectively. Put $h=h_{1}-h_{2}$ and then $n_{y}(h)=-1$ in
$[0,0.5]\times[0,0.5]$. With $N=8$, the coefficient matrix of $y=h(x)$ is
calculated as

$\begin{matrix}&\\\
\left(\kern-4.0pt\vphantom{\begin{matrix}12\end{matrix}}\right.\end{matrix}\begin{matrix}\scalebox{0.75}{$1$}&\scalebox{0.75}{$x$}&\scalebox{0.75}{$x^{2}$}&\scalebox{0.75}{$x^{3}$}&\scalebox{0.75}{$x^{4}$}&\scalebox{0.75}{$x^{5}$}&\scalebox{0.75}{$x^{6}$}&\scalebox{0.75}{$x^{7}$}&\\\
\hline\cr\\\\[-5.69054pt]
\scalebox{0.58}{$\rotatebox{0.0}{0.116319702320269}$}&\scalebox{0.58}{$\rotatebox{0.0}{0.713165167944210}$}&\scalebox{0.58}{$\rotatebox{0.0}{0.843601400792723}$}&\scalebox{0.58}{$\rotatebox{0.0}{-0.552538829658902}$}&\scalebox{0.58}{$\rotatebox{0.0}{-3.00696939156087}$}&\scalebox{0.58}{$\rotatebox{0.0}{19.9372907280740}$}&\scalebox{0.58}{$\rotatebox{0.0}{70.0568910222166}$}&\scalebox{0.58}{$\rotatebox{0.0}{-382.742430825347}$}\end{matrix}\begin{matrix}&\\\
\left.\vphantom{\begin{matrix}12\end{matrix}}\right).\end{matrix}$

In Figure 14, the values of $f_{1}(x,h(x),h_{2}(x,h(x))$ and
$f_{2}(x,h(x),h_{2}(x,h(x))$ are illustrated, where $y=h(x)$ and
$z=h_{2}(x,h(x))$ approximate $y=y(x)$ and $z=z(x)$ such that $f_{1}(x,y,z)=0$
and $f_{2}(x,y,z)=\mathbf{0}$, respectively. Furthermore, the approximated
curve $x\mapsto(h(x),h_{2}(x,h(x))$ is shown in Figure 15.

Figure 13: Surfaces of $f_{1}(x,y,z)=x^{2}+y^{2}+(x-1)z^{2}=0$,
$f_{2}(x,y,z)=x^{2}+2y^{2}+(y-1)z^{2}=0$ with $f(0,0,0)=(0,0)$ and a rectangle
$[0,0.5]\times[0,0.5]\times[0,3]$ which contains a corner $(0,0,0)$. The
adjacent netted surfaces illustrate the implicit function $f_{2}(x,y,z)=0$.

Figure 14: For $y=h(x)$ and $z=h_{2}(x,h(x))$, the values of
$f_{1}(x,h(x),h_{2}(x,h(x)))$ and $f_{2}(x,h(x),h_{2}(x,h(x)))$ are
illustrated.

Figure 15: In $(0,0,0)\in[0,0.5]\times[0,0.5]\times[0,3]$, two surfaces of
$f_{1}(x,y,z)=x^{2}+y^{2}+(x-1)z^{2}=0$ and
$f_{1}(x,y,z)=x^{2}+2y^{2}+(y-1)z^{2}=0$, and the curve of
$(h(x),h_{2}(x,h(x))$ are depicted.

######  Remark 5.1.

  * $(i)$

Throughout this article the partitions can be chosen differently on each axis.

  * $(ii)$

We close this article with the application of our formulation to a continuous
implicit function. For a continuous function
$f(x,y):\mathbb{R}^{n+1}\to\mathbb{R}$, suppose that there is an $R\times I$
and a unique continuous function $y=g(x):R\to I$ such that $f=0$. We do not
know whether the sequence of multivariate polynomials which are calculated in
Section 4, converges to the continuous implicit function. Their coefficient
matrices depend on the choice of partitions of $R$. Because of this, by
choosing the partitions suitably, the desired convergence can be expected.
More precisely, by using dyadic decomposition of partitions, the derived
multivariate polynomials can converge to the continuous implicit function in
weak-star topology. Indeed, put $\Gamma_{1}=\\{R\\}$ and let $\Gamma_{2}$ be a
collection of grid blocks $R_{2,\alpha}$ of $R$, which are caused by dyadic
decomposition on $R$. For a positive integer $N$, inductively, let
$\Gamma_{N+1}$ consist of all grid blocks $R_{2^{N+1},\alpha}$ of every
element of $\Gamma_{N}$ by dyadic decomposition. Now, construct a multivariate
polynomial $\tilde{g}_{N}$ whose coefficients are calculate by solving (25)
over $\Gamma_{N}$ according to (15). Let
$R_{2^{N^{\prime}},\alpha}\in\bigcup_{N}\Gamma_{N}$. For any $N\geq
N^{\prime}$, there are non-overlapping grid blocks in $\Gamma_{N}$ whose union
equals $R_{2^{N^{\prime}},\alpha}$, i.e.,
$R_{2^{N^{\prime}},\alpha}=\bigcup_{R_{2^{N},\beta}\subset
R_{2^{N^{\prime}},\alpha}}R_{2^{N},\beta}$. Furthermore, the grid points of
$\Gamma_{N^{\prime}}$ are contained in those of $\Gamma_{N}$. By (25) and
(14), we have

$\int_{R_{2^{N^{\prime}},\alpha}}\tilde{g}_{N}-g\,dx=\sum_{R_{2^{N},\beta}\subset
R_{2^{N^{\prime}},\alpha}}\int_{R_{2^{N},\beta}}\tilde{g}_{N}-g\,dx=0$ (36)

for every $N\geq N^{\prime}$. On the other hand, since the collection of all
finite linear combinations of characteristic functions which are supported in
dyadic grid blocks in $\bigcup_{N}\Gamma_{N}$ is dense in $L^{1}(R)$, the
multivariate polynomial $\tilde{g}_{N}$ converges to $g$ weak-star in
$L^{\infty}(R)$ as $N\to\infty$ by the Riesz representation theorem for the
Lebesgue spaces.

## References

  * [1] U. Dini. Lezioni di analisi infinitesimale, Vol 1, Pisa, 197–241, 1907.
  * [2] A. L. Dontchev and R. T. Rockafellar. Implicit functions and solution mappings: A view from variational analysis, 2nd edition, Springer, 2014.
  * [3] G. Folland. Real analysis: Modern techniques and their applications, 2nd edition, Wiley, 2007.
  * [4] S. G. Krantz and H. R. Parks. The implicit function theorems: History, theory, and applications, Birkhaüser, 2013.
  * [5] S. G. Krantz and H. R. Parks. A Primer of Real Analytic Functions, the 2nd edition, Birkhaüser, 2002.
  * [6] R. J. Mathar. Improved first estimates to the solution of Kepler’s Equation, arXiv:2108.03215, 2021.
  * [7] J. K. Moser. A new technique for the construction of solutions of nonlinear differential equations, Proceedings of the National Academy of Sciences of the United States of America 47, 1824–1831, 1961.
  * [8] J .F. Nash. The imbedding problem for Riemannian manifolds, Annals of Mathematics 63, 20–63, 1956.
  * [9] I. Newton. Mathematical papers of Isaac Newton, vol 2, edited by D.T. Whiteside, Cambridge University Press, 1968.
  * [10] V. Raposo-Pulido and J. Peláez. An efficient code to solve the Kepler equation. Elliptic case, Monthly Notices Royal Astro. Soc. 467, 1702–1713. 2017.
  * [11] D. J. Struik. A source book in mathematics, 1200–1800, Harvard University Press, Cambridge, Massachusetts, 1969.
  * [12] W. Rudin. Functional analysis, the 2nd edition, McGraw-Hill Science, 1991.