Solution of the Differential Equation $y^{(k)}=e^{ay}$, Special Values of Bell
Polynomials and $(k,a)$-Autonomous Coefficients
Ronald Orozco López
Departamento de Matemáticas
Universidad de los Andes
Bogotá, 111711
Colombia
<EMAIL_ADDRESS>

###### Abstract

In this paper special values of Bell polynomials are given by using the power
series solution of the equation $y^{(k)}=e^{ay}$. In addition, complete and
partial exponential autonomous functions, exponential autonomous polynomials,
autonomous polynomials and $(k,a)$-autonomous coefficients are defined.
Finally, we show the relationship between various numbers counting
combinatorial objects and binomial coefficients, Stirling numbers of second
kind and autonomous coefficients.

##  1 Introduction

It is a known fact that Bell polynomials are closely related to derivatives of
composition of functions. For example, Faa di Bruno [5], Foissy [6], and
Riordan [10] showed that Bell polynomials are a very useful tool in
mathematics to represent the $n$-th derivative of the composition of
functions. Also, Bernardini and Ricci [2], Yildiz et al. [12], Caley [3], and
Wang [13] showed the relationship between Bell polynomials and differential
equations. On the other hand, Orozco [9] studied the convergence of the
analytic solution of the autonomous differential equation $y^{(k)}=f(y)$ by
using Faa di Bruno’s formula. We can then look at differential equations as a
source for investigating special values of Bell polynomials.

In this paper we will focus on finding special values of Bell polynomials when
the vector field $f(x)$ of the autonomous differential equation $y^{(k)}=f(y)$
is the exponential function. We will not consider the convergence of the
solutions, but we will show that well known numbers such as reduced tangent
numbers, Bernoulli numbers, Euler zigzag numbers, Blasius numbers, among
others, can be constructed using Bell polynomials. In general, a special class
of numbers, which have not yet been studied, are constructed using Bell
polynomials. On the other hand, a new family of numbers, called
$(k,a)$-autonomous coefficients, is obtained for each value of $k$. Four
conjectures about these numbers are established.

This paper is divided as follows. We begin with a summary of results on
complete and partial Bell polynomials, which will be used to demonstrate the
main results presented here. Next, we introduce the complete and partial
exponential autonomous functions, the recurrence relations of these are
constructed using Bell polynomials, and some recurrence relations of solutions
of various initial value problems are given. In the fourth section the
$(k,a)$-autonomous coefficients are introduced. From these numbers we can
obtain the triangular numbers, the 8-sequence numbers of $[1,n]$ with 2
contiguous pairs, among others. We finish this work by studying the cases
$k=2,3,4$ for the autonomous differential equation $y^{(k)}=e^{ay}$.

##  2 Bell exponential polynomials

The following basic results can be found at Comtent [4], and Riordan [11].
Exponential Bell polynomials are used to encode information on the ways in
which a set can be partitioned, hence they are a very useful tool in
combinatorial analysis. Bell polynomials are obtained from the derivatives of
composite functions and are given by Faa Di Bruno’s formula [5]. Bell [1],
Gould [7] and Mihoubi [8] provided important results on these polynomials. We
start with the definition of partial Bell polynomials

######  Definition 1.

The exponential partial Bell polynomials are the polynomials

$B_{n,k}(x_{1},x_{n},...,x_{n-k+1})$

in the variables $x_{1},x_{2},...$ defined by the series expansion

$\mathop{\mbox{{exp}}}\nolimits\left(u\sum_{j=1}^{\infty}x_{j}\frac{t^{j}}{j!}\right)=1+\sum_{n=1}^{\infty}\frac{t^{n}}{n!}\sum_{k=1}^{n}u^{k}B_{n,k}(x_{1},x_{2},...,x_{n-k+1}).$
(1)

The following result gives the explicit way to calculate the partial Bell
polynomials

######  Theorem 2.

The partial or incomplete exponential Bell polynomials are given by

$B_{n,k}(x_{1},\ldots,x_{n-k+1})=\sum\frac{n!}{c_{1}!c_{2}!\cdots
c_{n-k+1}!}\left(\frac{x_{1}}{1!}\right)^{c_{1}}\cdots\left(\frac{x_{n-k+1}}{(n-k+1)!}\right)^{c_{n-k+1}}$

where the summation takes place over all integers
$c_{1},c_{2},...,c_{n-k+1}\geq 0$, such that

$\displaystyle c_{1}+2c_{2}+\cdots+(n-k+1)c_{n-k+1}$ $\displaystyle=$
$\displaystyle n$ $\displaystyle c_{1}+c_{2}+\cdots+c_{n-k+1}$
$\displaystyle=$ $\displaystyle k$

The following are special cases of partial Bell polynomials and will be very
useful for proving results in this paper

$\displaystyle B_{n,1}(x_{1},...,x_{n})$ $\displaystyle=x_{n},$ (2)
$\displaystyle B_{n,2}(x_{1},...,x_{n-1})$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{n-1}\binom{n}{k}x_{k}x_{n-k},$ (3)
$\displaystyle B_{n,n-a}(x_{1},...,x_{a+1})$
$\displaystyle=\sum_{j=a+1}^{2a}\frac{j!}{a!}\binom{n}{j}x_{1}^{n-j}B_{a,j-a}\left(\frac{x_{2}}{2},\cdots,\frac{x_{2(a+1)-j}}{2(a+1)-j}\right),$
$\displaystyle\ 1\leq a<n,$ (4) $\displaystyle B_{n,n}(x_{1})$
$\displaystyle=x_{1}^{n},$ (5) $\displaystyle B_{n,n-1}(x_{1},x_{2})$
$\displaystyle=\binom{n}{2}x_{1}^{n-2}x_{2},$ (6) $\displaystyle
B_{n,n-2}(x_{1},x_{2},x_{3})$
$\displaystyle=\binom{n}{3}x_{1}^{n-3}x_{3}+3\binom{n}{4}x_{1}^{n-4}x_{2}^{2},$
(7) $\displaystyle B_{n,n-3}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle=\binom{n}{4}x_{1}^{n-4}x_{4}+10\binom{n}{5}x_{1}^{n-5}x_{2}x_{3}+15\binom{n}{6}x_{1}^{n-6}x_{2}^{3},$
(8) $\displaystyle B_{n,n-4}(x_{1},x_{2},x_{3},x_{4},x_{5})$
$\displaystyle=\binom{n}{5}x_{1}^{n-5}x_{5}+5\binom{n}{6}x_{1}^{n-6}[3x_{2}x_{4}+2x_{3}^{2}]$
$\displaystyle+105\binom{n}{7}x_{1}^{n-7}x_{2}^{2}x_{3}+105\binom{n}{8}x_{1}^{n-8}x_{2}^{4}.$
(9)

Some values of partial Bell polynomials are

$\displaystyle B_{n,k}(0!,1!,...,(n-k)!)$
$\displaystyle=\genfrac{[}{]}{0.0pt}{}{n}{k}\ \ (\textit{Stirling number of
first kind}),$ $\displaystyle B_{n,k}(1!,...,(n-k)!)$
$\displaystyle=\binom{n-1}{k-1}\frac{n!}{k!}\ \ (\textit{Lah number}),$
$\displaystyle B_{n,k}(1,1,...,1)$
$\displaystyle=\genfrac{\\{}{\\}}{0.0pt}{}{n}{k}\ \ (\textit{Stirling number
of second kind}),$ $\displaystyle B_{n,k}(1,2,...,n-k+1)$
$\displaystyle=\binom{n}{k}k^{n-k}\ \ (\textit{Idempotent number}).$

Then we can see the beautiful relationship that exists between Bell
polynomials and numbers like the above.

On the other hand, the partial Bell polynomials can be efficiently computed by
means of the recurrence relation

$B_{n,k}(x_{1},...,x_{n-k+1})=\sum_{i=1}^{n-k+1}\binom{n-1}{i-1}x_{i}B_{n-i,k-1}(x_{1},...,x_{n-i-k+2}).$
(10)

The definition of complete Bell polynomials is as follows

######  Definition 3.

The sum

$B_{n}(x_{1},x_{2},...,x_{n})=\sum_{k=1}^{n}B_{n,k}(x_{1},x_{2},...,x_{n-k+1})$
(11)

is called $n$-th complete exponential Bell polynomials with exponential
generating function given by to make $u=1$ in (1)

$\mathop{\mbox{{exp}}}\nolimits\left(\sum_{m=1}^{\infty}x_{m}\frac{t^{m}}{m!}\right)=\sum_{n=0}^{\infty}B_{n}(x_{1},x_{2},...,x_{n})\frac{t^{n}}{n!}$
(12)

and $B_{0}=1$.

Some complete Bell polynomials are

$\displaystyle B_{1}(x_{1})$ $\displaystyle=x_{1},$ $\displaystyle
B_{2}(x_{1},x_{2})$ $\displaystyle=x_{1}^{2}+x_{2},$ $\displaystyle
B_{3}(x_{1},x_{2},x_{3})$ $\displaystyle=x_{1}^{3}+3x_{1}x_{2}+x_{3},$
$\displaystyle B_{4}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle=x_{1}^{4}+6x_{1}^{2}x_{2}+4x_{1}x_{3}+3x_{2}^{2}+x_{4},$
$\displaystyle B_{5}(x_{1},x_{2},x_{3},x_{4},x_{5})$
$\displaystyle=x_{1}^{5}+10x_{1}^{3}x_{2}+15x_{1}x_{2}^{2}+10x_{1}^{2}x_{3}+10x_{2}x_{3}+5x_{1}x_{4}+x_{5}.$

######  Theorem 4.

The complete Bell polynomials $B_{n}$ satisfy the identity

$B_{n+1}(x_{1},...,x_{n+1})=\sum_{i=0}^{n}\binom{n}{i}B_{n-i}(x_{1},...,x_{n-i})x_{i+1}.$
(13)

From this it follows that

$B_{2n+1}(0,x_{2},0,..,0,x_{2n+1})=0$ (14)

for all $n\geq 0$.

Another useful identity that Bell polynomials fulfill is as follows

$B_{n}(-x_{1},x_{2},-x_{3},...,(-1)^{n-1}x_{n})=(-1)^{n}B_{n}(x_{1},x_{2},x_{3},...,x_{n})$
(15)

##  3 Exponential autonomous functions

We will study the solution of the equation

$y^{(k)}=e^{ay}$ (16)

for any $a\in\mathbb{C}$. Making $y=u/a$ we obtain the equivalent equation

$u^{(k)}=ae^{u}$ (17)

Then without loss of generality we will focus on the equation (17). Now by
applying derivative to (17) we obtain another equation equivalent to (16)

$\displaystyle u^{(k+1)}$ $\displaystyle=$ $\displaystyle
ae^{u}u^{\prime}=u^{(k)}u^{\prime}$ (18)

Denote $(x_{1},x_{2},...,x_{k})$ the initial value problem $y(0)=x_{1}$,
$y^{\prime}(0)=x_{2}$,…,$y^{(k-1)}(0)=x_{k}$. In this section the general
solution and solutions with initial values $(x,0,0,...,0)$, $(x_{1}+k\ln
c,cx_{2},...,c^{k-1}x_{k})$, and $(x_{1},-x_{2},...,x_{2k-1},-x_{2k})$ of the
equation (17) are given. We will define the complete and partial exponential
autonomous functions and the exponential autonomous polynomials, which are the
coefficients of the power series solution of the equation (17). Moreover, we
will find special values of these functions by using Bell polynomials. We
begin with the following definition

######  Definition 5.

Take $a\in\mathbb{C}$. Suppose $\mathbf{x}=(x_{1},...,x_{k})$. Let
$f_{n}(\mathbf{x},a)$ denote the complete exponential autonomous functions of
order $k$, $k\geq 1$, recursively defined as

$\displaystyle f_{0}(\mathbf{x},a)$ $\displaystyle=x_{1},$ (19) $\displaystyle
f_{1}(\mathbf{x},a)$ $\displaystyle=x_{2},$ (20) $\displaystyle\vdots$
$\displaystyle f_{k-1}(\mathbf{x},a)$ $\displaystyle=x_{k},$ (21)
$\displaystyle f_{k}(\mathbf{x},a)$ $\displaystyle=ae^{x_{1}},$ (22)
$\displaystyle f_{n+k}(\mathbf{x},a)$
$\displaystyle=ae^{x_{1}}B_{n}(f_{1}(\mathbf{x},a),...,f_{n}(\mathbf{x},a)),\
\ n\geq 1,$ (23)

where $B_{n}(y_{1},...,y_{n})$ are the complete Bell polynomials. When
$x_{1}=0$, we define the exponential autonomous polynomials as
$q_{n}(x_{2},...,x_{k})=f_{n}(0,x_{2},...,x_{k})$, for $n\geq 1$.

When $a=1$ in the above definition, we will write
$f_{n}(\mathbf{x})=f_{n}(\mathbf{x},1)$. In this section we will restrict
ourselves to exponential autonomous functions and autonomous polynomials will
be dealt with in the next section.

The following are complete exponential autonomous functions for $k=1,2,3,4$.
They will be very useful in the next section:
When $k=1$, $f_{n}(x,a)=(n-1)!a^{n}e^{nx}$. When $k=2$

$\displaystyle f_{0}(x,y,a)$ $\displaystyle=x,$ $\displaystyle f_{1}(x,y,a)$
$\displaystyle=y,$ $\displaystyle f_{2}(x,y,a)$ $\displaystyle=ae^{x},$
$\displaystyle f_{3}(x,y,a)$ $\displaystyle=aye^{x},$ $\displaystyle
f_{4}(x,y,a)$ $\displaystyle=ae^{x}(ae^{x}+y^{2}),$ $\displaystyle
f_{5}(x,y,a)$ $\displaystyle=ae^{x}(4aye^{x}+y^{3}),$ $\displaystyle
f_{6}(x,y,a)$ $\displaystyle=ae^{x}(4a^{2}e^{2x}+11ay^{2}e^{x}+y^{4}).$

When $k=3$

$\displaystyle f_{0}(x,y,z,a)$ $\displaystyle=x,$ $\displaystyle
f_{1}(x,y,z,a)$ $\displaystyle=y,$ $\displaystyle f_{2}(x,y,z,a)$
$\displaystyle=z,$ $\displaystyle f_{3}(x,y,z,a)$ $\displaystyle=ae^{x},$
$\displaystyle f_{4}(x,y,z,a)$ $\displaystyle=aye^{x},$ $\displaystyle
f_{5}(x,y,z,a)$ $\displaystyle=ae^{x}(z+y^{2}),$ $\displaystyle
f_{6}(x,y,z,a)$ $\displaystyle=ae^{x}(ae^{x}+3yz+y^{3}),$ $\displaystyle
f_{7}(x,y,z,a)$ $\displaystyle=ae^{x}(5ye^{x}+3z^{2}+6y^{2}z+y^{4}).$

And finally, when $k=4$

$\displaystyle f_{0}(x,y,z,w,a)$ $\displaystyle=x,$ $\displaystyle
f_{1}(x,y,z,w,a)$ $\displaystyle=y,$ $\displaystyle f_{2}(x,y,z,w,a)$
$\displaystyle=z,$ $\displaystyle f_{3}(x,y,z,w,a)$ $\displaystyle=w,$
$\displaystyle f_{4}(x,y,z,w,a)$ $\displaystyle=ae^{x},$ $\displaystyle
f_{5}(x,y,z,w,a)$ $\displaystyle=aye^{x},$ $\displaystyle f_{6}(x,y,z,w,a)$
$\displaystyle=ae^{x}(z+y^{2}),$ $\displaystyle f_{7}(x,y,z,w,a)$
$\displaystyle=ae^{x}(w+3yz+y^{3}),$ $\displaystyle f_{8}(x,y,z,w,a)$
$\displaystyle=ae^{x}(ae^{x}+3z^{2}+4yw+6y^{2}z+y^{4}).$

In the following result we will show that the exponential generating function
of the complete exponential autonomous functions is solution of the equation
(17)

######  Theorem 6.

Let $\textbf{x}=(x_{1},...,x_{k})$. The series

$E_{k}(t,\textbf{x},a)=\sum_{n=0}^{\infty}f_{n}(\textbf{x},a)\frac{t^{n}}{n!}$
(24)

is solution of the differential equation (17).

###### Proof.

Taking derivative $k$ times with respect to $t$ the series
$E_{k}(t,\textbf{x},a)$, using the definition of the autonomous functions
$f_{n}(\textbf{x},a)$ and the equation (1), then

$\displaystyle\frac{\partial^{k}E_{k}(t,\textbf{x},a)}{\partial t^{k}}$
$\displaystyle=\sum_{n=0}^{\infty}f_{n+k}(\textbf{x},a)\frac{t^{n}}{n!}$
$\displaystyle=ae^{x_{1}}+\sum_{n=1}^{\infty}ae^{x_{1}}B_{n}(f_{1}(\mathbf{x},a),...,f_{n}(\mathbf{x},a))\frac{t^{n}}{n!}$
$\displaystyle=e^{ax_{1}}\left(1+\sum_{n=1}^{\infty}B_{n}(f_{1}(\mathbf{x},a),...,f_{n}(\mathbf{x},a))\frac{t^{n}}{n!}\right)$
$\displaystyle=ae^{x_{1}}e^{E_{k}(t,\textbf{x},a)-x_{1}}$
$\displaystyle=ae^{E_{k}(t,\textbf{x},a)}.$

∎

Now we define the partial exponential autonomous functions

######  Definition 7.

Let $g_{n,i}(\mathbf{x},a)$ denote the partial exponential autonomous
functions as

$g_{n,i}(\mathbf{x},a)=B_{n,i}(f_{1}(\mathbf{x},a),...,f_{n-i+1}(\mathbf{x},a))$
(25)

with $g_{0,0}(\mathbf{x},a)=1$, $g_{n,0}(\mathbf{x},a)=0$, for $n\geq 1$, and
$g_{0,i}(\mathbf{x},a)=0$ , for $i\geq 1$. Then

$f_{n+k}(\mathbf{x},a)=ae^{x_{1}}\sum_{i=1}^{n}g_{n,i}(\mathbf{x},a).$ (26)

In the following result we establish recurrence relations for the functions
$f_{n}(\mathbf{x},a)$ and $g_{n,i}(\mathbf{x},a)$. Many important results of
this paper will be proved using this theorem.

######  Theorem 8.

The autonomous functions $f_{n}(\textbf{x},a)$ and $g_{n,i}(\mathbf{x},a)$
fulfill the recurrence relations

$f_{n+k+1}(\textbf{x},a)=\sum_{i=0}^{n}\binom{n}{i}f_{n-i+k}(\textbf{x},a)f_{i+1}(\textbf{x},a)$
(27)

and

$g_{n,i}(\textbf{x},a)=\sum_{j=1}^{n-i+1}\binom{n-1}{j-1}f_{j}(\textbf{x},a)g_{n-j,i-1}(\textbf{x},a).$
(28)

###### Proof.

Making $y_{j}=f_{j}(\textbf{x},a)$ in (10) and (13) and multiplying these by
$ae^{x_{1}}$, we obtain the desired result. ∎

Now we will study the behavior of the functions $f_{n}(\mathbf{x},a)$
evaluated at $\mathbf{x}=(x_{1},0,0,...,0)$. From previous result we can
construct the first important sequence arising from the differential equation
(17)

######  Theorem 9.

The functions $f_{n}(\mathbf{x},a)$ take the following values at
$\mathbf{x}=(x_{1},0,...,0)$

  1. 1.

$f_{kn+1}(x_{1},0,...,0,a)=f_{kn+2}(x_{1},0,...,0,a)=\cdots=f_{kn+k-1}(x_{1},0,...,0,a)=0$,
$n\geq 0$,

  2. 2.

$f_{kn}(x_{1},0,...,0,a)=A_{n}^{(k)}(a)e^{nx_{1}}$, $n\geq 1$

where $A^{(k)}_{1}(a)=1$ and

$A^{(k)}_{n+2}(a)=\sum_{i=0}^{n}\binom{kn+k-1}{ki+k-1}A^{(k)}_{n-i+1}(a)A^{(k)}_{i+1}(a),\
\ $ (29)

$n\geq 0$, $k\geq 1$.

###### Proof.

Let $\mathbf{x}=(x_{1},0,...,0)$. Clearly, $f_{1}(\mathbf{x},a)=0$,
$f_{k+1}(\mathbf{x},a)=B_{1}(f_{1}(\mathbf{x},a))=f_{1}(\mathbf{x},a)=0$. Now
suppose it is true that $f_{ki+1}(\mathbf{x},a)=0$ for $2\leq i\leq n-1$. By
the theorem 8

$\displaystyle f_{kn+1}(\mathbf{x},a)$
$\displaystyle=f_{k(n-1)+k+1}(\mathbf{x},a)$
$\displaystyle=\sum_{i=0}^{k(n-1)}\binom{k(n-1)}{i}f_{kn-i}(\mathbf{x},a)f_{i+1}(\mathbf{x},a).$

Since the product $f_{kn-i}(\mathbf{x},a)f_{i+1}(\mathbf{x},a)$ contain the
functions $f_{kj+1}(\mathbf{x},a)$, then $f_{kn+1}(\mathbf{x},a)=0$ for all
$n$. Likewise it is proved for $f_{kn+j}(\mathbf{x},a)$, $j=2,...,k-1$. Now we
will prove $2$. We know that $f_{k}(\mathbf{x},a)=ae^{x_{1}}$,
$f_{2k}(\mathbf{x},a)=a^{2}e^{2x_{1}}$ and suppose that
$f_{kn}(\mathbf{x},a)=A_{n}^{(k)}(a)e^{nx_{1}}$. Then

$\displaystyle f_{kn+k}(\mathbf{x},a)$
$\displaystyle=f_{(kn-1)+k+1}(\mathbf{x},a)$
$\displaystyle=\sum_{i=0}^{kn-1}\binom{kn-1}{i}f_{kn-1-i+k}(\mathbf{x},a)f_{i+1}(\mathbf{x},a)$
$\displaystyle=\sum_{i=0}^{n-1}\binom{kn-1}{ki+1}f_{k(n-i)}(\mathbf{x},a)f_{ki+k}(\mathbf{x},a)$
$\displaystyle=\sum_{i=0}^{n-1}\binom{kn-1}{ki+1}A_{n-i}^{(k)}(a)e^{(n-i)x_{1}}A_{i+1}^{(k)}(a)e^{(i+1)x_{1}}$
$\displaystyle=e^{(n+1)x_{1}}\sum_{i=0}^{n-1}\binom{kn-1}{ki+1}A_{n-i}^{(k)}(a)A_{i+1}^{(k)}(a)$
$\displaystyle=e^{(n+1)x_{1}}A_{n+1}^{(k)}(a).$

∎

It is easy to show that $A_{n}^{(1)}(a)=(n-1)!a^{n}$ when $k=1$. We will use
the equation (29) to prove this result. Suppose it is true that
$A_{i}^{(1)}(a)=(i-1)!a^{i}$ for $i$ ranging between $1$ and $n+1$. We have

$\displaystyle A^{(1)}_{n+2}(a)$
$\displaystyle=\sum_{i=0}^{n}\binom{n}{i}A^{(1)}_{n-i+1}(a)A^{(1)}_{i+1}(a)$
$\displaystyle=\sum_{i=0}^{n}\binom{n}{i}(n-i)!a^{n-i+1}i!a^{i+1}$
$\displaystyle=a^{n+2}\sum_{i=0}^{n}\binom{n}{i}(n-i)!i!$
$\displaystyle=a^{n+2}n!(n+1)=a^{n+2}(n+1)!.$

We can extend the above result to all $k\geq 1$

######  Proposition 1.

For all $k\geq 1$ we have

$A_{n}^{(k)}(a)=a^{n}A_{n}^{(k)}(1)$ (30)

###### Proof.

Suppose by induction that $A_{i}^{(k)}(a)=a^{i}A_{i}^{(k)}(1)$ for all $i$
ranging between $1$ and $n+1$, then use the same steps as in the previous
proof. ∎

From the above proposition it follows that

$E_{k}(t,(x,0,...,0),a)=E_{k}(at,(x,0,...,0),1).$ (31)

Then without loss of generality it is sufficient to study the solution of (17)
with initial conditions $y(0)=x$,
$y^{\prime}(0)=y^{\prime}(0)=\cdots=y^{(k-1)}(0)=0$ and $a=1$ to generate the
sequence $A_{n}^{(k)}(1)$.

The following corollary of theorem 9 shows us that the numbers $A_{n}^{(k)}$
can be constructed using Bell polynomials

######  Corollary 1.

Numbers $A_{n}^{(k)}(a)$ fulfill the recurrence relation

$A_{n}^{(k)}(a)=B_{n-k}(\overbrace{0,...,0}^{k-1},A_{1}^{(k)}(a),...,\overbrace{0,...,0}^{k-1},A_{n-k}^{(k)}(a)),\
\ n\geq 1.$ (32)

In the following theorem we calculate some special values of the functions
$g_{n.i}(\mathbf{x},a)$

######  Theorem 10.

Let $\mathbf{x}=(x_{1},0,...,0)$. For all $a\in\mathbb{R}$ we have

  1. 1.

$g_{n,i}(\mathbf{x},a)=0$, si $k\not|n$.

  2. 2.

$g_{lk,1}(\mathbf{x},a)=A_{l}^{(k)}(a)e^{lx_{1}}$.

  3. 3.

$g_{lk,2}(\mathbf{x},a)=e^{lx_{1}}\sum_{j=1}^{l}\binom{kl-1}{kj-1}A_{j}^{(k)}(a)A_{l-j}^{(k)}(a)$.

  4. 4.

$g_{n,n}(\mathbf{x},a)=g_{n,n-1}(\mathbf{x},a)=g_{n,n-2}(\mathbf{x},a)=g_{n,n-3}(\mathbf{x},a)=g_{n,n-4}(\mathbf{x},a)=0$,
$k>1$.

###### Proof.

Suppose $k\not|n$ and $g_{n-j,i-1}(\mathbf{x},a)=0$ for all $k$ such that
$k\not|j$. Using theorem 8 and 9 we have that $f_{j}(\mathbf{x},a)=0$. This
proves 1. To prove 2 we have

$\displaystyle g_{lk,1}(\mathbf{x},a)$
$\displaystyle=\sum_{j=1}^{lk}\binom{lk-1}{j-1}f_{j}(\mathbf{x},a)g_{lk-j,0}(\mathbf{x},a)$
$\displaystyle=\binom{lk-1}{lk-1}f_{lk}(\mathbf{x},a)g_{0,0}(\mathbf{x},a)$
$\displaystyle=A_{l}^{(k)}(a)e^{lx_{1}}.$

On the other hand,

$\displaystyle g_{lk,2}(\mathbf{x},a)$
$\displaystyle=\sum_{j=1}^{lk-1}\binom{lk-1}{j-1}f_{j}(\mathbf{x},a)g_{lk-j,1}(\mathbf{x},a)$
$\displaystyle=\sum_{j=1}^{l}\binom{lk-1}{jk-1}f_{kj}(\mathbf{x},a)g_{lk-
kj,1}(\mathbf{x},a)$
$\displaystyle=\sum_{j=1}^{l}\binom{lk-1}{jk-1}A_{j}^{(k)}(a)e^{jx_{1}}A_{l-j}^{(k)}(a)e^{(l-j)x_{1}}$
$\displaystyle=e^{lx_{1}}\sum_{j=1}^{l}\binom{lk-1}{jk-1}A_{j}^{(k)}(a)A_{l-j}^{(k)}(a).$

Then this proves 3. To prove 4 we use the equations (5)-(9). ∎

We conclude this section with the following properties of the exponential
autonomous functions

######  Theorem 11.

For all $n\geq 1$, $k\geq 1$ and for all $a,c\in\mathbb{C}$ is fulfilled

$f_{n}(x_{1}+k\ln
c,cx_{2},...,c^{k-1}x_{k},a)=c^{n}f_{n}(x_{1},x_{2},...,x_{k},a)$ (33)

###### Proof.

Let $\mathbf{y}=(x_{1}+k\ln c,cx_{2},...,c^{k-1}x_{k})$ and
$\mathbf{x}=(x_{1},x_{2},...,x_{k})$. Suppose that the result is true for
$i\leq n$. Then

$\displaystyle f_{n+1}(\mathbf{y},a)$
$\displaystyle=f_{(n+1-k)+k}(\mathbf{y},a)$
$\displaystyle=ac^{k}e^{x_{1}}B_{n+1-k}(f_{1}(\mathbf{y},a),...,f_{n+1-k}(\mathbf{y},a))$
$\displaystyle=ac^{k}e^{x_{1}}B_{n+1-k}(cf_{1}(\mathbf{x},a),...,c^{n+1-k}f_{n+1-k}(\mathbf{x},a))$
$\displaystyle=ac^{k}c^{n+1-k}e^{x_{1}}B_{n}(f_{1}(\mathbf{x},a),...,f_{n}(\mathbf{x},a))$
$\displaystyle=c^{n+1}f_{n+1}(\mathbf{x},a).$

∎

The following is the corollary to the theorem 11 that allows to calculate the
solutions of (17) when the initial values are $(x_{1}+k\ln
c,cx_{2},...,c^{k-1}x_{k})$

######  Corollary 2.

$E_{k}(t,(x_{1}+k\ln c,cx_{2},...,c^{k-1}x_{k}),a)=k\ln
c+E_{k}(ct,(x_{1},x_{2},...,x_{k}),a)$ (34)

###### Proof.

Let $\mathbf{y}=(x_{1}+k\ln c,cx_{2},...,c^{k-1}x_{k})$ and
$\mathbf{x}=(x_{1},x_{2},...,x_{k})$. From the above theorem and the
definition of the function $E_{k}(t,\mathbf{x},a)$ we have

$\displaystyle E_{k}(t,\mathbf{y},a)$ $\displaystyle=x_{1}+k\ln
c+\sum_{n=1}^{\infty}f_{n}(\mathbf{y},a)\frac{t^{n}}{n!}$
$\displaystyle=x_{1}+k\ln
c+\sum_{n=1}^{\infty}c^{n}f_{n}(\mathbf{x},a)\frac{t^{n}}{n!}$
$\displaystyle=x_{1}+k\ln
c+\sum_{n=1}^{\infty}f_{n}(\mathbf{x},a)\frac{(ct)^{n}}{n!}$
$\displaystyle=k\ln c+E_{k}(ct,\mathbf{x},a).$

∎

Finally, we compute $f_{n}(\mathbf{x},a)$ when
$\mathbf{x}=(x_{1},-x_{2},...,x_{2k-1},-x_{2k})$

######  Theorem 12.

For all $n\geq 0$ and for all exponential autonomous functions of order $2k$
it is satisfied that

$f_{n}((x_{1},-x_{2},...,x_{2k-1},-x_{2k}),a)=(-1)^{n}f_{n}((x_{1},x_{2},...,x_{2k-1},x_{2k}),a)$
(35)

###### Proof.

Let $\mathbf{y}=(x_{1},-x_{2},...,x_{2k-1},-x_{2k})$, and let
$\mathbf{x}=(x_{1},x_{2},...,x_{2k-1},x_{2k})$. Suppose it is true for all
values less than or equal to $n$. Then

$\displaystyle f_{n+1}(\mathbf{y},a)$ $\displaystyle=$ $\displaystyle
f_{(n+1-2k)+2k}(\mathbf{y},a)$ $\displaystyle=$ $\displaystyle
ae^{x_{1}}B_{n+1-2k}(f_{1}(\mathbf{y},a),f_{2}(\mathbf{y},a),...,f_{n+1-2k}(\mathbf{y},a))$
$\displaystyle=$ $\displaystyle
ae^{x_{1}}B_{n+1-2k}(-f_{1}(\mathbf{x},a),f_{2}(\mathbf{x},a),...,(-1)^{n+1-2k}f_{n+1-2k}(\mathbf{x},a))$
$\displaystyle=$
$\displaystyle(-1)^{n+1-2k}ae^{x_{1}}B_{n+1-2k}(f_{1}(\mathbf{x},a),f_{2}(\mathbf{x},a),...,f_{n+1-2k}(\mathbf{x},a))$
$\displaystyle=$ $\displaystyle(-1)^{n+1}f_{n}(\mathbf{x},a).$

∎

Finally we have the corollary to theorem 12

######  Corollary 3.

$E_{2k}(-t,(x_{1},-x_{2},...,x_{2k-1},-x_{2k}),a)=E_{2k}(t,(x_{1},x_{2},...,x_{2k-1},x_{2k}),a).$
(36)

###### Proof.

From the above theorem and the definition of the function $E_{2k}(t,x,a)$ we
have

$\displaystyle E_{2k}(-t,\mathbf{y},a)$
$\displaystyle=\sum_{n=0}^{\infty}f_{n}(\mathbf{y},a)\frac{(-t)^{n}}{n!}$
$\displaystyle=\sum_{n=0}^{\infty}(-1)^{n}(-1)^{n}f_{n}(\mathbf{x},a)\frac{t^{n}}{n!}$
$\displaystyle=\sum_{n=0}^{\infty}f_{n}(\mathbf{x},a)\frac{t^{n}}{n!}$
$\displaystyle=E_{2k}(t,\mathbf{x},a)$

where $\mathbf{y}=(x_{1},-x_{2},...,x_{2k-1},-x_{2k})$ and
$\mathbf{x}=(x_{1},x_{2},...,x_{2k-1},x_{2k})$. ∎

##  4 (k,a)-autonomous coefficients

When $k=1$ we obtain the equation $y^{\prime}=ae^{y}$, which is the easiest to
solve for all $k$. Using the method of separation of variables we reach the
solution

$y(t)=-\ln(e^{-x}-at)$

with initial condition $y(0)=x$. On the other hand, by the theorem 6 the
solution in power series becomes

$\displaystyle E_{1}(t,x,a)$
$\displaystyle=x+\sum_{n=1}^{\infty}A_{n}^{(1)}(a)\frac{t^{n}}{n!}$
$\displaystyle=x+\sum_{n=1}^{\infty}(n-1)!a^{n}e^{nx}\frac{t^{n}}{n!}$
$\displaystyle=x+\sum_{n=1}^{\infty}\frac{(ae^{x})^{n}}{n}=x-\ln(1-ae^{x}t)$

Now we can use the results of the previous section to prove some results
already known. By the definition of complete exponential autonomous functions

$\displaystyle n!a^{n+1}e^{a(n+1)x}$
$\displaystyle=ae^{x}\sum_{i=1}^{n}B_{n,i}(0!a^{1}e^{x},1!a^{2}e^{2x},...,(n-i)!a^{n-i+1}e^{(n-i+1)x})$
$\displaystyle=ae^{x}a^{n}e^{nx}\sum_{i=1}^{n}B_{n,i}(0!,1!,...,(n-i)!)$
$\displaystyle=e^{a(n+1)x}a^{n+1}\sum_{i=1}^{n}B_{n,i}(0!,1!,...,(n-i)!)=e^{a(n+1)x}a^{n+1}\sum_{i=1}^{n}\genfrac{[}{]}{0.0pt}{}{n}{i}$

from which follows the result relating factorials and Stirling number of first
kind

$n!=\sum_{i=1}^{n}\genfrac{[}{]}{0.0pt}{}{n}{i}.$ (37)

Furthermore $g_{n,i}(\mathbf{x},1)=\genfrac{[}{]}{0.0pt}{}{n}{i}$ and by the
equation (28) we obtain the following finite-sum identity

$\displaystyle\genfrac{[}{]}{0.0pt}{}{n+1}{i+1}$
$\displaystyle=\sum_{j=1}^{n-i+1}\binom{n}{j-1}(j-1)!\genfrac{[}{]}{0.0pt}{}{n+1-j}{i}$
$\displaystyle=\sum_{j=n}^{i}\frac{n!}{j!}\genfrac{[}{]}{0.0pt}{}{j}{i}$
$\displaystyle=\sum_{j=i}^{n}\frac{n!}{j!}\genfrac{[}{]}{0.0pt}{}{j}{i}$
$\displaystyle=\sum_{j=0}^{n}\frac{n!}{j!}\genfrac{[}{]}{0.0pt}{}{j}{i}.$

On the other hand, from the equation (27) we obtain the trivial result

$\displaystyle(n+1)!$ $\displaystyle=\sum_{i=0}^{n}\binom{n}{i}(n-i)!i!$
$\displaystyle=\sum_{i=0}^{n}n!.$

The Stirling numbers of the first kind originally arose algebraically from the
expansion of the falling factorial

$(x)_{n}=x(x-1)(x-2)\cdots(x-n+1)$

and in polynomial form is as follows

$(x)_{n}=\sum_{i=0}^{n}(-1)^{n-k}\genfrac{[}{]}{0.0pt}{}{n}{i}x^{i}.$

Analogously, we want to define and study the coefficients of the expansion of
the autonomous exponential polynomials $q_{n}(\mathbf{x},a)$ with
$\mathbf{x}=(0,x,x,...,x)$. First we calculate the degree of
$q_{n}(\mathbf{x},a)$.

######  Proposition 2.

Let $\mathbf{x}=(0,x,x,...,x)$. Then the degree $q_{n}(\mathbf{x},a)$ of

$\textsf{gr}(q_{n}(\mathbf{x},a))=n-k,\ \ \ n\geq k.$ (38)

###### Proof.

By definition

$\displaystyle q_{n+k}(\mathbf{x},a)$
$\displaystyle=\sum_{i=1}^{n-1}a^{i}g_{n,i}(\mathbf{x},a)+a^{n}g_{n,n}(\mathbf{x},a)$
$\displaystyle=\sum_{i=1}^{n-1}a^{i}g_{n,i}(\mathbf{x},a)+a^{n}x_{1}^{n}.$

As $\textsf{gr}(g_{n,i}(\mathbf{x},a))\leq i$, then
$\textsf{gr}(q_{n+k}(\mathbf{x},a))=n$. ∎

We now define the autonomous polynomials and autonomous coefficients

######  Definition 13.

Let $A_{n}^{(k)}(x,a)=q_{n}(0,x,...,x,a)$ denote the autonomous polynomials of
degree $n-k$ for all $n\geq k$.

Using (23) we note that

$A_{n+k}^{(k)}(x,a)=aB_{n}(A_{1}^{(k)}(x,a),...,A_{n}^{(k)}(x,a))$ (39)

for all $n\geq 1$.

######  Definition 14.

We define the $(k,a)$-autonomous coefficients, denoted by
$\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{i}_{(k,a)}$, as the coefficients
of the autonomous polynomials $A_{n+k}^{(k)}(x,a)$, i.e.,

$A_{n+k}^{(k)}(x,a)=\sum_{i=0}^{n}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{i}_{(k,a)}x^{i}.$
(40)

Now we will give some values of the $(k,a)$-autonomous coefficients

######  Theorem 15.

Some values of the coefficients
$\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{i}_{(k,a)}$ are

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{0}_{(k,a)}$
$\displaystyle=\begin{cases}0,&\text{ si $k\not|n$;}\\\
a^{n/k}A_{n/k}^{(k)}(1),&\text{ si $k|n$}\end{cases}$ (41)
$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{0}{i}_{(k,a)}$
$\displaystyle=0,\text{ si $i\geq 1$},$ (42)
$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-l}_{(k,a)}$
$\displaystyle=a\genfrac{\\{}{\\}}{0.0pt}{}{n}{n-l},\ k>l+1,\ 0\leq l<n,$ (43)

###### Proof.

The equation (41) follows from theorem 9. By definition, $B_{0,i}=0$ for
$i\geq 1$. Then (42) is true. Finally, if $k>l+1$,

$\displaystyle B_{n,n-l}(A_{1}^{(k)}(x,a),...,A_{l+1}^{(k)}(x,a))$
$\displaystyle=B_{n,n-l}(x,...,x)$
$\displaystyle=\genfrac{\\{}{\\}}{0.0pt}{}{n}{n-l}x^{n-l}$

and from here follows (43). ∎

We will now show the relationship between the $(k,1)$-autonomous coefficients
and the binomial coefficients

######  Theorem 16.

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1}{0}_{(k,1)}$
$\displaystyle=\binom{n}{k-1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1-k}{0}_{(k,1)}$
$\displaystyle+\sum_{h=k+1}^{n}\binom{n}{h}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{0}_{(k,1)}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{h+1-k}{0}_{(k,1)},$
(44)

for $1\leq i\leq n-k+3$

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1}{i}_{(k,1)}$
$\displaystyle=\sum_{h=0}^{k-2}\binom{n}{h}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{i-1}_{(k,1)}$
$\displaystyle+\binom{n}{k-1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1-k}{i}_{(k,1)}+\binom{n}{k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-k}{i-1}_{(k,1)}$
$\displaystyle+\sum_{h=k+1}^{n}\binom{n}{h}\sum_{j+l=i}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{j}_{(k,1)}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{h+1-k}{l}_{(k,1)}$
(45)

and for $n-k+4\leq i\leq n+1$

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1}{i}_{(k,1)}$
$\displaystyle=\sum_{h=0}^{n-i+1}\binom{n}{h}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{i-1}_{(k,1)}$
$\displaystyle+\binom{n}{k-1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1-k}{i}_{(k,1)}+\binom{n}{k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-k}{i-1}_{(k,1)}$
$\displaystyle+\sum_{h=k+1}^{n}\binom{n}{h}\sum_{j+l=i}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{j}_{(k,1)}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{h+1-k}{l}_{(k,1)}$
(46)

###### Proof.

As

$\displaystyle A_{n+1+k}^{(k)}(x,1)$
$\displaystyle=\sum_{i=0}^{n}\binom{n}{i}A_{n-i+k}^{(k)}(x,1)A_{i+1}^{(k)}(x,1)$
$\displaystyle=\sum_{i=0}^{k-2}\binom{n}{i}A_{n-i+k}^{(k)}(x,1)x+\binom{n}{k-1}A_{n+1}^{(k)}(x,1)$
$\displaystyle+\binom{n}{k}A_{n}^{(k)}(x,1)x+\sum_{i=k+1}^{n}\binom{n}{i}A_{n-i+k}^{(k)}(x,1)A_{i+1}^{(k)}(x,1),$

then

$\displaystyle\sum_{i=0}^{n+1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1}{i}_{(k,1)}x^{i}$
$\displaystyle=\sum_{i=0}^{k-2}\binom{n}{i}\sum_{j=0}^{n-i}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-i}{j}_{(k,1)}x^{j+1}$
$\displaystyle+\binom{n}{k-1}\sum_{j=0}^{n+1-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1-k}{j}_{(k,1)}x^{j}+\binom{n}{k}\sum_{j=0}^{n-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-k}{j}_{(k,1)}x^{j+1}$
$\displaystyle+\sum_{i=k+1}^{n}\binom{n}{i}\left(\sum_{j=0}^{n-i}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-i}{j}_{(k,1)}x^{j}\sum_{l=0}^{i+1-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{i+1-k}{l}_{(k,1)}x^{l}\right).$

We multiply the two autonomous polynomials within the last sum

$\displaystyle\sum_{i=0}^{n+1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1}{i}_{(k,1)}x^{i}$
$\displaystyle=\sum_{i=0}^{k-2}\binom{n}{i}\sum_{j=0}^{n-i}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-i}{j}_{(k,1)}x^{j+1}$
$\displaystyle+\binom{n}{k-1}\sum_{j=0}^{n+1-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1-k}{j}_{(k,1)}x^{j}+\binom{n}{k}\sum_{j=0}^{n-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-k}{j}_{(k,1)}x^{j+1}$
$\displaystyle+\sum_{i=k+1}^{n}\binom{n}{i}\sum_{h=0}^{n+1-k}\left(\sum_{j+l=h}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-i}{j}_{(k,1)}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{i+1-k}{l}_{(k,1)}\right)x^{h}.$

Then by rearranging the first and fourth sums we obtain

$\displaystyle\sum_{i=0}^{n+1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1}{i}_{(k,1)}x^{i}$
$\displaystyle=\sum_{i=0}^{n-k+3}\left(\sum_{h=0}^{k-2}\binom{n}{h}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{i-1}_{(k,1)}\right)x^{i}$
$\displaystyle+\sum_{i=n-k+4}^{n+1}\left(\sum_{h=0}^{n+1-i}\binom{n}{h}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{i-1}_{(k,1)}\right)x^{i}$
$\displaystyle+\binom{n}{k-1}\sum_{i=0}^{n+1-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n+1-k}{i}_{(k,1)}x^{i}$
$\displaystyle+\binom{n}{k}\sum_{i=1}^{n-k+1}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-k}{i-1}_{(k,1)}x^{i}$
$\displaystyle+\sum_{i=0}^{n+1-k}\left(\sum_{h=k+1}^{n}\binom{n}{h}\sum_{j+l=i}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n-h}{j}_{(k,1)}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{h+1-k}{l}_{(k,1)}\right)x^{i}.$

For a suitable value of $i$ the desired results are attained. ∎

Finally, we show without proof the relationship between the $(k,1)$-autonomous
coefficients and the Stirling numbers of second kind.

######  Conjecture 1.

Suppose that $A_{1}^{(k)}(1,1)=\cdots=A_{k}^{(k)}(1,1)=1$. Then

$B_{n}(A_{1}^{(k)}(1,1),...,A_{n}^{(k)}(1,1))=\sum_{i=1}^{n}\genfrac{\\{}{\\}}{0.0pt}{}{n}{i}A_{i}^{(k)}(1,1),\
\ n\geq 1$ (47)

Then,

$A_{n+k}^{(k)}(1,1)=\sum_{i=1}^{n}\genfrac{\\{}{\\}}{0.0pt}{}{n}{i}A_{i}^{(k)}(1,1),\
\ n\geq 1$ (48)

and

$\sum_{i=0}^{n}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{i}_{(k,1)}=\sum_{j=1}^{k}\genfrac{\\{}{\\}}{0.0pt}{}{n}{j}+\sum_{j=k+1}^{n}\genfrac{\\{}{\\}}{0.0pt}{}{n}{j}\sum_{i=0}^{j-k}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{i}_{(k,1)}$
(49)

Equation (48) corresponds to the number of shifts left $k-1$ places under
Stirling transform.

##  5 Sequences related to the equation (17)

We conclude this article by showing sequences related to equation (17) for
values of $k=2,3,4$. Especially, we show that the numbers known as reduced
tangent numbers, Bernoulli numbers, Euler zigzag numbers, Eulerian numbers,
Blasius numbers, triangular numbers, number of shifts left 3 places under
Stirling transform, and number of 8-sequences of $[1,n]$ with 2 contiguous
pairs can be constructed using Bell polynomials, Stirling numbers of second
kind, binomial coefficient and autonomous coefficients.

###  5.1 Case k=2

The first case to be studied is

$y^{\prime\prime}=ae^{y}$ (50)

The equation (50) is equivalent to the equation
$y^{(3)}=y^{\prime\prime}y^{\prime}$, whose solution is

$y^{\prime}=\sqrt{2}\sqrt{c_{1}}\tan\left(\frac{1}{2}\sqrt{2}\sqrt{c_{1}}t+\frac{1}{2}\sqrt{2}\sqrt{c_{1}}c_{2}\right)$

and therefore

$y=\ln\left(\sec^{2}\left(\frac{1}{2}\sqrt{2}\sqrt{c_{1}}t+\frac{1}{2}\sqrt{2}\sqrt{c_{1}}c_{2}\right)\right)+c_{3}$
(51)

where $c_{1},c_{2}$ and $c_{3}$ are constants in $\mathbb{C}$. Since we want
$y(0)=x$, $y^{\prime}(0)=y$ and $y^{\prime}(0)=ae^{x}$, then

$\displaystyle\ln\left(\sec^{2}\left(\frac{1}{2}\sqrt{2}\sqrt{c_{1}}c_{2}\right)\right)+c_{3}$
$\displaystyle=x,$
$\displaystyle\sqrt{2}\sqrt{c_{1}}\tan\left(\frac{1}{2}\sqrt{2}\sqrt{c_{1}}c_{2}\right)$
$\displaystyle=y,$ $\displaystyle
c_{1}\sec^{2}\left(\frac{1}{2}\sqrt{2}\sqrt{c_{1}}c_{2}\right)$
$\displaystyle=ae^{x}.$

As result

$\displaystyle c_{1}$ $\displaystyle=ae^{x}-\frac{y^{2}}{2},$
$\displaystyle\frac{1}{2}\sqrt{2}\sqrt{c_{1}}c_{2}$
$\displaystyle=\arctan\left(\frac{y}{\sqrt{2ae^{x}-y^{2}}}\right),$
$\displaystyle c_{3}$
$\displaystyle=x-\ln\left(1+\frac{y^{2}}{2ae^{x}-y^{2}}\right).$

Thus, the function

$\displaystyle E_{2}(t,(x,y),a)$
$\displaystyle=x+\ln\sec^{2}\left(\frac{\sqrt{2ae^{x}-y^{2}}t}{2}+\arctan\left(\frac{y}{\sqrt{2ae^{x}-y^{2}}}\right)\right)$
$\displaystyle-\ln\left(1+\frac{y^{2}}{2ae^{x}-y^{2}}\right)$ (52)

is the solution of the equation (50) with initial value $y(0)=x$,
$y^{\prime}(0)=y$.

The following is a list of particular solutions of (50) which are obtained
from the equation (5.1)

$\displaystyle E_{2}(t,(x,0),a)$
$\displaystyle=x+\ln\left(\sec^{2}\left(\frac{\sqrt{a}e^{x/2}t}{\sqrt{2}}\right)\right),\
\ a>0,$ (53) $\displaystyle E_{2}(t,(x,0),-a)$
$\displaystyle=x+\ln\left(\operatorname{sech}^{2}\left(\frac{\sqrt{a}e^{x/2}t}{\sqrt{2}}\right)\right),\
\ a>0,$ (54) $\displaystyle E_{2}(t,(0,y),a)$
$\displaystyle=\ln\sec^{2}\left(\frac{\sqrt{2a-y^{2}}t}{2}+\arctan\left(\frac{y}{\sqrt{2a-y^{2}}}\right)\right)$
$\displaystyle-\ln\left(1+\frac{y^{2}}{2a-y^{2}}\right),\ \ a>0,$ (55)
$\displaystyle E_{2}(t,(0,y),-a)$
$\displaystyle=\ln\operatorname{sech}^{2}\left(\frac{\sqrt{2a+y^{2}}t}{2}+\operatorname{arctanh}\left(\frac{y}{\sqrt{2a+y^{2}}}\right)\right)$
$\displaystyle-\ln\left(1-\frac{y^{2}}{2a+y^{2}}\right),\ \ a>0.$ (56)

We now show the relationship between reduced tangent numbers (A002105 in OEIS)
and Bell polynomials and binomial coefficients

######  Theorem 17.

Let

$(T_{n})_{n\geq 1}=(1,1,4,34,496,\ldots)$ (57)

denote the sequence of reduced tangent numbers. Then

  1. 1.

$A_{n}^{(2)}(a)=a^{n}T_{n}$.

  2. 2.

$T_{n}=B_{n}(0,T_{1},\ldots,0,T_{n-1})$, $n\geq 2$.

  3. 3.

$(-1)^{n}T_{n}=B_{n}(0,-T_{1},\ldots,0,(-1)^{n-1}T_{n-1})$, $n\geq 2$.

  4. 4.

$T_{n+2}=\sum_{i=0}^{n}\binom{2n+1}{2i+1}T_{n-i+2}T_{i+1}$, $n\geq 0$.

###### Proof.

Another way to write (53) is

$E_{2}(t,(x,0),a)=x+\sqrt{2}\int_{0}^{\sqrt{a}e^{x/2}t}\tan\left(\frac{u}{\sqrt{2}}\right)du.$

Then

$\displaystyle E_{2}(t,(x,0),a)$
$\displaystyle=x+\int_{0}^{\sqrt{a}e^{x/2}t}\sum_{n=1}^{\infty}T_{n}\frac{u^{2n-1}}{(2n-1)!}$
$\displaystyle=x+\sum_{n=1}^{\infty}T_{n}\frac{(\sqrt{a}e^{x/2}t)^{2n}}{(2n)!}.$

By comparison, $f_{2n}(x,0,a)=a^{n}T_{n}e^{nx}$. Thus follows 1, 2, and 3. The
recurrence relation 4 follows from equation (29). ∎

In general, the solution (53) is the generating function of the sequence

$(a^{n}T_{n})_{n\geq 1}=(a,a^{2},4a^{3},34a^{4},496a^{5},\ldots).$

On the other hand, it is known that $T_{n}=\frac{2^{n}(2^{2n}-1)|b_{2n}|}{n}$,
where the $b_{2n}$ are the Bernoulli numbers (A000367, A002445 in OEIS). Then
the theorem 17 provides a relation between Bell polynomials and Bernoulli
numbers, that is

$\frac{2^{n}(2^{2n}-1)|b_{2n}|}{n}=B_{n}\left(0,6|b_{2}|,0,30|b_{4}|,\ldots,0,\frac{2^{n-1}(2^{2n-2}-1)|b_{2n-2}|}{n-1}\right)$
(58)

We now show the relationship between Euler zigzag numbers (A000111 in OEIS)
and Bell polynomials, binomial coefficients, and Stirling numbers of second
kind

######  Theorem 18.

Suppose

$(A_{n})_{n\geq 0}=(1,1,1,2,5,16,61,272,\ldots)$ (59)

the sequence of Euler zigzag numbers. Then

  1. 1.

$A_{n+1}=B_{n}(A_{0},\ldots,A_{n-1})$, $n\geq 1$.

  2. 2.

$(-1)^{n}A_{n+1}=B_{n}(-A_{0},A_{1},\ldots,(-1)^{n-1}A_{n-1})$, $n\geq 1$.

  3. 3.

$A_{n+2}=\sum_{i=0}^{n}\binom{n}{i}A_{n-i+1}A_{i}$, $\neq 0$.

  4. 4.

$A_{n+2}=\sum_{i=1}^{n}\genfrac{\\{}{\\}}{0.0pt}{}{n}{i}A_{i}$, $n\geq 1$.

  5. 5.

$A_{2n+2}=T_{n}+\sum_{i=1}^{n}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{2n}{2i}_{(2,1)}$.

  6. 6.

$A_{2n+3}=\sum_{i=1}^{n}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{2n+1}{2i+1}_{(2,1)}$

###### Proof.

From equation (55),

$\displaystyle E_{2}(t,(0,1),1)$
$\displaystyle=\ln\left(\sec^{2}\left(\frac{t}{2}+\frac{\pi}{4}\right)\right)-\ln(2)$
$\displaystyle=\ln(\sec^{2}(t)+\sec(t)\tan(t)).$

Then

$\displaystyle E_{2}(t,(0,1),1)$
$\displaystyle=\int_{0}^{t}(\sec(u)+\tan(u))du$
$\displaystyle=\int_{0}^{t}\left(1+\sum_{n=1}^{\infty}A_{n}\frac{u^{n}}{n!}\right)du$
$\displaystyle=t+\sum_{n=1}^{\infty}A_{n}\int_{0}^{t}\frac{u^{n}}{n!}du$
$\displaystyle=t+\sum_{n=1}^{\infty}A_{n}\frac{t^{n+1}}{(n+1)!}$
$\displaystyle=\sum_{n=1}^{\infty}A_{n-1}\frac{t^{n}}{n!}.$

We apply the equation (23) to obtain 1. By corollary 3 it follows that

$E_{2}(t,(0,-1),1)=E(-t,(0,1),1)=\ln(\sec^{2}(t)-\sec(t)\tan(t)).$ (60)

From equation (15) follows 2. Formula 3 follows from equation (27). From
equation (48) follows 4. The identities 5 and 6 follow because the Euler
zigzag numbers are obtained when $x=1$ in $A_{n}^{(2)}(x,1)$. ∎

When $k=2$ the exponential autonomous polynomials and the autonomous
polynomials match. Some autonomous polynomials of the equation (50) are

$\displaystyle q_{1}(y,a)$ $\displaystyle=$ $\displaystyle y,$ $\displaystyle
q_{2}(y,a)$ $\displaystyle=$ $\displaystyle a,$ $\displaystyle q_{3}(y,a)$
$\displaystyle=$ $\displaystyle ay,$ $\displaystyle q_{4}(y,a)$
$\displaystyle=$ $\displaystyle a(a+y^{2}),$ $\displaystyle q_{5}(y,a)$
$\displaystyle=$ $\displaystyle a(4ay+y^{3}),$ $\displaystyle q_{6}(y,a)$
$\displaystyle=$ $\displaystyle a(4a^{2}+11ay^{2}+y^{4}),$ $\displaystyle
q_{7}(y,a)$ $\displaystyle=$ $\displaystyle a(34a^{2}y+26ay^{3}+y^{5}),$
$\displaystyle q_{8}(y,a)$ $\displaystyle=$ $\displaystyle
a(34a^{3}+180a^{2}y^{2}+57ay^{4}+y^{6}).$

From the above we obtain the first $(2,a)$-autonomous coefficients

$n$ $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
0 | $a$ |  |  |  |  |  |
1 | 0 | $a$ |  |  |  |  |
2 | $a^{2}$ | 0 | $a$ |  |  |  |
3 | 0 | $4a^{2}$ | 0 | $a$ |  |  |
4 | $4a^{3}$ | 0 | $11a^{2}$ | 0 | $a$ |  |
5 | 0 | $34a^{3}$ | 0 | $26a^{2}$ | 0 | $a$ |
6 | $34a^{4}$ | 0 | $180a^{3}$ | 0 | $57a^{2}$ | 0 | $a$
Table 1: $(2,a)$-autonomous coefficients

######  Theorem 19.

Some values of $(2,a)$-autonomous coefficients are

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{2n}{2i+1}_{(2,a)}$
$\displaystyle=$
$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{2n+1}{2i}_{(2,a)}=0$
(61)

for all $i$.

###### Proof.

The equation (61) follows from theorem 9. ∎

######  Conjecture 2.

$\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-2}_{(2,a)}=a^{2}(2^{n}-n-1)$
(62)

The sequence

$\displaystyle 2^{n}-n-1$
$\displaystyle=(0,0,1,4,11,26,57,120,247,502,1013,2036,4083,8178,16369,$
$\displaystyle 32752,65519,131054,262125,524268,1048555,2097130,\ldots)$

is known as Eulerian numbers (A000295 in OEIS).

###  5.2 Case k=3

When $k=3$ we obtain the equation

$y^{(3)}=ae^{y}$ (63)

Solving (63) with initial conditions $(0,0,x)$ and $a=-1$ we get the solution
of Blasius equation

$u^{(3)}+u^{\prime\prime}u=0.$ (64)

The Blasius equation [14] describes the velocity profile of the fluid in the
boundary layer which forms when fluid flows along a flat plate. Using the
theorem 9 and the corollary 1 we reach the following result on Blasius numbers
(A018893 in OEIS)

######  Theorem 20.

Let

$(\mathrm{b}_{n})_{n\geq 1}=(1,1,11,375,27.897,\ldots)$ (65)

denote the sequence of Blasius numbers. Then

  1. 1.

$\mathrm{b}_{n}=B_{n}(0,0,\mathrm{b}_{1},\ldots,0,0,\mathrm{b}_{n-1})$, $n\geq
2$.

  2. 2.

$\mathrm{b}_{n+2}=\sum_{i=0}^{n}\binom{3n+2}{3i+2}\mathrm{b}_{n-i+1}\mathrm{b}_{i+1}$,
$n\geq 0$.

On the other hand, the autonomous polynomials for the equation (63) are

$\displaystyle A_{3}^{(3)}(x,a)$ $\displaystyle=a,$ $\displaystyle
A_{4}^{(3)}(x,a)$ $\displaystyle=ax,$ $\displaystyle A_{5}^{(3)}(x,a)$
$\displaystyle=a(x+x^{2}),$ $\displaystyle A_{6}^{(3)}(x,a)$
$\displaystyle=a(a+3x^{2}+x^{3}),$ $\displaystyle A_{7}^{(3)}(x,a)$
$\displaystyle=a(5ax+3x^{2}+6x^{3}+x^{4}),$ $\displaystyle A_{8}^{(3)}(x,a)$
$\displaystyle=a(11ax+16ax^{2}+15x^{3}+10x^{4}+x^{5}),$ $\displaystyle
A_{9}^{(3)}(x,a)$
$\displaystyle=a(11a^{2}+84ax^{2}+(42a+15)x^{3}+45x^{4}+15x^{5}+x^{6}),$
$\displaystyle A_{10}^{(3)}(x,a)$
$\displaystyle=a(117a^{2}x+129ax^{2}+384ax^{3}+(99a+105)x^{4}+105x^{5}+21x^{6}+x^{7})$

and from here we obtain the following table with the first $(3,a)$-autonomous
coefficients

$n$ $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---|---
0 | $a$ |  |  |  |  |  |  |
1 | 0 | $a$ |  |  |  |  |  |
2 | 0 | $a$ | $a$ |  |  |  |  |
3 | $a^{2}$ | 0 | $3a$ | $a$ |  |  |  |
4 | 0 | $5a^{2}$ | $3a$ | $6a$ | $a$ |  |  |
5 | 0 | $11a^{2}$ | $16a^{2}$ | $15a$ | $10a$ | $a$ |  |
6 | $11a^{3}$ | 0 | $84a^{2}$ | $42a^{2}+15a$ | $45a$ | $15a$ | $a$ |
7 | 0 | $117a^{3}$ | $129a^{2}$ | $384a^{2}$ | $99a^{2}+105a$ | $105a$ | $21a$ | $a$
Table 2: $(3,a)$-autonomous coefficients

######  Theorem 21.

Some values of $(3,a)$-autonomous coefficients are

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n}_{(3,a)}$
$\displaystyle=a,$
$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-1}_{(3,a)}$
$\displaystyle=a\binom{n}{2}.$

###### Proof.

The results follow from theorem 15 with $l=0,1$ and by keeping in mind that
$\genfrac{\\{}{\\}}{0.0pt}{}{n}{n-1}=\binom{n}{2}$. ∎

######  Conjecture 3.

$\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-2}_{(3,a)}=a\binom{\binom{n}{2}}{2}$
(66)

The numbers $\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-2}_{(3,1)}$ are
the triangular numbers

$(0,0,3,15,45,105,210,378,630,990,1485,\ldots)$ (67)

(A050534 in OEIS).

Finally, by the equations (15), (39) and (48) we have

######  Theorem 22.

Let

$(\mathrm{e}_{n})_{n\geq 1}=(A_{n}^{(3)}(1,1))_{n\geq
1}=(1,1,1,1,2,5,15,53,213,\ldots).$ (68)

denote the number of shifts 3 places left under exponentiation (A007548 in
OEIS). Then

  1. 1.

$\mathrm{e}_{n+3}=B_{n}(\mathrm{e}_{1},\ldots,\mathrm{e}_{n})$, $n\geq 1$.

  2. 2.

$(-1)^{n}\mathrm{e}_{n+3}=B_{n}(-\mathrm{e}_{1},\mathrm{e}_{2},\ldots,(-1)^{n-1}\mathrm{e}_{n})$,
$n\geq 1$.

  3. 3.

$\mathrm{e}_{n+3}=\sum_{i=1}^{n}\genfrac{\\{}{\\}}{0.0pt}{}{n}{i}\mathrm{e}_{i}$,
$n\geq 1$.

  4. 4.

$\mathrm{d}_{3n}=\mathrm{b}_{n}+\sum_{i=1}^{3n}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{3n}{i}_{(3,1)}$.

  5. 5.

$\mathrm{d}_{3n+j}=\sum_{i=1}^{3n+j}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{3n+j}{i}_{(3,1)}$,
$j=1,2$.

###  5.3 Case k=4

The equation to be studied is

$y^{(4)}=ae^{y}.$ (69)

This equation is not commonly studied in the literature. Here we show the
relation of this equation with the number of shifts left 3 places under
Stirling transform, and also the relation with the numbers $A_{n}^{(4)}(1)$. A
list of exponential autonomous polynomials of the equation (69) is as follows:

$\displaystyle q_{1}(y,z,w,a)$ $\displaystyle=y,$ $\displaystyle
q_{2}(y,z,w,a)$ $\displaystyle=z,$ $\displaystyle q_{3}(y,z,w,a)$
$\displaystyle=w,$ $\displaystyle q_{4}(y,z,w,a)$ $\displaystyle=a,$
$\displaystyle q_{5}(y,z,w,a)$ $\displaystyle=ay,$ $\displaystyle
q_{6}(y,z,w,a)$ $\displaystyle=a(z+y^{2}),$ $\displaystyle q_{7}(y,z,w,a)$
$\displaystyle=a(w+3yz+y^{3}),$ $\displaystyle q_{8}(y,z,w,a)$
$\displaystyle=a(a+3z^{2}+4yw+6y^{2}z+y^{4}).$

From the equation (29) we calculate the first numbers $A_{n}^{(4)}(1)$,

$\displaystyle q_{4}(0,0,0,1)$ $\displaystyle=A_{1}^{(4)}(1)=1,$
$\displaystyle q_{8}(0,0,0,1)$
$\displaystyle=A_{2}^{(4)}(1)=\binom{3}{3}A_{1}^{(4)}A_{1}^{(4)}=1,$
$\displaystyle q_{12}(0,0,0,1)$
$\displaystyle=A_{3}^{(4)}(1)=\binom{7}{3}A_{2}^{(4)}A_{1}^{(4)}+\binom{7}{7}A_{1}^{(4)}A_{2}^{(4)}=35,$
$\displaystyle q_{16}(0,0,0,1)$
$\displaystyle=A_{4}^{(4)}(1)=\binom{11}{3}A_{3}^{(4)}A_{1}^{(4)}+\binom{11}{7}A_{2}^{(4)}A_{2}^{(4)}+\binom{11}{11}A_{1}^{(4)}A_{3}^{(4)}=6140.$

Following theorem 9, corollary 1 and equation (15) we have the following
recurrence relations for the numbers $A_{n}^{(4)}(1)$

######  Theorem 23.

Let $(\mathrm{c}_{n})_{n\geq 1}=(A_{n}^{(4)}(1))_{n\geq
1}=(1,1,35,6140,\ldots)$. Then

  1. 1.

$\mathrm{c}_{n}=B_{n}(0,0,0,\mathrm{c}_{1},\ldots,0,0,0,\mathrm{c}_{n-1})$,
$n\geq 2$.

  2. 2.

$\mathrm{c}_{n+2}=\sum_{i=0}^{n}\binom{4n+3}{4i+3}\mathrm{c}_{n-i+1}\mathrm{c}_{i+1}$,
$n\geq 0$.

  3. 3.

$(-1)^{n}c_{n}=B_{n}(0,0,0,-\mathrm{c}_{1},\ldots,0,0,0,(-1)^{n}\mathrm{c}_{n-1})$

The autonomous polynomials associated with the equation (69) are

$\displaystyle A_{1}^{(4)}(x,a)$
$\displaystyle=A_{2}^{(4)}(x,a)=A_{3}^{(4)}(x,a)=x,$ $\displaystyle
A_{4}^{(4)}(x,a)$ $\displaystyle=a,$ $\displaystyle A_{5}^{(4)}(x,a)$
$\displaystyle=ax,$ $\displaystyle A_{6}^{(4)}(x,a)$
$\displaystyle=a(x+x^{2}),$ $\displaystyle A_{7}^{(4)}(x,a)$
$\displaystyle=a(x+3x^{2}+x^{3}),$ $\displaystyle A_{8}^{(4)}(x,a)$
$\displaystyle=a(a+7x^{2}+6x^{3}+x^{4}),$ $\displaystyle A_{9}^{(4)}(x,a)$
$\displaystyle=a(6ax+10x^{2}+25x^{3}+10x^{4}+x^{5}),$ $\displaystyle
A_{10}^{(4)}(x,a)$
$\displaystyle=a(16ax+32ax^{2}+75x^{3}+65x^{4}+15x^{5}+x^{6}),$ $\displaystyle
A_{11}^{(4)}(x,a)$
$\displaystyle=a(36ax+136ax^{2}+(64a+175)x^{3}+315x^{4}+140x^{5}+21x^{6}+x^{7}).$

We now derive recurrence relations of the numbers $A_{n}^{(4)}(1,1)$ using the
equations (15), (39), and (48).

######  Theorem 24.

Suppose

$(\mathrm{d}_{n})_{n\geq 1}=(A_{n}^{(4)}(1,1))_{n\geq
1}=(1,1,1,1,1,2,5,15,53,222,1115,6698,\ldots)$ (70)

the number of shifts left 3 places under Stirling transform (A336020 in OEIS).
Then

  1. 1.

$\mathrm{d}_{n+4}=B_{n}(\mathrm{d}_{1},\ldots,\mathrm{d}_{n})$, $n\geq 1$.

  2. 2.

$(-1)^{n}\mathrm{d}_{n+4}=B_{n}(-\mathrm{d}_{1},\mathrm{d}_{2}...,(-1)^{n-1}\mathrm{d}_{n})$,
$n\geq 1$.

  3. 3.

$\mathrm{d}_{n+4}=\sum_{i=1}^{n}\genfrac{\\{}{\\}}{0.0pt}{}{n}{i}\mathrm{d}_{i}$,
$n\geq 1$.

  4. 4.

$\mathrm{d}_{4n}=\mathrm{c}_{n}+\sum_{i=1}^{4n}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{4n}{i}_{(4,1)}$.

  5. 5.

$\mathrm{d}_{4n+j}=\sum_{i=1}^{4n+j}\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{4n+j}{i}_{(4,1)}$,
$j=1,2,3$.

The following is a table of the first $(4,a)$-autonomous coefficients

$n$ $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---|---
0 | $a$ |  |  |  |  |  |  |
1 | 0 | $a$ |  |  |  |  |  |
2 | 0 | $a$ | $a$ |  |  |  |  |
3 | 0 | $a$ | $3a$ | $a$ |  |  |  |
4 | $a^{2}$ | 0 | $7a$ | $6a$ | $a$ |  |  |
5 | 0 | $6a^{2}$ | $10a$ | $25a$ | $10a$ | $a$ |  |
6 | 0 | $16a^{2}$ | $32a^{2}$ | $75a$ | $65a$ | $15a$ | $a$ |
7 | 0 | $36a^{2}$ | $136a^{2}$ | $64a+175$ | $315a$ | $140a$ | $21a$ | $a$
Table 3: $(4,a)$-autonomous coefficients

######  Theorem 25.

Some values of $(4,a)$-autonomous coefficients are

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n}_{(4,a)}$
$\displaystyle=a$ (71)
$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-1}_{(4,a)}$
$\displaystyle=a\binom{n}{2}$ (72)
$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-2}_{(4,a)}$
$\displaystyle=a\genfrac{\\{}{\\}}{0.0pt}{}{n+2}{n}$ (73)

###### Proof.

The equations (71)-(73) arise from theorem 15 with $l=0,1,2$. ∎

######  Conjecture 4.

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-3}_{(4,a)}$
$\displaystyle=$ $\displaystyle\frac{5a}{2}(n-1)\binom{n}{5},\ n\geq 5$ (74)

The sequence

$\displaystyle\genfrac{\llbracket}{\rrbracket}{0.0pt}{}{n}{n-3}_{(4,1)}$
$\displaystyle=$
$\displaystyle(10,75,315,980,2520,5670,11550,21780,38610,65065,$
$\displaystyle 105105,163800,247520,364140,523260,736440,1017450,$
$\displaystyle 1382535,1850695,2443980,3187800,4111250,5247450,\ldots)$

counts the number of 8-sequences of $[1,n]$ with 2 contiguous pairs, (A027778
in the OEIS).

## References

  * [1] E. T. Bell, Exponential polynomials, Ann. of Math.  35(2)(1934), 258–277.
  * [2] A. Bernardini and P. E. Ricci, Bell Polinomials and Differential Equations of Freud-Type Polynomials, Math. Comput. Model. 36(2002), 1115-1119.
  * [3] A. Caley, On the theory of analytical forms called trees, Philos. Mag. 13(19)(1857),4-9.
  * [4] L. Comtent, Analyse Combinatoire, Tome Premier, Presses Universitaires de France (1970).
  * [5] C. Faà di Bruno, Sullo sviluppo delle funzioni, Annali di Scienze Matematiche e Fisiche 6(1855), 479–480.
  * [6] L. Foissy, Faá di Bruno subalgebras of the Hopf of planar trees from combinatorial Dyson-Schwinger equations, Adv. Math. 218(1)(2008),136-162.
  * [7] H. W. Gould and J. Quaintance, A Linear Binomial Rrecurrence and the Bell Numbers and Polynomials, Appl. Anal. Discrete Math. 1(2007), 371-385.
  * [8] M. Mihoubi, M. S. Maamra, Touchard polynomials, partial Bell polynomials and polynomials of binomial type, J. Integer Seq., 14(2011) 11.3.1.
  * [9] R. Orozco, Solución Analítica de la Ecuación Diferencial Autónoma de orden k, Revista De Matemática Teoría y Aplicaciones vol(23) fasc.1: 63-84.
  * [10] J. Riordan, Derivatives of composite functions, Bull. Amer. Math. Soc. 52(8)(1946), 664–667.
  * [11] J. Riordan, An Introduction to Combinatorial Analysis, John Wiley & Sons, Inc., New York (1958).
  * [12] G. Yıldız, G. Tınaztepe, and M. Sezer, Bell Polynomial Approach for the Solutions of Fredholm Integro-Differential Equations with Variable Coefficients, CMES-Computer Modeling in Engineering & Sciences, 123(3)(2020), 973–993.
  * [13] Y. Wang and Y. Chen, Bell polynomials approach for two higher-order KdV-type equations in fluids, Nonlinear Anal. Real World Appl., 31(2016), 533-551.
  * [14] H. Blasius, Grenzschichten in Flussigkeiten mit kleiner Reibung, Zeitschr. Math. Phys. 56(1908), 1-37.

2020 Mathematics Subject Classification: Primary 34A34; Secondary 11B37,
11B68, 11B73, 11B83.

_Keywords:_ Bell’s polynomial, autonomous function, autonomous polynomial,
autonomous coefficient, Stirling number.

(Concerned with sequences A000111, A000295, A000367, A002105, A002445,
A007548, A018893, A027778, A050534, and A336020.)


# Oscillations in Mixed-Feedback Systems

Amritam Das<EMAIL_ADDRESS>Thomas Chaffey<EMAIL_ADDRESS>Rodolphe Sepulchre
<EMAIL_ADDRESS>University of Cambridge, Department of Engineering,
Trumpington Street, Cambridge CB2 1PZ

###### Abstract

A new method is presented for the analysis of limit cycle oscillations in
mixed-feedback systems. The calculation of the limit cycle is reformulated as
the zero finding of a mixed-monotone relation, that is, of the difference of
two maximally monotone relations. The problem can then be solved efficiently
by borrowing existing algorithms that minimize the difference of two convex
functions. The potential of the method is illustrated on the classical Van der
Pol oscillator.

###### keywords:

Limit cycle analysis, describing function method, maximal monotonicity, mixed
monotonicity, difference of convex functions.

††thanks: Corresponding author: Amritam Das

, , ,

##  1 Introduction

Mixed-feedback amplification is a ubiquitous mechanism in engineering and
biology to generate and control robust oscillations [1]. The feedback
oscillation results from a balance between (local) positive feedback, a
mechanism that destabilises an otherwise stable equilibrium, and (global)
negative feedback, a mechanism that guarantees boundedness of the solutions.
The physical interpretation is that the oscillation results from an energy
balance between locally active and globally dissipative elements. This energy
balance is conveniently expressed in the language of dissipativity theory [2],
[3].

The mixed-feedback mechanism has long been acknowledged in simple models like
the Van der Pol oscillator, where the desired energy balance is achieved by
inserting a negative resistance element in a passive RLC circuit [4]. The same
mechanism is at the core of the biophysical Hodgkin-Huxley model of neuronal
excitability [5], where the oscillation is controlled by the balance between
two distinct ion channels, one of them providing negative conductance (see
[6]).

Despite being intuitive, the mathematical analysis of oscillations in mixed-
feedback systems still remains challenging beyond the phase portrait analysis
of two-dimensional models. Describing function analysis remains the tool of
choice when the oscillation is nearly harmonic (see [7], Chapter 5). The
limitations of this approximate tool are well understood, and it is, for
instance, known to be of limited use in the analysis of relaxation
oscillations. Successful alternative methods have been developed for relay
feedback systems [8]. Relay feedback systems can be regarded as mixed-feedback
systems, where the hysteresis of the relay models the positive feedback loop.
However, a general theory of oscillations in mixed-feedback systems is still
lacking.

The present paper explores an algorithmic angle of attack, grounded in the
theory of maximal monotone relations. The key proposal is to model a mixed-
feedback system as a mixed feedback loop of monotone relations. The negative
feedback loop preserves monotonicity, whereas the positive feedback loop
locally destroys monotonicity. In recent work, we have explored maximal
monotonicity to algorithmically compute the input-output solutions of monotone
relations [9]. We follow the same approach here, but extend the algorithm from
monotone to mixed-monotone relations. This extension is not new in the field
of optimization, and efficient algorithms have been proposed for the
minimization of difference of convex functions [10], [11]. Such algorithms are
directly applicable to the question of the present paper. We illustrate the
success of that bridge on the classical model of the Van der Pol oscillator.

The rest of the paper is organized as follows. After some preliminary
notations and definitions in Section II, Section III presents the input-output
configuration of mixed-feedback systems. In Section IV, the problem of
computing the output of a mixed-feedback system for a given input is
formulated as a mixed-monotone inclusion problem. Followed by a summary of the
recent result on computing a monotone feedback system’s output in Section V,
Section VI presents the algorithmic solution of the mixed-monotone inclusion
problem by borrowing theories of optimization involving the difference between
two convex functions. In Section VII, the algorithm is illustrated on the Van
der Pol oscillator. Finally, Section VIII presents some discussions on the
developed algorithm and directions of future research.

##  2 Preliminaries

####  2.0.1 Notations

The gradient of a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is denoted
by $\nabla f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. The space of
$T$-periodic, square-summable sequences in $\mathbb{R}^{n}$ is denoted by
$l_{2,T}$. For any $u_{1},u_{2}\in l_{2,T}$, the inner product is denoted by
$\langle u_{1},u_{2}\rangle:=\sum_{n=1}^{T}u_{1}^{\top}(n)u_{2}(n)$ and the
norm of $u\in l_{2,T}$ is denoted by $\mid\mid
u\mid\mid:=\sqrt{\sum_{n=1}^{T}u^{\top}(n)u(n)}$.

####  2.0.2 Convex Function

A function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is convex if for all
$x_{1},x_{2}\in\mathbb{R}^{n}$, $\beta\in(0,1)$, the following inequality
holds:

$f(\beta x_{1}+(1-\beta)x_{2})\leq\beta f(x_{1})+(1-\beta)f(x_{2}).$

Strictness in the above inequality implies strict convexity.

####  2.0.3 Relation

A relation on a space $X$ is a subset $S\subseteq X\times X$.

  * •

We write $y\in S(u)$ to denote $(u,y)\in S$.

  * •

The usual operations on functions can be extended to relations:

$\displaystyle S^{-1}:=$ $\displaystyle\\{(u,y)\mid y\in S(u)\\}$
$\displaystyle S+R:=$ $\displaystyle\\{(u,y+z)\mid(y,u)\in S(u),(z,u)\in R\\}$
$\displaystyle SR:=$ $\displaystyle\\{(x,z)\mid\exists\ y\text{ such that
}(x,y)\in R,(y,z)\in S\\}.$

  * •

The inverse $(\cdot)^{-1}$, in this paper, is always associated with the
relational inverse unless stated otherwise. Note that, unlike the functional
inverse, the relational inverse always exists , and if the functional inverse
does exists, the two inverses coincide. However, given a relation $S$, in
general, $SS^{-1}\neq I$ where $I:=\\{(x,x)\mid x\in X\\}$.

####  2.0.4 Monotone Relation

Let $\mathcal{H}$ be a Hilbert space, equipped with an inner product
$\langle\cdot,\cdot\rangle:\mathcal{H}\times\mathcal{H}\to\mathbb{R}$ and a
norm $\mid\mid\cdot\mid\mid$. A relation $R:\mathcal{H}\to\mathcal{H}$ is
called monotone if, for all $u_{1},u_{2}\in\mathcal{H}$,

$\displaystyle\Big{\langle}u_{1}-u_{2},R(u_{1})-R(u_{2})\Big{\rangle}\geq 0.$
(1)

  * •

Note that this definition refers to monotonicity in the operator theoretic
sense, and this is distinct from the notion of monotonicity in the sense of
partial order preservation by a state-space system (cf. [12]).

  * •

In classical input-output literature, montonicity on $L_{2}$ is known as
incremental positivity. For causal relations, this is equivalent to
incremental passivity. See [13] for more details.

####  2.0.5 Maximal Monotone Relation

A monotone relation $R:\mathcal{H}\to\mathcal{H}$ is called _maximal_ if its
graph $\\{u,y\,|\,y=R(u)\\}$ is not properly contained within the graph of any
other monotone relation.

  * •

Any continuous monotone relation is maximal.

  * •

The (sub)gradient of a convex and proper function is a monotone operator and
if the same function is also closed then its gradient (subgradient) is maximal
monotone [14](pp. 6, 28).

####  2.0.6 Lipschitz Relation

A relation $S:\mathcal{H}\rightarrow\mathcal{H}$ has a Lipschitz constant of
$L>0$ if, for all $(u,w),(v,y)\in S$,

$\displaystyle\mid\mid u-v\mid\mid\leq L\mid\mid w-y\mid\mid$ (2)

  * •

If $L<1$, $S$ is called a contraction map.

####  2.0.7 Coercive and Coccoercive Relation

The relation $R:\mathcal{H}\to\mathcal{H}$ is said to be $\alpha$-coercive if,
for all $u_{1},u_{2}\in\mathcal{H}$, there exists an $\alpha>0$ such that

$\displaystyle\Big{\langle}u_{1}-u_{2},R(u_{1})-R(u_{2})\Big{\rangle}\geq\alpha\mid\mid
u_{1}-u_{2}\mid\mid^{2}.$ (3)

The relation $R$ is called $\beta$-cocoercive for all
$u_{1},u_{2}\in\mathcal{H}$, there exists a $\beta>0$ such that

$\displaystyle\Big{\langle}u_{1}-u_{2},R(u_{1})-R(u_{2})\Big{\rangle}\geq\beta\mid\mid
R(u_{1})-R(u_{2})\mid\mid^{2}.$ (4)

  * •

For causal relations, coercivity is equivalent to input-strict passivity and
cocoercivity is equivalent to output-strict passivity.

  * •

If $A$ is $\alpha$-coercive, then $A^{-1}$ is $\frac{1}{\alpha}$-cocoercive.
If $B$ is $\beta$-cocoercive, then it follows from the Cauchy-Schwarz
inequality that $B$ is $\frac{1}{\beta}$-Lipschitz (for details, see [14],
section 2.2).

##  3 Input-Output Relation of Mixed-Feedback Systems

A mixed-feedback system is defined as a mixed feedback loop around monotone
relations. Its input-output mapping is depicted as a signal-flow graph in
Figure 1. Such a mixed-feedback system maps inputs signals
$u\in\mathcal{H}_{u}$ to output signals $y\in\mathcal{H}_{y}$. At this stage,
the input and output spaces $\mathcal{H}_{u}$ and $\mathcal{H}_{y}$ are
considered to be Hilbert spaces. For example, in order to compute periodic
solutions of mixed feedback systems, the input and output spaces can be chosen
as the spaces of $T$-periodic, square-summable sequences $l_{2,T}$. We
consider $H:\mathcal{H}_{u}\rightarrow\mathcal{H}_{y}$ and
$E_{1},E_{2}:\mathcal{H}_{y}\rightarrow\mathcal{H}_{u}$ to be maximal monotone
relations. Moreover, in this paper, we assume that $H$ is Linear Time
Invariant (LTI). Hence, for a given $H$ (e.g. a transfer function) one may
easily verify the passivity condition and consequently infer its maximal
monotonicity (see [15], [9] for more details).

Figure 1: Input-output configuration of mixed-feeback system

We note that the forward relation can even be chosen as the identity without
any loss of generality. The chosen representation is however preferred as
physical examples of mixed-feedback systems typically define an input-output
relation between the two terminals of a port, such as currents and voltages or
forces and velocities. For instance, see figure 3 for an illustration of the
Van der Pol electrical circuit.

The main focus of this paper is devoted to the question of algorithmically
determining the existence of a limit cycle solution in the mixed feedback
system represented in figure 1. In such case, we consider zero input or
periodic input $u$ and look for a periodic output $y$ as a solution of the
feedback system.

The most classical approach to answer this question would be to perform
describing function analysis (c.f. [16], [17] and the references therein).
Describing function analysis provides a graphical and analytic, nevertheless
_approximate_ , method to predict limit cycles of nonlinear systems. Both the
merits and limitations of describing function analysis are well understood
[7]. The aim of this paper is to complement the describing function analysis
with an algorithmic approach as rooted in the theory of convex optimization.

##  4 Problem Formulation

The key proposal in this paper is to reformulate the existence of a limit
cycle as the zero finding of a specific relation.

To this end, the input and output spaces are chosen as the spaces of
$T$-periodic, square-summable sequences $l_{2,T}$. We rewrite the input-output
relation of mixed-feedback system in figure 1 as

$\displaystyle u\in H^{-1}y+E_{1}(y)-E_{2}(y).$ (5)

Therefore, for a given input $u^{\star}\in l_{2,T}$, the problem of
determining the output from (5) amounts to computing $y=y^{\star}\in l_{2,T}$
such that

$\displaystyle 0\in A(y^{\star})-B(y^{\star}),$ with, $\displaystyle
A(y):=H^{-1}y+E_{1}(y)-u^{\star},\ B(y):=E_{2}(y).$ (6)

Here, both $A$ and $B$ are maximal monotone relations. Since, $H$ is maximal
monotone, $H^{-1}$ is also maximal monotone, as maximal monotonicity is
preserved by the relational inversion. The main observation from (4) is that
the mixed-feedback structure translates into the difference of two maximal
monotone relations. Therefore, finding the output $y$ requires computing the
zero of a _mixed-motonone_ relation, termed as _mixed-monotone zero finding_.

The remainder of this paper is dedicated to make use of efficient algorithms
that are developed in the field of optimization and exploit them to solve the
mixed-monotone zero finding problem.

##  5 Algorithmic Analysis of Monotone Feedback Systems

The property of maximal monotonicity first arose in the study of networks of
nonlinear resistors [18], [19], and has grown to become a fundamental property
in the theory of optimization [20], [21], [14]. Early connections between
maximal monotonicity and passive LTI systems are found in the literature on
nonsmooth dynamical systems [22], [23], and fixed point methods are used to
compute periodic solutions in nonsmooth Lur’e systems [24]. Recently, In [9],
it is shown that the periodic outputs of periodically forced maximal monotone
input-output relations, built from port connections of basic circuit elements,
can be computed using fixed point methods. The main idea of [9] is as follows.

Given a maximal monotone relation $R:l_{2,T}\to l_{2,T}$, we can compute a
solution $x\in l_{2,T}$ to the problem

$0\in R(x),$ (7)

by solving for a fixed point $x_{f}\in l_{2,T}$ of a related relation $F$:

$\displaystyle x_{f}=F(x_{f}).$ (8)

With a suitable choice of $F$, $x_{f}$ is a solution to (7). The fixed point
iteration, or Picard iteration, is given by the update rule

$\displaystyle x^{k+1}=F(x^{k}).$ (9)

If $F$ is a contractive or averaged relation, the Picard iteration converges
to a fixed point of $F$, by the Banach fixed point theorem [25] or the
Krasnosel’skiĭ-Mann theorem [26, Theorem 5.14], respectively. Given an
relation $R$, there is a wide range of choices for the fixed point relation
$F$, depending on the particular properties of $R$ (see [14], [21]).

In (5), assume that the positive feedback is not present (i.e. $E_{2}(y)=0$
for all $y$). Therefore, the input-output relation is now given by a maximal
monotone relation, i.e. $u=H^{-1}y+E_{1}(y)$. Hence, finding the output
$y^{\star}\in l_{2,T}$ given a particular input $u^{\star}\in l_{2,T}$
requires finding the zero of a monotone relation as follows:

$\displaystyle 0\in H^{-1}y^{\star}+E_{1}(y^{\star})-u^{\star}.$ (10)

Note that $u^{\star}$ is an offset, which does not affect monotonicity, and
maximal monotonicity is preserved under relational summation [14]. Hence, (10)
can be solved using a fixed point iteration (see [14], Chapter 2). In [9],
this approach is used to compute the periodic output of a periodically forced
port circuit, where $u^{\star}$ represents an input current or voltage.

##  6 Connection between Mixed-Monotone Zero-Finding Problem and DC Programing

The theory of maximal monotone relations has a close connection with the
theory of convex optimization. In fact, for minimizing a closed, convex, and
proper function $F(x)$, the first-order optimality condition is to find a zero
of its (sub)gradient, $\nabla F(x)$, which is a maximal monotone relation
[14]. In turn, if the objective is the difference of two convex functions, say
$F(x)-G(x)$ where $G(x)$ is an another closed, convex, and proper function,
then the first-order optimality condition is to find a zero of $\nabla
F(x)-\nabla G(x)$ which is the difference of two maximal monotone relations.
This correspondence suggests that methods for optimizing the difference of two
convex functions can be utilised to find a zero of the difference of two
maximal monotone relations, which are not necessarily the (sub)gradients of
closed, convex, and proper functions.

###  6.1 Difference of Convex (DC) Programing

In literature, optimizing the difference of convex functions is known as
Difference of Convex (DC) programing.

######  Definition 1.

(Unconstrained DC Programing) Let $F:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and
$G:\mathbb{R}^{n}\rightarrow\mathbb{R}$ are closed, strictly convex, and
proper functions. A Difference of Convex (DC) programing amounts to finding a
local minimum of $F-G$.

DC programing is a natural way to make use of convex programing in non-convex
optimization. The reader is referred to [10] and the references therein for an
overview of DC programing and its potential applications. Although based on
heuristics, it has been used with success in a number of structured problems
that include minimizing quadratic function with indefinite weight, polynomial
functions etc. In fact, any function whose second partial derivatives are
continuous everywhere can be expressed as the DC functions and, hence, its
minimization can be cast as a DC program. Minimization of DC functions may
also include constraints that are rewritten in terms of inequalities on DC
functions.

###  6.2 Disciplined DC Programing

Particulary relevant to the present paper is the disciplined DC programing
approach recommended in [27]. The reasoning is as follows: if $G$ is an affine
function in definition 1, then DC programing is equivalent to convex
programing and can be solved efficiently. Hence, one way to solve a DC program
and find a local optimum value is to iteratively linearize the concave
function $-G$, and solve a locally convex problem. The corresponding iterative
procedure is known as Convex-Concave Procedure (CCP) [28] and can be
implemented using the following algorithm.

Algorithm 1 : CCP for Unconstrained DC Program

1:Initialize: $i:=0$

2:Given: An initial guess $x_{0}$, $\epsilon_{1}$

3:for $i=i+1$ do

4: Linearize $G(x)$:  $\hat{G}(x;x_{i})=G(x_{i})+\nabla
G(x_{i})^{\top}(x-x_{i})$

5: Locally solve convex optimization problem:
$x_{i+1}=\arg\min\limits_{x}F(x)-\hat{G}(x;x_{i})$

6:end for when $\frac{\max\mid x_{i+1}-x_{i}\mid}{\max\mid
x_{i}\mid}<\epsilon_{1}$ with $x^{\star}=x_{i+1}$

7:return $x^{\star}$

In Algorithm 1, it has been proven that all the iterates are feasible and the
objective value locally converges [29].

###  6.3 Solving Mixed-Monotone Inclusion Problems Using DC Programing

Disciplined DC programing provides a clear path to an iterative algorithm that
computes a zero of the difference of the monotone operator $A(x)-B(x)$. Here,
linearizing the concave term at each iteration in Algorithm 1 simply
translates into evaluating the non-monotone term at each iteration, since
$A(x)-B(a)$ is monotone for each constant $a$. This leads to the following
iteration:

Algorithm 2 : Mixed-Monotone Zero Finding

1:Initialize: $i:=0$

2:Given: An initial guess $x_{0}$, $\epsilon_{1}$

3:for $i=i+1$ do

4: Solve for $x_{i+1}$ given:  $0\in A(x_{i+1})-B(x_{i})$

5:end for when $\frac{\max\mid x_{i+1}-x_{i}\mid}{\max\mid
x_{i}\mid}<\epsilon_{1}$ with $x^{\star}=x_{i+1}$

6:return $x^{\star}$

######  Remark 2.

(Reasoning Behind Algorithm 2)

If $x^{\star}$ is a root of $A(x)-B(x)$, then $A(x)-B(x^{\star})$ clearly has
$x=x^{\star}$ as a root. The point $x^{\star}$ is therefore a fixed point of
the scheme above. The converse statement also holds true.

Algorithm 2 provides an iterative approach for solving the mixed-monotone zero
finding problem by solving a sequence of monotone zero-finding problems. This
means that each iteration can be solved efficiently by using a fixed point
iteration (c.f. [14], Chapter 2).

###  6.4 Local Convergence of Mixed-Monotone Zero Finding Algorithm

In Algorithm 2, the mixed monotone update rule, which solves $0\in
A(x_{i+1})-B(x_{i})$ can be rewritten as

$x_{i+1}\in A^{-1}B(x_{i})$ (11)

Suppose that $A$ is $\alpha$-coercive and $B$ is $\beta$-cocoercive. Then
$A^{-1}B$ is $\frac{\alpha}{\beta}$-Lipschitz. If $\frac{\alpha}{\beta}<1$,
the update rule is a contraction mapping, and the iteration converges to a
unique fixed point. Moreover, this fixed point is a solution to the original
problem.

We can generalize the convergnce result to the case where the properties hold
only locally, allowing the result to be applicable to the problems with
multiple solutions.

To this end, we use the following version of Banach fixed-point theorem.

######  Lemma 3.

Suppose $\mathcal{D}$ is a $T$-invariant, closed subset of a Hilbert space
where $T:\mathcal{D}\rightarrow\mathcal{D}$ is a contraction map. Then $T$
admits a unique fixed-point $x^{*}\in\mathcal{D}$ (i.e. $T(x^{*})=x^{*}$).
Furthermore, starting with an arbitrary element $x_{0}\in\mathcal{D}$ there
exists a sequence $\\{x_{i}\\}$ given by $x_{i+1}=T(x_{i})$ for $i\geq 1$ such
that $x_{i}\rightarrow x^{*}$ when $i\rightarrow\infty$.

The invariance of $\mathcal{D}$ under $T$ ensures that the sequence of
iterates remains in $\mathcal{D}$. Moreover, the closedness of $\mathcal{D}$
and the completeness of a Hilbert space mean that $\mathcal{D}$ is complete,
therefore contains the limit of the Cauchy sequence. Using these properties,
the proof is identical to the standard Banach fixed-point theorem (e.g. see
[30]).

######  Theorem 4.

(Local Convergence of Algorithm 2)

Given a Hilbert space $\mathcal{H}$ let $A:\mathcal{H}\rightarrow\mathcal{H}$,
$B:\mathcal{H}\rightarrow\mathcal{H}$ and $\mathcal{D}\subseteq\mathcal{H}$
such that

  1. 1.

$\mathcal{D}$ is closed,

  2. 2.

$A^{-1}B(\mathcal{D})\subseteq\mathcal{D}$,

  3. 3.

$A$ is $\alpha$-coercive on $\mathcal{D}$,

  4. 4.

$B$ is $\beta$-cocoercive on $\mathcal{D}$,

  5. 5.

$\frac{\alpha}{\beta}<1$.

Then $A^{-1}B$ admits a unique fixed-point $x^{*}\in\mathcal{D}$ (i.e.
$A^{-1}B(x^{*})=x^{*}$). Furthermore, starting with an arbitrary element
$x_{0}\in\mathcal{D}$ there exists a sequence $\\{x_{i}\\}$ given by
$x_{i+1}=A^{-1}B(x_{i})$ for $i\geq 1$ such that $x_{i}\rightarrow x^{*}$ when
$i\rightarrow\infty$.

The proof follows directly from Lemma 3, noting that the conditions of Theorem
4 ensure that $A^{-1}B$ is $\frac{\alpha}{\beta}$-Lipschitz on $\mathcal{D}$.

######  Remark 5.

According to Theorem 4, existence of the closed subset
$\mathcal{D}\subseteq\mathcal{H}$ guarentess local existence and uniqueness of
a fixed point in the iterative scheme $0\in A(x_{i+1})-B(x_{i})$, hence the
existence and uniqueness of a $x^{\star}$ such that $0\in
A(x^{\star})-B(x^{\star})$. However, $A-B$ may have none, one or more than one
zeros, making the selection of $\mathcal{D}$ problem-specific.

###  6.5 Illustration: Double-Well Potential Function

As an example, consider the minimization problem of the double-well potential
function (also considered in [29])

$\displaystyle\arg\min\limits_{x\in\mathbb{R}}\
\frac{x^{4}}{12}-\frac{x^{2}}{2}.$ (12)

According to the first-order optimality condition, its critical point can be
determined by solving the following mixed monotone zero finding problem.

$\displaystyle 0\in\frac{x^{3}}{3}-x.$ (13)

Apart from the trivial solution $x=0$, there are two distinct zeros that are
located on the left and right half plane equidistant from the origin
respectively.

Using Algorithm 2, at every iteration $i>0$ and given $x_{i}$, one has to
determine $x_{i+1}$ by solving $0=\frac{x_{i+1}^{3}}{3}-x_{i}$ that yields
$x_{i+1}=\sqrt[3]{3x_{i}}$.

To make the iterative scheme $x_{i+1}=\sqrt[3]{3x_{i}}$ locally convergent,
one has to satisfy the conditions of Theorem 4 by choosing a closed subset
$\mathcal{D}$ from which the initial guess $x_{0}$ is selcted. Note that,
cubic function is $\alpha$-coercive on $\mathbb{R}\setminus\\{0\\}$ where the
value of $\alpha$ is the slope of its tangent and it varies based on the
domain. On the other hand, linear function is $\beta$-cocoercive on entire
$\mathbb{R}$ where the value of $\beta$ is the slope. Therefore, we select two
individual closed intervals $\mathcal{D}^{+}\subseteq\mathbb{R}^{+}$ and
$\mathcal{D}^{-}\subseteq\mathbb{R}^{-}$. The condition
$\frac{\alpha}{\beta}<1$ is achieved by keeping the lower interval of
$\mathcal{D}^{+}$ and $\mathcal{D}^{-}$ small.

Once intialized with $x_{0}=0$, the algorithm converges to the trivial
solution $x=0$. On the othr hand, based on whether the intial guess $x_{0}$ is
chosen to be $\mathcal{D}^{+}$ or $\mathcal{D}^{-}$, the algorithm always
converges to either the positive or the negative fixed point respectively (in
this case, $x=\pm 1.7321$). Therefore, in case of mutiple solutions, Algorithm
2 still converges, however, the convergence is local and depends on the intial
guess. Grphical illustration of the algorithm is depicted in Figure 2.

Figure 2: Geometrical enactment of Algorithm 2 for finding zeros of
$\frac{x^{3}}{3}-x$ (given in  ) that essentially computes the points of
intersection between $y=x$ and $y=\frac{x^{3}}{3}$ (given in ). From two
distinct intial guess $x_{0}=0.3,-0.3$, evolution of each iterate is depcited
in  and  respectively.

##  7 Finding a Periodic Solution to Mixed-Feedback System

We aim to use the Algorithm 2 to compute periodic solution of a mixed-feedback
system as shown in figure 1 casted as a mixed monotone zero finding problem.
in other words, for a given input $u^{\star}\in l_{2,T}$, we aim to compute
the output $y^{\star}\in l_{2,T}$ such that

$\displaystyle 0\in A(y^{\star})-B(y^{\star}),$ with, $\displaystyle
A(y):=H^{-1}y+E_{1}(y)-u^{\star},\ B(y):=E_{2}(y).$ (14)

According to Alogorithm 2, finding $y^{\star}$ from (7) requires to start from
an intial guess $y_{0}$ that belongs to a preselected set
$\mathcal{D}\subseteq l_{2,T}$ and iteratively compute the zero of the maximal
monotone (locally) relation $A(y)-B(y_{i})$ for every iteration $i>0$ untill
the stopping criterion is satisfied.

###  7.1 Computing the zero of a locally monotone relation

At every $i^{\text{th}}$ iteration, the maximal monotone relation
$R(y;y_{i}):=A(y)-B(y_{i})$ often admits a natural splitting such as
$R=F_{1}+F_{2}$ where $F_{1}$ and $F_{2}$ are maximal monotone operators. In
such a case, the Douglas-Rachford splitting algorithm (c.f. [31], [32]) is
applied, and the computation for the $F_{1}$ and $F_{2}$ components are
separated. Given $y^{j}$, the Douglas-Rachford algorithm is given by the
following iteration

$\displaystyle w^{j+1/2}=$
$\displaystyle\operatorname{res}_{F_{1},\lambda}(y^{j})$ $\displaystyle
z^{j+1/2}=$ $\displaystyle 2w^{j+1/2}-y^{j}$ $\displaystyle w^{j+1}=$
$\displaystyle\operatorname{res}_{F_{2},\lambda}(z^{j+1/2})$ $\displaystyle
y^{j+1}=$ $\displaystyle y^{j}+w^{j+1}-w^{j+1/2},$ (15)

where $\operatorname{res}_{A,\lambda}$ denotes the resolvent of $A$, given by
$(I+\lambda A)^{-1}$, for parameter $\lambda>0$. The estimates $y$ and $w$
both converge to a fixed point of this iteration which is a zero of
$R=F_{1}+F_{2}$, as both $F_{1}$ and $F_{2}$ are maximal monotone.

###  7.2 Complete Algorithm

The complete algorithm employed to determine the solution to the mixed-
monotone inclusion problem (7) is as follows:

Algorithm 3 : Computing Solution to Mixed-Feedback Systems

1:Given: $H,E_{1},E_{2},u^{\star}$

2:Initialize: $i:=0$

3:Given: An initial guess $y_{0}$, $\lambda$, $\epsilon_{1}$, $\epsilon_{2}$

4:for $i=i+1$ do

5: Define $F_{1},F_{2}$ :  $\text{ such that }R=F_{2}+F_{2}\text{ and
}R(y,y_{i}):=A(y)-B(y_{i})$

6: Initialize $j=0$, given an initial guess $y^{0}=y_{i}$

7: for $j=j+1$ do $\displaystyle w^{j+1/2}=$
$\displaystyle\operatorname{res}_{F_{1},\lambda}(y^{j})$ $\displaystyle
z^{j+1/2}=$ $\displaystyle 2w^{j+1/2}-y^{j}$ $\displaystyle w^{j+1}=$
$\displaystyle\operatorname{res}_{F_{2},\lambda}(z^{j+1/2})$ $\displaystyle
y^{j+1}=$ $\displaystyle y^{j}+w^{j+1}-w^{j+1/2}$

8: end for when $\frac{\max\mid y^{j+1}-y^{j}\mid}{\max\mid
y^{j}\mid}<\epsilon_{2}$ with  $y_{\star}=y^{j+1}$

9: return $y_{\star}$

10: Update iteratively:  $y_{i+1}=y_{\star}$

11:end for when $\frac{\max\mid y_{i+1}-y_{i}\mid}{\max\mid
y_{i}\mid}<\epsilon_{1}$ with $y^{\star}=y_{i+1}$

12:return $y^{\star}$

##  8 Illustration: Van Der Pol Oscillator

The differential equation governing the Van der Pol Oscillator is

$\ddot{x}+K(x^{2}-1)\dot{x}+x=0,$ (16)

where, for time $t$, $x(t)\in\mathbb{R}$ is the state-variable and $K\geq 0$
is a constant parameter. The limit cycle solution to (16) and how it changes
as a function of the parameter $K$ is well understood. In particular, when
$K\rightarrow 0$, the limit cycle solution approaches a pure sinusoid with
frequency $1\text{rad}/\text{sec}$. This limiting case is well captured by
harmonic analysis. On the other hand, when $K\rightarrow\infty$, the limit
cycle solution approaches a square-wave solution typical of relaxation
oscillators [33]. It is well known that the prediction of the describing
function analysis is a periodic solution that has a natural frequency of
$1\text{rad}/\text{sec}$ and amplitude of $2$, regardless of the value of $K$
(see [7], Chapter 5). As a result, the describing function analysis fails to
capture the change in the periodic solution for larger values of $K$.

Here, the periodic solutions to the Van der Pol oscillator model is determined
by the Algorithm 3 and involves the following considerations.

####  8.0.1 Reformulation as a Mixed-Feedback System

The decomposition of Van der Pol model as a mixed feedback system is shown in
Figure 3. The positive feedback is naturally identified as the negative
damping term $-K\dot{x}$ in (16).

Figure 3: The mixed-feeback block diagram of the Van der Pol oscillator model.
Here, $s(\cdot):=\frac{\text{d}}{\text{dt}}(\cdot)$

As a result, computing a periodic solutio to the Van der Pol oscillator model
amounts to computing the output $x^{\star}\in l_{2,T}$ such that

$\displaystyle 0\in A(x^{\star})-B(x^{\star}),$ with, $\displaystyle
A(x):=\Bigg{(}\frac{s^{2}+1}{s}\Bigg{)}x+K\frac{x^{3}}{3},\ B(x):=Kx.$ (17)

####  8.0.2 Computation of the Algorithm

According to Algorithm 3, one must iteratively find a zero of
$R(x,x_{i}):=A(x)-B(x_{i})$ starting from an initial guess $x_{0}$ for each
$i>0$ till the stopping criterion is met. Here, $R(x,x_{i})$ admits a natural
splitting such that $R=F_{1}+F_{2}$ where

$\displaystyle F_{1}(x):=$ $\displaystyle\Bigg{(}\frac{s^{2}+1}{s}\Bigg{)}x,$
$\displaystyle F_{2}(x,x_{i}):=$ $\displaystyle K\frac{x^{3}}{3}-Kx_{i}.$ (18)

Hence, to compute a zero of $A(x)-B(x_{i})$ the Douglas-Rachford splitting
algorithm is applied. Starting from an intial guess $x^{0}=x_{i}$, for each
$j>0$ the resolvent of $F_{1}$ and $F_{2}$ are computed at the point $x^{j}$
and (7.1) is computed iteratively. Since $F_{1}(x)$ is linear in $x$, for a
given $x^{j}$, the computation of
$w^{j+1/2}=\operatorname{res}_{F_{1},\lambda}(x^{j})$ amounts to solve
$w^{j+1/2}$ such that

$\displaystyle\Bigg{(}\lambda s^{2}+s+\lambda\Bigg{)}w^{j+1/2}=sx^{j}.$ (19)

To numerically evaluate the differential relations, the central difference
scheme is employed, where for a function $f(t)$,
$[s(f)](t):=\frac{\text{d}}{\text{d}t}f(t)\approx\frac{f(t+h)-f(t-h)}{2h}$ and
$[s^{2}(f)](t):=\frac{\text{d}^{2}}{\text{d}t^{2}}f(t)\approx\frac{f(x+h)-2f(x)+f(x-h)}{h^{2}}$
for a given $h>0$. Moreover, to retain the periodicity, two conditions,
$f(0-h)=f(T)$ and $f(T+h)=f(0)$, are enforced. As a result, determining
$w^{j+1/2}$ in (19) amounts to solving a set of linear equations.

On the other hand, finding the resolvent of
$F_{2}(x):=K\frac{x^{3}}{3}-Kx_{i}$ amounts to numerically finding the root of
a cubic polynomial function with coefficients $(-Kx_{i},0,0,K)$.

####  8.0.3 Simulation Setting

Algorithm 3 is implemented to compute the periodic solution of (16) for three
different values of $K\in\\{0.0002,1.5,10\\}$, capturing three distinct
regimes of oscillatory behaviour. It is expected that with increasing values
of $K$ the periodic solution will change from pure sinusoids to square-waves,
leading to the behavior of a relaxation oscillator.

To obtain a periodic signal from the simulation, for each value of
$K\in\\{0.0002,1.5,10\\}$, the time-interval is chosen to be $[0,nT_{K}]$,
$n=1,2,\cdots$, where $T_{K}$ is the fundamental period of solution for a
given value of $K$. Since the accurate value of the period is unknown for
higher value of $K$, a rough approximation of the period is provided using M.
Cartwright’s semi-empirical formula (c.f. [34]). The time-interval is
partitioned in $N$ equidistant points and the distance between two adjacent
time points is denoted by $h$ where $h=\frac{nT_{K}}{N}$. For all simulations,
we use $n=1$ and $N=5000$. The remaining user-defined paramaters are given in
Table 1. Moreover, the choice of initial guess is motivated by the
illustration of double-well potential function (in subsection 6.5) since the
addition of linear operator preserves the coercivity and cocoercivity.

Table 1: User-defined parameters for simulation Case | $\lambda$ | $\epsilon_{1}=\epsilon_{2}$ | Initial Guess
---|---|---|---
$K=0.0002$ | $0.05$ | $0.01$ | Ramp with slope 1
$K=1.5$ | $0.05$ | $0.01$ | Ramp with slope 1
$K=10$ | $0.01$ | $0.01$ | Ramp with slope 1

####  8.0.4 Simulation Results

The outputs of the algorithm provide the periodic solutions of (16) and are
depicted in figure 4-6 along with the required number of iterations till the
solution converges. As expected, the algorithm demonstrates the effective
change in the oscillatory solution of the Van der Pol osciollator as the value
of $K$ in (16) increases from $0.0002$ till $10$. Evidently, the amplitude and
period of the solution changes with different values of $K$, which is a
significant distinction from the conclusion drawn from the describing function
analysis for the Van der Pol oscillator.

Figure 4: Periodic solution to the Van der Pol Oscillator for $K=0.0002$. The
required number of iteration $i$ is $10$. Figure 5: Periodic solution to the
Van der Pol Oscillator for $K=1.5$. The required number of iteration $i$ is
$8$. Figure 6: Periodic solution to the Van der Pol Oscillator for $K=10$. The
required number of iteration $i$ is $8$.

##  9 Discussion and Future Work

This paper borrows an effective algorithm of DC programing to compute the
limit cycle periodic solution of a general class of mixed-feedback systems.
Solutions of the feedback system are formulated as zeros of an relation. The
mixed nature of the feedback system translates into the mixed monotonicity of
the relation. For a given input, computing the output trajectory is formulated
as finding zero of the difference between two maximal monotone relations and
an iterative algorithm is presented that solves a zero finding problem for a
sequence of monotone relations. The application of the algorithm to the
classical Van der Pol oscillator for different parameter regimes suggests the
efficiency and generality of the proposed approach.

Due to the natural splitting between the feedforward passive operator and the
feedback operators, the presented algorithm is scalable for the passive
operator of arbitrarily large order. For instance, it is directly applicable
to all mixed-feedback models considered in [2], [3] regardless of the
dimension of the LTI passive operator.

As the operators are evaluated on the space $l_{2,T}$ of $T$-periodic, square-
summable signals, the algorithm requires some knowledge about the period.
However, the knowledge of the period does not need to be accurate. The period
is easily adapted over the course of the algorithm based on the zero crossings
of the iterates.

We anticipate that the proposed method is general and offers an attractive
complement to describing function analysis. We also envision natural
extensions to the spatio-temporal mixed feeedback systems such as the PDE
version of the Fitzhug-Nagumo model [35], or standard models of Lateral
Inhibition [36]. Ultimately, we hope to exploit the algorithmic framework of
this paper for the design and analysis of multi-scale mixed-feedback systems
such as those encountered in neuroscience and neuromorphic engineering (see
[1] for more details).

The research leading to these results has received funding from the European
Research Council under the Advanced ERC Grant Agreement Switchlet no. 670645.

## References

  * [1] R. Sepulchre, A Franci, and G. Drion.  Control across scales by positive and negative feedback.  Annual Reviews in Control, Robotics, and Autonomous systems, 2019 (in press).
  * [2] Rodolphe Sepulchre and Guy-Bart Stan.  Feedback mechanisms for global oscillations in lure systems.  Systems & Control Letters, 54(8):809–818, 2005.
  * [3] G. Stan and R. Sepulchre.  Analysis of interconnected oscillators by dissipativity theory.  IEEE Transactions on Automatic Control, 52(2):256–270, 2007.
  * [4] João Carlos Ferreira de Almeida Casaleiro, Luís Augusto Bica Gomes Oliveira, and Igor M. Filanovsky.  Quadrature RC Oscillators: The van der Pol Approach, pages 35–51.  Springer International Publishing, Cham, 2019.
  * [5] A. L. Hodgkin and A. F. Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  The Journal of Physiology, 117(4):500–544, 1952.
  * [6] Rodolphe Sepulchre, Guillaume Drion, and Alessio Franci.  Emerging Applications of Control and Systems Theory: A Festschrift in Honor of Mathukumalli Vidyasagar, pages 269–280.  Springer International Publishing, Cham, 2018.
  * [7] Jean-Jacques E Slotine, Weiping Li, et al.  Applied nonlinear control, volume 199.  Prentice hall Englewood Cliffs, NJ, 1991.
  * [8] Karl J. ÅStröM.  Oscillations in systems with relay feedback.  In K. J. Åström, G. C. Goodwin, and P. R. Kumar, editors, Adaptive Control, Filtering, and Signal Processing, pages 1–25, New York, NY, 1995. Springer New York.
  * [9] Thomas Chaffey and Rodolphe Sepulchre.  Monotone Circuits.  arXiv:2012.11533 [cs, eess, math], December 2020.
  * [10] R. Horst and N. V. Thoai.  Dc programming: Overview.  Journal of Optimization Theory and Applications, 103(1):1–43, 1999\.
  * [11] Thai Quynh Phong.  An algorithm for solving general d.c. programming problems.  Operations Research Letters, 15(2):73–79, 1994.
  * [12] D. Angeli and E.D. Sontag.  Monotone control systems.  IEEE Transactions on Automatic Control, 48(10):1684–1698, 2003\.
  * [13] C. Desoer and M. Vidyasagar.  Feedback Systems, chapter 6, pages 168–227.
  * [14] Ernest K. Ryu and Wotao Yin.  Large-Scale Convex Optimization via Monotone Operators.  Draft edition, 2020.
  * [15] M. K. Camlibel and A. J. van der Schaft.  Incrementally port-hamiltonian systems.  In 52nd IEEE Conference on Decision and Control, pages 2538–2543, 2013.
  * [16] Jay C. Hsu and Andrew U Meyer.  Modern control principles and applications [by] Jay C. Hsu [and] Andrew U. Meyer.  McGraw-Hill New York, 1968.
  * [17] O. R. Fendrich.  Describing functions and limit cycles.  IEEE Transactions on Automatic Control, 37(4):486–487, 1992.
  * [18] George J. Minty.  Monotone networks.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 257(1289):194–212, September 1960.
  * [19] George J. Minty.  On the maximal domain of a “monotone” function.  The Michigan Mathematical Journal, 8(2):135–137, 1961.
  * [20] R. Tyrrell Rockafellar.  Monotone Operators and the Proximal Point Algorithm.  SIAM Journal on Control and Optimization, 14(5):877–898, August 1976.
  * [21] Neal Parikh and Stephen Boyd.  Proximal Algorithms.  Foundations and Trends in Optimization, 1(3):123–231, 2013.
  * [22] M. K. Camlibel and J. M. Schumacher.  Linear passive systems and maximal monotone mappings.  Mathematical Programming, 157(2):397–420, 2016.
  * [23] Bernard Brogliato.  Absolute stability and the Lagrange–Dirichlet theorem with monotone multivalued mappings.  Systems & Control Letters, 51(5):343–353, 2004.
  * [24] W.P.M.H. Heemels, V. Sessa, F. Vasca, and M.K. Camlibel.  Computation of periodic solutions in maximal monotone dynamical systems with guaranteed consistency.  Nonlinear Analysis: Hybrid Systems, 24:100–114, 2017.
  * [25] Stefan Banach.  Sur les opérations dans les ensembles abstraits et leur application aux équations intégrales.  Fundamenta Mathematicae, 3:133–181, 1922.
  * [26] Heinz H. Bauschke and Patrick L. Combettes.  Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  CMS Books in Mathematics. Springer New York, New York, NY, 2011.
  * [27] Xinyue Shen, Steven Diamond, Yuantao Gu, and Stephen Boyd.  Disciplined convex-concave programming, 2016.
  * [28] Alan Yuille and Anand Rangarajan.  The concave-convex procedure.  Neural Computation, 15:915–936, 04 2003.
  * [29] T. Lipp and Stephen P. Boyd.  Variations and extension of the convex–concave procedure.  Optimization and Engineering, 17:263–287, 2016.
  * [30] P. Agarwal, M. Jleli, and B. Samet.  Fixed Point Theory in Metric Spaces: Recent Advances and Applications, pages 1–23.  Springer Singapore, Singapore, 2018.
  * [31] Jim Douglas and H. H. Rachford.  On the Numerical Solution of Heat Conduction Problems in Two and Three Space Variables.  Transactions of the American Mathematical Society, 82(2):421–439, 1956.
  * [32] P. L. Lions and B. Mercier.  Splitting Algorithms for the Sum of Two Nonlinear Operators.  SIAM Journal on Numerical Analysis, 16(6):964–979, 1979.
  * [33] Balth. van der Pol Jun. D.Sc.  Lxxxviii. on “relaxation-oscillations”.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):978–992, 1926.
  * [34] Jean-Marc Ginoux and Christophe Letellier.  Van der pol and the history of relaxation oscillations: Toward the emergence of a concept.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(2):023120, 2012.
  * [35] Félix A. Miranda-Villatoro and Rodolphe Sepulchre.  Differential dissipativity analysis of reaction–diffusion systems.  Systems & Control Letters, 148:104858, 2021.
  * [36] S. Amari.  Dynamics of pattern formation in lateral-inhibition type neural fields.  Biological Cybernetics, 27:77–87, 1977.