Dynamics of Some Perturbed Morse-Type Oscillators: Simulations and Applications

Abstract

The purpose of this paper is to investigate some Morse-type oscillators. In its original form, it is a model for describing the vibrations of a diatomic molecule. The Morse potential generalizes the harmonic oscillator by introducing deviations from the classical theoretical model. In the present study, we perturbed the Morse differential equation by several periodic terms based on the $\mathrm{cosine}$ function and by a damping term. The frequency is driven by different coefficients. The size of the deviations is controlled by another constant. We provide two modifications w.r.t. the damping term. The Melnikov approach is applied as an indicator of the possible chaotic opportunities. We also propose a novel approach for stochastic control of the perturbations. It is based on the assumption that the coefficients of the periodic terms are the probabilities of underlying distribution. As a result, the dynamics are driven by its characteristic function. Several applications are considered. We demonstrate some specialized modules for investigating the dynamics of the proposed models, along with the synthesis of radiating antenna patterns.

1. Introduction

A number of authors have devoted their research to the classical Morse [1] oscillator. The reader can find detailed information about the dynamics of the Morse oscillator—analytical expressions for trajectories, action–angle variables, and chaos dynamics—in the magnificent studies [1,2,3,4,5,6]. The driven Morse oscillator is an equation frequently used in theoretical chemistry to describe the photodissociate of molecules [5]. The equation is given by

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}}\hfill & =& y\hfill \\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}\hfill & =& -a(e^{-x}-e^{-2x})+ϵg_{1}cos(\omega t)\hfill \end{array}\right.$$

(1)

with initial conditions $x(t_0)=x_0$, $y(t_0)=y_0$, where $a,g_1,\omega >0$ and $0\le ϵ<1$. From the Hamiltonian, taking $H(x,y)$ as a fixed value h, we can see that $y=\pm \sqrt{2h-a(-2e^{-x}+e^{-2x})}$. The homoclinic orbit is the solution for which $h=0$ (see Figure 1). For $ϵ=0$, $(x,y)=(\infty ,0)$ is a nonhyperbolic fixed point of (1) that is connected to itself by a homoclinic orbit. We would like to apply Melnikov’s theory [7] to (1) in order to see if (1) has horseshoes; however, the fixed point having the homoclinic orbit is nonhyperbolic. Therefore, the classical theory does not immediately apply. Schecter [8] has extended Melnikov’s method so that it applies to nonhyperbolic fixed points. The following transformation of variables is known:

$$x=-2logu;y=v;{\displaystyle \frac{\mathrm{d}s}{\mathrm{d}t}}=-{\displaystyle \frac{u}{2}}$$

(2)

(see, for example, Ref. [3] (p. 722)), where the resulting equation with new variables is

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}s}}\hfill & =& v\hfill \\ {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}s}}\hfill & =& a(2u-2u^3)+ϵ\left(-{\displaystyle \frac{2g_1}{u}}cos(\omega t(s))\right),\hfill \end{array}\right.$$

(3)

where $t(s)=-2\int {\displaystyle \frac{\mathrm{d}s}{u(s)}}$.

For $ϵ=0$, Equation (3) has a hyperbolic fixed point at the origin that is connected to itself by a homoclinic orbit (see Figure 2).

$$\begin{array}{c}{u}_0(s,s_0)=\sqrt{2}\mathrm{sech}(\sqrt{2a}(s-s_0)),\hfill \\ v_0(s,s_0)=-2\sqrt{a}\mathrm{sech}(\sqrt{2a}(s-s_0))\tanh (\sqrt{2a}(s-s_0)).\hfill \end{array}$$

(4)

This problem was originally solved by Bruhn [6]. The transformation (2) is known as a “McGehee transformation” [9]. In [10], the author investigated the following model:

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}\hfill & =& y\hfill \\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}\hfill & =& -a(e^{-x}-e^{-2x})-ϵ\left(Ay-g_{1}cos(\omega t)\right).\hfill \end{array}\right.$$

In our article, we have used Melnikov’s approach [7] on a novel modified Morse type oscillator to identify the existence of chaotic behavior into it. We developed specific modules for the dynamics research of these hypothetical oscillators. The obtained outcomes can be useful as an integral part of application for numerical computations—for more specifics, see [11,12,13,14]. Some interesting results can be found in [15,16,17].

To introduce a specific structure for the perturbations, we assume that their components are governed by a discrete stochastic distribution via its probabilities. As a particular examples, we investigated several examples based on binomial, $\beta$-binomial, and geometric distributions.

The plan of the paper is as follows. We state our models in Section 2 and investigate them in light of Melnikov’s approach. In addition, we discuss the considered models driving the perturbations through a stochastic distribution. Some possible applications are provided in Section 3. We conclude in Section 4. For specialists working in the field of antenna array theory, we provide in Appendix A an explicit presentation of the Melnikov function for what we tentatively call Model B in this particular case.

2. Some New Models — Some Simulations

2.1. Model A

We consider the following new modified model of the form:

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}}\hfill & =& y\hfill \\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}\hfill & =& -a(e^{-x}-e^{-2x})-ϵ\left(Ay-{\displaystyle \sum _{i=1}^N}{g}_{i}cos(i\omega t)\right),\hfill \end{array}\right.$$

(5)

with initial conditions $x(t_0)=x_0$, $y(t_0)=y_0$, where $a,A,g_i,\omega >0$, $0\le ϵ<1$ and $N\ge 2$ is integer.

Example 1.

For a given $N=2$, $A=0.2$, $\omega =0.1$, $g_1=0.1$, $g_2=0.5$, $a=0.1$, $ϵ=0.0085$, the simulations of system (5) for $x_0=0.3$, $y_0=0.1$ are depicted in Figure 3.

Example 2.

For a given $N=3$, $A=0.2$, $\omega =0.1$, $g_1=0.1$, $g_2=0.5$, $g_3=0.9$, $a=0.1$, $ϵ=0.0085$, the simulations of system (5) for $x_0=0.3$, $y_0=0.1$ are depicted in Figure 4.

Example 3.

For a given $N=4$, $A=0.5$, $\omega =0.08$, $g_1=0.07$, $g_2=0.4$, $g_3=0.81$, $g_4=0.4$, $a=0.1$, $ϵ=0.0085$, the simulations of system (5) for $x_0=0.3$, $y_0=0.1$ are depicted in Figure 5.

2.1.1. Considerations in Light of Melnikov’s Approach

The Melnikov function provides a measure of the leading order distance between the stable and unstable manifolds when $ϵ\ne 0$ and can be used to tell where the stable and unstable manifolds intersect transversely.

We can prove the following proposition.

Proposition 1.

The roots of Melnikov function $M(t_0)$ are given as solutions of the nonlinear equation

$$M(t_0)=2\pi \left(A\sqrt{a}+{\displaystyle \sum _{i=1}^N}{g}_{i}sin(i\omega t_0)e^{{\displaystyle \frac{-i\omega }{2\sqrt{a}}}}\right)=0.$$

(6)

Proof. 

We see that $t(s,s_0)=t_0-{\displaystyle \frac{1}{2\sqrt{a}}}\sinh (\sqrt{2a}(s-s_0))$. Here, $t_0$ is the integration constant that corresponds to the initial time. The integral of Melnikov is given by

$$M(t_0,{\psi }_0)={\displaystyle {\int }_{-\infty }^{\infty }}\left({\displaystyle \frac{2A{v}_0^2}{u_0}}-{\displaystyle \frac{2v_0}{u}}{\displaystyle \sum _{i=1}^N}{g}_{i}cos(i\omega t(s,s_0))\right)\mathrm{d}s.$$

where the functions $u_0$ and $v_0$ are defined in Equation (4). After the substitution

$$s^{\prime }=\sinh (\sqrt{2a}(s-s_0))$$

we find

$$M(s^{\prime })={\displaystyle {\int }_{-\infty }^{\infty }}\left(4A\sqrt{a}{\displaystyle \frac{{(s^{\prime })}^2}{{(1+{(s^{\prime })}^2)}^2}}+{\displaystyle \frac{2s^{\prime }}{1+{(s^{\prime })}^2}}{\displaystyle \sum _{i=1}^N}{g}_{i}cos(i\omega t(s^{\prime }))\right)\mathrm{d}{s}^{\prime },$$

where $t(s^{\prime })=t_0-{\displaystyle \frac{s^{\prime }}{2\sqrt{a}}}$.

Evidently,

$${\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle \frac{{(s^{\prime })}^2}{{(1+{(s^{\prime })}^2)}^2}}\mathrm{d}{s}^{\prime }={\displaystyle \frac{\pi }{2}}$$

$${\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle \frac{s^{\prime }}{1+{(s^{\prime })}^2}}cos({\displaystyle \frac{N\omega s^{\prime }}{2\sqrt{a}}})\mathrm{d}{s}^{\prime }=0$$

$${\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle \frac{s^{\prime }}{1+{(s^{\prime })}^2}}sin({\displaystyle \frac{N\omega s^{\prime }}{2\sqrt{a}}})\mathrm{d}{s}^{\prime }=\pi e^{{\displaystyle \frac{-N\omega }{2\sqrt{a}}}}.$$

From

$$\begin{array}{ccc}cos(N\omega t(s^{\prime }))\hfill & =& cos(N\omega (t_0-{\displaystyle \frac{s^{\prime }}{2\sqrt{a}}}))\hfill \\ & =& cos(N\omega t_0)cos({\displaystyle \frac{N\omega s^{\prime }}{2\sqrt{a}}})+sin(N\omega t_0)sin({\displaystyle \frac{N\omega s^{\prime }}{2\sqrt{a}}})\hfill \end{array}$$

we have

$$M(t_0)=2\pi \left(A\sqrt{a}+{\displaystyle \sum _{i=1}^N}{g}_{i}sin(i\omega t_0)e^{{\displaystyle \frac{-i\omega }{2\sqrt{a}}}}\right).$$

□

This completes the proof of Proposition 1.

Note. We will explicitly note that in the proof of Proposition 1 we have used methodological studies detailed in works [1,10].

From a numerical point of view, the task of finding the roots of $M(t_0)$ are more interesting given that the parameters appearing in the proposed differential model are subject to a number of restrictions from the field of theoretical chemistry.

If $M(t_0)=0$ and ${\displaystyle \frac{\mathrm{d}M(t_0)}{\mathrm{d}{t}_0}}\ne 0$ for some $t_0$ and some sets of parameters, then chaos occurs.

Clearly, in the case $N=1$, the Melnikov function has ”simple zeros”, indicating the existence of homoclinic orbit and the associated Smale horseshoe type chaos. For details, see [1].

For example, fixing $N=2$, the equation $M(t_0)=0$ is depicted in Figure 6.

For $A=1.01,\omega =0.4,a=0.5,g_1=0.59,g_2=0.95$, the roots are $12.4817$ and $14.3304$ (see Figure 6a).

For $A=1.01,\omega =0.4,a=0.5,g_1=0.32,g_2=0.95$, the root with multiplicity two is $13.453$ (Figure 6b).

For $A=1.01,\omega =0.1,a=0.5,g_1=0.1,g_2=0.5$$M(t_0)=0$, there are no roots (see Figure 6c).

Fixing $N=3$, $A=0.1,\omega =0.5,a=0.5,g_1=0.59,g_2=0.95,g_3=0.4$, the equation $M(t_0)=0$ is depicted in Figure 7.

Fixing $N=4$, $A=0.1,\omega =0.5,a=0.5,g_1=0.59,g_2=0.95,g_3=0.4$, $g_4=0.9$, the equation $M(t_0)=0$ is depicted in Figure 8.

From Proposition 1, for any fixed N, the reader may formulate the Melnikov criterion for the appearance of the intersection between the perturbed and unperturbed separatrixes.

2.1.2. Probability Distributions Driving the Perturbation

Suppose now that the coefficients $g_i$ that drive the perturbation in oscillator (5) are the probabilities of some discrete distribution supported on the set $\left\{1,2,\dots ,N\right\}$. We can assume this without loss of generality after scaling ${\displaystyle \sum _{k=1}^N}{g}_k$ for $g_k>0$. Let $\xi$ be a random variable distributed under this law and $\psi \left(·\right)$ be its characteristic function. Some most popular choices are the discrete uniform, binomial, $\beta$-binomial, and hypergeometric distributions. For some similar constructions, we refer to [14]. In addition, if we obtain an infinite but countable set of $g_k$’s imposing the condition ${\displaystyle \sum _{k=1}^{\infty }}{g}_k<\infty$, then we can use a discrete distribution stated on the infinite support $\left\{1,2,3,\dots \right\}$—for example, the geometric distribution.

Having in mind the well-known complex presentation of the cos-function

$$cos\left(x\right)={\displaystyle \frac{e^{ix}+e^{-ix}}{2}},$$

we rewrite oscillator dynamics (5) as

$$\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}\hfill & =& {\displaystyle -a\left(e^{-x}-e^{-2x}\right)-ϵ\left(Ay-\sum _{k=1}^N{g}_{k}cos\left(k\omega t\right)\right)}\hfill \\ & =& {\displaystyle -a\left(e^{-x}-e^{-2x}\right)-ϵ\left(Ay-\sum _{k=1}^N{g}_k\frac{e^{ik\omega t}+e^{-ik\omega t}}{2}\right)}\hfill \\ & =& {\displaystyle -a\left(e^{-x}-e^{-2x}\right)-ϵ\left(Ay-\frac{1}{2}\left(\psi \left(\omega t\right)+\psi \left(-\omega t\right)\right)\right)}.\hfill \end{array}$$

Let the scaling coefficient be denoted by $K={\displaystyle \sum _{k=1}^N}{g}_k$. We present four examples in Figure 9 based in binomial, $\beta$-binomial, and geometric distributions. Their characteristic functions are

$$\begin{array}{cc}{\displaystyle {\psi }^{\mathrm{binomial}}\left(x\right)}\hfill & ={\displaystyle {\left(1-p+p{e}^{ix}\right)}^{n-1}},\hfill \\ {\displaystyle {\psi }^{\beta \text{-}\mathrm{binomial}}\left(x\right)}\hfill & ={\displaystyle {\phantom{.}}_2{F}_1\left(-n+1,\alpha ,\alpha +\beta ,1-e^{ix}\right)},\hfill \\ {\displaystyle {\varphi }^{\mathrm{geometric}}\left(x\right)}\hfill & ={\displaystyle \frac{p{e}^{ix}}{1-\left(1-p\right)e^{ix}}}.\hfill \end{array}$$

The used parameters for the provided examples are as follows:

• Binomial distribution: $n=50$, $p=0.8$, $a=0.1$, $A=10$, $\omega =0.01$, $ϵ=0.1$, $K=200$;
• $\beta$-Binomial distribution: $n=10$, $\alpha =0.2$, $\beta =0.25$, $a=10,$$A=1$, $\omega =0.01$, $ϵ=0.1$, $K=200$;
• $\beta$-Binomial distribution: $n=20$, $\alpha =0.2$, $\beta =0.25$, $a=10,$$A=1$, $\omega =0.01$, $ϵ=0.1$, $K=200$;
• Geometric distribution: $p=2$, $a=10$, $A=0.1$, $\omega =0.01$, $ϵ=0.1$, $K=200$.

Note that the $\beta$-binomial examples differ only by the number of trials—we consider $n=10$ or $n=20$. We can see in Figure 9c–f the similarities and differences in the dynamics they generate.

Let us consider Proposition 1. Using the complex form of the sin-function,

$$sin\left(x\right)={\displaystyle \frac{e^{ix}-e^{-ix}}{2i}},$$

(7)

we see that the Melnikov function (6) can be rewritten as

$$\begin{array}{cc}{\displaystyle M\left(t_0\right)}\hfill & ={\displaystyle 2\pi \left(A\sqrt{a}+\sum _{k=1}^N{g}_{k}sin\left(k\omega t_0\right)e^{\frac{-k\omega }{2\sqrt{a}}}\right)}\hfill \\ & {\displaystyle =2\pi \left(A\sqrt{a}+\sum _{k=1}^N{g}_k{e}^{\frac{-k\omega }{2\sqrt{a}}}\frac{e^{ik\omega t_0}-e^{-ik\omega t_0}}{2i}\right)}\hfill \\ & {\displaystyle =2\pi \left(A\sqrt{a}+\frac{\psi \left(\omega \left(t_0+\frac{i}{2\sqrt{a}}\right)\right)-\psi \left(\omega \left(-t_0+\frac{i}{2\sqrt{a}}\right)\right)}{2i}\right)}.\hfill \end{array}$$

Thus, Equation (6) can be rewritten as

$$\psi \left(\omega \left(t_0+{\displaystyle \frac{i}{\sqrt{2a}}}\right)\right)-\psi \left(\omega \left(-t_0+{\displaystyle \frac{i}{\sqrt{2a}}}\right)\right)=-2iA\sqrt{a}.$$

(8)

2.2. The Model B

We consider the following new modified model of the form:

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}}\hfill & =& y\hfill \\ {\displaystyle \frac{\mathrm{d}y(y)}{\mathrm{d}t}}\hfill & =& -a(e^{-x}-e^{-2x})-ϵ\left({Ay\left|y\right|}^{p-1}-{\displaystyle \sum _{i=1}^N}{g}_{i}cos(i\omega t)\right),\hfill \end{array}\right.$$

(9)

with initial conditions $x(t_0)=x_0$, $y(t_0)=y_0$, where $a,A,g_i,\omega >0$, $0\le ϵ<1$, $N\ge 2$ is integer, and $p=2,4,6,8,\dots$ is the damping exponent.

We will look at some interesting simulations of model (9):

Example 4.

For given $N=4$, $p=4$, $A=0.5$, $\omega =0.08$, $g_1=0.07$, $g_2=0.4$, $g_3=0.81$, $g_{4}0.4$, $a=0.1$, $ϵ=0.0085$, the simulations of system (9) for $x_0=0.2$, $y_0=0.15$ are depicted in Figure 10.

Example 5.

For given $N=6$, $p=6$, $A=0.5$, $\omega =0.5$, $g_1=0.91$, $g_2=0.4$, $g_3=0.91$, $g_4=0.14$, $g_5=0.4$, $g_6=0.6$, $a=0.1$, $ϵ=0.0085$, the simulations of system (9) for $x_0=0.2$, $y_0=0.15$ are depicted in Figure 11.

Example 6.

For given $N=8$, $p=8$, $A=0.5$, $\omega =0.5$, $g_1=0.9$, $g_2=0.4$, $g_3=0.8$, $g_4=0.16$, $g_5=0.38$, $g_6=0.59$, $g_7=0.32$, $g_8=0.42$, $a=0.1$, $ϵ=0.0085$, the simulations of system (9) for $x_0=0.2$, $y_0=0.15$ are depicted in Figure 12.

Considerations in Light of Melnikov’s Approach

We can prove the following propositions:

Proposition 2.

For fixed $p=2$, the roots of Melnikov function $M(t_0)$ are given as solutions of the nonlinear equation

$$M(t_0)=2aA+\pi {\displaystyle \sum _{i=1}^N}{g}_{i}sin(i\omega t_0)e^{{\displaystyle \frac{-i\omega }{2\sqrt{a}}}}=0.$$

Examples. For fixed $p=2$, $N=2$:

(i)

For $A=0.1, \omega =0.8, a=0.5, g_1=0.19, g_2=0.35$, the equation $M(t_0)=0$ is depicted in Figure 13.

(ii)

For $A=0.7, \omega =0.7, a=0.5, g_1=0.09, g_2=0.167$, the equation $M(t_0)=0$ has a root with multiplicity two (see Figure 14).

(iii)

For $A=0.7, \omega =0.2, a=0.5, g_1=0.1, g_2=0.15$, the equation $M(t_0)=0$ has no roots (see Figure 15).

Proposition 3.

For fixed $p=4$, the roots of Melnikov function $M(t_0)$ are given as solutions of the nonlinear equation

$$M(t_0)={\displaystyle \frac{4}{3}}{a}^{2}A+\pi {\displaystyle \sum _{i=1}^N}{g}_{i}sin(i\omega t_0)e^{{\displaystyle \frac{-i\omega }{2\sqrt{a}}}}=0.$$

For example, for $p=4$, $N=4$, $A=0.1,\omega =0.8,a=0.2,g_1=0.19,g_2=0.2$, $g_3=0.3,g_4=0.21$, the equation $M(t_0)=0$ is depicted in Figure 16.

Proposition 4.

For fixed $p=6$, the roots of Melnikov function $M(t_0)$ are given as solutions of the nonlinear equation

$$M(t_0)={\displaystyle \frac{16}{15}}{a}^{3}A+\pi {\displaystyle \sum _{i=1}^N}{g}_{i}sin(i\omega t_0)e^{{\displaystyle \frac{-i\omega }{2\sqrt{a}}}}=0.$$

For example, for $p=6$, $N=6$, $A=0.1,\omega =0.8,a=0.55,g_1=0.19,g_2=0.2$, $g_3=0.3,g_4=0.111$, the equation $M(t_0)=0$ is depicted in Figure 17.

The proof of Propositions 2–4 follows the ideas given in this article and will be omitted.

3. Possible Applications

Let us now focus on $M(t)$.

We define the hypothetical normalized antenna factor as follows: $M^{*}(\theta )={\displaystyle \frac{1}{D}}|M(Kcos\theta$ + $k_1)|$, where $\theta$ is the azimuth angle; $K=kd;k={\displaystyle \frac{2\pi }{\lambda }};\lambda$ is the wave length; d is the distance between emitters; and $k_1$ is the phase difference. For more details and related literature, see [18,19].

Following the classical research by Dolph–Chebyshev, we will call $M^{*}(\theta )$ the Melnikov antenna array (see also [13,14]).

Model A.

(case I) For fixed $N=1$, $A=0.4,\omega =0.689,g_1=0.9,K=6.24,k_1=0.001$, $a=0.096$ (from Proposition 1), see Figure 18.

(case II) For fixed $N=2$, $A=0.09,\omega =0.689,g_1=0.9,g_2=1.31,K=6.24$, $k_1=0.001,a=0.096$ (from Proposition 1), see Figure 19.

Model B.

For fixed $N=4$, $p=4$, $A=0.0001,\omega =0.2,g_1=0.08,g_2=0.06,g_3=0.1$, $g_4=0.05,K=15.5,k_1=0,a=0.5$ from Proposition 2, the Melnikov function is depicted in Figure 20. The Melnikov antenna array is depicted in Figure 21.

For fixed $N=6$, $p=6$, $A=0.01,\omega =0.18,g_1=0.08,g_2=0.06,g_3=0.01$, $g_4=0.11,g_5=0.04,g_6=0.07,K=14.95,k_1=0,a=0.01$ from Proposition 3, the Melnikov function is depicted in Figure 22. The Melnikov antenna array is depicted in Figure 23.

Of course, after serious consideration by specialists working in this scientific direction, the hypothetical Melnikov array proposed by us can be seen as a supplement to the Array Antenna Theory.

4. Conclusions

In this article, we propose a new modified Morse-type oscillator. Investigations in light of Melnikov’s approach is considered. We also demonstrate some specialized modules for investigating the dynamics of these oscillators. This will be included as an integral part of a planned much more general Web-based application for scientific computing.

From a numerical point of view, the task of finding the roots of $M(t_0)$ are more interesting given that the parameters appearing in the proposed differential model are subject to a number of restrictions from the field of theoretical chemistry. Numerical methods for solving nonlinear equations can be found in [20].