<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\documentclass[12pt,reqno]{article}

\usepackage[usenames]{color}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{amscd}
\usepackage{graphicx}

\usepackage[colorlinks=true,
linkcolor=webgreen,
filecolor=webbrown,
citecolor=webgreen]{hyperref}

\definecolor{webgreen}{rgb}{0,.5,0}
\definecolor{webbrown}{rgb}{.6,0,0}

\usepackage{color}
\usepackage{fullpage}
\usepackage{float}

\usepackage{psfig}
\usepackage{graphics}
\usepackage{latexsym}
\usepackage{epsf}
\usepackage{breakurl}

\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{.1in}
\setlength{\evensidemargin}{.1in}
\setlength{\topmargin}{-.1in}
\setlength{\textheight}{8.4in}

\newcommand{\seqnum}[1]{\href{https://oeis.org/#1}{\underline{#1}}}

\begin{document}

\begin{center}
\epsfxsize=4in
\leavevmode\epsffile{logo129.eps}
\end{center}

\theoremstyle{plain}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{conjecture}[theorem]{Conjecture}

\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}

\newcommand{\gseq}{\left ( \gamma^{n}_{k}\right )_{n \geq 0,~ 0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}}
\newcommand{\useq}{\left ( u^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}}
\newcommand{\ds}{\displaystyle}

\begin{center}
\vskip 1cm{\LARGE\bf
When Numerical Analysis Crosses Paths with \\
\vskip .1in
Catalan and Generalized Motzkin Numbers}
\vskip 1cm \large
Paul Eloe and Catherine Kublik \\
Department of Mathematics\\
University of Dayton\\
300 College Park\\
Dayton, OH 45469\\
USA \\
\href{mailto:peloe1@udayton.edu}{\tt peloe1@udayton.edu}\\
\href{mailto:ckublik1@udayton.edu}{\tt ckublik1@udayton.edu}\\
\end{center}

\vskip .2 in

\begin{abstract}
We study a linear doubly indexed sequence that contains the Catalan numbers and relates to a class of generalized Motzkin numbers.  We obtain a closed form formula, a generating function and a nonlinear recursion relation for this sequence. We show that a finite difference scheme with compact stencil applied to a nonlinear differential operator acting on the Euclidean distance function is exact, and exploit this exactness to produce the nonlinear recursion relation. In particular, the nonlinear recurrence relation is obtained by using standard error analysis techniques from numerical analysis. This work shows a connection between numerical analysis and number theory, and illustrates an interesting occurrence of the Catalan and generalized Motzkin numbers in a context a priori void of combinatorial objects.
 \end{abstract}

\section{Introduction}

In this paper, we study the doubly indexed sequence $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ defined by the recursion relation and initial conditions
\begin{equation} \label{gamma_linear_recursion}
\begin{cases}
\gamma^{n}_{0} = 1, &amp; \forall ~ n \geq 0; \\
\gamma^{2n-1}_{n} = 0, &amp; \forall ~ n \geq 1; \\
\gamma^{n}_{k} = \left ( 1 + \frac{2k}{n+2} \right ) \gamma^{n-1}_{k} + \frac{2k-n-1}{n+2} \gamma^{n-1}_{k-1}, &amp; \forall~ n \geq 2, 1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor.
\end{cases}
\end{equation}
We observe that the recursion relation is linear in $\gamma^{n}_{k}$ and that $\gseq$ is uniquely determined by \eqref{gamma_linear_recursion}.

In this paper, we show that the sequence $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ is related to the well-known Catalan and Motzkin numbers via its closed form formula, and show that $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ also satisfies a nonlinear recursion relation. We obtain the nonlinear recursion relation by using a finite difference scheme which is exact when applied to a specific differential operator acting on the Euclidean distance function. 


In the western world, the first appearance of the Catalan numbers dates back to Euler around 1751. At that time Euler was considering the problem of counting triangulations of convex $n+2$-gons. With the help of Segner and Goldbach, Euler showed that for $n \geq 0$, the number of triangulations of a convex $n+2$-gon is
$$
C_n = \frac{(2n)!}{n!(n+1)!} = \frac{1}{n+1} \binom{2n}{n}.
$$
A major element in the proof was the recursion relation (also known as Segner's relation)
\begin{equation} \label{recursion_segner}
C_{n+1} = \sum_{k=0}^{n} C_{k} C_{n-k},~~~~ n \geq 0,
\end{equation}
satisfied by $C_n$, the number of triangulations of an $n+2$-gon. Segner proved the relation via a combinatorial argument.
They also discovered the generating function for the Catalan numbers
$$
A(x) = \sum_{n=0}^{\infty} C_n x^n = \frac{1 - \sqrt{1 - 4x}}{x},~~~~ |x| &lt; \frac{1}{4}.
$$
The Catalan number $C_n$ is named after Catalan, a French and Belgian mathematician who studied this sequence and was the first to obtain what are now standard formulas
$$
C_n = \frac{(2n)!}{n!(n+1)!} = \binom{2n}{n} - \binom{2n}{n-1}, ~~~ n \geq 0.
$$
Many mathematicians have studied the Catalan numbers and there is now a very large literature on that topic. For example, in addition to the triangulation problem of Euler, there are now more than two hundred known combinatorial problems that are counted by the Catalan numbers. See Stanley \cite{stanley15} for this enumeration. There is also a very large literature on generalizations of the Catalan numbers. Some of the well-known generalizations of the Catalan numbers include the Fuss-Catalan numbers \cite{fuss1791, heubach_li_mansour08, hilton_pedersen91, kim08, pak99} and their generalizations \cite{armstrong_rhoades_williams13, gould64}. We finish by a remark. In 1730, Ming Antu, a Chinese scientist and mathematician wrote a book that included a number of trigonometric identities, some which involved the Catalan numbers. This was observed by J. J. Luo \cite{larcombe99, luo88} in 1988. This is the earliest known appearance of what we now call Catalan numbers.

The Motzkin numbers are connected to the Catalan numbers by
 $$C_{n+1}=\sum_{k=0}^{n}\binom{n}{k}M_{k}, ~~~~~ n \geq 0,$$
 and are defined as
 \begin{equation} \label{motzkin_number}
 M_{n}=\sum_{k=0}^{\lfloor \frac{n}{2}\rfloor}\binom{n}{2k}C_{k},
 \end{equation}
 for $n \geq 0$. The Motzkin numbers first appeared in the seminal article by Motzkin \cite{motzkin48} in which Motzkin enumerated the number of partitions of $2n$ points on a circle into $n$ pairs without crossings. Later, Donaghey and Shapiro \cite{donaghey_shapiro77} presented fourteen combinatorial objects enumerated by the Motzkin numbers. We also refer the reader to Aigner \cite{aigner98} for some of the properties of the Motzkin numbers and Bernhart \cite{bernhart99} for their relation with the Catalan numbers. 
 In related papers, Dattoli and co-authors \cite{artioli_dattoli_licciardi_pagnutti17, blasiak_dattoli_horzela_penson_zhukovsky08} have employed hybrid polynomials to study Motzkin numbers and their associated forms. In 2014, Sun \cite{sun14} introduced a generalized Motzkin number $M_n(a,b)$ depending on two additional integer parameters $a$ and $b$, defined as
\begin{equation} \label{gen_motzkin}
M_n(a,b) = \sum_{k=0}^{\left \lfloor \frac{n}{2} \right \rfloor} \binom{n}{2k} C_k b^k a^{n-2k},
\end{equation}
for $n \geq 0$ and $a,b \in \mathbb Z$. In particular, Sun showed that $M_n(a,b)$ satisfies the linear recursion relation
\begin{equation} \label{gen_motzkin_recurrence}
\begin{cases}
M_0(a,b) = 1,  &amp; \\
M_1(a,b) = a, &amp; \\
(n+3) M_{n+1}(a,b) = b(2n+3) M_n(a,b) - (a^2-4b) n M_{n-1}(a,b), &amp;  n \geq 1,
\end{cases}
\end{equation}
which enabled him to obtain the generating function
\begin{equation} \label{motzkin_genfn}
\mathcal{M}_{a,b}(x) = \sum_{n=0}^{\infty} M_{n}(a,b) x^n = \frac{1-ax-\sqrt{(1-ax)^2 - 4bx^2}}{2bx^2}.
\end{equation}
We refer the reader to Wang et al.~\cite{wangzhang15} and  Zhao et al.~\cite{zhaoqi17} for some of the properties of these generalized Motztkin numbers. 

Throughout the literature, the development of the Catalan numbers, the Motzkin numbers and their generalizations were primarily motivated by the process of enumerating combinatorial objects.  The motivation in this paper is different, to our knowledge new, and represents the primary contribution of this manuscript. We apply a finite difference scheme with compact stencil to a nonlinear differential operator acting on a function of two variables.  When the function of two variables is the Euclidean distance function, the application of the finite difference algorithm is exact and the sequence \eqref{gamma_linear_recursion} is generated in a corresponding error analysis.  To our knowledge, the observation of Catalan and Motzkin-like numbers in numerical analysis is new.

The paper is organized as follows. In Section~\ref{preliminary} we give results on rising polynomials and Catalan numbers that are used later on in the paper. In Section~\ref{distfn_gamma} we establish the relationship between the two-dimensional Euclidean distance function and the sequence $\gseq$. We obtain the closed form formula (exhibiting the Catalan numbers) for the sequence $\gseq$ in Section~\ref{closedform} and derive its nonlinear recursion relation in Section~\ref{nlr_section}. We finish in Section~\ref{generatingfn} by presenting the generating function for $\gseq$, and conclude in Section~\ref{conclusion}.

\section{Preliminaries} \label{preliminary}

In this section we discuss two results that are needed to preserve the completeness of the paper.  The first result is the binomial theorem for special polynomials, called rising polynomials or Pochhammer symbols. The second result is the recursion relation for the Catalan numbers.

\begin{definition}[Rising polynomial]
For $t, n \in \mathbb R$, we define the rising polynomial $t^{\overline{n}}$ as $$t^{\overline{n}}=\frac{\Gamma (t+n)}{\Gamma (t)},$$ with the convention that if $t$ is a nonpositive integer then $t^{\overline{n}}=0$, and if $n=0$ then  $t^{\overline{0}}=1$. If $n$ is a positive integer then $$t^{\overline{n}}=t(t+1) \cdots (t+n-1).$$
  \end{definition}
We note that $t^{\overline{n}}$ is also known as the Pochhammer symbol, denoted $(t)_{n}$. Also, just like the regular polynomial $t^n$, the rising polynomial $t^{\overline{n}}$ satisfies the binomial theorem.

\begin{theorem}[Binomial theorem for rising polynomials] \label{binomthm_risingpoly}
Assume $n$ denotes a nonnegative integer and let $a,b$ be real.  Then
\begin{equation}\label{binomial theorem}
(a+b)^{\overline{n}}=\sum_{k=0}^{n}\binom{n}{k}a^{\overline{k}}b^{\overline{n-k}}.
\end{equation}
\end{theorem}

\begin{proof}
The proof is done by induction.

\begin{enumerate}

\item Consider the base case $n=0$. Then $(a+b)^{\overline{0}} = 1$ and $\ds \sum_{k=0}^{0}\binom{n}{k}a^{\overline{k}}b^{\overline{n-k}} = \binom{0}{0}a^{\overline{0}}b^{\overline{0}} =1$.


\item Assume $n$ is a nonnegative integer and assume \eqref{binomial theorem} is true.  Then
\begin{align*}
(a+b)^{\overline{n+1}}&amp;=(a+b+n)(a+b)^{\overline{n}}=(a+b+n)\sum_{k=0}^{n}\binom{n}{k}a^{\overline{k}}b^{\overline{n-k}}\\
&amp;=\sum_{k=0}^{n}\binom{n}{k}(a+k)a^{\overline{k}}b^{\overline{n-k}}+\sum_{k=0}^{n}\binom{n}{k}a^{\overline{k}}(b+n-k)b^{\overline{n-k}}\\
&amp;=\sum_{k=0}^{n}\binom{n}{k}a^{\overline{k+1}}b^{\overline{n-k}}+\sum_{k=0}^{n}\binom{n}{k}a^{\overline{k}}b^{\overline{n+1-k}}\\
&amp;=a^{\overline{n+1}}+\sum_{k=0}^{n-1}\binom{n}{k}a^{\overline{k+1}}b^{\overline{n-k}}+\sum_{k=1}^{n}\binom{n}{k}a^{\overline{k}}b^{\overline{n+1-k}}+b^{\overline{n+1}}\\
&amp;=a^{\overline{n+1}}+\sum_{k=1}^{n}\binom{n}{k-1}a^{\overline{k}}b^{\overline{n+1-k}}+\sum_{k=1}^{n}\binom{n}{k}a^{\overline{k}}b^{\overline{n+1-k}}+b^{\overline{n+1}}\\
&amp;=a^{\overline{n+1}}+\sum_{k=1}^{n}\left (\binom{n}{k-1}+\binom{n}{k} \right )a^{\overline{k}}b^{\overline{n+1-k}}+b^{\overline{n+1}}\\
&amp;=\sum_{k=0}^{n+1}\binom{n+1}{k}a^{\overline{k}}b^{\overline{n+1-k}}.
\end{align*}
Thus \eqref{binomial theorem} is true for $n+1$ and therefore for all $n \geq 0$.
\end{enumerate}
\end{proof}

\begin{theorem} \label{segner_relation}
For $n \geq 0$, the $n$-th Catalan number is defined as
$$
C_n = \frac{(2n)!}{n! (n+1)!}.
$$
Then $\left ( C_n \right )_{n \geq 0}$ can be defined recursively as
\begin{equation}\label{catalan}
\begin{cases}
C_0 = 1, &amp; \\
C_{n+1} = \sum\limits_{k=0}^{n} C_k C_{n-k}, &amp; n \geq 0.
\end{cases}
\end{equation}
\end{theorem}
\begin{proof}
We verify \eqref{catalan} by using the generating function for the Catalan numbers.
Assume $(v_{n})_{n\geq 0}$ denotes a sequence of real numbers satisfying~\eqref{segner_relation}, namely
\begin{equation}\label{Segner}
v_{0}=1, \quad v_{n+1}=\sum_{k=0}^{n}v_{k}v_{n-k}.
\end{equation}
Note that $(v_{n})_{n\geq0}$ is uniquely determined by the recursion \eqref{Segner}.

  Let
$$p(x)=\sum_{n=0}^{\infty}v_{n}x^{n}$$ be the generating function for the sequence  $(v_{n})_{n\geq 0}$.  Employ \eqref{Segner} to obtain
$$xp^{2}(x)=x \left (\sum_{n=0}^{\infty}v_{n}x^{n} \right )^{2}=x\sum_{n=0}^{\infty}\left ( \sum_{k=0}^{n}v_{k}v_{n-k}\right ) x^{n}=\sum_{n=0}^{\infty}v_{n+1}x^{n+1}=p(x)-1.$$
Thus $p$ satisfies the quadratic equation
$$
xp^2(x) - p(x) + 1 = 0.
$$
In particular,
$$p(x)=\frac{1-\sqrt{1-4x}}{2x}=\sum_{n=0}^{\infty}\frac{(2n)!}{n!(n+1)!}x^{n}, \quad |x|&lt; \frac{1}{4}.$$
Thus, since $\left (v_n \right )_{n\geq 0}$ is unique, we obtain that
$$
v_n = \frac{(2n)!}{n!(n+1)!} = C_n,
$$
and thus \eqref{catalan} holds.
\end{proof}



\section{Distance function and the sequence \texorpdfstring{$\gseq$}{g}}  \label{distfn_gamma}

In this section we show that the sequence $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ naturally arises in the derivatives of the two-dimensional Euclidean distance function to a fixed point. Indeed, let $(x_0,y_0) \in \mathbb R^2$. We define $d: \mathbb R^2 \rightarrow \mathbb R^+$ to be the Euclidean distance from any point $\mathbf{x} = (x,y) \in \mathbb R^2$ to the point $(x_0,y_0)$, namely
$$
d(\mathbf{x}) = d(x,y) = \sqrt{(x-x_0)^2 + (y-y_0)^2}.
$$
Without loss of generality, we take $x_0 = y_0 = 0$ and consider the distance function to the origin
\begin{equation} \label{distancefn}
d(\mathbf{x}) = d(x,y) = \sqrt{x^2 + y^2}.
\end{equation}

We define the quantities
\begin{equation}  \label{def_X}
X := x^2,
\end{equation}
\begin{equation} \label{def_Y}
Y := y^2,
\end{equation}
and
\begin{equation} \label{def_D}
D := X+Y.
\end{equation}
It follows that the distance $d$ defined in \eqref{distancefn} can be written as
\begin{equation}
d = D^{\frac{1}{2}},
\end{equation}
and that
\begin{equation} \label{deriv_D}
D_x = 2x, \mbox{ }  D_y = 2y , ~ \\
 \mbox{ and } D_{xx} =D_{yy} = 2. \\
\end{equation}

\begin{definition} \label{define_P}
For $m \in \mathbb N_{0}$, we define the homogeneous polynomials of degree $m$ in $X$ and $Y$, $P^{2m+2}_{m}$ and $P^{2m+3}_{m}$, as
$$
P^{2m+2}_{m}(X,Y) = \sum_{k=0}^{m} \gamma^{2m}_{k} Y^{k} X^{m-k} ,
$$
and
$$
P^{2m+3}_{m}(X,Y) = \sum_{k=0}^{m} \gamma_{k}^{2m+1} Y^{k} X^{m-k},
$$
where $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ is defined in \eqref{gamma_linear_recursion}.
\end{definition}
Note that the two polynomials $P^{2m+2}_{m}$ and $P^{2m+3}_{m}$ can be defined by the single equation
$$
P^{n+2}_{\left \lfloor \frac{n}{2} \right \rfloor}(X,Y) = \sum_{k=0}^{\left \lfloor \frac{n}{2} \right \rfloor } \gamma^{n}_{k} Y^kX^{\left \lfloor \frac{n}{2} \right \rfloor  -k}.
$$
for $n \geq 0$ and $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ as defined in \eqref{gamma_linear_recursion}.

The following theorem establishes the connection between the sequence $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ and the higher order derivatives of the distance function $d$ defined in \eqref{distancefn}.

\begin{theorem} \label{formula_deriv_d}
Let $d$ be the distance function defined in \eqref{distancefn}. Then for $m \geq 0$, the partial derivatives of $d$ with respect to $x$ can be written as
\begin{equation} \label{deriv_d_even}
\frac{\partial^{(2m+2)}d }{\partial x^{2m+2}} = \frac{(2m+2)!}{2} Y D^{-\frac{4m+3}{2}} P^{2m+2}_{m}(X,Y),
\end{equation}
\begin{equation} \label{deriv_d_odd}
\frac{\partial^{(2m+3)} d }{\partial x^{2m+3}}  = -\frac{(2m+3)!}{2} x Y D^{-\frac{4m+5}{2}} P^{2m+3}_{m}(X,Y),
\end{equation}
where $P^{2m+2}_{m}$ and $P^{2m+3}_{m}$ are defined in Definition~\ref{define_P} and $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ is defined in \eqref{gamma_linear_recursion}.

Similarly, the partial derivatives of $d$ with respect to $y$ can be written as
$$
\frac{\partial^{(2m+2)}d }{\partial y^{2m+2}} = \frac{(2m+2)!}{2} X D^{-\frac{4m+3}{2}} P^{2m+2}_{m}(Y,X),
$$
and
$$
\frac{\partial^{(2m+3)} d }{\partial y^{2m+3}}  = -\frac{(2m+3)!}{2} y X D^{-\frac{4m+5}{2}} P^{2m+3}_{m}(Y,X).
$$
\end{theorem}

\begin{proof} The proof is done by induction.

First, we remark that since $d$ is symmetric in $x$ and $y$ (and hence in $X$ and $Y$), the derivatives of $d$ with respect to $y$ can be obtained from the derivatives of $d$ with respect to $x$ by simply interchanging the roles of $x$ and $y$ (and hence $X$ and $Y$). Thus we only prove the formulas for the partial derivatives of $d$ with respect to $x$.
\begin{enumerate}
\item Consider the base case $m=0$.
 \begin{align*}
 d_{xx}  &amp; = \left (D^{\frac{1}{2}} \right )_{xx} = \frac{1}{2} \left (D^{-\frac{1}{2}} D_x \right )_x = \frac{1}{2} \left (-\frac{1}{2} D^{-\frac{3}{2}} \underbrace{D_x^2}_{=4X} + D^{-\frac{1}{2}} \underbrace{D_{xx}}_{=2} \right ) \\
 &amp; =  D^{-\frac{3}{2}} \left (D - X \right ) =  Y D^{-\frac{3}{2}}  = \frac{(2)!}{2} Y D^{-\frac{3}{2}} \gamma^{0}_{0}X^{1-1},
 \end{align*}
 where $\gamma^{0}_{0} = 1$.
\begin{align*}
d_{xxx} &amp; = \underbrace{Y_x}_{=0} D^{-\frac{3}{2}} -\frac{3}{2} Y D^{-\frac{5}{2}} \underbrace{D_x}_{=2x} =  -3x Y D^{-\frac{5}{2}} \\
&amp; = -\frac{3!}{2} x Y D^{-\frac{4(1)+1}{2}}\gamma^{1}_{0} X^{1-1}, \\
\end{align*}
with $ \gamma^{1}_{0} =1$.

\item For the induction hypothesis, assume \eqref{deriv_d_even} and \eqref{deriv_d_odd} both hold for $m \geq 0$. Then for $m+1\geq 0$ we have
\begin{align*}
\frac{\partial^{(2m+4)} d }{\partial x^{2m+4}} &amp; = \frac{\partial }{\partial x} \left (\frac{\partial^{(2m+3)} d }{\partial x^{2m+3}}  \right ) \\
&amp; = -\frac{(2m+3)!}{2} Y \left (  x D^{-\frac{4m+5}{2}} P^{2m+3}_{m} \right )_x \\
&amp; = -\frac{(2m+3)!}{2} Y \left (  \frac{\partial }{\partial x} \left (x D^{-\frac{4m+5}{2}} \right ) P^{2m+3}_{m} + x D^{-\frac{4m+5}{2}} \left ( P^{2m+3}_{m}\right )_x \right ). 
\end{align*}
By using Equations~\eqref{def_X}, \eqref{def_Y}, \eqref{def_D}, and \eqref{deriv_D} we can calculate the product rule above as follows:
\begin{align*}
&amp;  \frac{\partial}{\partial x} \left (x D^{-\frac{4m+5}{2}} \right ) P^{2m+3}_{m} + x D^{-\frac{4m+5}{2}} \left ( P^{2m+3}_{m}\right )_x  \\
&amp; = D^{-\frac{4m+7}{2}} \left (D - (4m+5)x^2 \right )P^{2m+3}_{m} \\
&amp; + x D^{-\frac{4m+7}{2}} D \sum\limits_{k=0}^{m} \gamma^{2m+1}_k (m-k) Y^k X^{m-k-1} X_x \\
&amp; = D^{-\frac{4m+7}{2}} \left (X+Y - (4m+5)X \right )P^{2m+3}_{m} \\
&amp; + 2X (X+Y) D^{-\frac{4m+7}{2}} \sum\limits_{k=0}^{m} \gamma^{2m+1}_k (m-k) Y^k X^{m-k-1}\\
&amp; = D^{-\frac{4m+7}{2}}  \mathcal{E}, 
\end{align*}
where
$$
\mathcal{E} =  (Y-4(m+1)X) P^{2m+3}_{m} + 2X(X+Y) \sum_{k=0}^{m} \gamma^{2m+1}_{k} (m-k) Y^k X^{m-1-k}.
$$
Note that
\begin{align*}
\mathcal{E}  &amp;  =  \sum_{k=0}^{m} \left ( 1 + 2(m-k) \right) \gamma^{2m+1}_{k} Y^{k+1} X^{m-k}   \\
&amp; +  \sum_{k=0}^{m} \left ( 2(m-k) - 4(m+1) \right )\gamma^{2m+1}_{k}  Y^k X^{m+1-k}  \\
&amp; = \sum_{k=1}^{m+1} \left ( 2m+ 3 - 2k \right )  \gamma^{2m+1}_{k-1} Y^{k} X^{m+1-k} -  2 \sum_{k=0}^{m} \left ( k+2+m  \right )\gamma^{2m+1}_{k}  Y^k X^{m+1-k} \\
&amp; = \gamma^{2m+1}_{m} Y^{m+1} - 2 (m+2) \gamma^{2m+1}_{0} X^{m+1} \\
&amp; + \sum_{k=1}^{m} \left ( \left ( 2m+ 3 - 2k \right )  \gamma^{2m+1}_{k-1} - 2 \left ( k+2+m  \right )\gamma^{2m+1}_{k}  \right )Y^k X^{m+1-k}.
\end{align*}
Thus we can write
\begin{align*}
 \frac{\partial^{(2m+4)} d }{\partial x^{2m+4}}  &amp; = - \frac{(2m+3)!}{2} Y D^{-\frac{4m+7}{2}} \mathcal{E} \\
 &amp; =  \frac{(2m+4)!}{2} Y D^{-\frac{4m+7}{2}} \frac{\mathcal{-E}}{2m+4},
\end{align*}
where
\begin{align*}
\frac{\mathcal{-E}}{2m+4} &amp; = \underbrace{\gamma^{2m+1}_{0}}_{=1} X^{m+1} + \sum_{k=1}^{m} \left ( 2 \frac{k+2+m}{2m+4} \gamma^{2m+1}_{k} - \frac{2m+ 3 - 2k}{2m+4}  \gamma^{2m+1}_{k-1}  \right )Y^k X^{m+1-k} \\
&amp; - \frac{\gamma^{2m+1}_{m}}{2m+4}Y^{m+1} \\
&amp; = X^{m+1} +  \sum_{k=1}^{m} \left ( \left ( 1 + \frac{2k}{2m+4} \right ) \gamma^{2m+1}_{k} + \frac{2k-2m-3}{2m+4}  \gamma^{2m+1}_{k-1}  \right )Y^k X^{m+1-k} \\
&amp; - \frac{\gamma^{2m+1}_{m}}{2m+4}Y^{m+1} \\
&amp; = \sum_{k=0}^{m+1} \gamma^{2m+2}_{k} Y^k X^{m+1-k},
\end{align*}
where for $m \geq 0$ we have
$$
\begin{cases}
\gamma^{2m+2}_{0} = 1, &amp; \\
\gamma^{2m+2}_{m+1} = -\frac{1}{2m+4} \gamma^{2m+1}_{m}, &amp; \\
\gamma^{2m+2}_{k} =   \left ( 1 + \frac{2k}{2m+4} \right ) \gamma^{2m+1}_{k} + \frac{2k-2m-3}{2m+4}  \gamma^{2m+1}_{k-1}, &amp; \mbox{ for } 1 \leq k \leq m.\\
\end{cases}
$$
By definition, $\gamma^{2m+1}_{m+1} = 0$ for any $m \geq 0$; thus
\begin{equation} \label{recursion_relation_even}
\gamma^{2m+2}_{k} =   \left ( 1 + \frac{2k}{2m+4} \right ) \gamma^{2m+1}_{k} + \frac{2k-2m-3}{2m+4}  \gamma^{2m+1}_{k-1}  \mbox{ for } 1 \leq k \leq m,
\end{equation}
is still true for $k=m+1$ and reduces to
$$
\gamma^{2m+2}_{m+1} = -\frac{1}{2m+4} \gamma^{2m+1}_{m}.
$$
Recognizing \eqref{recursion_relation_even} as the recursion relation for the sequence  $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ for an even index $n$, it follows that \eqref{deriv_d_even} is true for $m+1$ and thus is true for all $m \geq 0$.

Similarly, we can show Equation~\eqref{deriv_d_odd} for $m+1 \geq 0$, namely
$$
\frac{\partial^{(2m+5)} d }{\partial x^{2m+5}} = \frac{(2m+5)!}{2} Y xD^{-\frac{4m+9}{2}} \sum_{k=0}^{m+1} \gamma^{2m+3}_{k} Y^k X^{m+1-k}
$$
with for $m \geq 0$,
\begin{equation} \label{recursion_relation_odd}
\begin{cases}
\gamma^{2m+3}_{0} = 1, &amp; \\
\gamma^{2m+3}_{k} = \left ( (1 + \frac{2k}{2m+3}) \gamma^{2m+2}_{k} + \frac{2k-2m-4}{2m+3} \gamma^{2m+2}_{k-1} \right ), &amp; \mbox{for } 1 \leq k \leq m+1,
\end{cases}
\end{equation}
Recognizing the recursion relation for the sequence $\left ( \gamma^{n}_{k}\right )_{n \geq 0, ~0 \leq  k \leq \left \lfloor \frac{n}{2} \right \rfloor}$ for an odd index $n$, it follows that \eqref{deriv_d_odd} is true for $m+1$ and thus is true for all $m \geq 0$.
\end{enumerate}
\end{proof}

We will exploit the relationship between the sequence $\gseq$ and the derivatives of the distance function to prove that $\gseq$ also satisfies a nonlinear recursion relation. This will be done by employing a finite difference scheme on $\nabla \cdot \left (d \nabla d \right )$, where $d$ is the Euclidean distance function defined in \eqref{distancefn}.

\section{Closed form formula for \texorpdfstring{$\gseq$}{g}} \label{closedform}

In this section, we obtain a closed form formula for the general term of the sequence $\gseq$.

\begin{theorem} \label{closedform_thm}
Let $\gseq$ be defined by the recursion relation and initial conditions given in \eqref{gamma_linear_recursion}. Then the general term of the sequence $\gseq$ is
$$
\gamma^{n}_{k} = \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k,
$$
where $\ds C_k = \frac{(2k)!}{k! (k+1)!} = \frac{1}{k+1} \binom{2k}{k}$ for $k \geq 0$, is the $k$-th Catalan number.
\end{theorem}

\begin{remark} Observe that $n \geq 0, 0 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor$ is equivalent to $ k \geq 0, n \geq 2k$. Thus, the binomial coefficient $\ds \binom{n}{2k}$ represents the number of ways $2k$ objects (and $2k \leq n$) can be taken from a set of $n$ objects. In the case $2k&gt;n$, we define $\ds \binom{n}{2k}$ to be $0$. \end{remark}

\begin{remark} Note that the term $\ds \binom{n}{2k} C_k$ can be recognized as the summand of the Motzkin number $M_n$ defined in \eqref{motzkin_number}. \end{remark}

\begin{proof} Define $$u^{n}_{k}: = \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k, ~~~n \geq 0,  0 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor.$$ We will show that the sequence $\useq$ coincides with the sequence $\gseq$. \\

\begin{enumerate}

\item Let's first calculate $u^{n}_{0}$ for all $n \geq 0$ and $u^{2n-1}_{n}$ for all $n \geq 1$.

For $n \geq 0$, we have
$$
 u^{n}_{0} = \left ( -\frac{1}{4} \right )^0 \binom{n}{0} C_0 = 1.
$$

For $n \geq 1$, we have
$$
u^{2n-1}_{n} = \left ( -\frac{1}{4} \right )^n \underbrace{\binom{2n-1}{2n}}_{=0} C_n = 0.
$$
Thus we have
$$
\begin{cases}
u^{n}_{0} = 1, &amp; \forall ~ n \geq 0; \\
u^{2n-1}_{n} = 0, &amp; \forall ~ n \geq 1. \\
\end{cases}
$$

\item Suppose $n \geq 2$ and $1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor$. Then we have
\begin{align*}
&amp; \left ( 1 + \frac{2k}{n+2} \right ) u^{n-1}_{k} + \frac{2k-n-1}{n+2} u^{n-1}_{k-1} \\
&amp; = \left ( -\frac{1}{4} \right )^k \left ( \frac{n+2+2k}{n+2}\binom{n-1}{2k} C_k + 4 \frac{n+1-2k}{n+2} \binom{n-1}{2k-2} C_{k-1}  \right).
\end{align*}
Using the definition of the Catalan numbers we can write
$$
C_{k-1} = \frac{k+1}{2(2k-1)}C_k, ~ k \geq 1.
$$
Similarly, using the definition of the binomial coefficients, we can write
$$
\binom{n-1}{2k} = \frac{n-2k}{n} \binom{n}{2k}, ~ n \geq 2, 1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor,
$$
and
$$
\binom{n-1}{2k-2} = \frac{2k(2k-1)}{n(n+1-2k)} \binom{n}{2k}, ~n \geq 2, 1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor.
$$
Consequently,
\begin{align*}
&amp; \left ( 1 + \frac{2k}{n+2} \right ) u^{n-1}_{k} + \frac{2k-n-1}{n+2} u^{n-1}_{k-1}  \\
&amp; = \left ( -\frac{1}{4} \right )^k  \left ( \frac{n+2+2k}{n+2} \frac{n-2k}{n} \binom{n}{2k} C_k + 4 \frac{n+1-2k}{n+2} \frac{2k(2k-1)}{n(n+1-2k)} \binom{n}{2k} \frac{k+1}{2(2k-1)}C_k  \right) \\
&amp; = \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k  \underbrace{\frac{(n+2-2k)(n-2k) + 4 k(k+1)}{n(n+2)}}_{=1} \\
&amp; = \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k  = u^{n}_{k}.
\end{align*}

\begin{remark}The calculations above include the case $n=2m$ and $k=m$. Note that in this case, we have $\ds \binom{n-1}{2k} = \binom{2m-1}{2m} = 0$, but the calculations above still hold in that case since
\begin{align*}
 \left ( 1 + \frac{2k}{n+2} \right ) u^{n-1}_{k} + \frac{2k-n-1}{n+2} u^{n-1}_{k-1}  &amp;= \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k \frac{4k(k+1)}{n(n+2)}\\
&amp; =\left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k\frac{4m(m+1)}{2m(2m+2)} \\
&amp; = \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k.
\end{align*}
\end{remark}

This shows that the sequence $\useq$ satisfies
$$
\begin{cases}
u^{n}_{0} = 1, &amp; \forall ~ n \geq 0; \\
u^{2n-1}_{n} = 0, &amp; \forall ~ n \geq 1; \\
u^{n}_{k} = \left ( 1 + \frac{2k}{n+2} \right ) \gamma^{n-1}_{k} + \frac{2k-n-1}{n+2} \gamma^{n-1}_{k-1}, &amp; \forall~ n \geq 2, 1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor,
\end{cases}
$$
which is the recursion and initial conditions for $\gseq$ defined by \eqref{gamma_linear_recursion}. Since  \eqref{gamma_linear_recursion} defines a unique sequence, it follows that the sequence $\useq$ is the same as the sequence $\gseq$. Thus for all $n \geq 0$ and $0 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor$,
$$
\gamma^{n}_{k} = \left ( -\frac{1}{4} \right )^k \binom{n}{2k} C_k.
$$
\end{enumerate}
\end{proof}

\section{Nonlinear recursion relation for \texorpdfstring{$\gseq$}{g}} \label{nlr_section}

In this section we obtain a nonlinear recursion relation for the sequence $\gseq$ by using its relationship with the partial derivatives of the two-dimensional distance function. Specifically, we introduce a numerical scheme and apply it to $\nabla \cdot \left (d \nabla d \right )$,  where $d$ is the distance function defined in~\eqref{distancefn}.  We will show that the scheme is exact on $\left ( \nabla \cdot \left (d \nabla d \right )\right )(x,y)$ (as long as $(x,y) \neq (0,0)$) which will enable us to obtain the nonlinear recursion relation for $\gseq$.

\subsection{Two-dimensional finite difference scheme}

In this section we introduce a finite difference scheme for discretizing the differential operator $\nabla \cdot \left (f \nabla f \right )$ where $f$ is a $C^\infty$ function on $\mathbb R^2$ except possibly at finitely many points in $\mathbb R^2$. Let $\alpha&lt;\beta$ and assume $(0,0) \notin [\alpha,\beta]^2 $. Assume $N \in \mathbb N$. We consider the uniform discretization of $[\alpha,\beta]^2$, where $h = \frac{\beta-\alpha}{N}$ is the grid size in each of the coordinate directions. Then for a grid point $(x_i,y_j)$, $0\leq i, j\leq N$, and a $C^\infty$ function $f$ defined on  $[\alpha,\beta]^2$, we define the forward difference operator in $x$ as
$$
D^{x}_{+} \left ( f_{i,j} \right ) := \frac{f(x_{i+1},y_j) - f(x_i,y_j)}{h} = \frac{f(x_i+h,y_j) - f(x_i,y_j)}{h},
$$
and the forward difference operator in $y$ as
$$
D^{y}_{+} \left ( f_{i,j} \right ) := \frac{f(x_{i},y_{j+1}) - f(x_i,y_j)}{h} = \frac{f(x_i,y_j+h) - f(x_i,y_j)}{h}.
$$
Similarly, we define the backward difference operator in $x$ as
$$
D^{x}_{-} \left ( f_{i,j} \right ) := \frac{f(x_i,y_j) - f(x_{i-1},y_j)}{h} = \frac{f(x_i,y_j) - f(x_i-h,y_j)}{h},
$$
and the backward difference operator in $y$ as
$$
D^{x}_{-} \left (f_{i,j} \right ) := \frac{f(x_i,y_j) - f(x_i,y_{j-1})}{h} = \frac{f(x_i,y_j) - f(x_i,y_j-h)}{h}.
$$

We now define
$$
f_{i-\frac{1}{2},j} := \frac{1}{2} \left ( f_{i-1,j} +f_{i,j} \right ),
$$
and
$$
f_{i,j-\frac{1}{2}} := \frac{1}{2} \left ( f_{i,j} +f_{i,j-1} \right ),
$$
and construct the numerical scheme
\begin{equation}\label{finite_diff_scheme}
D_{+}^{x} \left ( f_{i-\frac{1}{2},j} D_{-}^{x} (f_{i,j}) \right ) + D_{+}^{y} \left ( f_{i,j-\frac{1}{2}} D_{-}^{y} (f_{i,j}) \right ).
\end{equation}
This scheme \eqref{finite_diff_scheme} is a discretized version of $\nabla \cdot \left ( f \nabla f \right )$. We note that this scheme uses a compact stencil in the sense that the stencil used is the smallest stencil on which we can construct a second order convergent discretization scheme for the differential operator  $\nabla \cdot \left ( f \nabla f \right )$. Indeed, for a smooth function $f$, we can show that this numerical scheme is second order. However, with the choice $f=d$, this discretization coincides with its non-discretized version $\nabla \cdot \left (d \nabla d \right )$. Now in the case of the distance function, the scalar $\left (\nabla \cdot \left (d \nabla d \right )\right)(x,y)$ can be calculated to be $2$ at every point $(x,y) \neq (0,0)$.

\begin{theorem} \label{exact_scheme}
Let $d$ be the distance function defined in \eqref{distancefn}. Choose $[\alpha,\beta]^2$ such that $(0,0) \notin [\alpha,\beta]^2$ and discretize $[\alpha,\beta]^2$ uniformly with grid spacing $h = \frac{\beta-\alpha}{N}$ in each coordinate direction. Then
$$
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right ) + D_{+}^{y} \left (d_{i,j-\frac{1}{2}} D_{-}^{y} (d_{i,j}) \right ) = 2.
$$
\end{theorem}

\begin{proof} We prove that $D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right )=1$. Since the distance function $d$ and the numerical scheme are symmetric in $x$ and $y$, we can obtain $D_{+}^{y} \left (d_{i,j-\frac{1}{2}} D_{-}^{y} (d_{i,j}) \right )=1$ by interchanging the roles of $x$ and $y$.


\begin{align*}
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right ) &amp; = D_{+}^{x} \left ( \frac{d_{i,j}+d_{i-1,j}}{2} \frac{d_{i,j} - d_{i-1,j}}{h} \right )\\
&amp; = \frac{1}{2h} D_{+}^{x} \left ( d_{i,j}^2 - d_{i-1,j}^2 \right ) \\
&amp; = \frac{1}{2h^2} \left ( d_{i+1,j}^2 - d_{i,j}^2 - d_{i,j}^2 + d_{i-1,j}^2 \right ) \\
&amp; = \frac{d_{i+1,j}^2 - 2d_{i,j}^2 + d_{i-1,j}^2}{2h^2}.
\end{align*}
Now recall that $$d_{i,j}^2 = x_{i}^2 + y_{j}^2,$$ $$d_{i+1,j}^2 = x_{i+1}^2 + y_{j}^2 = (x_{i}+h)^2 + y_{j}^2$$ and $$d_{i-1,j}^2 = x_{i-1}^2 + y_{j}^2 = (x_{i}-h)^2 + y_{j}^2.$$ Thus,
\begin{align*}
d_{i+1,j}^2 - 2d_{i,j}^2 + d_{i-1,j}^2 &amp; = (x_{i}+h)^2 + y_{j}^2 - 2x_{i}^2 -2 y_{j}^2 + (x_{i}-h)^2 + y_{j}^2  \\
&amp; = x_i^2 - 2hx_i + h^2 - 2 x_i^2 + x_i^2 - 2hx_i+h^2 \\
&amp; = 2h^2.
\end{align*}
Consequently, we obtain
$$
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right ) = \frac{d_{i+1,j}^2 - 2d_{i,j}^2 + d_{i-1,j}^2}{2h^2} = 1.
$$
\end{proof}

\begin{remark} Since the proof of Theorem~\ref{exact_scheme} is done coordinate-wise, it follows that the result still holds if we use a different grid spacing in each coordinate direction. In other words, as long as the discretization is uniform in each coordinate direction, it is not necessary to take the same grid spacing. \end{remark}

\subsection{Nonlinear recursion relation}

In this section, we obtain the nonlinear recursion relation for $\gseq$.

\begin{theorem} \label{nonlinear_thm}
Let $\gseq$ be defined by the recursion relation and initial conditions given in \eqref{gamma_linear_recursion}. Then the sequence satisfies the nonlinear recursion relation
\begin{equation} \label{gamma_nonlinear_recursion}
\begin{cases}
\gamma^{n}_{0} = 1, &amp; \forall ~ n \geq 0; \\
\gamma^{2n-1}_{n} = 0, &amp; \forall ~ n \geq 1; \\
\gamma^{n}_{k} = \gamma^{n-1}_{k} - \frac{1}{4} \sum\limits_{i=0}^{k-1} \sum\limits_{j=2i}^{n + 2(i-k)} \gamma^{j}_{i} \gamma^{n-2-j}_{k-i-1}, &amp; \forall~ n \geq 2, 1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor.

\end{cases}
\end{equation}
\end{theorem}

For the sake of clarity, we shall first prove several lemmas, which will then be employed to obtain the result of Theorem~\ref{nonlinear_thm}.

Let's first start with a brief discussion of our method. The basic idea is very simple and comes from standard numerical analysis techniques. In numerical analysis, it is important to obtain the convergence rate of a specific numerical scheme, namely we would like to know how fast the finite difference approximation converges to the true solution as the grid size $h$ decreases to $0$. In other words, we want to obtain a relation of the type
\begin{equation} \label{error_num_scheme}
\mbox{Finite Difference Scheme} = \mbox{True Solution} + O(h^p),
\end{equation}
where $h$ is the grid spacing in each coordinate direction (assumed to be the same in each coordinate direction here) and $p \in \mathbb N$. $p$ is called the order of convergence of the scheme. The higher $p$ is, the more accurate the scheme is, since it converges faster to the true solution.  To obtain such a relation, the standard procedure in numerical analysis is to use Taylor expansion for the finite difference scheme.

We do the same with our finite difference scheme and note that in our case \eqref{error_num_scheme} is valid for all $p \in \mathbb N$. Using a Taylor expansion, we write the finite difference scheme as a power series in $h$, knowing that all the coefficients are zero except for the constant term (since the scheme is exact by Theorem~\ref{exact_scheme}). Using the general term for the sequence of coefficients and setting it to zero will give us the nonlinear recursion.

We note that it suffices to consider $D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right )$ to obtain the nonlinear recursion relation. Thus from now on, we only consider the finite difference in $x$
$$
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right ).
$$
since an analysis of $D_{+}^{y} \left ( f_{i,j-\frac{1}{2}} D_{-}^{y} (f_{i,j}) \right )$ generates the same nonlinear recursion relation.

For the sake of clarity, we introduce some notation.
\begin{definition} \label{def_pq}
For $d$ defined in \eqref{distancefn}, we define the sequences of functions $\left ( p_{k}(x,y) \right )_{k\geq0}$ and $\left (q_{k}(x,y) \right )_{k\geq0}$ as
$$
p_{k}(x,y) : = \begin{cases} d, &amp; \mbox{ if } k=0; \\ \frac{(-1)^k}{2k!} \frac{\partial ^{(k)} d }{\partial x^k}, &amp; \mbox{ if } k \in \mathbb N, \end{cases}
$$
and
$$
q_{k}(x,y) := \frac{(-1)^{k}}{(k+1)!} \frac{\partial ^{(k+1)}  d}{\partial x^{k+1}},  \mbox{ for } k \in \mathbb N_{0}.
$$
\end{definition}
Note that we dropped the dependence on $x$ and $y$ in $d$ above. From now on, if the dependence on $x$ and $y$ is clear from context we will not write it.

\begin{definition} \label{define_S}
We define, for an integer $i \geq 0$,
$$
S_{i}:=\sum_{\ell = 0}^{i} p_{\ell} q_{i-\ell},
$$
where $\left ( p_{k}(x,y) \right )_{k\geq0}$ and $\left (q_{k}(x,y) \right )_{k\geq0}$ are defined in Definition~\ref{def_pq}.
\end{definition}

We are now ready to state our first lemma.
\begin{lemma} \label{taylor_scheme}
Let $d$ be the distance function defined in~\eqref{distancefn}. Then we have
\begin{equation}
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right ) = \sum_{n=0}^{\infty} a_n h^n,
\end{equation}
where
\begin{equation} \label{def_an}
a_n = \sum_{k=0}^{n}   \frac{1}{(n-k+1)!}\frac{\partial ^{(n-k+1)} S_{k}}{\partial x^{n-k+1}} ,
\end{equation}
with $S_k$ defined in Definition \ref{define_S}.
\end{lemma}

\begin{proof}
Let $(x_i,y_j)$ be any point in $[\alpha,\beta]^2$ such that $(0,0) \notin [\alpha,\beta]^2$. Since $d$ is $C^{\infty}$ in a neighborhood of $(x_i,y_j)\neq (0,0)$, we can expand $D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right )$ with Taylor series expansions. For simplicity in the exposition, we drop the argument $(x_i,y_j)$ from the distance function $d$. It is thus assumed that the distance function $d$ and any of its derivatives are evaluated at the grid node $(x_i,y_j)$.


\begin{align*}
D_{-}^{x}(d_{i,j})  &amp;= \frac{d_{i,j} - d_{i-1,j}}{h} = \frac{d(x_{i},y_{j}) - d(x_{i}-h,y_{j})}{h} \\
&amp; = \frac{1}{h} \left (d(x_{i},y_{j}) - \left ( d(x_{i},y_{j}) + \sum_{k=1}^{\infty} \frac{1}{k!} \frac{\partial ^{(k)} d}{\partial x^{k}} (-1)^k h^k \right )  \right )  \\
&amp; = \sum_{k=1}^{\infty} \frac{(-1)^{k-1} h^{k-1} }{k!} \frac{\partial ^{(k)} d}{\partial x^{k}}  \\
&amp; = \sum_{k=0}^{\infty} \frac{(-1)^{k}}{(k+1)!} \frac{\partial ^{(k+1)} d}{\partial x^{k+1}}   h^{k}\\
&amp; = \sum_{k=0}^{\infty}  q_{k} h^{k},
\end{align*}
by Definition~\ref{def_pq}.
Now
\begin{align*}
d_{i-\frac{1}{2},j}  &amp; = \frac{1}{2} \left ( d_{i,j} + d_{i-1,j} \right ) \\
&amp; = \frac{1}{2} \left ( d_{i,j} + d_{i,j} + \sum_{k=1}^{\infty} \frac{(-1)^k h^k}{k!} \frac{\partial ^{(k)} d}{\partial x^{k}} (x_{i},y_{j}) \right ) \\
&amp; = d_{i,j} + \frac{1}{2}  \sum_{k=1}^{\infty} \frac{(-1)^k h^k}{k!} \frac{\partial ^{(k)} d}{\partial x^{k}}\\
&amp; =  \sum_{k=0}^{\infty} p_{k} h^k,
\end{align*}
by Definition~\ref{def_pq} again.
Applying the Cauchy product, we obtain
$$
d_{i-\frac{1}{2},j} D_{-}^{x}(d)  = \sum_{k=0}^{\infty} \left ( \sum_{\ell=0}^{k} p_{\ell} q_{k-\ell}\right ) h^k = \sum_{k=0}^{\infty} S_{k} h^k,
$$
using Definition~\ref{define_S}.
Now since the Taylor series expansion of the forward difference $D_{+}^{x}$ of a function $f$ at $(x_i,j_j)$ is
$$
\sum_{m=0}^{\infty} \frac{1}{(m+1)!} \frac{\partial ^{(m+1)} f}{\partial x^{m+1}} h^{m},
$$
we can write
$$
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x}(d)  \right ) = \sum_{m=0}^{\infty} \frac{1}{(m+1)!} \sum_{k=0}^{\infty}  \frac{\partial ^{(m+1)}S_k }{\partial x^{m+1}}  h^{k+m}.
$$
Changing variables in the above sum, we obtain
\begin{align*}
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d) \right )  &amp; = \sum_{m+k=n, n \geq 0}  \frac{1}{(m+1)!} \frac{\partial ^{(m+1)} S_k }{\partial x^{m+1}} h^{n} \\
&amp; = \sum_{n=0}^{\infty} \underbrace{\sum_{k=0}^{n}   \frac{1}{(n-k+1)!}\frac{\partial ^{(n-k+1)}S_k }{\partial x^{n-k+1}} }_{=a_n} h^n \\
&amp; = \sum_{n=0}^{\infty} a_n h^n.
\end{align*}
\end{proof}

\begin{lemma} \label{lemma_pq}
We have
$$
p_{2m} =
\begin{cases}
D^{\frac{1}{2}}, &amp; \mbox{ if } m=0; \\
\frac{1}{4} Y D^{-\frac{4m-1}{2}} P_{m-1}^{2m}(X,Y), &amp; \mbox{ if } m \geq 1.
\end{cases}
$$

$$
p_{2m+1} =
\begin{cases}
-\frac{1}{2} D^{-\frac{1}{2}} x, &amp; \mbox{ if } m=0; \\
\frac{1}{4} x Y D^{-\frac{4m+1}{2}} P^{2m+1}_{m-1}(X,Y), &amp; \mbox{ if } m \geq 1.
\end{cases}
$$

$$
q_{2m} =
\begin{cases}
D^{-\frac{1}{2}} x, &amp; \mbox{ if } m=0; \\
-\frac{1}{2} x Y D^{-\frac{4m+1}{2}} P^{2m+1}_{m-1}(X,Y), &amp; \mbox{ if } m \geq 1.
\end{cases}
$$

$$
q_{2m-1} =
-\frac{1}{2}Y D^{-\frac{4m-1}{2}} P^{2m}_{m-1}(X,Y), \mbox{ if } m \geq 1.
$$
\end{lemma}

\begin{proof} These relations are consequences of combining Definition~\ref{def_pq} and the result of Theorem~\ref{formula_deriv_d}.

\end{proof}

\begin{lemma}\label{lemma_product_pq} For $p,q$ defined in Definition~\ref{def_pq}, we have
$$
p_{2k} q_{2m-2k} =
\begin{cases}
-\frac{1}{2} x Y D^{-2m} P^{2m+1}_{m-1}, &amp; k=0, m \geq 1; \\
-\frac{1}{2} x Y D^{-2m} \left ( \frac{1}{4} Y P^{2k}_{k-1} P^{2m-2k+1}_{m-k-1}\right ), &amp; k \geq 1, m \geq k+1.
\end{cases}
$$

$$
p_{2k+1} q_{2m-2k-1} =
\begin{cases}
-\frac{1}{2} x Y D^{-2m} \left ( -\frac{1}{2} P^{2m}_{m-1}\right ), &amp; k=0, m \geq 1; \\
 -\frac{1}{2} x Y D^{-2m} \left (\frac{1}{4} Y P^{2k+1}_{k-1} P^{2m-2k}_{m-k-1} \right ), &amp;  k \geq 1, m \geq k+1.
\end{cases}
$$

$$
p_{2k} q_{2m+1-2k} =
\begin{cases}
-\frac{1}{2} Y D^{-(2m+1)} P^{2m+2}_{m}, &amp; k=0, m \geq 0; \\
-\frac{1}{2} Y D^{-(2m+1)} \left ( \frac{1}{4} Y P^{2k}_{k-1} P^{2(m-k+1)}_{m-k} \right ), &amp; k \geq 1, m \geq k.
\end{cases}
$$
$$
p_{2k+1} q_{2m-2k} =
\begin{cases}
- \frac{1}{2} Y D^{-(2m+1)} \left (-\frac{1}{2} X P^{2m+1}_{m-1}  \right ), &amp; k=0, m \geq 1;  \\
-\frac{1}{2} Y D^{-(2m+1)} \left ( \frac{1}{4} XY P^{2k+1}_{k-1} P^{2(m-k)+1}_{m-k-1} \right ), &amp;  k \geq 1, m \geq k+1.
\end{cases}
$$


\end{lemma}

\begin{proof} These relations are consequences of Lemma~\ref{lemma_pq}.

\end{proof}


\begin{lemma} \label{symmetry_pq}
For $\ell \geq 1$,$k \geq \ell$, we have
\begin{equation}
p_{\ell} q_{k-\ell} = p_{k-\ell+1}q_{\ell-1}.
\end{equation}
\end{lemma}

\begin{proof}
Suppose $\ell \geq 1$,$k \geq \ell$. Then
\begin{align*}
p_{\ell} q_{k-\ell} &amp; = \frac{(-1)^\ell}{2\ell!} \frac{\partial^{(\ell)} d}{\partial x^\ell} \frac{(-1)^{(k-\ell)}}{(k-\ell+1)!} \frac{\partial^{(k-\ell+1)} d}{\partial x^{k-\ell+1}} \\
&amp; = \frac{(-1)^k}{2\ell! (k-\ell+1)!} \frac{\partial^{(\ell)} d}{\partial x^\ell} \frac{\partial^{(k-\ell+1)} d}{\partial x^{k-\ell+1}}.
\end{align*}
Thus
\begin{align*}
p_{k-\ell+1}q_{\ell-1} &amp; = \frac{(-1)^{k-\ell+1}}{2(k-\ell+1)!} \frac{\partial^{(k-\ell+1)}d}{\partial x^{k-\ell+1}} \frac{(-1)^{\ell-1}}{\ell!} \frac{\partial ^{\ell} d}{\partial x^\ell} \\
&amp; = \frac{(-1)^k}{2\ell! (k-\ell+1)!} \frac{\partial ^{\ell} d}{\partial x^\ell} \frac{\partial^{(k-\ell+1)}d}{\partial x^{k-\ell+1}} \\
&amp; = p_{\ell} q_{k-\ell}.
\end{align*}
\end{proof}

\begin{lemma} \label{S_with_gamma}
For $m \geq 0$, $$S_{2m+2} =  \left ( -\frac{1}{2} x Y D^{-(2m+2)}\right ) \mathcal A_{m},$$ where $\mathcal A_{m}$ is a homogeneous polynomial of degree $m$ in $X$ and $Y$. Similarly for $m \geq 0$, $$S_{2m+3} = \left ( -\frac{1}{2} Y D^{-(2m+3)} \right ) \mathcal B_{m+1},$$ where $\mathcal B_{m+1}$ is a homogeneous polynomial of degree $m+1$ in $X$ and $Y$. In addition, we have the following:
\begin{itemize}

\item The coefficient of $X^{m}$ in $\mathcal A_m$ is
$$
\gamma^{2m+1}_{0} - \gamma^{2m}_{0}.
$$

\item The coefficient of $Y^n X^{m-n}$ in  $\mathcal A_m$, for $1 \leq n \leq m$  is
$$
\gamma^{2m+1}_{n} - \gamma^{2m}_{n} + \frac{1}{4} \sum_{k=0}^{n-1} \sum_{j = 2k}^{2m+1+2(k-n)} \gamma^{j}_{k} \gamma^{2m-1-j}_{n-k-1}.
$$

\item The coefficient of $X^{m+1}$ in $\mathcal B_{m+1}$ is
$$
\gamma^{2m+2}_{0} - \gamma^{2m+1}_{0}.
$$

\item The coefficient of $Y^{n}X^{m+1-n}$ in $\mathcal B_{m+1}$, for $1 \leq n \leq m$ is
$$
 \gamma^{2m+2}_{n} - \gamma^{2m+1}_{n} + \frac{1}{4}  \sum_{k=0}^{n-1} \sum_{j = 2k}^{2m+2+2(k-n)} \gamma^{j}_{k} \gamma^{2m-j}_{n-k-1}.
$$

\item The coefficient of $Y^{m+1}$ in $\mathcal B_{m+1}$ is
$$
 \gamma^{2m+2}_{m+1} + \frac{1}{4} \sum_{k=0}^{m} \gamma^{2k}_{k} \gamma^{2(m-k)}_{m-k}.
$$
\end{itemize}

\end{lemma}

\begin{proof}

Let $m \geq 0$. We shall show the proof for $S_{2m+3}$. The argument is similar (and slightly simpler) for $S_{2m+2}$. 

Using Definition~\ref{define_S}, the symmetry stated in Lemma~\ref{symmetry_pq} and the result of Lemma~\ref{lemma_product_pq}, we can write
\begin{align*}
S_{2m+3} &amp; = p_{0} q_{2m+3} + 2 p_{1}q_{2m+2} + \sum_{\ell=2}^{2m+2} p_{\ell} q_{2m+3-\ell} \\
&amp; = p_{0} q_{2m+3} + 2 p_{1}q_{2m+2} + \sum_{j=1}^{m+1} p_{2j} q_{2m+3-2j} + \sum_{j=1}^{m} p_{2j+1} q_{2m+2-2j}  \\
&amp; = \left (\frac{-1}{2} Y D^{-(2m+3)} \right ) \mathcal B_{m+1},
\end{align*}
where
\begin{align*}
\mathcal B_{m+1} &amp;= \left ( P^{2m+4}_{m+1} - X P^{2m+3}_{m} + \frac{Y}{4} \sum_{j=1}^{m+1} P^{2j}_{j-1} P^{2m-2j+4}_{m+1-j} + \frac{XY}{4} \sum_{j=1}^{m} P^{2j+1}_{j-1} P^{2m-2j+3}_{m-j}\right )\\
&amp; = \sum_{n=0}^{m+1} \gamma^{2m+2}_{n} Y^{n} X^{m+1-n} - X \sum_{n=0}^{m} \gamma^{2m+1}_{n} Y^n X^{m-n} \\
&amp; + \frac{Y}{4} \sum_{j=1}^{m+1} \left ( \sum_{k=0}^{j-1} \gamma^{2j-2}_{k} Y^k X^{j-1-k} \right ) \left ( \sum_{i=0}^{m+1-j} \gamma^{2m-2j+2}_{i} Y^i X^{m+1-j-i} \right ) \\
&amp; + \frac{XY}{4} \sum_{j=1}^{m} \left ( \sum_{k=0}^{j-1} \gamma^{2j-1}_{k} Y^k X^{j-1-k} \right ) \left ( \sum_{i=0}^{m-j} \gamma^{2m-2j+1}_{i} Y^i X^{m-j-i} \right ).
\end{align*}
Thus
\begin{align*}
\mathcal B_{m+1} &amp; = \sum_{n=0}^{m+1} \gamma^{2m+2}_{n} Y^{n} X^{m+1-n} - X \sum_{n=0}^{m} \gamma^{2m+1}_{n} Y^n X^{m-n} \\
&amp; + \frac{Y}{4} \sum_{j=1}^{m+1} \left ( \sum_{n=0}^{j-1+m+1-j} \left (  \sum_{k = \max(0,n-m+j-1)}^{\min(n,j-1)} \gamma^{2j-2}_{k} \gamma^{2m-2j+2}_{n-k} X^{j-1-k} X^{m+1-j-n+k} \right )Y^n \right ) \\
&amp; + \frac{XY}{4} \sum_{j=1}^{m} \left ( \sum_{n=0}^{j-1+m-j} \left (  \sum_{k=\max(0,n-m+j)}^{\min(n,j-1)} \gamma^{2j-1}_{k} \gamma^{2m-2j+1}_{n-k} X^{j-1-k} X^{m-j-n+k}\right ) Y^n \right ) \\
&amp; = \sum_{n=0}^{m+1} \gamma^{2m+2}_{n} Y^{n} X^{m+1-n} - X \sum_{n=0}^{m} \gamma^{2m+1}_{n} Y^n X^{m-n} \\
&amp; + \frac{1}{4} \sum_{j=1}^{m+1} \sum_{n=0}^{m} \left ( \sum_{k = \max(0,n-m+j-1)}^{\min(n,j-1)} \gamma^{2j-2}_{k} \gamma^{2m-(2j-2)}_{n-k}  \right ) Y^{n+1} X^{m-n} \\
&amp; + \frac{1}{4} \sum_{j=1}^{m} \sum_{n=0}^{m-1} \left (  \sum_{k = \max(0,n-m+j)}^{\min(n,j-1)} \gamma^{2j-1}_{k} \gamma^{2m-(2j-1)}_{n-k} \right ) Y^{n+1} X^{m-n} \\
&amp; = \sum_{n=0}^{m+1} \gamma^{2m+2}_{n} Y^{n} X^{m-n+1} -  \sum_{n=0}^{m} \gamma^{2m+1}_{n} Y^n X^{m-n+1} \\
&amp; + \frac{1}{4}  \sum_{n=1}^{m+1} \sum_{j=1}^{m+1} \left ( \sum_{k = \max(0,n-m+j-2)}^{\min(n-1,j-1)} \gamma^{2j-2}_{k} \gamma^{2m-(2j-2)}_{n-k-1}  \right ) Y^{n} X^{m-n+1}  \\
&amp; + \frac{1}{4} \sum_{n=1}^{m} \sum_{j=1}^{m} \left (  \sum_{k = \max(0,n-m+j-1)}^{\min(n-1,j-1)} \gamma^{2j-1}_{k} \gamma^{2m-(2j-1)}_{n-k-1} \right ) Y^{n} X^{m-n+1}.
\end{align*}

Consequently, the coefficient of $X^{m+1}$ in $\mathcal B_{m+1}$ is
$$
\gamma^{2m+2}_{0} - \gamma^{2m+1}_{0}.
$$
Also, the coefficient of $Y^{m+1}$ in $\mathcal B_{m+1}$ $(n=m+1)$ is
\begin{align*}
&amp; \gamma^{2m+2}_{m+1} + \frac{1}{4} \sum_{j=1}^{m+1} \sum_{k=\max(0,j-1)=j-1}^{\min(m,j-1)=j-1} \gamma^{2j-2}_{k} \gamma^{2m-(2j-2)}_{m-k} \\
&amp; = \gamma^{2m+2}_{m+1} + \frac{1}{4}  \sum_{j=1}^{m+1} \gamma^{2j-2}_{j-1} \gamma^{2m-(2j-2)}_{m-j+1} \\
&amp; = \gamma^{2m+2}_{m+1} + \frac{1}{4}  \sum_{j=0}^{m} \gamma^{2j}_{j} \gamma^{2(m-j)}_{m-j}.
\end{align*}

Now for $1 \leq n \leq m$, the coefficient of $Y^n X^{m+1-n}$ is
\begin{align*}
 \gamma^{2m+2}_{n}  - \gamma^{2m+1}_{n} &amp; +\frac{1}{4} \sum_{j=1}^{m+1} \left ( \sum_{k = \max(0,n-m+j-2)}^{\min(n-1,j-1)} \gamma^{2j-2}_{k} \gamma^{2m-(2j-2)}_{n-k-1}  \right )\\
  &amp; + \frac{1}{4}\sum_{j=1}^{m} \left (  \sum_{k = \max(0,n-m+j-1)}^{\min(n-1,j-1)} \gamma^{2j-1}_{k} \gamma^{2m-(2j-1)}_{n-k-1} \right ) \\
 = \gamma^{2m+2}_{n}  - \gamma^{2m+1}_{n} &amp; + \frac{1}{4} \sum_{r=0}^{2m} \sum_{k = \max(0,n-m-1+\frac{r}{2})}^{\min(n-1,\frac{r}{2})} \gamma^{r}_{k} \gamma^{2m-r}_{n-k-1} \\
 =   \gamma^{2m+2}_{n}  - \gamma^{2m+1}_{n} &amp; + \frac{1}{4} \sum_{k=0}^{n-1} \sum_{r = 2k}^{2m+2+2(k-n)}  \gamma^{r}_{k} \gamma^{2m-r}_{n-k-1}.
\end{align*}
\end{proof}

\begin{lemma} \label{Sn_zero}
Let $\left (S_n\right)_{n \geq 0}$ be defined as in Definition~\ref{define_S}. Then
\begin{equation}
\begin{cases}
S_0 = x, &amp; \\
S_1 = -\frac{1}{2,} &amp; \\
S_n = 0, &amp; n \geq 2.
\end{cases}
\end{equation}
\end{lemma}

\begin{proof} We start with $S_0$ and $S_1$ and obtain the result by a direct calculation. We then show that $S_n = 0$ for $n \geq 2$ by induction.

We calculate $S_0$ using the definitions of $\left ( p_{k}(x,y) \right )_{k\geq 0}$ and $\left (q_{k}(x,y) \right )_{k\geq0}$ given in Definition~\ref{def_pq} and obtain
\begin{equation} \label{S0}
S_0 = p_{0}q_0 = d \frac{\partial d}{\partial x} = D^{\frac{1}{2}}\frac{x}{D^{\frac{1}{2}}} = x.
\end{equation}

For $S_1$ we use the results of Lemmas~\ref{lemma_pq} and \ref{lemma_product_pq} to write
\begin{align} \label{S1}
S_1 &amp; = p_0q_1 + p_1q_0 = -\frac{1}{2}YD^{-1}\underbrace{P^{2}_0}_{=\gamma^{0}_{0}=1} - \frac{1}{2}D^{-\frac{1}{2}}x D^{-\frac{1}{2}}x \nonumber\\
&amp; = -\frac{1}{2}D^{-1} \underbrace{\left (Y+X \right )}_{=D}  = -\frac{1}{2}.
\end{align}

The remainder of the proof focuses on $S_n$, for $n \geq 2$.
Combining the results of Theorem~\ref{exact_scheme} and Lemma~\ref{taylor_scheme}, we can write
$$
D_{+}^{x} \left ( d_{i-\frac{1}{2},j} D_{-}^{x} (d_{i,j}) \right ) = \sum_{n=0}^{\infty} a_n h^n = 2,
$$
where $(a_n)_{n \geq 0}$ is given in Equation~\eqref{def_an}. It follows that
\begin{equation} \label{an_zero}
\begin{cases}
a_0 = 2, &amp; \\
a_n =  \sum\limits_{k=0}^{n} \frac{1}{(n-k+1)!}\frac{\partial ^{(n-k+1)} S_{k}}{\partial x^{n-k+1}} = 0, &amp; n \geq 1.
\end{cases}
\end{equation}
We now show by induction that Equation~\eqref{an_zero} implies $S_n = 0$, for all $n \geq 2$.

\begin{enumerate}

\item First consider the base cases $n=2$ and $n=3$.

\begin{itemize}
\item For $n=2$ we have
$$
a_2 = \sum_{k=0}^{2} \frac{1}{(3-k)!}\frac{\partial ^{(3-k)} S_{k}}{\partial x^{3-k}} =  \frac{1}{6} \frac{\partial ^{3} S_{0}}{\partial x^{3}} +  \frac{1}{2} \frac{\partial ^{2} S_{1}}{\partial x^{2}} +  \frac{\partial  S_{2}}{\partial x}.
$$
Since $S_0$ is a degree $1$ polynomial in $x$ and $S_1$ is constant by \eqref{S0} and \eqref{S1}, it follows that $\ds \frac{\partial ^{3} S_{0}}{\partial x^{3}}~ =~ \frac{\partial ^{2} S_{1}}{\partial x^{2}} ~= ~0$. Consequently
$$
a_2 = \frac{\partial  S_{2}}{\partial x}.
$$
Additionally, since $a_2 = 0$ by Equation~\eqref{an_zero}, it follows that
$$
\frac{\partial  S_{2}}{\partial x} = 0.
$$
Hence $S_2$ is constant in $x$. Furthermore, using Lemma~\ref{S_with_gamma} we know that
$$
S_2 = -\frac{1}{2}xYD^{-2} \mathcal{A},
$$
where $\mathcal{A}$ is constant. Note that $-\frac{1}{2}xYD^{-2}$ is a known non-constant and non-zero rational function of $x$ where the numerator is a polynomial of degree $1$ in $x$ and the denominator a polynomial of degree $4$ in $x$. Since $S_2$ is independent of $x$, it follows that necessarily $\mathcal A = 0$.
Thus
$$
S_2 = 0.
$$

\item For $n=3$ we have
$$
a_3 = \sum_{k=0}^{3} \frac{1}{(4-k)!}\frac{\partial ^{(4-k)} S_{k}}{\partial x^{4-k}} = \frac{1}{24}   \frac{\partial ^{4} S_{0}}{\partial x^{4}} +  \frac{1}{6} \frac{\partial ^{3} S_{1}}{\partial x^{3}} + \frac{1}{2} \frac{\partial ^{2} S_{2}}{\partial x^{2}} + \frac{\partial  S_{3}}{\partial x} = 0.
$$
Since $S_1$ is a polynomial of degree $1$, $S_0$ is constant and $S_2=0$, we have
$$
a_3 = \frac{\partial  S_{3}}{\partial x} = 0.
$$
Thus $S_3$ is constant in $x$. Furthermore, using Lemma~\ref{S_with_gamma} we know that
$$
S_3 = -\frac{1}{2} Y D^{-3} \mathcal{B},
$$
where $\mathcal{B}$ is a polynomial of degree $2$ in $x$. Again we observe that $-\frac{1}{2} Y D^{-3}$ is a known non-constant and non-zero rational function of $x$ where the numerator is constant in $x$ and the denominator a polynomial of degree $6$ in $x$. Since $S_3$ is independent of $x$ it follows that necessarily $\mathcal{B} = 0$. Hence
$$
S_3 = 0.
$$

\end{itemize}

\item For the induction hypothesis, we let $n \geq 2$ and assume that $S_k = 0$ for all $2 \leq k \leq n$. Then we have
\begin{align*}
a_{n+1} &amp; = \sum\limits_{k=0}^{n+1} \frac{1}{(n-k+2)!}\frac{\partial ^{(n-k+2)} S_{k}}{\partial x^{n-k+2}}  \\
&amp; = \frac{1}{(n+2)!}\frac{\partial ^{n+2} S_{0}}{\partial x^{n+2}}  + \frac{1}{(n+1)!}\frac{\partial ^{n+1} S_{1}}{\partial x^{n+1}} + \frac{1}{n!}\frac{\partial ^{n} S_{2}}{\partial x^{n}} + \cdots +\frac{\partial^2  S_{n}}{\partial x^2} +  \frac{\partial  S_{n+1}}{\partial x}.
\end{align*}

Since $n \geq 2$, it follows that $n+2\geq 4$ and $n+1\geq 3$. Consequently, since $S_0$ is a degree $1$ polynomial in $x$ and $S_1$ is constant, we have $$\ds \frac{\partial ^{n+2} S_{0}}{\partial x^{n+2}}  = \frac{\partial ^{n+1} S_{1}}{\partial x^{n+1}} = 0.$$ Additionally, using the induction hypothesis and Equation~\eqref{an_zero}, we obtain
$$
a_{n+1} = \frac{\partial  S_{n+1}}{\partial x}=0.
$$
Thus $S_{n+1}$ is constant in $x$. Using the result of Lemma~\ref{S_with_gamma}, we see that $S_{n+1}$ can be expressed as a product of a fixed non-constant and non-zero rational function of $x$ and a polynomial in $x$. (The expression of the rational function and the degree of the polynomial will vary depending on whether $n$ is even or odd).
In any case, since $S_{n+1}$ is independent of $x$ it follows that the polynomial in $S_{n+1}$ needs to be zero. Hence
$$
S_{n+1} = 0.
$$

Since $S_2=S_3=0$ it follows that $S_n = 0$ for all $n \geq 2$.
\end{enumerate}
\end{proof}

\begin{proof} Proof of Theorem~\ref{nonlinear_thm}. 

Let $\gseq$ be defined by~\eqref{gamma_linear_recursion}. Combining Lemmas~\ref{S_with_gamma} and \ref{Sn_zero}, we obtain for $m \geq 0$
$$
S_{2m+2} = S_{2m+3} =0.
$$
Thus
$$
S_{2m+2} = \left ( -\frac{1}{2} x Y D^{-(2m+2)}\right ) \mathcal A_{m} = 0,
$$
and $\mathcal A_m = 0$ for $m \geq 0$ (proved in the proof of Lemma~\ref{Sn_zero}).

Similarly,
$$
S_{2m+3} = \left ( -\frac{1}{2} Y D^{-(2m+3)} \right ) \mathcal B_{m+1} = 0,
$$
and $\mathcal{B}_{m+1}$ = 0 for $m \geq 0$ (proved in the proof of Lemma~\ref{Sn_zero}).
Since $\mathcal{A}_m$ and $\mathcal{B}_{m+1}$ are polynomials, it follows that all their coefficients are zero. Hence by Lemma~\ref{S_with_gamma} we obtain
$$
\begin{cases}
\gamma^{2m+1}_{0} - \gamma^{2m}_{0} = 0, &amp; \mbox{ for } m \geq 0;\\
\gamma^{2m+2}_{0} - \gamma^{2m+1}_{0} = 0, &amp; \mbox{ for } m \geq 0; \\
  \gamma^{2m+2}_{m+1} + \frac{1}{4} \sum\limits_{k=0}^{m} \gamma^{2k}_{k} \gamma^{2(m-k)}_{m-k}=0, &amp; \mbox{ for } m \geq 0;\\
\gamma^{2m+1}_{n} - \gamma^{2m}_{n} + \frac{1}{4} \sum\limits_{k=0}^{n-1} \sum\limits_{j = 2k}^{2m+1+2(k-n)} \gamma^{j}_{k} \gamma^{2m-1-j}_{n-k-1}=0, &amp; \mbox{ for } m \geq 1, 1 \leq n \leq m;\\
 \gamma^{2m+2}_{n} - \gamma^{2m+1}_{n} + \frac{1}{4}  \sum\limits_{k=0}^{n-1} \sum\limits_{j = 2k}^{2m+2+2(k-n)} \gamma^{j}_{k} \gamma^{2m-j}_{n-k-1} =0, &amp; \mbox{ for } m \geq 1, 1 \leq n \leq m.\\
\end{cases}
$$
We remark that we could also have obtained the first two equations by using $\gamma^{n}_{0} = 1$ for all $n \geq 0$. Observe also that the last three equations correspond to the single equation
\begin{equation} \label{recursion_gamma}
\gamma^{p}_{n} -\gamma^{p-1}_{n} + \frac{1}{4} \sum_{k=0}^{n-1} \sum_{j = 2k}^{p+2(k-n)} \gamma^{j}_{k} \gamma^{p-2-j}_{n-k-1} =0,
\end{equation}
where $p \geq 2$ and $1 \leq n \leq \left \lfloor \frac{p}{2} \right \rfloor$. The last equation corresponds to \eqref{recursion_gamma} for $p=2m+2$ and $1 \leq k \leq \left \lfloor \frac{p}{2} \right \rfloor-1$. The equation above it corresponds to \eqref{recursion_gamma} for $p=2m+1$ and $1 \leq k \leq \left \lfloor \frac{p}{2} \right \rfloor$, and finally the third equation above corresponds to \eqref{recursion_gamma}  for $p=2m+2$ and $k=m+1$.

Thus, by rewriting~\eqref{recursion_gamma} as
$$
\gamma^{p}_{n}  = \gamma^{p-1}_{n} - \frac{1}{4} \sum_{k=0}^{n-1} \sum_{j = 2k}^{p+2(k-n)} \gamma^{j}_{k} \gamma^{p-2-j}_{n-k-1},
$$
for $p \geq 2$ and $1 \leq n \leq \left \lfloor \frac{p}{2} \right \rfloor$,  we obtain the nonlinear recursion relation stated in Theorem~\ref{nonlinear_thm}.

\end{proof}


\begin{remark} The nonlinear recursion can also be verified by using the closed form formula for the sequence $\gseq$
$$
\gamma^{n}_{k} = \left  (-\frac{1}{4} \right )^k \binom{n}{2k} C_k,
$$
for $n \geq 0$ and $0 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor$.
In this case, for $n \geq 2$ and $1 \leq k \leq \left \lfloor \frac{n}{2} \right \rfloor$, we can write
$$
 - \frac{1}{4} \sum\limits_{i=0}^{k-1} \sum\limits_{j=2i}^{n + 2(i-k)} \gamma^{j}_{i} \gamma^{n-2-j}_{k-i-1}  = \left ( -\frac{1}{4} \right )^{k} \sum\limits_{i=0}^{k-1} C_i C_{k-1-i} \sum\limits_{j=2i}^{n - 2k + 2i}  \binom{j}{2i} \binom{n-2-j}{2k-2-2i}.
$$
Changing variables in the inner sum, we can write
$$
 \sum\limits_{j=2i}^{n - 2k + 2i}  \binom{j}{2i} \binom{n-2-j}{2k-2-2i}  =  \sum\limits_{j=0}^{n - 2k}  \binom{j+2i}{2i} \binom{n-2-(j+2i)}{2k-2-2i}.
$$
Using Theorem~\ref{binomthm_risingpoly} (binomial theorem for rising polynomials) we can rewrite the inner sum as
\begin{align*}
\sum\limits_{j=0}^{n - 2k}  \binom{j+2i}{2i} &amp;\binom{n-2-(j+2i)}{2k-2-2i} \\&amp; = \frac{1}{(n-2k)!} \sum\limits_{j=0}^{n - 2k}  \frac{(n-2k)!}{j! (n-2k-j)!} \frac{(j+2i)!}{(2i)!} \frac{(n-2-(j+2i))!}{(2k-2i-2)!}\\
&amp; =  \frac{1}{(n-2k)!} \sum_{j=0}^{n-2k} \binom{n-2k}{j} (2i+1)^{\overline{j}} (2k-2i-1)^{\overline{n-k-j}} \\
&amp; = \frac{1}{(n-2k)!}  (2k)^{\overline{n-2k}} \\
&amp; =  \frac{(n-1)!}{(2k)!(n-2k)!} \\
&amp; = \binom{n-1}{2k-1}.
\end{align*}
Thus, the inner sum is simply  the binomial coefficient $ \ds \binom{n-1}{2k-1}$. Consequently, we have
$$
 - \frac{1}{4} \sum\limits_{i=0}^{k-1} \sum\limits_{j=2i}^{n + 2(i-k)} \gamma^{j}_{i} \gamma^{n-2-j}_{k-i-1} = \left ( -\frac{1}{4} \right )^{k}  \binom{n-1}{2k-1} \sum\limits_{i=0}^{k-1} C_i C_{k-1-i}.
 $$
We now employ the recursion relation for the Catalan numbers \eqref{recursion_segner} to obtain
$$
- \frac{1}{4} \sum\limits_{i=0}^{k-1} \sum\limits_{j=2i}^{n + 2(i-k)} \gamma^{j}_{i} \gamma^{n-2-j}_{k-i-1} =  \left ( -\frac{1}{4} \right )^{k}  \binom{n-1}{2k-1} C_k.
$$
Finally, since $\ds \binom{n-1}{2k-1} = \binom{n}{2k} - \binom{n-1}{2k}$, we have
\begin{align*}
- \frac{1}{4} \sum\limits_{i=0}^{k-1} \sum\limits_{j=2i}^{n + 2(i-k)} \gamma^{j}_{i} \gamma^{n-2-j}_{k-i-1} &amp; = \left ( -\frac{1}{4} \right )^{k}  C_k \left (\binom{n}{2k} - \binom{n-1}{2k} \right )  \\
&amp; = \gamma^{n}_{k} - \gamma^{n-1}_{k}.
\end{align*}
\end{remark}

\section{Generating function} \label{generatingfn}

In this section we calculate the generating function for the sequence $\gseq$. Interestingly, the calculations involve the generalized Motzkin numbers \eqref{gen_motzkin} and their generating function \eqref{motzkin_genfn}.  We note that the generalized Motzkin numbers $M_n(a,b)$ introduced by Sun \cite{sun14} were only studied for $a,b \in \mathbb Z$. Therefore, the generating function \eqref{motzkin_genfn} was shown to exist for $a,b \in \mathbb Z$. We first point out that the generalized Motzkin numbers as defined in \eqref{gen_motzkin} are well defined for $a,b \in \mathbb R$, as is their generating function \eqref{motzkin_genfn}. Secondly, we can show that the generalized Motzkin numbers \eqref{gen_motzkin} and the generating function \eqref{motzkin_genfn} both imply the recurrence relation \eqref{gen_motzkin_recurrence}. Since there is a unique solution to the recurrence relation \eqref{gen_motzkin_recurrence}, it follows that the generating function for the generalized Motzkin numbers for $a,b \in \mathbb R$ is the generating function \eqref{motzkin_genfn} obtained by Sun. Finally, we remark that if $x$ satisfies
\begin{equation} \label{conv_genfn_genmotzkin}
|x| &lt; \min \left ( \frac{1}{2|a|}, \frac{1}{\left |a + 2 \sqrt{|b|} \right |},   \frac{1}{\left |a - 2 \sqrt{|b|} \right |}\right ),
\end{equation}
then \eqref{motzkin_genfn} holds and its series is absolutely convergent.

\begin{theorem} Let $\gseq$ be defined by the recursion relation and initial conditions given in \eqref{gamma_linear_recursion}.  Then its generating function is
$$
F_{\gamma}(x,y) = \sum_{n=0}^{\infty} \sum_{k=0}^{\left \lfloor \frac{n}{2} \right \rfloor} \gamma^{n}_{k} x^k y^{n-2k} = 2 \frac{y-1 + \sqrt{(1-y)^2 + x}}{x}.
$$
\end{theorem}

\begin{proof} The proof uses the generalized Motzkin numbers \eqref{gen_motzkin} and their generating function~\eqref{motzkin_genfn}.

\begin{align*}
F_{\gamma}(x,y) &amp; = \sum_{n=0}^{\infty} \sum_{k=0}^{\left \lfloor \frac{n}{2} \right \rfloor} \gamma^{n}_{k} x^k y^{n-2k} \\
&amp; = \sum_{n=0}^{\infty} \sum_{k=0}^{\left \lfloor \frac{n}{2} \right \rfloor} \binom{n}{2k} C_k \left ( \frac{-x}{4} \right )^k y^{n-2k}  \\
&amp; = \sum_{n=0}^{\infty} M_{n}\left (y,-\frac{x}{4}\right). \\
\end{align*}
This series is the generating function for the Motzkin numbers $\ds M_{n}\left (y,-\frac{x}{4} \right )$ at $1$. Thus, using \eqref{conv_genfn_genmotzkin}, we see that if we select $x$ and $y$ such that
$$
\begin{cases}
|y| &lt; \frac{1}{2}, &amp; \\
\left | y + \sqrt{|x|}\right| &lt; 1, &amp; \\
\left | y - \sqrt{|x|}\right| &lt; 1, &amp; \\
\end{cases}
$$
then the series $\ds \sum_{n=0}^{\infty} M_{n}\left (y,-\frac{x}{4}\right)$ converges absolutely to $\ds \mathcal{M}_{y,-\frac{x}{4}}(1)$. Thus we have
\begin{align*}
F_{\gamma}(x,y) &amp; = \sum_{n=0}^{\infty} M_{n}\left (y,-\frac{x}{4}\right)\\
&amp; = \mathcal{M}_{y,-\frac{x}{4}}(1) \\
&amp; = \frac{1-y - \sqrt{(1-y)^2 - 4 \left ( -\frac{x}{4} \right )}}{2 \frac{-x}{4}} \\
&amp; =  2 \frac{y-1 + \sqrt{(1-y)^2 + x}}{x}.
\end{align*}

\end{proof}
\section{Conclusion} \label{conclusion} In this paper, we studied a doubly indexed sequence related to the Catalan and Motzkin numbers. We showed that this sequence satisfies both a linear and a nonlinear recursion relation. The nonlinear recursion relation was obtained by applying an exact finite difference scheme with compact stencil to the differential operator $\left ( \nabla \cdot \left (d \nabla d \right )\right )(x,y)$, where $d$ is the two-dimensional Euclidean distance function. We also calculated the generating function associated with this doubly indexed sequence. This work is, to our knowledge, the first instance of the usage of a finite difference scheme to obtain identities about a sequence, and another demonstration of the ubiquity of the Catalan numbers.


\begin{thebibliography}{10}

\bibitem{aigner98}
Martin Aigner, Motzkin numbers, {\em European J. Combin.} {\bf 19} (1998),
  663--675.

\bibitem{armstrong_rhoades_williams13}
Drew Armstrong, Brendon Rhoades, and Nathan Williams, Rational Catalan
  combinatorics: the associahedron, {\em Discrete Math. Theoret. Comput. Sci.}
  (2013), 933--944.

\bibitem{artioli_dattoli_licciardi_pagnutti17}
M.~Artioli, G.~Dattoli, S.~Licciardi, and S.~Pagnutti, Motzkin numbers: an
  operational point of view, preprint, 2017.  Available at
  \url{https://arxiv.org/abs/1703.07262}.

\bibitem{bernhart99}
F.~R. Bernhart, Catalan, {M}otzkin, and {R}iordan numbers, {\em Discrete Math.}
  {\bf 204} (1999), 73--112.

\bibitem{blasiak_dattoli_horzela_penson_zhukovsky08}
P.~Blasiak, G.~Dattoli, A.~Horzela, K. A. Penson, and K.~Zhukovsky, 
\href{https://cs.uwaterloo.ca/journals/JIS/VOL11/Penson/penson131.html}{Motzkin
  numbers, central trinomial coefficients and hybrid polynomials}, {\em J.
  Integer Sequences} {\bf 11} (2008), Article 08.1.1.

\bibitem{donaghey_shapiro77}
R.~Donaghey and L. W. Shapiro, Motzkin numbers, {\em J. Comb. Theory Ser. A}
  {\bf 23} (1977), 291--301.

\bibitem{fuss1791}
Nikolao Fuss, Solutio quaestionis, quot modis polygonum $n$ laterum in polygona
  $m$ laterum, per diagonales resolvi quaeat, {\em Novi Acta Academiae.
  Scientiarum imperialis Petropolitanae} {\bf 9} (1791), 243--251 (in Latin).

\bibitem{gould64}
H.~Gould, Binomial coefficients, the bracket functions, and compositions with
  relatively prime summands, {\em Fibonacci Quart.} {\bf 4} (1964), 241--259.

\bibitem{heubach_li_mansour08}
S.~Heubach, N.~Y. Li, and T.~Mansour, Staircase tilings and $k$-{C}atalan
  structures, {\em Discrete Math.} {\bf 308} (2008), 5954--5964.

\bibitem{hilton_pedersen91}
P.~Hilton and J.~Pedersen, Catalan numbers, their generalization, and their
  uses, {\em Math. Intelligencer} {\bf 13} (1991), 64--75.

\bibitem{kim08}
D.~Kim, On the $(n,k)$-th {C}atalan numbers, {\em Commun. Korean Math. Soc.}
  {\bf 23} (2008), 349--356.

\bibitem{larcombe99}
P.~J. Larcombe, The $18$-th century {C}hinese discovery of the {C}atalan
  numbers, {\em Math. Spectrum} {\bf 32} (1999), 5--7.

\bibitem{luo88}
J.~J. Luo, Antu {M}ing, the first inventor of {C}atalan numbers of the world,
  {\em Neimenggu Daxue Xuebao} {\bf 19} (1988), 239--245 (in Chinese).

\bibitem{motzkin48}
T.~Motzkin, Relations between hypersurface cross ratios, and a combinatorial
  formula for partitions of a polygon, for permanent preponderance, and for
  non-associative products, {\em Bull. Amer. Math. Soc.} {\bf 54} (1948),
  352--360.

\bibitem{pak99}
I.~Pak, Reduced decompositions of permutations in terms of star transpositions,
  generalized {C}atalan numbers and $k$-ary trees, {\em Discrete Math.} {\bf
  204} (1999), 329--335.

\bibitem{stanley15}
R.~Stanley, {\em Catalan Numbers}, Cambridge University Press, 2015.

\bibitem{sun14}
Z.-W. Sun, Congruences involving generalized central trinomial coefficients,
  {\em Sci. China Math.} {\bf 57} (2014), 1375--1400.

\bibitem{wangzhang15}
Y.~Wang and Z.~Zhang, 
\href{https://cs.uwaterloo.ca/journals/JIS/VOL18/Wang/wang21.html}{Combinatorics of generalized {M}otzkin numbers}, {\em J.
  Integer Sequences} {\bf 18} (2015) Article 15.2.4.

\bibitem{zhaoqi17}
J.~Zhao and F.~Qi, Two explicit formulas for the generalized {M}otzkin numbers,
  {\em J. Inequal. Appl.} {\bf 19} (2017).

\end{thebibliography}

\bigskip
\hrule
\bigskip

\noindent 2000 {\it Mathematics Subject Classification}: Primary 11B37;
Secondary 39A14, 65Q30,  05A15.

\noindent \emph{Keywords:}
Catalan number, Motzkin number, generalized Motzkin number, recurrence
relation, explicit formula, finite difference scheme, generating function.

\bigskip
\hrule
\bigskip

\noindent (Concerned with sequences
\seqnum{A000108} and \seqnum{A001006}.)

\bigskip
\hrule
\bigskip

\vspace*{+.1in}
\noindent
Received May 11 2018;
revised version received August 27 2018.
Published in {\it Journal of Integer Sequences}, November 24 2018.

\bigskip
\hrule
\bigskip

\noindent
Return to
\htmladdnormallink{Journal of Integer Sequences home page}{http://www.cs.uwaterloo.ca/journals/JIS/}.
\vskip .1in


\end{document}

</pre></body></html>