Center of Marine Numerical Experiment, State Key Laboratory of Ocean
Engineering, Shanghai 200240, China
School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong
University, Shanghai 200240, China

# Avoiding small denominator problems
by means of the homotopy analysis method

Shijun Liao<EMAIL_ADDRESS>(Shijun Liao)

###### Abstract

The so-called “small denominator problem” was a fundamental problem of
dynamics, as pointed out by Poincaré. Small denominators appear most commonly
in perturbative theory. The Duffing equation is the simplest example of a non-
integrable system exhibiting all problems due to small denominators. In this
paper, using the forced Duffing equation as an example, we illustrate that the
famous “small denominator problems” never appear if a non-perturbative
approach based on the homotopy analysis method (HAM), namely “the method of
directly defining inverse mapping” (MDDiM), is used. The HAM-based MDDiM
provides us great freedom to directly define the inverse operator of an
undetermined linear operator so that all small denominators can be completely
avoided and besides the convergent series of multiple limit-cycles of the
forced Duffing equation with high nonlinearity are successfully obtained. So,
from the viewpoint of the HAM, the famous “small denominator problems” are
only artifacts of perturbation methods. Therefore, completely abandoning
perturbation methods but using the HAM-based MDDiM, one would be never
troubled by “small denominators”. The HAM-based MDDiM has general meanings in
mathematics and thus can be used to attack many open problems related to the
so-called “small denominators”.

###### Key Words.:

Small denominator problem, Duffing equation, limit cycle, homotopy analysis
method (HAM), MDDiM.

41A58, 34C25

##  1 Origin of “small denominator problem”

Poincaré Poincare1957  pointed out that the so-called “small denominator
problem” was “the fundamental problem of dynamics”. The small denominator was
first mentioned by Delaunay Delaunay1867  in his 900 pages book about
celestial motions using perturbation method. Poincaré Poincare1957  first
recognized that, when small denominator appears, the coefficients of
perturbation series may grow too large too often, threatening the convergence
of the series. As pointed out by Pérez Rodrigo2011AMS , “small denominators
are found most commonly in the perturbative theory”. It often appears when
perturbation methods are used to solve problems in classical and celestial
mechanics Arnold1963 , fluid mechanics Iooss2005ARMA ; Iooss2015 , and so on
Marmi1988JPA ; Chierchia1996CMP .

What is the origin of the so-called “small denominator problem”? As pointed
out by Giorgilli Giorgilli1998 , the Duffing equation Alex2006ND  “is perhaps
the simplest example of a non-integrable system exhibiting all problems due to
the small denominators”. So, without loss of generality, let us focus on the
forced Duffing equation

${\cal
N}[u(t)]=u^{\prime\prime}(t)+2\;\xi\;u^{\prime}(t)+u(t)+\beta\;u^{3}(t)-\alpha\cos(\Omega
t)=0,$ (1)

where $\cal N$ is a nonlinear operator, the prime denotes the differentiation
with respect to the time $t$, $\alpha$ and $\Omega$ is the amplitude and
frequency of the external force $F=\alpha\cos(\Omega t)$, $\xi>0$ is the
resistance coefficient, and $\beta>0$ is a physical parameter related to
nonlinearity, respectively.

As pointed out by Kartashova Kartashova2011 , “physical classification of PDEs
is based not on the form of equations, but on the form of solutions”. So, let
us consider here the stationary periodic limit-cycle of $u(t)$ as
$t\to+\infty$ of the forced Duffing equation (1), which can be expressed in
the form:

$\displaystyle u(t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{+\infty}\Big{\\{}a_{n}\cos(\omega_{n}t)+b_{n}\sin(\omega_{n}t)\Big{\\}},$
(2)

where $a_{n}$, $b_{n}$ are constants and

$\omega_{n}=(2n-1)\Omega,\hskip 14.22636ptn\geq 1.$ (3)

This is mainly because the common solution

$A\exp(-\xi t)\cos(t)+B\exp(-\xi t)\sin(t)$

of the linear equation

$u^{\prime\prime}(t)+2\xi u^{\prime}(t)+u(t)=0$

tends to zero as $t\to+\infty$ for arbitrary constants $A$ and $B$, and thus
disappear in the so-called “solution-expression” (2) of the limit-cycle.

Let us first show how perturbation technique Kevorkian1996 ; Rand1987  can
bring the so-called small denominators into the above-mentioned problem. Let
$\beta$ is a small parameter and assume that $u(t)$ can be expanded in such a
series

$u(t)=u_{0}(t)+\sum_{n=1}^{+\infty}u_{n}(t)\;\beta^{n}.$ (4)

Substituting it into (1) and equating the like power of $\beta$, we have the
perturbation equations at different orders of $\beta$:

$\displaystyle\beta^{0}$ $\displaystyle:$ $\displaystyle
u_{0}^{\prime\prime}(t)+2\xi u_{0}^{\prime}(t)+u_{0}(t)=\alpha\cos(\Omega t),$
(5) $\displaystyle\beta^{1}$ $\displaystyle:$ $\displaystyle
u_{1}^{\prime\prime}(t)+2\xi u_{1}^{\prime}(t)+u_{1}(t)=-u^{3}_{0}(t),$ (6)
$\displaystyle\beta^{2}$ $\displaystyle:$ $\displaystyle
u_{2}^{\prime\prime}(t)+2\xi
u_{2}^{\prime}(t)+u_{2}(t)=-3u^{2}_{0}(t)u_{1}(t),$ (7)
$\displaystyle\beta^{3}$ $\displaystyle:$ $\displaystyle
u_{3}^{\prime\prime}(t)+2\xi
u_{3}^{\prime}(t)+u_{3}(t)=-3u^{2}_{0}(t)u_{2}(t)-3u_{0}(t)u^{2}_{1}(t),$
$\displaystyle\hskip 85.35826pt\vdots$

The above perturbation equations have the unique linear operator

${\cal L}_{0}[u(t)]=u^{\prime\prime}(t)+2\xi u^{\prime}(t)+u(t),$ (9)

whose inverse operator ${\cal L}_{0}^{-1}$ reads

$\displaystyle{\cal L}_{0}^{-1}\left[\cos(\omega t)\right]$ $\displaystyle=$
$\displaystyle\frac{(1-\omega^{2})\cos(\omega t)+2\;\xi\;\omega\sin(\omega
t)}{(1-\omega^{2})^{2}+4\xi^{2}\omega^{2}},$ (10) $\displaystyle{\cal
L}_{0}^{-1}\left[\sin(\omega t)\right]$ $\displaystyle=$
$\displaystyle\frac{(1-\omega^{2})\sin(\omega t)-2\;\xi\;\omega\cos(\omega
t)}{(1-\omega^{2})^{2}+4\xi^{2}\omega^{2}},$ (11)

where $\omega$ is a frequency. Note that the denominator

$(1-\omega^{2})^{2}+4\xi^{2}\omega^{2}$

becomes rather small when $\omega\to 1$ and $\xi\to 0$ so that small
denominators appear. This is the origin of the so-called “small denominator
problem”. Thus, for Eq. (1), the “small denominator problem” occurs when $\xi$
is small and $\Omega=1/(2n-1)$, i.e. $\omega_{n}=(2n-1)\Omega=1$, for any a
positive integer $n\geq 1$. Without loss of generality, let us focus on here
the fixed values $\alpha=1$ and $\Omega=1/3$, say, $\omega_{2}=3\Omega=1$, but
investigate different values of $\xi$ and $\beta$. In this case, the so-called
“resonance” occurs when $\xi=0$ and the “near resonance” occurs for a small
$\xi$, corresponding to the “small denominator problem”.

In the frame of the perturbation approach, the unique initial guess is given
according to (5) and (10), say,

$u_{0}(t)=\frac{(1-\Omega^{2})\cos(\Omega t)+2\;\xi\;\Omega\sin(\Omega
t)}{(1-\Omega^{2})^{2}+4\xi^{2}\Omega^{2}}.$ (12)

Note that small denominator appears when $\Omega\to 1$ and $\xi\to 0$ for this
perturbation initial approximation! This is the reason why we choose
$\Omega=1/3=\omega_{1}$ in this paper, otherwise the perturbation method fails
at the very beginning.

Let us first consider the perturbation method in case of $\Omega=1/3$ and
$\alpha=1$. It is found that, when $\xi=1/100$, say, there exists the so-
called “small denominator problem” for the terms $\cos(\omega_{2}t)$ and
$\sin(\omega_{2}t)$, where $\omega_{2}=3\Omega=1$, the perturbation series is
divergent even for a small value $\beta=0.012$, corresponding to a rather weak
nonlinearity. Thus, when the so-called “small denominator problem” appears,
the perturbation method is indeed invalid in practice. Besides, it is found
that, when $\xi=1$, the perturbation series is also divergent for $\beta\geq
0.55$, indicating that the perturbation series is divergent for high
nonlinearity even if the “small denominator problem” does not occur!
Therefore, the perturbation approach indeed does not work for high
nonlinearity and/or when the so-called “small denominator problem” occurs. As
mentioned by Arnol’d Arnold1963 , there often exist two difficulties with
perturbation method in many classical and celestial problems: (1) the
appearance of small denominator , and (2) the divergence of solution series.

The above-mentioned perturbation approach has the following disadvantages:

  1. (1)

small/large physical parameters should exist;

  2. (2)

there is no freedom to choose its linear operator;

  3. (3)

there is no freedom to choose its initial approximation,

  4. (4)

it is convergent only for weak nonlinearity.

These limitations of perturbation methods are well-known. They are the origin
of the so-called “small denominator problem” and the divergence of solution
series.

Note that, due to some historic reasons, the so-called “small denominator
problem” has very close relationships with perturbation methods. Indeed,
“small denominators are found most commonly in the perturbative theory”, as
pointed out by Pérez Rodrigo2011AMS . Is perturbation method the only way to
solve these problems? What happens if we completely abandon perturbation
methods ?

In this paper, we use the forced Duffing equation (1) to illustrate that the
so-called “small denominator problem” can never appear if we completely
abandon perturbation techniques but use a non-perturbative technique, namely
the homotopy analysis method (HAM) Liao1992PhD ; Liao2003Book ; Liao2012Book ;
Liao1995IJNLM ; Liao1997IJNLM ; Liao1999IJNLM ; Liao2004AMC ; Liao2007SAM ;
Liao2009CNSNS ; Liao2010CNSNS ; Liao2016JFM ; Liao2020SCPMA ; Zhong2018JFM ,
which is based on the basic concept “homotopy” in topology and can overcome
all restrictions of perturbation methods. A new HAM-based approach is
proposed, which provides us great freedom to directly define the inverse
operator of an auxiliary linear operator so that all small denominators can be
completely avoided. Convergent series of multiple limit-cycles of the Duffing
equation are successfully obtained, although the directly defined inverse
operators might be beyond the traditional mathematical theories. Thus, from
the viewpoint of the HAM, the famous “small denominator problems” are only
artifacts of perturbation methods. Therefore, completely abandoning
perturbation methods, one would be never troubled by small denominators, as
illustrated below in this paper.

##  2 Basic ideas of the HAM

Can we avoid the famous “small denominator problem” in a systematic way? The
answer is yes. Here, we give an approach based on the homotopy analysis method
(HAM) Liao1992PhD ; Liao2003Book ; Liao2012Book , which can completely avoid
the “small denominator problem”.

“Small denominators are found most commonly in the perturbative theory”, as
pointed out by Pérez Rodrigo2011AMS . So, in order to avoid “small denominator
problem”, we must abandon perturbation methods completely. The homotopy
analysis method (HAM) was proposed by Liao in 1992 in his dissertation
Liao1992PhD . Based on the basic concept homotopy in topology Hilton1953 ,
i.e. a continuous deformation, the HAM Liao2003Book ; Liao2012Book ;
Liao1995IJNLM ; Liao1997IJNLM ; Liao1999IJNLM ; Liao2004AMC ; Liao2007SAM ;
Liao2009CNSNS ; Liao2010CNSNS ; Liao2016JFM ; Liao2020SCPMA ; Zhong2018JFM
has the following advantages:

  1. (a)

unlike perturbation techniques, the HAM works even if there exist no
small/large physical parameters;

  2. (b)

the HAM provides great freedom to choose an auxiliary linear operator;

  3. (c)

the HAM provides great freedom to choose an initial guess;

  4. (d)

different from other approximation methods, the HAM can guarantee the
convergence of solution series even for highly nonlinear problem.

The HAM has been broadly used and its above-mentioned advantages have been
verified and confirmed in thousands of articles by scientists and engineers
all over the world Zhu2006QF ; Bouremel2007CNSNS ; Nassar2011 ;
Kimiaeifar2011CMA ; Ghotbi2011 ; Duarte2015CSF ; Sardanyes2015 ; VanGorder2017
; Pfeffer2017 ; Cullen2019JCP ; Sultana2019EPJP ; Massa2020 ; Kaur2022JMAA ;
Masjedi2022AMM ; Botton2022AMM . In this paper, we use the forced Duffing
equation (1) as an example to illustrate how to completely avoid the “small
denominator problem” by means of a HAM-based approach.

First, let us briefly describe the basic ideas of the HAM using (1) as an
example. Let

$\displaystyle{\cal
S}=\sum_{n=1}^{+\infty}\Big{[}A_{n}\cos(\omega_{n}t)+B_{n}\sin(\omega_{n}t)\Big{]}$
(13)

denote a vector space, where $\omega_{n}$ is defined by (3) and $A_{n},B_{n}$
are arbitrary constants. Let $u(t)\in{\cal S}$, ${\cal L}$ denote an auxiliary
linear operator with the property ${\cal L}[0]=0$, which we have great freedom
to choose, $q\in[0,1]$ be a homotopy parameter, $u_{0}(t)\in{\cal S}$ be an
initial guess of $u(t)$, $c_{0}$ be a constant having no physical meanings,
namely “the convergence-control parameter”, whose value will be determined
later, respectively. Then, we construct a continuous deformation
$\phi(t,q)\in{\cal S}$ from the initial guess $u_{0}\in{\cal S}$ to the
solution $u(t)\in{\cal S}$ of the forced Duffing equation (1), governed by the
so-called zeroth-order deformation equation

$(1-q){\cal L}\left[\phi(t,q)-u_{0}(t)\right]=c_{0}\;q\;{\cal
N}[\phi(t,q)],\hskip 28.45274ptq\in[0,1],$ (14)

where the nonlinear operator ${\cal N}[u]$ is defined by Eq. (1). When $q=0$,
due to the property ${\cal L}[0]=0$ of the auxiliary linear operator, we have
the solution

$\phi(t,0)=u_{0}(t)$ (15)

of Eq. (14). When $q=1$, Eq. (14) is exactly the same as the original equation
(1), thus we have the solution

$\phi(t,1)=u(t),$ (16)

where $u(t)$ is the solution (limiting cycle) of the original Duffing equation
(1). So, as $q$ increases from 0 to 1, $\phi(t,q)$ deforms continuously from
the initial guess $u_{0}(t)$ to the solution $u(t)$ of the original nonlinear
equation (1), since both of $u_{0}(t)\in{\cal S}$ and $u(t)\in{\cal S}$ can be
expressed by the so-called “solution expression” (2). Then, expanding
$\phi(t,q)$ in a power series of $q$, we have according to (15) the homotopy-
series

$\phi(t,q)=u_{0}(t)+\sum_{k=1}^{+\infty}u_{k}(t)\;q^{k}.$ (17)

Note that $\phi(t,q)\in{\cal S}$ is dependent upon the so-called convergence-
control parameter $c_{0}$, which has no physical meanings. Therefore,
$u_{k}(t)\in{\cal S}$ $(k\geq 1)$ in (17) is also dependent upon $c_{0}$ so
that the convergence radius of the series (17) is determined by $c_{0}$. Here,
it should be emphasized that we have great freedom to choose the auxiliary
linear operator $\cal L$, the initial guess $u_{0}$ and the convergence-
control parameter $c_{0}$. This is the key point of the HAM. Assuming that the
auxiliary linear operator $\cal L$, the initial guess $u_{0}$ and the
convergence-control parameter $c_{0}$ are so properly chosen that the
Maclaurin series (17) is convergent at $q=1$, we have due to (16) the
homotopy-series solution

$u(t)=u_{0}(t)+\sum_{k=1}^{+\infty}u_{k}(t).$ (18)

So, even for given auxiliary linear operator $\cal L$ and initial guess
$u_{0}$, the convergence-control parameter $c_{0}$ provides us an additional
way to guarantee the convergence of the solution series, which overcomes the
limitations of perturbation methods mentioned above, as illustrated below in
this paper and other publications Zhu2006QF ; Bouremel2007CNSNS ; Nassar2011 ;
Kimiaeifar2011CMA ; Ghotbi2011 ; Duarte2015CSF ; Sardanyes2015 ; VanGorder2017
; Pfeffer2017 ; Cullen2019JCP ; Sultana2019EPJP ; Massa2020 ; Kaur2022JMAA ;
Masjedi2022AMM ; Botton2022AMM .

Substituting the power series (17) into the zeroth-order deformation equation
(14) and equating the like-power of $q$, we have the high-order deformation
equation

${\cal L}\left[u_{m}(t)-\chi_{m}\;u_{m-1}(t)\right]=c_{0}\;R_{m-1}(t),\hskip
28.45274ptm\geq 1,$ (19)

where

$\displaystyle R_{k}(t)$ $\displaystyle=$
$\displaystyle\left.\frac{1}{k!}\frac{d^{k}{\cal
N}[\phi(t,q)]}{dq^{k}}\right|_{q=0}$ (20)

in general and

$\displaystyle R_{k}(t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l}u^{\prime\prime}_{0}+2\xi
u^{\prime}_{0}+u_{0}+\beta u_{0}^{3}-\alpha\cos(\Omega t),\hskip
28.45274pt\mbox{when $k=0$},\\\ \\\ u^{\prime\prime}_{k}+2\xi
u^{\prime}_{k}+u_{k}+\beta\sum\limits_{i=0}^{k}\sum\limits_{j=0}^{i}u_{k-i}u_{i-j}u_{j},\hskip
21.33955pt\mbox{when $k\geq 1$},\end{array}\right.$ (24)

for the forced Duffing equation (1) considered here, and besides

$\chi_{m}=\left\\{\begin{array}[]{l}0,\hskip 28.45274pt\mbox{when $m=1$},\\\
1,\hskip 28.45274pt\mbox{when $m>1$}.\end{array}\right.$ (25)

The general solution of the linear $m$th-order deformation equation (19) reads

$u_{m}(t)=\chi_{m}\;u_{m-1}(t)+c_{0}\;{\cal
L}^{-1}\Big{[}R_{m-1}(t)\Big{]}+\sum_{n=1}^{\mu}A_{m,n}\psi_{n}(t),\hskip
28.45274ptm\geq 1,$ (26)

where $\mu$ is a positive integer, ${\cal L}^{-1}$ is the inverse operator of
$\cal L$, $A_{m,n}$ is an arbitary constant, and $\psi_{n}(t)\in{\cal S}$ is a
base function satisfying

${\cal L}\Big{[}\sum_{n=1}^{\mu}A_{m,n}\psi_{n}(t)\Big{]}=0,\hskip
14.22636pt1\leq n\leq\mu.$ (27)

In other words, we have

$\ker[{\cal L}]=\sum_{n=1}^{\mu}A_{m,n}\psi_{n}(t),$ (28)

say, the kernel of the auxiliary linear operator $\cal L$ is a vector space in
dimension $\mu$. Note that the linear part (9) of the original Duffing
equation (1) is a second-order differential equation, whose kernel is a vector
space in two dimension. However, we have great freedom to choose the auxiliary
linear operator $\cal L$ and its kernel $\ker[{\cal L}]$, as mentioned below,
which might be a breakthrough in nonlinear differential equations.

The $M$th-order HAM approximation is given by

$u^{*}\approx u_{0}(t)+\sum_{k=1}^{M}u_{k}(t).$ (29)

Since the HAM provides us great freedom to choose the initial guess $u_{0}$,
we can further use the above $M$th-order approximation as a new initial guess,
say, $u_{0}=u^{*}$, to gain another $M$th-order approximation, and so on. This
provides us the $M$th-order iteration approach of the HAM. Note that, for the
HAM iteration approach, in order to avoid the exponential increment in growth
of the terms in the solution expression (2), we eliminate the terms of
$\cos(\omega_{n}t)$, $\sin(\omega_{n}t)$ whose coefficients are less than a
small value, such as $10^{-20}$ for the forced Duffing equation (1) considered
in this paper.

It should be emphasized once again that, different from all other
approximation methods (including perturbation techniques), the homotopy
analysis method (HAM) can guarantee the convergence of solution series by
means of choosing a proper value of the so-called “convergence-control
parameter” $c_{0}$. This is the fundamental difference of the HAM from all
other approaches! The optimal value of the “convergence-control parameter”
$c_{0}$ is determined by the minimum of the residual error square

${\cal E}=\int_{0}^{T}\Big{(}{\cal
N}[u(t)]\Big{)}^{2}dt\approx\frac{1}{K+1}\sum_{j=0}^{K}\Big{(}{\cal
N}[u(j\Delta t)]\Big{)}^{2},$ (30)

where $\Delta t=T/(K+1)$ is a time-step for numerical simulation, $T$ is the
period of the limiting cycle for the considered problem, $K>0$ is a large
enough integer, $u(t)$ is an approximation of limiting cycle of the original
equation (1), ${\cal N}$ is the nonlinear operator defined by (1),
respectively.

What is the relationship between perturbation method and the HAM? Generally
speaking, perturbation approach is often a special case of the HAM, if we
choose the perturbation initial approximation as $u_{0}$, the original linear
operator as the auxiliary linear operator, i.e. ${\cal L}={\cal L}_{0}$, and
besides $c_{0}=-1$. For example, when we choose (9) as the auxiliary linear
operator $\cal L$, (12) as the initial approximation $u_{0}$, and besides set
$c_{0}=-1$, the $k$th-order deformation equation (19) is exactly the same as
the $k$th-order perturbation equation mentioned in § 1. Therefore, the
perturbation approach can be indeed regarded as a special case of the homotopy
analysis method! However, the perturbation approach corresponds to only one
choice, but there exist many other much better choices in the frame of the
HAM, which can avoid the “small denominator problems” completely, as
illustrated below.

In summary, the above-mentioned HAM has the following characteristics:

  1. (A)

the homotopy-series (17) is expanded in the homotopy parameter $q\in[0,1]$
that has no physical meanings at all. So, the HAM has nothing to do with any
small/large physical parameters: it works no matter whether small/large
physical parameters exist or not;

  2. (B)

the HAM provides us great freedom to choose its auxiliary linear operator;

  3. (C)

the HAM provides us great freedom to choose its initial approximation;

  4. (D)

the so-called convergence-control parameter $c_{0}$ has no physical meanings
but can guarantee the convergence of the solution series even for high
nonlinearity, as illustrated below and verified in many related publications
Zhu2006QF ; Bouremel2007CNSNS ; Nassar2011 ; Kimiaeifar2011CMA ; Ghotbi2011 ;
Duarte2015CSF ; Sardanyes2015 ; VanGorder2017 ; Pfeffer2017 ; Cullen2019JCP ;
Sultana2019EPJP ; Massa2020 ; Kaur2022JMAA ; Masjedi2022AMM ; Botton2022AMM .

Thus, the HAM can indeed overcome all limitations and restrictions of
perturbation methods.

##  3 How to avoid “small denominator problem”

As mentioned above, the HAM provides us great freedom to choose the auxiliary
linear operator $\cal L$ and the initial guess $u_{0}$: it is such kind of
freedom that provides us possibility to avoid the so-called “small denominator
problem”, as described below.

###  3.1 Choice of the auxiliary linear operator

Obviously, for the forced Duffing equation (1), the origin of the “small
denominator problem” is mainly due to the original linear operator (9), which
is unique from the viewpoint of perturbation theory. So, in order to avoid
“small denominator problem”, we must abandon (9) thoroughly. Different from
other approximation techniques, the HAM provides us great freedom to choose an
auxiliary linear operator $\cal L$, as illustrated by Liao and Tan Liao2007SAM
and Liao and Zhao Liao2016NA . In most applications of the HAM, one often
chooses a proper auxiliary linear operator $\cal L$ to gain the solution of
the high-order deformation equation (19). However, the freedom is so large
that we can here directly define its inverse operator ${\cal L}^{-1}$ and the
kernel of ${\cal L}$ in (26). In fact, based on the HAM, Liao and Zhao
Liao2016NA  proposed the so-called “method of directly defining inverse
mapping”, i.e. the MDDiM, which has been successfully applied to solve many
types of nonlinear equations KV2018NA ; KV2018AMC ; KV2019CNA ;
thesis2019Dewasurendra ; KV2021AMNS ; KV2021CAM ; KV2022CMA ; Nave2018JBD ;
Nave2018JMM ; Nave2019JMM ; Nave2020BS ; Munjam2022IJAE ; Munjam2023IJAE .
According to the solution expression (2) and the definition (24) of
$R_{m}(t)$, the right-hand side of the high-order deformation equation (19)
contains terms $\cos(\omega_{n}t)$ and $\sin(\omega_{n}t)$, where
$\omega_{n}=(2n-1)\Omega$. Thus, we directly define here its inverse operator

${\cal
L}^{-1}\Big{[}A\cos(\omega_{n}t)+B\sin(\omega_{n}t)\Big{]}=\frac{A\cos(\omega_{n}t)+B\sin(\omega_{n}t)}{\lambda^{2}-\omega_{n}^{2}},\hskip
14.22636pt\left|\lambda^{2}-\omega_{n}^{2}\right|>\delta$ (31)

and its kernel

${\cal
L}\Big{[}A^{\prime}\cos(\omega_{n}t)+B^{\prime}\sin(\omega_{n}t)\Big{]}=0,\hskip
14.22636pt\left|\lambda^{2}-\omega_{n}^{2}\right|\leq\delta,$ (32)

for arbitrary constants $A,B,A^{\prime},B^{\prime}$, where we have great
freedom to choose the two parameters $\lambda>0$ and $\delta\geq 0$. Note that
both of $\cal L$ and ${\cal L}^{-1}$ are linear, say,

${\cal L}\Big{[}A\cos(\omega_{n}t)+B\sin(\omega_{n}t)\Big{]}=A\;{\cal
L}\Big{[}\cos(\omega_{n}t)\Big{]}+B\;{\cal L}\Big{[}\sin(\omega_{n}t)\Big{]},$
(33) ${\cal
L}^{-1}\Big{[}A^{\prime}\cos(\omega_{n}t)+B^{\prime}\sin(\omega_{n}t)\Big{]}=A^{\prime}\;{\cal
L}^{-1}\Big{[}\cos(\omega_{n}t)\Big{]}+B^{\prime}\;{\cal
L}^{-1}\Big{[}\sin(\omega_{n}t)\Big{]},$ (34)

for arbitrary constants $A,B,A^{\prime}$ and $B^{\prime}$. The above
definitions are complete, according to the so-called “solution-expression”
(2).

Let

$W_{\lambda,\delta}=\Big{\\{}\omega_{n}:\left|\lambda^{2}-\omega_{n}^{2}\right|\leq\delta\Big{\\}}$
(35)

denote a set containing all frequencies $\omega_{n}$ that satisfies
$\left|\lambda^{2}-\omega_{n}^{2}\right|\leq\delta$ for a given pair of
$\lambda$ and $\delta$, where $\omega_{n}=(2n-1)\Omega$ is the frequency
defined by (3). Then, its inverse operator (31) and the kernel (32) of the
auxiliary linear operator $\cal L$ can be rewritten by

${\cal
L}^{-1}\Big{[}A\cos(\omega_{n}t)+B\sin(\omega_{n}t)\Big{]}=\frac{A\cos(\omega_{n}t)+B\sin(\omega_{n}t)}{\lambda^{2}-\omega_{n}^{2}},\hskip
14.22636pt\omega_{n}\notin W_{\lambda,\delta}$ (36)

and

${\cal
L}\Big{[}A^{\prime}\cos(\omega_{n}t)+B^{\prime}\sin(\omega_{n}t)\Big{]}=0,\hskip
14.22636pt\omega_{n}\in W_{\lambda,\delta},$ (37)

respectively, where $A,B,A^{\prime}$ and $B^{\prime}$ are arbitrary constants.
Assume that $W_{\lambda,\delta}$ contains $\kappa$ frequencies. Then, the
kernel of the auxiliary linear operator ${\cal L}$ defined by (31) and (32) is
a vector space in dimension $\mu=2\kappa$, say,

$\ker\left[{\cal L}\right]=\sum_{\omega_{n}\in
W_{\lambda,\delta}}\Big{[}A_{n}\cos(\omega_{n}t)+B_{n}\sin(\omega_{n}t)\Big{]}.$
(38)

for arbitrary constants $A_{n}$ and $B_{n}$. Note that, in the frame of the
HAM, we have great freedom to choose the dimension $\mu$ that is determined by
$\lambda$ and $\delta$, as mentioned below. This is completely different from
the traditional mathematical theory for a second-order differential equation.
Thus, the general solution of the $m$th-order deformation equation (19) reads

$\displaystyle u_{m}(t)$ $\displaystyle=$
$\displaystyle\chi_{m}\;u_{m-1}(t)+c_{0}\;{\cal
L}^{-1}\Big{[}R_{m-1}(t)\Big{]}$ (39) $\displaystyle+$
$\displaystyle\sum_{\omega_{n}\in
W_{\lambda,\delta}}\Big{[}A_{m,n}\cos(\omega_{n}t)+B_{m,n}\sin(\omega_{n}t)\Big{]},$

where $A_{m,n}$ and $B_{m,n}$ are unknown constants. Note that ${\cal
L}^{-1}\Big{[}R_{m-1}(t)\Big{]}$ can be directly obtained using the inverse
operator defined by (31) or (36), and the term $u_{m-1}(t)$ is known. The
unknown constants $A_{m,n}$ and $B_{m,n}$ are determined via $R_{m}(t)$ in the
way described below.

According to the definition (31) or (36) of the inverse operator ${\cal
L}^{-1}$, $R_{m}(t)$ can not contain the terms $\cos(\omega_{n}t)$ and
$\sin(\omega_{n}t)$, where $\omega_{n}\in W_{\lambda,\delta}$, since there are
no definitions on them. Substituting $u_{m}$ defined by (39) into $R_{m}(t)$
defined by (24), we have

$R_{m}(t)=\sum\Big{[}Q_{m,n}\cos(\omega_{n}t)+S_{m,n}\sin(\omega_{n}t)\Big{]}.$
(40)

To avoid the appearance of $\cos(\omega_{n}t)$ and $\sin(\omega_{n}t)$ terms
in the above expression, where $\omega_{n}\in W_{\lambda,\delta}$, we had to
enforce the following coefficients to be zero:

$Q_{m,n}=0,\hskip 14.22636ptS_{m,n}=0,\hskip 28.45274pt\omega_{n}\in
W_{\lambda,\delta},$ (41)

which give us $\mu=2\kappa$ linear algebraic equations that determine the
$\mu=2\kappa$ unknown coefficients $A_{m,n}$ and $B_{m,n}$ of $u_{m}$ defined
by (39). In this way, we successively gain the solution $u_{m}(t)$ of the
$m$th-order deformation equation (19), where $m=1,2,3$ and so on, without any
small denominators.

According to the definition (31) or (36) of the linear inverse operator ${\cal
L}^{-1}$, we can always choose a proper pair of $\lambda$ and $\delta$ so as
to avoid the so-called “small denominator problem”, as illustrated below in §
4. Note that it is the HAM Liao2007SAM ; Liao2016NA  that provides us such
kind of great freedom.

###  3.2 Choice of the initial guess

The perturbation method provides us the unique initial guess (12), which
unfortunately contains the small denominator
$(1-\omega^{2})^{2}+4\xi^{2}\omega^{2}$ when the frequency $\omega$ is close
to 1 and $\xi$ is small. It is well-known that perturbation approaches become
invalid when small denominators appear. In addition, in the frame of
perturbation techniques, there is no freedom to choose initial guess. So, we
had to abandon the initial guess (12) of the perturbation method.

In the frame of the HAM, for the limit cycle of the Duffing equation (1), all
approximations should be in the form of (2), called the “solution expression”.
So, the initial guess should agree with the “solution expression” (2).
Besides, it should contain at least one or two lowest frequencies, such as
$\cos(\omega_{1}t),\sin(\omega_{1}t)$, $\cos(\omega_{2}t)$ and
$\sin(\omega_{2}t)$ for the forced Duffing equation (1). In addition, it
should contain the kernel (38) of the auxiliary linear operator $\cal L$, too.
For the forced Duffing equation (1), considering the solution-expression (2)
and the kernel (38) of the auxiliary linear operator $\cal L$, we choose the
initial guess in the form

$\displaystyle u_{0}(t)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{\gamma}\Big{[}a_{0,j}\cos(\omega_{j}t)+b_{0,j}\sin(\omega_{j}t)\Big{]}$
(42) $\displaystyle+$ $\displaystyle\sum_{\omega_{n}\in
W_{\lambda,\delta}}\Big{[}A_{0,n}\cos(\omega_{n}t)+B_{0,n}\sin(\omega_{n}t)\Big{]},$

where $\gamma\geq 1$ is an integer, and $a_{0,j},b_{0,j},A_{0,n},B_{0,n}$ are
unknown constants. Note that the same terms in the above expression should be
combined. All of these unknown constants in (42) are determined in the way
mentioned below.

According to (24), $R_{0}(t)={\cal N}[u_{0}(t)]$ denotes the residual error of
the forced Duffing equation (1) for the initial guess $u_{0}$. So, using the
initial guess (42), we have

$R_{0}=\sum\Big{[}Q_{0,n}\cos(\omega_{n}t)+S_{0,n}\sin(\omega_{n}t)\Big{]}.$
(43)

To avoid the appearance of the terms $\cos(\omega_{n}t),\sin(\omega_{n}t)$ in
the above expression, where $\omega_{n}\in W_{\lambda,\delta}$ as defined by
(35), we had to enforce the following coefficients to be zero:

$Q_{0,n}=0,\hskip 14.22636ptS_{0,n}=0,\hskip 28.45274pt\omega_{n}\in
W_{\lambda,\delta},$ (44)

which provides us $\mu=2\kappa$ algebraic equations. Besides, if necessary, we
had better enforce the disappearance of the base functions with the lowest
frequencies, such as $\cos(\omega_{1}t),\sin(\omega_{1}t)$,
$\cos(\omega_{2}t)$ and $\sin(\omega_{2}t)$, in the above expression of
$R_{0}(t)$. In this way, all unknown constants in the initial guess (42) could
be gained. Note that it is a set of nonlinear algebraic equations with a few
unknowns, which can be solved by means of some well-known symbolic computation
software, such as the commends FindRoot and NSolve of mathematica.

Assume that, in the iteration approach of the HAM mentioned in § 2, we have a
known approximation $u^{*}(t)$. Then, we choose the initial guess

$u_{0}(t)=u^{*}(t)+\sum_{\omega_{n}\in
W_{\lambda,\delta}}\Big{[}A_{0,n}\cos(\omega_{n}t)+B_{0,n}\sin(\omega_{n}t)\Big{]},$
(45)

where the unknown constants $A_{0,n},B_{0,n}$ are determined by enforcing the
disappearance of the kernel terms of ${\cal L}$ in $R_{0}(t)$ in the similar
way as mentioned above.

In this way, we can avoid the “small denominator problem” in the initial guess
$u_{0}(t)$, which occurs for the perturbation initial guess (12) when
$\Omega\to 1$ and $\xi\to 0$. More importantly, a set of nonlinear algebraic
equations often has multiple solutions, which might lead to multiple solutions
of the limit-cycle for the forced Duffing equation (1), as described below.

Finally, we should emphasize that it is the HAM that provides us great freedom
to choose the initial guess $u_{0}$.

##  4 Some examples

In this section, let us use the forced Duffing equation (1) to illustrate the
validity and novelty of the HAM approach mentioned in § 2 and § 3. Without
loss of generality, let us consider the case of $\alpha=1$, $\Omega=1/3$ but
various values of $\beta$ and $\xi$. Note that $\beta$ is a measurement of the
nonlinearity of the forced Duffing equation (1): the larger the value of
$\beta$, the higher the nonlinearity of the forced Duffing equation (1).

Since $\Omega=1/3$ is fixed, we always have $\omega_{2}=3\Omega=1$. Thus, from
the viewpoint of perturbation techniques, the so-called “small denominator
problem” happens when $\xi$ is small, such as $10^{-4}\leq\xi\leq 10^{-2}$, so
that the perturbation method fails, as mentioned in § 1. However, we
illustrate here that such kind of “small denominator problem” never appears in
the frame of the HAM approach, as long as we properly choose a pair of
$\lambda$ and $\delta$.

Since $\Omega=1/3$ and $\omega_{n}=(2n-1)\Omega$, we have

$\omega_{1}=\frac{1}{3},\;\omega_{2}=1,\;\omega_{3}=\frac{5}{3},\;\omega_{4}=\frac{7}{3},\;\omega_{5}=3,\;\omega_{6}=\frac{11}{3},\;\cdots$
(46)

in this paper. In theory, there are an infinite number of ways to choose the
values of $\lambda$ and $\delta$. In this section, we just consider the
following two cases:

  1. (a)

$\lambda=\sqrt{2}$ and $\delta=0$;

  2. (b)

$\lambda=\omega_{1}$ and
$\delta=\Big{|}\omega_{1}^{2}-\omega^{2}_{\kappa}\Big{|}$ with $\kappa\geq 1$.

All of them can completely avoid the so-called “small denominator problem”.

###  4.1 In case of $\lambda=\sqrt{2}$ and $\delta=0$

Table 1: Residual error square of the forced Duffing equation (1) for $u(t)$
at different order of approximations in case of $\alpha=1,\beta=1,\Omega=1/3$
and different values of $\xi$, given by the HAM approach in the case of
$\lambda=\sqrt{2}$ and $\delta=0$ described in § 4.1.

Order of $\xi=0$ $\xi=10^{-4}$ $\xi=0.01$ $\xi=0.1$ approximation
$c_{0}=-9/10$ $c_{0}=-9/10$ $c_{0}=-9/10$ $c_{0}=-8/10$ 0 1.2E-3 1.2E-3 1.2E-3
4.9E-4 1 4.1E-4 4.1E-4 4.1E-4 4.4E-4 3 5.3E-6 5.3E-6 5.2E-6 2.5E-5 5 2.4E-7
2.4E-7 2.5E-7 2.1E-6 10 9.5E-11 9.5E-11 1.0E-10 4.6E-9 15 4.0E-14 4.0E-14
4.7E-14 1.4E-11 20 4.2E-17 4.2E-17 4.8E-17 5.1E-14 25 5.7E-20 5.7E-20 6.7E-20
2.1E-16 30 4.1E-23 4.1E-23 5.3E-23 9.1E-19

---

Figure 1: Residual error square of the forced Duffing equation (1) for $u(t)$
at different orders of approximation versus the convergence-control parameter
$c_{0}$ in case of $\alpha=1,\beta=1,\Omega=1/3$ with different values of
$\xi$, given by the HAM approach in the case of $\lambda=\sqrt{2}$ and
$\delta=0$ described in § 4.1. Solid line: 1st-order approximation; Dashed-
line: 2nd-order approximation; Dash-dotted line: 3rd-order approximation. (a)
$\xi=10^{-4}$ using the initial guess (53); (b) $\xi=0.1$.

---

Figure 2: Convergent series solutions of limit-cycle of the forced Duffing
equation (1) in cases of $\alpha=1,\beta=1,\Omega=1/3$ and different values of
$\xi$. Symbols: 10th-order approximation; Solid line: 30th-order
approximation. (a) $\xi=10^{-4}$ using $c_{0}=-0.9$ and the initial guess
(53); (b) $\xi=1/10$ using $c_{0}=-0.8$.

---

Figure 3: Residual error square of the forced Duffing equation (1) versus the
times of iteration in case of $\alpha=1,\beta=1$ and $\Omega=1/3$, given by
the HAM iteration approach in the case of $\lambda=\sqrt{2}$ and $\delta=0$
described in § 4.1. Cycle: the 1st-order HAM iteration; Square: the 2nd-order
HAM iteration; Delta: the 3rd-order HAM iteration. (a) $\xi=0.0001$ using
$c_{0}=-0.9$ and the initial guess (53); (b) $\xi=0.1$ using $c_{0}=-0.8$.

In case of $\lambda=\sqrt{2}$, according to the definitions (31) and (3), we
have the denominator

$\lambda^{2}-\omega_{n}^{2}=2-\frac{(2n-1)^{2}}{9}=\frac{17+4n-4n^{2}}{9}.$
(47)

Setting $|\lambda^{2}-\omega_{n}^{2}|=\delta=0$, we have the two solutions
$n\approx 2.62132$ and $n\approx-1.62132$, which however are not positive
integers. Thus, nothing belongs to the set $W_{\lambda,\delta}$ when we choose
$\lambda=\sqrt{2}$ and $\delta=0$, say, its kernel of the corresponding
auxiliary linear operator $\cal L$ has nothing, say, $\ker[{\cal
L}]=\emptyset$, corresponding to $W_{\lambda,\delta}=\emptyset$. Besides, the
corresponding denominators $(\lambda^{2}-\omega_{n}^{2})$ read

$\frac{17}{9},\;\;1,\;\;-\frac{7}{9},\;\;-\frac{31}{9},\;\;-7,\;\;-\frac{103}{9},\cdots$

which are far away from zero, and therefore all denominators are not small so
that the “small denominator problem” does not appear at all for arbitrary
values of $\alpha,\beta$ and $\xi$.

In fact, in the case of $\lambda=\sqrt{2}$ and $\delta=0$, the definitions
(31) and (32) are equivalent to such an auxiliary linear operator

${\cal L}[u]=u^{\prime\prime}+2\;u,$ (48)

with the property

${\cal
L}\Big{[}A\cos\left(\sqrt{2}\;t\right)+B\sin\left(\sqrt{2}\;t\right)\Big{]}=0,$
(49)

whose inverse operator reads

${\cal
L}^{-1}\Big{[}A^{\prime}\cos(\omega_{n}t)+B^{\prime}\sin(\omega_{n}t)\Big{]}=\frac{A^{\prime}\cos(\omega_{n}t)+B^{\prime}\sin(\omega_{n}t)}{2-\omega_{n}^{2}},\hskip
14.22636ptn\geq 1,$ (50)

where $A,B,A^{\prime},B^{\prime}$ are arbitrary constants. Note that (48) has
no relationship with the original linear operator ${\cal L}_{0}$ defined by
(9).

In general, letting $n\geq 1$ be an integer, one can choose
$\lambda=2n\Omega$, which is far away from all $\omega_{n}=(2n-1)\Omega$,
since $\lambda-\omega_{n}=\omega_{n+1}-\lambda=\Omega$. So, in theory there
are an infinite number of ways to choose a proper $\lambda$ (and $\delta=0$)
so that the so-called “small denominator problem” never appears for the forced
Duffing equation (1)!

Note that the solution of the $m$th-order deformation equation (19) does not
contain the terms $\cos\left(\sqrt{2}\;t\right)$ and
$\sin\left(\sqrt{2}\;t\right)$, since they do not agree with the solution
expression (2) and thus must be disappeared. This agrees with the conclusion
that the kernel of the corresponding auxiliary linear operator $\cal L$ is an
empty set, say, $\ker[{\cal L}]=\emptyset$. Thus, the solution of the $m$th-
order deformation equation (19) reads

$u_{m}(t)=\chi_{m}\;u_{m-1}(t)+c_{0}\;{\cal
L}^{-1}\left[R_{m-1}(t)\right],\hskip 28.45274ptm\geq 1.$ (51)

According to (42), since $\ker[{\cal L}]=\emptyset$, we choose the initial
guess in the form

$u_{0}(t)=\sum_{n=1}^{2}\left[a_{0,n}\cos(\omega_{n}t)+b_{0,n}\sin(\omega_{n}t)\right],$
(52)

where the four unknown constants $a_{0,n},b_{0,n}$ $(n=1,2)$ are determined by
enforcing the disappearance of the terms $\cos(\omega_{1}t)$,
$\sin(\omega_{1}t)$, $\cos(\omega_{2}t)$ and $\sin(\omega_{2}t)$ in
$R_{0}(t)$. For details, please refer to § 3.2.

Without loss of generality, let us first consider here the case of
$\alpha=1,\beta=1,\Omega=1/3$ and $\xi=10^{-4}$, corresponding to the “small
denominator problem” from the viewpoint of perturbation method. Following the
method described in § 3.2, we have its corresponding initial guess

$\displaystyle u_{0}(t)$ $\displaystyle=$ $\displaystyle
0.775251\cos(\omega_{1}t)-0.127485\cos(\omega_{2}t)$ (53) $\displaystyle+$
$\displaystyle 4.98191\times 10^{-5}\sin(\omega_{1}t)-5.24821\times
10^{-5}\sin(\omega_{2}t),$

which is the unique real solution of the related set of nonlinear algebraic
equations.

Although all physical parameters are given, we always have one unknown
parameter, i.e. the convergence-control parameter $c_{0}$, which has no
physical means but can guarantee the convergence of solution series given by
the HAM. To choose an optimal value of $c_{0}$, we check the residual error
squares of the first several orders of approximation, defined by (30), as
shown in Figure 1, which give us the optimal value $c_{0}\approx-0.9$ for the
considered case. It is found that, using $c_{0}=-0.9$, the corresponding
solution series indeed converge very quickly: the residual error square
decreases about 20 orders of magnitude at the 30th-order of approximation,
say, from $1.2\times 10^{-3}$ at the very beginning to $4.1\times 10^{-23}$,
as shown in Table 1. In fact, the 10th-order approximation already agrees
quite well with the 30th-order approximation, as shown in Figure 2. Similarly,
we also check the validity of our HAM approach in the case of $\xi=0$, 0.01
and $0.1$, respectively, with $\alpha=1,\Omega=1/3$ and $\beta=1$. In all of
these cases, the solution series converge rather quickly, as shown in Table 1.
It is found that, given fixed values of $\alpha,\beta$ and $\Omega$, the
series solutions for small values of $0\leq\xi\leq 0.01$ are almost the same,
and the solution series converges almost at the same rate, as shown in Table
1. This is reasonable in physics, because the small resistance coefficient
$\xi$ has a very small influence on the limit-cycle. All of these illustrate
the validity of the HAM approach mentioned above.

As mentioned in § 3.2, one $M$th-order HAM approximation can be used as a new
initial guess to gain a better approximation, and so on. As shown in Figure 3,
in case of $\alpha=1,\beta=1$, $\Omega=1/3$ and $\xi=0.0001$, using the
initial guess (53), we gain the convergent series solution by means of the
first, second and third-order HAM iteration approach, and the corresponding
residual error square decreases rather quickly: from $1.2\times 10^{-3}$ at
the beginning to $10^{-23}$ at the 16 iterations for the 2nd-order formula, or
at the 11 iterations for the 3rd-order formula, respectively. It is found that
the higher the order of the HAM iteration approach, the faster the solution
series converges. Similarly, in case of $\xi=0.1$, the iteration also
converges rather quickly, as shown in Figure 3. So, choosing an optimal
convergence-control parameter $c_{0}$ and using HAM iteration approach, we can
quickly gain the convergent series solution of the limiting cycle of the
forced Duffing equation (1) by means of the HAM approach described in § 2 and
§ 3. This illustrates the validity of our HAM iteration approach.

As shown in Figure 3, the residual error square stops decreasing at the order
of magnitude $10^{-23}$. This is mainly because, in the solution expression
(2) of the limiting cycle $u(t)$, we delete all terms $a_{n}\cos(\omega_{n}t)$
and $b_{n}\sin(\omega_{n}t)$ when $|a_{n}|<\epsilon$ and $|b_{n}|<\epsilon$,
where we choose $\epsilon=10^{-20}$ in this paper. This is specially necessary
for the HAM iteration approach, otherwise the number of the base-functions
increases exponentially so that the HAM iteration approach can not work. It is
found that, when a smaller value of $\epsilon$ such as $\epsilon=10^{-30}$ is
used, the residual error square stops decreasing at a much smaller level.
However, for the problem considered in this paper, $\epsilon=10^{-20}$ is
small enough for the cases under consideration.

Note that, from the viewpoint of perturbation method, the “small denominator
problem” occurs when $\alpha=1,\beta=1$, $\Omega=1/3$ and $0.0001\leq\xi\leq
0.01$. However, as shown in Table 1, the series solutions given by our HAM
approach converge almost in the same rate in case of $0\leq\xi\leq 0.01$.
Therefore, our HAM approach indeed can avoid the small denominators. In other
words, from the viewpoint of the HAM approach (in the case of $\lambda=2$ and
$\delta=0$), the so-called “small denominator problem” does not really exist
at all: they are just the artifacts of perturbation methods.

###  4.2 In case of $\lambda=\omega_{1}$ and
$\delta=\Big{|}\omega_{1}^{2}-\omega^{2}_{\kappa}\Big{|}$ with $\kappa\geq 1$

Table 2: Residual error square of $u(t)$ at different order of approximations
of the forced Duffing equation (1) in case of $\alpha=1,\beta=1,\Omega=1/3$
and different values of $\xi$, given by the HAM approach described in § 4.2
when $\lambda=\omega_{1}^{2}$ and $\delta=\omega_{2}^{2}-\omega_{1}^{2}$.

Order of $\xi=0$ $\xi=10^{-4}$ $\xi=0.01$ $\xi=0.1$ approximation $c_{0}=-3/2$
$c_{0}=-3/2$ $c_{0}=-3/2$ $c_{0}=-8/5$ 0 1.2E-3 1.2E-3 1.2E-3 1.1E-3 1 4.0E-4
4.0E-4 4.0E-4 4.1E-4 3 1.2E-5 1.2E-5 1.2E-5 1.7E-5 5 5.2E-7 5.2E-7 5.2E-7
1.3E-6 10 1.7E-10 1.7E-10 1.7E-10 1.4E-9 15 1.3E-13 1.3E-13 1.3E-13 2.0E-12 20
6.4E-17 6.4E-17 6.4E-17 3.6E-15 25 1.2E-19 1.2E-19 1.2E-19 7.6E-18 30 7.9E-23
7.9E-23 7.9E-23 1.7E-20

---

Figure 4: Residual error squares of approximations of $u(t)$ at different-
orders versus the convergence-control parameter $c_{0}$ in case of
$\alpha=1,\beta=1,\Omega=1/3$ with different values of $\xi$, given by the HAM
approach described in § 4.2 when $\lambda=\omega_{1}$ and
$\delta=\omega_{2}^{2}-\omega_{1}^{2}$. Solid line: 2nd-order approximation;
Dashed-line: 3rd-order approximation; Dash-dotted line: 4th-order
approximation. (a) $\xi=10^{-4}$ using the initial guess (53); (b) $\xi=0.1$.

Figure 5: Residual error square of the forced Duffing equation (1) versus the
times of iteration in case of $\alpha=1,\beta=1,\Omega=1/3$, given by the HAM
approach described in § 4.2 when $\lambda=\omega_{1}$ and
$\delta=\omega_{2}^{2}-\omega_{1}^{2}$. Cycle: the 2nd-order HAM iteration;
Square: the 3rd-order HAM iteration; Delta: the 4th-order HAM iteration. (a)
$\xi=0.0001$ using the initial guess (53) and $c_{0}=-3/2$; (b) $\xi=0.1$
using $c_{0}=-8/5$.

Let us further consider the case $\lambda=\omega_{1}$ and
$\delta=\Big{|}\omega_{1}^{2}-\omega^{2}_{\kappa}\Big{|}=\omega^{2}_{\kappa}-\omega_{1}^{2}$,
where $\kappa\geq 1$ is an integer. According to (32),
$|\lambda^{2}-\omega_{n}^{2}|\leq\delta$ leads to the following equation

$\omega_{n}^{2}-\omega_{1}^{2}\leq\omega^{2}_{\kappa}-\omega_{1}^{2},\hskip
28.45274ptn\geq 1,$ (54)

which holds for $1\leq n\leq\kappa$. Thus, the corresponding set
$W_{\lambda,\delta}$ has $\kappa$ members, say,

$W_{\lambda,\delta}=\Big{\\{}\omega_{1},\omega_{2},\cdots,\omega_{\kappa}\Big{\\}}.$
(55)

In this case, (31) and (32) are equivalent to the following definitions:

${\cal
L}^{-1}\Big{[}A\cos(\omega_{n}t)+B\sin(\omega_{n}t)\Big{]}=\frac{A\cos(\omega_{n}t)+B\sin(\omega_{n}t)}{\omega_{1}^{2}-\omega_{n}^{2}},\hskip
14.22636ptn>\kappa,$ (56)

and

${\cal
L}\left\\{\sum_{n=1}^{\kappa}\Big{[}A^{\prime}_{m,n}\cos(\omega_{n}t)+B^{\prime}_{m,n}\sin(\omega_{n}t)\Big{]}\right\\}=0,$
(57)

for arbitrary constants $A,B,A^{\prime}_{m,n},B^{\prime}_{m,n}$, where we have
great freedom to choose the value of $\kappa$. According to (57), the kernel
of the corresponding linear operator is a vector space with $\mu=2\kappa$
dimension, say,

$\ker[{\cal
L}]=\sum_{n=1}^{\kappa}\Big{[}A^{\prime}_{m,n}\cos(\omega_{n}t)+B^{\prime}_{m,n}\sin(\omega_{n}t)\Big{]}.$
(58)

Thus, in this case, the solution of the $m$th-order deformation equation (19)
reads

$u_{m}(t)=\chi_{m}\;u_{m-1}(t)+c_{0}\;{\cal
L}^{-1}\Big{[}R_{m-1}(t)\Big{]}+\sum_{n=1}^{\kappa}\Big{[}A^{\prime}_{m,n}\cos(\omega_{n}t)+B^{\prime}_{m,n}\sin(\omega_{n}t)\Big{]},$
(59)

where $A^{\prime}_{m,n}$ and $B^{\prime}_{m,n}$ are $2\kappa$ unknown
constants, which are determined by enforcing the coefficients of
$\cos(\omega_{n}t)$ and $\sin(\omega_{n}t)$ in $R_{m}(t)$ being zero, where
$1\leq n\leq\kappa$.

In this case, we choose the initial guess in the form

$u_{0}(t)=\sum_{n=1}^{\max\\{2,\kappa\\}}\Big{[}a_{0,n}\cos(\omega_{n}t)+b_{0,n}\sin(\omega_{n}t)\Big{]},$
(60)

where the unknown coefficients $a_{0,n}$ and $b_{0,n}$ are determined by
enforcing the coefficients of $\cos(\omega_{n}t)$ and $\sin(\omega_{n}t)$ in
$R_{0}(t)$ being zero, where $1\leq n\leq\max\\{2,\kappa\\}$.

It should be emphasized here that $\kappa$ can be greater than 1, since we
have great freedom to choose its value! For example, when $\kappa=2$,
$\ker[{\cal L}]$, i.e. the kernel of the corresponding auxiliary linear
operator ${\cal L}$, is a vector space of 4 dimension. When $\kappa=3$,
$\ker[{\cal L}]$ is a vector space of 6 dimension! Note that the forced
Duffing equation (1) is just a 2nd-order nonlinear differential equation. In
the frame of the perturbation method, the forced Duffing equation (1) is
transferred into an infinite number of 2nd-order linear differential
equations, as shown in § 1. Note also that, according to the traditional
mathemtical theory, the kernel of a second-order linear differential operator
is a vector space of 2 dimension only. Thus, in case of $\kappa>1$, our HAM
approach is beyond the traditional mathematical theory about differential
equations. This also indicates the novelty of our HAM approach in mathematics.

####  4.2.1 Results when $\kappa=2$

In this case, we have $\delta=\omega_{2}^{2}-\omega_{1}^{2}$ and that the
kernel of the corresponding auxiliary linear operator $\cal L$ is a vector
space of 4 dimension, say,

$\ker[{\cal
L}]=\sum_{n=1}^{2}\left[A_{n}\cos(\omega_{n}t)+B_{n}\sin(\omega_{n}t)\right],$
(61)

where $A_{n}$ and $B_{n}$ are arbitrary constants.

Without loss of generality, let us consider the case of
$\alpha=1,\beta=1,\Omega=1/3$ and $\xi=10^{-4}$. According to (60), we can
choose the same initial guess as (53), since they have the same physical
parameters. Note that, unlike all other approximation methods, the HAM
contains the so-called “convergence-control parameter” $c_{0}$, which has no
physical meanings but can guarantee the convergence of the solution series. As
shown in Figure 4, the optimal convergence-control parameter is about
$c_{0}\approx-1.5$. Using $c_{0}=-3/2$ and the initial guess (53), the
corresponding series solution converges very quickly, from $1.2\times 10^{-3}$
at the beginning to $7.9\times 10^{-23}$ at the 30th-order of approximation,
about 20 orders of magnitude less, as shown in Table 2. Similarly, the
solution series converge very quickly for $\xi=0,0.01$ and $0.1$, as shown in
Table 2. It is found that the corresponding limit-cycles given by
$\lambda=\omega_{1}$ and $\delta=\omega_{2}^{2}-\omega_{1}^{2}$ in case of
$\xi=0.0001$ and $\xi=0.1$ are exactly the same as (a) and (b) in Figure 2
given by $\lambda=2$ and $\delta=0$, respectively. In addition, the 2nd, 3rd
and 4th-order HAM iteration formulas also give convergent series solutions
rather quickly, as shown in Figure 5. It is found again that, the higher the
order of iteration formula, the faster the solution series converges.

All of these confirm the validity of the HAM approach described in § 4.2. It
is important that, in the case of $\lambda=\omega_{1}$ and
$\delta=\omega_{2}^{2}-\omega_{1}^{2}$, our HAM approach has nothing to do
with the so-called “small denominators”! In other words, the “small
denominator problem” never appears from the viewpoint of the HAM.

Note that, according to traditional mathematical theories, a linear
differential operator $\cal L$, whose kernel is the same as the vector space
of 4 dimension defined by (61), should correspond to the 4th-order
differential equation

${\cal
L}[u]=u^{(4)}+\left(\omega_{1}^{2}+\omega_{2}^{2}\right)u^{\prime\prime}+\omega_{1}^{2}\omega_{2}^{2}u=0,$
(62)

whose inverse operator reads

$\displaystyle{\cal L}^{-1}\Big{[}A\cos(\omega t)+B\sin(\omega
t)\Big{]}=\frac{A\cos(\omega t)+B\sin(\omega
t)}{(\omega^{2}-\omega_{1}^{2})(\omega^{2}-\omega_{2}^{2})},\hskip
14.22636pt\omega\notin\left\\{\omega_{1},\omega_{2}\right\\}$ (63)

for arbitrary constants $A$ and $B$. However, the above expression is
obviously different from our inverse operator (56) that looks like one for a
2nd-order linear differential equation! In fact, we even do not know how to
explicitly express the corresponding auxiliary linear operator $\cal L$ when
$\kappa=2$, say, $\lambda=\omega_{1}$ and
$\delta=\omega_{2}^{2}-\omega_{1}^{2}$ in the HAM approach described in § 4.2,
but fortunately it is unnecessary to know it in the frame of the HAM! The most
important fact is that our HAM-based approach is valid and the corresponding
solution series of the limiting cycle converge quickly, as shown in Table 2
and Figures 4 & 5. This verifies the validity and novelty of our HAM approach
mentioned above.

####  4.2.2 Results when $\kappa=3$

In this case we have $\delta=\omega_{3}^{2}-\omega_{1}^{2}$ and the kernel of
the auxiliary linear operator $\cal L$ is a vector space of 6 dimensions, say,

$\ker[{\cal
L}]=\sum_{n=1}^{3}\left[A_{n}\cos(\omega_{n}t)+B_{n}\sin(\omega_{n}t)\right],$
(64)

where $A_{n}$ and $B_{n}$ are arbitrary constants.

Besides, according to (60), the initial guess should be in the form

$u_{0}(t)=\sum_{n=1}^{3}\Big{[}a_{0,n}\cos(\omega_{n}t)+b_{0,n}\sin(\omega_{n}t)\Big{]}.$
(65)

Substituting it into the forced Duffing equation (1) and enforcing the
coefficients of the terms $\cos(\omega_{n}t)$ and $\sin(\omega_{n}t)$
$(n=1,2,3)$ being zero, we have a set of six nonlinear algebraic equations,
whose real solutions determine the six unknown constants in (65). It is
interesting that the set of these six nonlinear algebraic equations has
multiple real solutions (complex solutions have no physical meanings here) in
many cases, for example, such as $\alpha=1,\xi=0,\Omega=1/3$ but a large value
of $\beta$, i.e. $\beta=5$:

$\displaystyle a_{0,1}=0.333781,a_{0,2}=0.107352,a_{0,3}=-0.509166;$ (66)
$\displaystyle a_{0,1}=0.526136,a_{0,2}=-0.11976,a_{0,3}=0.181401;$ (67)
$\displaystyle a_{0,1}=0.482000,a_{0,2}=0.0264671,a_{0,3}=-0.200815;$ (68)

with $b_{0,1}=b_{0,2}=b_{0,3}=0$, respectively, corresponding to the three
initial guesses in the form (65).

Using the initial guess (66) and the corresponding optimal convergence-control
parameter $c_{0}=-1$, we gain a convergent series solution, shown as (a) in
Figure 6, by means of the 2nd-order HAM iteration: the residual error square
of the forced Duffing equation (1) decreases from 0.11 at the beginning to
$4.9\times 10^{-22}$ at the 30th iteration.

Using the initial guess (67) and the corresponding optimal convergence-control
parameter $c_{0}=-3/2$, we gain the convergent series solution, shown as (b)
in Figure 6, by means of the 2nd-order HAM iteration: the residual error
square of the forced Duffing equation (1) decreases from 0.11 at the beginning
to $2.1\times 10^{-30}$ at the 20th iteration.

Using the initial guess (68) and the corresponding optimal convergence-control
parameter $c_{0}=-1$, we gain the convergent series solution, shown as (c) in
Figure 6, by means of the 2nd-order HAM iteration: the residual error square
of the forced Duffing equation (1) decreases from 0.11 at the beginning to
$7.7\times 10^{-30}$ at the 30th iteration.

It is interesting that we have the three initial guesses (66) - (68) in the
case of $\kappa=3$, which give us three different limit-cycles, as shown in
(a), (b) and (c) of Figure 6, respectively. It should be emphasized that, it
is the HAM that provides us such kind of great freedom to choose the initial
guess. Note also that, in the frame of the perturbation method, there exists
the unique initial guess (12) only, and thus in theory it is impossible to
find these multiple limit-cycles by the perturbation method111In fact, as
mentioned in § 1, the unique perturbation series diverges even when
$\beta=0.012$ in case of $\Omega=1/3$ and $\alpha=1$, corresponding to a very
weak nonlinearity. This illustrates the advantages and novelty of the HAM
beyond perturbation.

Note that the forced Duffing equation (1) contains the nonlinear term $\beta
u^{3}$. So, the larger the value of $\beta$, the higher the nonlinearity of
the Duffing equation. As mentioned in § 1, the perturbation approach is
invalid even for $\beta\geq 0.012$ and small $\xi$. However, using our HAM
approach in a similar way, we can gain convergent series solution even in the
cases with rather high nonlinearity, such as $\alpha=1,\Omega=1/3,\xi=0$ and
$10\leq\beta\leq 40$, as shown in Figure 7. It is found that, when $\kappa=3$,
using the approach mentioned in § 3.2, there exist only one initial guess in
the form (65) for each $\beta\in[10,40]$, say,

$\displaystyle\beta=10:a_{0,1}=0.450482,a_{0,2}=-0.0931627,a_{0,3}=0.0779936;$
(69)
$\displaystyle\beta=20:a_{0,1}=0.371919,a_{0,2}=-0.0737289,a_{0,3}=0.0477389;$
(70)
$\displaystyle\beta=30:a_{0,1}=0.330788,a_{0,2}=-0.0646278,a_{0,3}=0.0380275;$
(71)
$\displaystyle\beta=40:a_{0,1}=0.303900,a_{0,2}=-0.0589144,a_{0,3}=0.0328625,$
(72)

where $b_{0,n}=0$ for $n=1,2,3$, corresponding to the four initial guesses in
the form (65). In all of these cases, the convergent series solutions are
obtained by means of the 2nd-order HAM iteration using a proper convergence-
control parameter $c_{0}$, as shown in Figure 7. Thus, our HAM approach is
indeed valid for high nonlinearity. Besides, it should be emphasized that the
“small denominator problem” never appears in all cases. This illustrates the
validity of the HAM approach for high nonlinearity in the case of
$\lambda=\omega_{1}$ and $\delta=\omega_{3}^{2}-\omega_{1}^{2}$ and its
advantages beyond perturbation.

---

Figure 6: Multiple limit-cycles of the forced Duffing equation (1) in cases of
$\alpha=1,\Omega=1/3,\xi=0$ and $\beta=5$, given by the 2nd-order HAM
iteration described in § 4.2 when $\lambda=\omega_{1}$ and
$\delta=\omega_{\kappa}^{2}-\omega_{1}^{2}$. (a) using the initial guess (66),
$c_{0}=-1$ and $\kappa=3$; (b) using the initial guess (67), $c_{0}=-3/2$ and
$\kappa=3$; (c) using the initial guess (68), $c_{0}=-1$ and $\kappa=3$; (d)
using the initial guess (80), $c_{0}=-1$ and $\kappa=4$.

---

Figure 7: Multiple limit-cycles of the forced Duffing equation (1) in cases of
$\alpha=1,\Omega=1/3,\xi=0$ and some large values of $\beta$, given by the
2nd-order HAM iteration described in § 4.2 when $\lambda=\omega_{1}$ and
$\delta=\omega_{3}^{2}-\omega_{1}^{2}$, corresponding to $\kappa=3$. (a)
$\beta=10$ using the unique initial guess (69) and $c_{0}=-3/2$; (b)
$\beta=20$ using the unique initial guess (70) and $c_{0}=-1/2$; (c)
$\beta=30$ using the unique initial guess (71) and $c_{0}=-1/5$; (d)
$\beta=40$ using the unique initial guess (72) and $c_{0}=-1/25$.

Note that, according to traditional mathematical theories, a linear
differential operator $\cal L$, whose kernel is the same as the vector space
of 6 dimension defined by (64), should correspond to the 6th-order
differential equation

$\displaystyle{\cal L}[u]$ $\displaystyle=$ $\displaystyle
u^{(6)}+\left(\omega_{1}^{2}+\omega_{2}^{2}+\omega_{3}^{2}\right)u^{(4)}$ (73)
$\displaystyle+$
$\displaystyle\left(\omega_{1}^{2}\omega_{2}^{2}+\omega_{1}^{2}\omega_{3}^{2}+\omega_{2}^{2}\omega_{3}^{2}\right)u^{\prime\prime}+\omega_{1}^{2}\omega_{2}^{2}\omega_{3}^{2}\;u=0,$

whose inverse operator reads

$\displaystyle{\cal L}^{-1}\Big{[}A\cos(\omega t)+B\sin(\omega t)\Big{]}$ (74)
$\displaystyle=$ $\displaystyle\frac{A\cos(\omega t)+B\sin(\omega
t)}{(\omega^{2}-\omega_{1}^{2})(\omega^{2}-\omega_{2}^{2})(\omega^{2}-\omega_{3}^{2})},\hskip
14.22636pt\omega\notin\left\\{\omega_{1},\omega_{2},\omega_{3}\right\\}$

for arbitrary constants $A$ and $B$. However, the above expression is
obviously different from our inverse operator (56) that looks like one for a
2nd-order linear differential equation! In fact, we even do not know how to
explicitly express the corresponding auxiliary linear operator $\cal L$ when
$\kappa=3$, say, $\lambda=\omega_{1}$ and
$\delta=\omega_{3}^{2}-\omega_{1}^{2}$, but fortunately it is unnecessary to
know it in the frame of the HAM. The most important fact is that our HAM-based
approach is valid and the corresponding solution series converge quickly, as
mentioned above, which verifies the validity and novelty of our HAM approach.

####  4.2.3 Results given by $\kappa=4$

---

Figure 8: Multiple limit-cycles of the Duffing equation (1) in cases of
$\alpha=1,\Omega=1/3,\xi=0$ and $\beta=40$, given by the HAM iteration
approach described in § 4.2 when $\lambda=\omega_{1}$ and
$\delta=\omega_{4}^{2}-\omega_{1}^{2}$. (a) using the initial guess (81) and
$c_{0}=-3/2$; (b) using the initial guess (82) and $c_{0}=-4/5$; (c) using the
initial guess (83) and $c_{0}=-4/5$; (d) using the initial guess (84) and
$c_{0}=-1$.

In this case, we have $\lambda=\omega_{1}$ and
$\delta=\omega_{4}^{2}-\omega_{1}^{2}$ so that the kernel of the corresponding
auxiliary linear operator $\cal L$ is a vector space of 8 dimensions, say,

$\ker[{\cal
L}]=\sum_{n=1}^{4}\left[A_{n}\cos(\omega_{n}t)+B_{n}\sin(\omega_{n}t)\right],$
(75)

where $A_{n}$ and $B_{n}$ are arbitrary constants.

According to (60), the initial guess should be in the form

$u_{0}(t)=\sum_{n=1}^{4}\Big{[}a_{0,n}\cos(\omega_{n}t)+b_{0,n}\sin(\omega_{n}t)\Big{]}.$
(76)

Substituting it into the forced Duffing equation (1) and enforcing the
coefficients of the terms $\cos(\omega_{n}t)$ and $\sin(\omega_{n}t)$
$(n=1,2,3,4)$ being zero, we have a set of eight nonlinear algebraic
equations, whose real solutions determine the eight unknown constants in (76).
It is found that the set of these eight nonlinear algebraic equations has four
real solutions in the case of $\alpha=1,\Omega=1/3$, $\xi=0$ and $\beta=5$:

$\displaystyle
a_{0,1}=0.261766,a_{0,2}=0.0766644,a_{0,3}=-0.563565,a_{0,4}=-0.0899237;\hskip
42.67912pt$ (77) $\displaystyle
a_{0,1}=0.524251,a_{0,2}=-0.119311,a_{0,3}=0.150425,a_{0,4}=0.0487805;$ (78)
$\displaystyle
a_{0,1}=0.420975,a_{0,2}=0.0508348,a_{0,3}=-0.278947,a_{0,4}=-0.0849981;$ (79)
$\displaystyle
a_{0,1}=0.104111,a_{0,2}=0.000412556,a_{0,3}=-0.00628524,a_{0,4}=1.07865,$
(80)

where $b_{0,n}=0$ for $n=1,2,3,4$, corresponding to the four initial guesses
in the form (76).

Using the initial guess (77) and the convergence-control parameter
$c_{0}=-2/3$, we gain a convergent series solution of the limiting cycle,
which is exactly the same as (a) in Figure 6, by means of the 2nd-order HAM
iteration: the residual error square of the forced Duffing equation (1)
decreases from 0.14 at the beginning to $3.6\times 10^{-30}$ at the 30th
iteration.

Using the initial guess (78) and the convergence-control parameter $c_{0}=-1$,
we gain a convergent series solution of the limiting cycle, which is exactly
the same as (b) in Figure 6, by means of the 2nd-order HAM iteration: the
residual error square of the forced Duffing equation (1) decreases from 0.021
at the beginning to $8.4\times 10^{-30}$ at the 15th iteration.

Using the initial guess (79) and $c_{0}=-3/2$, we gain a convergent series
solution of the limiting cycle, which is exactly the same as (c) in Figure 6,
by means of the 2nd-order HAM iteration: the residual error square of the
forced Duffing equation (1) decreases from 0.022 at the beginning to
$6.6\times 10^{-30}$ at the 20th iteration.

Using the initial guess (80) and $c_{0}=-1$, we gain a new convergent series
solution, shown as (d) in Figure 6, by means of the 2nd-order HAM iteration:
the residual error square of the forced Duffing equation (1) decreases from
1.48 at the beginning to $5.6\times 10^{-16}$ at the 30th iteration. It is
interesting that this is a new solution, which is however not found by means
of $\lambda=\omega_{1}$ and $\kappa=3$. This is mainly because, when
$\kappa=4$, we should solve two more nonlinear algebraic equations to gain the
initial guess than the case of $\kappa=3$. This leads to one more initial
guess that gives one more limit-cycle by means of the HAM approach described
in this paper. It seems that, in the frame of the HAM described in § 2 and §
3, the larger the value of $\delta$, the greater the possibility to find
multiple solutions (if they indeed exist). This further shows the validity and
novelty of our HAM approach beyond perturbation.

As shown in Figure 7, for a given $\beta\in[10,40]$ (corresponding to high
nonlinearity), only one limit-cycle is found by the HAM approach using
$\lambda=\omega_{1}$ and $\delta=\omega_{3}^{2}-\omega_{1}^{2}$, corresponding
to $\kappa=3$. Do multiple limit-cycles exist in high nonlinearity, say, for a
large $\beta$? Without loss of generality, let us consider the case of
$\alpha=1,\Omega=1/3,\xi=0$ and $\beta=40$. It is found that, when $\kappa=4$,
we have the four corresponding initial guesses in the form of (76):

$\displaystyle
a_{0,1}=0.131123,a_{0,2}=-0.00645067,a_{0,3}=-0.0312269,a_{0,4}=0.336815;\hskip
42.67912pt$ (81) $\displaystyle
a_{0,1}=0.287691,a_{0,2}=-0.0231836,a_{0,3}=-0.0268927,a_{0,4}=0.0697927;$
(82) $\displaystyle
a_{0,1}=0.296537,a_{0,2}=-0.0580039,a_{0,3}=0.0566195,a_{0,4}=-0.0810271;$
(83) $\displaystyle
a_{0,1}=0.261690,a_{0,2}=-0.0276672,a_{0,3}=0.0714539,a_{0,4}=-0.156464;$ (84)

where $b_{0,n}=0$ $(n=1,2,3,4)$.

Using the initial guess (81) and $c_{0}=-3/2$, we gain a convergent series
solution, shown as (a) in Figure 8, by means of the 5th-order HAM iteration:
the residual error square of the forced Duffing equation (1) decreases from
0.26 at the beginning to $10^{-21}$ at the 10th iterations.

Using the initial guess (82) and $c_{0}=-4/5$, we gain a convergent series
solution, shown as (b) in Figure 8, by means of the 5th-order HAM iteration:
the residual error square of the forced Duffing equation (1) decreases from
0.016 at the beginning to $2.7\times 10^{-17}$ at the 10th iterations. This
limit-cycle is exactly the same as (d) in Figure 7, the only one limit-cycle
when $\beta=40$ given by means of $\lambda=\omega_{1}$ and
$\delta=\omega_{3}^{2}-\omega_{1}^{2}$, corresponding to $\kappa=3$.

Using the initial guess (83) and $c_{0}=-4/5$, we gain a convergent series
solution, shown as (c) in Figure 8, by means of the 5th-order HAM iteration:
the residual error square of the forced Duffing equation (1) decreases from
0.013 at the beginning to $8.1\times 10^{-17}$ at the 10th iterations.

Using the initial guess (84) and $c_{0}=-1$, we gain a convergent series
solution, shown as (d) in Figure 8, by means of the 5th-order HAM iteration:
the residual error square of the forced Duffing equation (1) decreases from
0.044 at the beginning to $6.4\times 10^{-21}$ at the 10th iterations.

Note that when $\kappa=3$, i.e. $\lambda=\omega_{1}$ and
$\delta=\omega_{3}^{2}-\omega_{1}^{2}$, we found only one limit-cycle for
$\beta\in[10,40]$. However, when $\kappa=4$, i.e. $\lambda=\omega_{1}$ and
$\delta=\omega_{4}^{2}-\omega_{1}^{2}$, we successfully gain four limit-cycles
in the case of $\beta=40$, corresponding to a very high nonlinearity. It seems
that, the larger the value of $\kappa$, say, the larger of $\delta$, more
limit-cycles of the forced Duffing equation (1) could be found. Note that
$\beta=40$ corresponds to a high nonlinearity: this verifies the validity of
our HAM approach for high nonlinearity. This is one of advantages of the HAM,
which has been proved in many articles (for example, please refer to Zhong and
Liao Zhong2018JFM ). In summary, all of these results illustrate the validity
and novelty of our HAM approach described in § 2 and § 3. Note that the so-
called “small denominator problem” never appears for the forced Duffing
equation (1) by means of the HAM approach.

Note that, according to traditional mathematical theories, a linear
differential operator $\cal L$, whose kernel is the same as the vector space
of 8 dimension defined by (75), should correspond to the 8th-order
differential equation

$\displaystyle{\cal L}[u]$ $\displaystyle=$ $\displaystyle
u^{(8)}+\left(\omega_{1}^{2}+\omega_{2}^{2}+\omega_{3}^{2}+\omega_{4}^{2}\right)u^{(6)}$
(85) $\displaystyle+$
$\displaystyle\left(\omega_{1}^{2}\omega_{2}^{2}+\omega_{1}^{2}\omega_{3}^{2}+\omega_{1}^{2}\omega_{4}^{2}+\omega_{2}^{2}\omega_{3}^{2}++\omega_{2}^{2}\omega_{4}^{2}++\omega_{3}^{2}\omega_{4}^{2}\right)u^{(4)}$
$\displaystyle+$
$\displaystyle\left(\omega_{1}^{2}\omega_{2}^{2}\omega_{3}^{2}+\omega_{1}^{2}\omega_{2}^{2}\omega_{4}^{2}+\omega_{1}^{2}\omega_{3}^{2}\omega_{4}^{2}+\omega_{2}^{2}\omega_{3}^{2}\omega_{4}^{2}\right)u^{\prime\prime}$
$\displaystyle+$
$\displaystyle\omega_{1}^{2}\omega_{2}^{2}\omega_{3}^{2}\omega_{4}^{2}u=0,$

whose inverse operator reads

$\displaystyle{\cal L}^{-1}\Big{[}A\cos(\omega t)+B\sin(\omega t)\Big{]}$ (86)
$\displaystyle=$ $\displaystyle\frac{A\cos(\omega t)+B\sin(\omega
t)}{(\omega^{2}-\omega_{1}^{2})(\omega^{2}-\omega_{2}^{2})(\omega^{2}-\omega_{3}^{2})(\omega^{2}-\omega_{4}^{2})},\hskip
14.22636pt\omega\notin\left\\{\omega_{1},\omega_{2},\omega_{3},\omega_{4}\right\\}\hskip
28.45274pt$

for arbitrary constants $A$ and $B$. However, the above expression is
obviously different from our inverse operator (56) that looks like one for a
2nd-order linear differential equation! In fact, we even do not know how to
explicitly express the corresponding auxiliary linear operator $\cal L$ when
$\kappa=4$, say, $\lambda=\omega_{1}$ and
$\delta=\omega_{4}^{2}-\omega_{1}^{2}$, but fortunately it is unnecessary to
know it in the frame of the HAM. The most important fact is that our HAM-based
approach is valid and the corresponding solution series of the limiting cycles
converge quickly, as mentioned above, which verifies the validity and novelty
of our HAM approach mentioned in § 2 and § 3.

##  5 Discussions and concluding remarks

First of all, when perturbation method is used to solve the forced Duffing
equation (1), the so-called “small denominator problem” is unavoidable when
$\omega\to 1$ and $\xi\to 0$, which leads to the divergence of the
perturbative series even for rather small $\beta$, corresponding to a very
weak nonlinearity. However, for the HAM approach described in § 2 and §3, such
kind of small denominators never appear for arbitrary values of physical
parameters $\alpha,\beta,\omega$ and $\xi$ so that the so-called “small
denominator problem” never occurs! Note that even in the case of large
$\beta$, corresponding to high nonlinearity, multiple limit-cycles are
successfully found by our HAM approach. All of these illustrate the validity
and novelty of the HAM approach. Thus, from the viewpoint of the HAM approach
described in this paper, the so-called “small denominator problem” does not
really exist! This suggests that whether or not the so-called “small
denominator problem” really exists should highly depend on the used method: it
indeed exists for perturbation methods, but not for the HAM! Thus, the origin
of the so-called “small denominator problem” comes from the limitations and
restrictions of perturbation method as a methodology. In other words, the
“small denominator problem” is only an artifact of perturbation method. Thus,
abandoning perturbation method but using the HAM, we can completely avoid the
“small denominator problem”. Note that the “small denominator problem” has
been regarded as a huge obstacle for many open problems in science. So, the
HAM provides us a new way to attack them.

Secondly, unlike all other approximation techniques (including perturbation
methods), we can directly define the inverse operator ${\cal L}^{-1}$ of an
undetermined linear operator in the frame of the HAM so as to easily gain the
solutions of the linear high-order equations. It should be emphasized that it
is the HAM that provides us such kind of great freedom Liao2007SAM ;
Liao2016NA . Using such kind of freedom, the so-called “small denominator
problem” can be completely avoided, as illustrated in this paper. Note that,
according to traditional mathematical theories, a linear differential operator
$\cal L$, whose kernel is a vector space of 4 dimension defined by (61),
should correspond to the 4th-order linear differential equation (62), whose
inverse operator should be expressed by (63). Similarly, a linear differential
operator $\cal L$, whose kernel is a vector space of 6 dimension defined by
(64), should correspond to the 6th-order linear differential equation (73),
whose inverse operator should be expressed by (74). In addition, a linear
differential operator $\cal L$, whose kernel is a vector space of 8 dimension
defined by (75), should correspond to the 8th-order linear differential
equation (85), whose inverse operator should be expressed by (86). However,
when $\lambda=\omega_{1}$ and $\delta=\omega^{2}_{\kappa}-\omega^{2}_{1}$,
although its kernel defined by (32) is a vector space of 4, 6, 8 dimension for
$\kappa=2,3,4$, respectively, its inverse operator defined by (31) always
looks like that of a second-order linear operator whose kernel should be a
vector space of 2 dimension according to the traditional mathematical theorms!
Obviously, the inverse operator (31), which we directly define in the frame of
the HAM, is quite different from (63), (74) and (86). Not that we even do not
know how to explicitly express its corresponding auxiliary linear operator
$\cal L$. Fortunately, it is unnecessary to know the undetermined linear
operator $\cal L$ in the frame of the HAM. Thus, to the best of author’s
knowledge, the auxiliary linear operator defined by (31) and (32) is
fundamentally different from all known traditional linear operators. Note
that, in the previous applications of the HAM Liao1992PhD ; Liao2003Book ;
Liao2012Book ; Liao1995IJNLM ; Liao1997IJNLM ; Liao1999IJNLM ; Liao2004AMC ;
Liao2007SAM ; Liao2009CNSNS ; Liao2010CNSNS ; Liao2016JFM ; Liao2020SCPMA ;
Zhong2018JFM ; Zhu2006QF ; Bouremel2007CNSNS ; Nassar2011 ; Kimiaeifar2011CMA
; Ghotbi2011 ; Duarte2015CSF ; Sardanyes2015 ; VanGorder2017 ; Pfeffer2017 ;
Cullen2019JCP ; Sultana2019EPJP ; Massa2020 ; Kaur2022JMAA ; Masjedi2022AMM ;
Botton2022AMM , one mostly chooses a proper linear auxiliary operator $\cal L$
and then find its corresponding inverse operator ${\cal L}^{-1}$ so as to
solve the high-order equations. However, in this paper, we directly define the
inverse operator ${\cal L}^{-1}$ but do not care about the explicit expression
of the corresponding auxiliary linear operator $\cal L$ at all. This might be
a breakthrough in the field of differential equations. It further illustrates
the novelty and great potential of the so-called “method of directly defining
inverse mapping” (MDDiM), which was proposed by Liao and Zhao Liao2016NA  in
the frame of the HAM and has been successfully applied to solve many types of
nonlinear equations KV2018NA ; KV2018AMC ; KV2019CNA ; thesis2019Dewasurendra
; KV2021AMNS ; KV2021CAM ; KV2022CMA ; Nave2018JBD ; Nave2018JMM ; Nave2019JMM
; Nave2020BS ; Munjam2022IJAE ; Munjam2023IJAE .

Thirdly, unlike perturbation techniques, the HAM provides us great freedom to
choose initial guesses. Using such kind of freedom, we can gain the multiple
limit-cycles of the forced Duffing equation (1) by means of the HAM. Note
that, when $\lambda=\omega_{1}$ and
$\delta=\omega^{2}_{\kappa}-\omega^{2}_{1}$ for $\kappa\geq 2$, the larger the
value of $\delta$ in the definitions (31) and (32), the greater the
probability to find more limit-cycles of the forced Duffing equation (1). In
contrast, perturbation method provides only one initial guess (12) and thus at
most one limit-cycle. Thus, this illustrates the novelty of our HAM approach
and its advantages beyond perturbation.

Unlike all other approximation techniques (including perturbation methods),
the HAM contains the so-called convergence-control parameter $c_{0}$, which
provides a simple way to guarantee the convergence of solution series even
when the nonlinearity is very high, as illustrated in this paper and also in
other publications about the HAM Zhu2006QF ; Bouremel2007CNSNS ; Nassar2011 ;
Kimiaeifar2011CMA ; Ghotbi2011 ; Duarte2015CSF ; Sardanyes2015 ; VanGorder2017
; Pfeffer2017 ; Cullen2019JCP ; Sultana2019EPJP ; Massa2020 ; Kaur2022JMAA ;
Masjedi2022AMM ; Botton2022AMM . This guarantees that the HAM-based approach
is generally valid for high nonlinearity.

As pointed out by Giorgilli Giorgilli1998 , Duffing equation “is perhaps the
simplest example of a non-integrable system exhibiting all problems due to the
small denominators”. So, although the forced Duffing equation (1) is used here
as an example to illustrate the validity and novelty of the HAM approach and
its advantages beyond perturbation, most conclusions mentioned above have
general meanings.

What will happen if the homotopy analysis method (HAM) instead of perturbation
method is first proposed by an intelligent being on a planet in the universe?
Certainly, using the HAM-based approach mentioned in this paper, this kind of
intelligent being should have no ideas of “small denominator problem” at all!
Thus, the famous “small denominator problem” does not really exist and should
be an artifact of perturbation method. Therefore, completely abandoning
perturbation methods but using the HAM-based MDDiM, we can thoroughly avoid
“small denominator problems” and besides could attack many open problems
related to small denominators. In addition, we illustrated here that a
nonlinear differential equation can be solved by directly defining a proper
inverse operator of an undetermined linear operator. Hopefully, this fact
might lead to a breakthrough in the field of differential equations.

In summary, completely abandoning perturbation methods but using the HAM-based
MDDiM, one would be never troubled by “small denominator problems”!

## Acknowledgements

This work is partly supported by National Natural Science Foundation of China
(Approval No. 12272230) and Shanghai Pilot Program for Basic Research -
Shanghai Jiao Tong University ( No. 21TQ1400202).

## Data Availability Statements

The author confirms that the data supporting the findings of this study are
available within the article.

## Author contributions

S.L. is the unique author of this article.

## Declaration of competing interest

The author declares that he has no competing financial interests.

## References

  * (1) H. Poincaré, Les éthodes Nouvelles de la Mécanique Céleste. Tome II. Méthodes de MM. Newcomb, Gyldén, Lindstedt et Bohlin, Dover Publications, New York, 1957.
  * (2) C. E. Delaunay, Théorie du Mouvement de la Lune, Gauthier-Villars, Paris, 1867.
  * (3) R. A. Pérez, A brief but historic article of Siegel, Notices of the American Mathematical Society 58 (2011) 558 – 566.
  * (4) V. I. Arnol’d, Small denominators and problems of stability of motion in classical and celestial mechanics, Russ. Math. Surv. 18 (1963) 85.  http://iopscience.iop.org/0036-0279/18/6/R02.
  * (5) G. Iooss, P. I. Plotnikov, J. F. Toland, Standing waves on an infinitely deep perfect fluid under gravity, Arch. Rational Mech. Anal. 177 (2005) 367–478.
  * (6) G. Iooss, Small divisor problems in fluid mechanics, J. Dyn. Diff. Equat. 27 (2015) 787 – 802.
  * (7) S. Marmi, A method for accurate stability bounds in a small denominator problem, J. Phys. A: Math. Gen. 21 (1988) L961.  doi:10.1088/0305-4470/21/20/001.
  * (8) L. Chierchia, C. Falcolini, Compensations in small divisor problems, Commun. Math. Phys. 175 (1996) 135 –160.
  * (9) A. Giorgilli, Small denominators and exponential stability: From Poincaré to the present time, Seminario Mat. e. Fis. di Milano 68 (1998) 19.
  * (10) A. Elias-Zúniga, A general solution of the Duffing equation, Nonlinear Dyn. 45 (2006) 227 – 235.
  * (11) E. Kartashova, Nonlinear Resonance Analysis: Theory, Computation, Applications, Cambridge University Press, 2011.
  * (12) J. K. Kevorkian, J. D. Cole, Multiple Scale and Singular Perturbation Methods, Springer New York, 1996.
  * (13) R. H. Rand, D. Armbruster, Perturbation Methods, Bifurcation Theory and Computer Algebra, 1st Edition, Springer, 1987.
  * (14) S. Liao, The Proposed Homotopy Analysis Technique for the Solution of Nonlinear Problems, Ph.D. thesis, Shanghai Jiaotong University, Shanghai, China (1992).
  * (15) S. Liao, Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman & Hall/CRC, 2003.
  * (16) S. Liao, Homotopy Analysis Method in Nonlinear Differential Equations, Springer-Verlag, 2012.
  * (17) S. Liao, An approximate solution technique not depending on small parameters: a special example, International Journal of Non-Linear Mechanic 30 (3) (1995) 371–380.  doi:10.1016/0020-7462(94)00054-E.
  * (18) S. Liao, A kind of approximate solution technique which does not depend upon small parameters – II, International Journal of Non-Linear Mechanics 32 (5) (1997) 815–822.  doi:10.1016/S0020-7462(96)00101-1.
  * (19) S. Liao, An explicit, totally analytic approximate solution for Blasius’ viscous flow problems, International Journal of Non-Linear Mechanics 34 (4) (1999) 759–778.  doi:10.1016/S0020-7462(98)00056-0.
  * (20) S. Liao, On the homotopy analysis method for nonlinear problem, Applied Mathematics and Computation 147 (2004) 499 – 513.  doi:10.1016/S0096-3003(02)00790-7.
  * (21) S. Liao, Y. Tan, A general approach to obtain series solutions of nonlinear differential equations, Studies in Applied Mathematics 119 (2007) 297 – 354.  doi:10.1111/j.1467-9590.2007.00387.x.
  * (22) S. Liao, Notes on the homotopy analysis method: some definitions and theorems, Communications in Nonlinear Science and Numerical Simulation 14 (4) (2009) 983–997.  doi:10.1016/j.cnsns.2008.04.013.
  * (23) S. Liao, An optimal homotopy-analysis approach for strongly nonlinear differential equations, Communications in Nonlinear Science and Numerical Simulation 15 (2010) 2003–2016.  doi:10.1016/j.cnsns.2009.09.002.
  * (24) S. Liao, D. Xu, M. Stiassnie, On the steady-state nearly resonant waves, J. Fluid Mech. 794 (2016) 175–199.  doi:10.1017/jfm.2016.162.
  * (25) S. Liao, A new non-perturbative approach in quantum mechanics for time-independent Schrödinger equations, Science China Physics, Mechanics & Astronomy 63 (3) (2020) 234612.  doi:10.1007/s11433-019-9430-4.
  * (26) X. Zhong, S. Liao, On the limiting Stokes wave of extreme height in arbitrary water depth, J. Fluid Mech. 843 (2018) 653–679.  doi:10.1017/jfm.2018.171.
  * (27) P. J. Hilton, An Introduction to Homotopy Theory, Cambridge University Press, 1953\.
  * (28) S. Zhu, An exact and explicit solution for the valuation of American put options, Quantitative Finance 6 (3) (2006) 229 – 242.  doi:10.1080/14697680600699811.
  * (29) Y. Bouremel, Explicit series solution for the Glauert-jet problem by means of the homotopy analysis method, Communications in Nonlinear Science and Numerical Simulation 12 (2007) 714 – 724.  doi:10.1016/j.cnsns.2005.07.001.
  * (30) C. J. Nassar, J. F. Revelli, R. J. Bowman, Application of the homotopy analysis method to the Poisson-Boltzmann equation for semiconductor devices, Commun. Nonlinear Sci. Numer. Simul. 16 (6) (2011) 2501.  doi:10.1016/j.cnsns.2010.09.015.
  * (31) A. Kimiaeifar, E. Lund, O. Thomsen, J. Sorensen, Application of the homotopy analysis method to determine the analytical limit state functions and reliability index for large deflection of a cantilever beam subjected to static co-planar loading, Computers and Mathematics with Applications 62 (2011) 4646 – 4655.  doi:10.1016/j.camwa.2011.10.050.
  * (32) A. R. Ghotbi, M. Omidvar, A. Barari, Infiltration in unsaturated soils - An analytical approach, Computers and Geotechnics 38 (2011) 777–782.  doi:10.1016/j.compgeo.2011.05.007.
  * (33) J. Duarte, C. Januário, N. Martins, Explicit series solution for a glucose-induced electrical activity model of pancreatic beta-cells, Chaos, Solitons $\&$ Fractals 76 (2015) 1–9.  doi:10.1016/j.chaos.2015.02.029.
  * (34) J. Sardanyés, C. Rodrigues, C. Januário, N. Martins, G. Gil-Gómez, J. Duarte, Activation of effector immune cells promotes tumor stochastic extinction: A homotopy analysis approach, Appl. Math. Comput. 252 (2015) 484.  doi:10.1016/j.amc.2014.12.005.
  * (35) R. A. Van Gorder, On the utility of the homotopy analysis method for non-analytic and global solutions to nonlinear differential equations, Numerical Algorithms 76 (1) (2017) 151–162.  doi:10.1007/s11075-016-0248-y.
  * (36) T. Pfeffer, L. Pollet, A stochastic root finding approach: the homotopy analysis method applied to Dyson - Schwinger equations, New J. Phys. 19 (2017) 043005.  doi:10.1088/1367-2630/aa631f.
  * (37) A. C. Cullen, S. R. Clarke, A fast, spectrally accurate homotopy based numerical method for solving nonlinear differential equations, Journal of Computational Physics 385 (2019) 106–118.  doi:10.1016/j.jcp.2019.01.057.
  * (38) J. Sultana, Obtaining analytical approximations to black hole solutions in higher-derivative gravity using the homotopy analysis method, European Physical Journal Plus 134 (111).  doi:10.1140/epjp/i2019-12616-0.
  * (39) L. Massa, J. A. Schetz, Hypersonic heat flux reconstruction with distributed temperature sensors, Journal of Thermophysics and Heat Transfer 34 (2) (2020) 331 – 346.  doi:10.2514/1.T5778.
  * (40) G. Kaur, R. Singh, H. Briesen, Approximate solutions of aggregation and breakage population balance equations, J. Math. Anal. Appl. 512 (2022) 126166\.  doi:10.1016/j.jmaa.2022.126166.
  * (41) P. K. Masjedi, P. M. Weaver, Analytical solution for arbitrary large deflection of geometrically exact beams using the homotopy analysis method, Applied Mathematical Modelling 103 (2022) 516 – 542.  doi:10.1016/j.apm.2021.10.037.
  * (42) E. d. Botton, J. B. Greenberg, A. Arad, D. Katoshevski, V. Vaikuntanathan, M. Ibach, B. Weigand, An investigation of grouping of two falling dissimilar droplets using the homotopy analysis method, Applied Mathematical Modelling 104 (2022) 486 – 498.  doi:10.1016/j.apm.2021.12.001.
  * (43) S. Liao, Y. Zhao, On the method of directly defining inverse mapping for nonlinear differential equations, Numer. Algor. 72 (2016) 989 – 1020.  doi:10.1007/s11075-015-0077-4.
  * (44) M. Baxter, M. Dewasurendra, K. Vajravelu, A method of directly defining the inverse mapping for solutions of coupled systems of nonlinear differential equations, Numer. Algor. 77 (2018) 1199 – 1211.  doi:10.1007/s11075-017-0359-0.
  * (45) M. Dewasurendra, M. Baxter, K. Vajravelu, A method of directly defining the inverse mapping for solutions of non-linear coupled systems arising in convection heat transfer in a second grade fluid, Applied Mathematics and Computation 339 (2018) 758 – 767.  doi:10.1016/j.amc.2018.07.015.
  * (46) K. Vajravelu, M. Dewasurendra, Y. Zhang, A method of directly defining the inverse mapping (MDDiM) for solutions of non-linear coupled systems arising in SIR and SIS epidemic models, Communications in Numerical Analysis 2019 (2) (2019) 64–77.  doi:10.5899/2019/cna-00385.
  * (47) M. Dewasurendra, Semi-analytical solutions of non-linear differential equations arising in science and engineering, Ph.D. thesis, University of Central Florida, Florida, USA (2019).  https://stars.library.ucf.edu/etd/6476.
  * (48) M. Dewasurendra, Y. Zhang, N. Boyette, I. Islam, K. Vajravelu, A method of directly defining the inverse mapping for a HIV infection of CD4 + T-cells model, Applied Mathematics and Nonlinear Sciences 6 (2) (2021) 469 – 482.  doi:10.2478/AMNS.2020.2.00035.
  * (49) C. W. Sahabandu, D. Karunarathna, P. Sewvandi, Z. A. M. S. Juman, M. Dewasurendra, K. Vajravelu, A method of directly defining the inverse mapping for a nonlinear partial differential equation and for systems of nonlinear partial differential equations, Computational and Applied Mathematics 40 (2021) 234.  doi:10.1007/s40314-021-01627-y.
  * (50) C. Sahabandu, M. Dewasurendra, Z. Juman, K. Vajravelu, A. J. Chamkha, Semi-analytical method for propagation of harmonic waves in nonlinear magneto-thermo-elasticity, Computers and Mathematics with Applications 105 (2022) 107–111.  doi:10.1016/j.camwa.2021.11.020.
  * (51) O. Nave, M. Elbaz, Combination of singularly perturbed vector field method and method of directly defining the inverse mapping applied to complex ODE system prostate cancer model, Journal of Biological Dynamics 12 (1) (2018) 961 – 986.  doi:10.1080/17513758.2018.1541104.
  * (52) O. Nave, A new method to find the base functions for the method of directly defining the inverse mapping (MDDiM), Journal of Multiscale Modelling 9 (4) (2018) 1850008.  doi:10.1142/S1756973718500087.
  * (53) O. Nave, A semi-analytical method for solving problems on the role of prey taxis in a biological control-mathematical model, Journal of Multiscale Modelling 10 (2) (2019) 1850009.  doi:10.1142/S1756973718500099.
  * (54) O. Nave, M. Elbaz, S. B. Mendrazitsky, Analysis of a breast cancer mathematical model by a new method to find an optimal protocol for HER2 - positive cancer, BioSystems 197 (2020) 104191.  doi:10.1016/j.biosystems.2020.104191.
  * (55) K. Gangadhar, M. V. S. Rao, S. Kumar, S. Sharma, S. R. Munjam, A series-form solution of the coupled nonlinear equations by the method of directly defined inverse mapping and SRM, International Journal of Ambient Energy 43 (1) (2022) 1345–1354.  doi:10.1080/01430750.2019.1688679.
  * (56) S. R. Munjam, K. Gangadhar, R. Seshadri, M. Rajeswar, Novel technique MDDiM solutions of MHD flow and radiative Prandtl - Eyring fluid over a stretching sheet with convective heating, International Journal of Ambient Energy(published online).  doi:10.1080/01430750.2021.1922498.