A class of discontinuous systems exhibit perturbed period doubling bifurcation

1.   Introduction

In the field of physics, there has been significant attention given to disturbances affecting a nonlinear dynamic system that approaches a period doubling bifurcation (PDB). This objective has been studied in such a way that the systems whose equations are either discrete or continuous. For instance, in Josephson-junction amplifiers, the basic amplifying mechanism depends on the PDB ^[1]. Recently, it has been shown that the PDB has a significant effect on subharmonics and ultra-harmonic emissions of bubble oscillations ^[2]. Further, the PDB is a type of bifurcation commonly observed in many systems corresponding to applications from different fields of science, such as biological, chemical, electrical, hydrodynamical, and optical systems, for instances see ^[3,4,5].

The differential equations with discontinuous vector fields may show complicated nonlinear phenomena, such as, hidden and nonclassical various types of bifurcations and chaos, which may occur in practical mechanical models ^[6,7]. Therefore, the exploration of novel mathematical approaches for analyzing discontinuous systems is intriguing, yet challenging. This is because traditional methods like Lyapunov-Schmidt reduction and center manifold theory are no longer directly suitable for addressing sudden alterations in the governing vector fields. There has been a lot of progress in improving the traditional methods of bifurcations to study and classify singularities, attractors, and hidden chaotic attractors in discontinuous systems, see ^{[6,8,9,10,11]}. Researchers are primarily focused on understanding the division of periodic and quasi-periodic orbits, as well as the count and arrangement of limit cycles, sliding behaviour, grazing bifurcation, and chaos in discontinuous systems defined in two or many regions separated by single or multiple switching manifolds, see ^{[12,13,14,15,16,17]}. For example, in ^[12], a smooth system with nonsmooth perturbations was shown to develop a flat cone. The Poincaré map, the Melnikov-like function, and averaging theory have been developed as theoretically important tools for characterizing the bifurcation to a solution of the perturbed discontinuous system ^[18,19,20]. The Melnikov function, derived from the expansion of a Poincaré map, has been employed to examine the existence of crossing periodic orbits in discontinuous systems, as demonstrated in previous studies ^{[9,21,22,23,24]}. The Melnikov functions have recently been developed for nonlinear ship rolling systems that use delayed feedback controls and are excited by random waves ^[25,26].

The main contribution of this work is to develop a perturbation technique, which measures the effect of the nonsmooth perturbation on the PDB of the unperturbed system. To achieve this objective, the unperturbed system is supposed to be close to PDB, such that when the bifurcation parameter $\alpha$ passes through $\tilde{\alpha}$, the response changes from a period one to a period two limit cycle. We turn next to examine initially the necessary conditions for the existence of PDB in an unperturbed system by the so-called generalized Poincaré map. In this context, we derive the generalized determination of the Poincaré map associated to perturbed systems subjected to nonsmooth transitions. This map, computed here, is characterized by the composition of sub-maps that take into account the interaction with the discontinuity surface and its associated time intersection functions. We then introduce a function involving an unknown variable, where the zeros pinpoint the PDB that remains after the nonsmooth perturbation. Finally, we present a simple and novel discontinuous system that demonstrates how a nonsmooth perturbation causes a period doubling orbit to emerge from a flat periodic orbit. Another instance is given to demonstrate the formation of a perturbed period doubling orbit near an unperturbed one.

2.   Period doubling bifurcation

In a continuous nonlinear dynamical system, the slow variation of the bifurcation parameters can lead to instability, resulting in qualitative changes in the system's long-term behaviour. For instance, one of the most abundant complex and interesting behaviour is chaos. A system is classified as chaotic if it shows sensitivity to initial conditions and exhibits an infinite number of unstable periodic trajectories with varying periods ^[27,28]. One of the main ways chaos happens in parameter space is through a series of PDB, also known as a period doubling cascade. PDB occurs when the bifurcation parameter $\alpha$ passes through $\tilde{\alpha}$, an attracting periodic cycle becomes repelling and spawns an attracting orbit having twice the period, see Figure 1. As the bifurcation parameter is further modified, this period doubling may become unstable and give birth to a sequence of period doubling cycles.

Before stating our results for non-smooth differential models, we need to define and evaluate the bifurcation process for the period one to period two transition for such map. The goal of the following theorem is to outline a set of conditions that ensure the presence of PDB in a given map. Based on the results of a classical Period doubling bifurcation theorem, for one-dimensional smooth maps (i.e., $d = 1$, see in [29, Theorem 3.5.1, p. 158]) and for the result for $d > 1$, see ^[16], we present the following theorem to define the criteria for the existence of PDB in a continuous map.

Theorem 1. Assume that the mapping $\mathcal{\tilde{P}}(\xi, \alpha): \mathbb{R}^d\times \mathbb{R} \longrightarrow \mathbb{R}^d$ is smooth (i.e., it is uniformly continuous $\mathcal{\tilde{P}}, \in \mathcal{C}^r(\mathbb{R}^d), \; r\geq1)$. Then, $\mathcal{\tilde{P}}$ has a PDB at $(\tilde{\xi}, \tilde{\alpha})$ if $\mathcal{\tilde{P}}$ satisfies the following conditions:

● The point $(\tilde{\xi}, \tilde{\alpha})$ is a critical (fixed) point of $\mathcal{\tilde{P}}$, i.e., $\mathcal{\tilde{P}}(\tilde{\xi}, \tilde{\alpha}) = \tilde{\xi}$.

● There exist $U = [u_1, u_2]\times X\subset \mathbb{R}^d$ and $\theta_i :[u_1, u_2] \rightarrow int X$ such that for $\alpha \in [\tilde{\alpha}, \alpha_2]$ there holds

where $\theta_1(\alpha), \theta_2(\alpha)$, and $\xi(\alpha)$ are fixed points of $\mathcal{\tilde{P}}^{(2)}(\xi, \alpha)$.

● The Jacobian matrix $D_\xi \mathcal{\tilde{P}}(\xi, \alpha)|_{(\tilde{\xi}, \tilde{\alpha})}$ has eigenvalues $\lambda_c(\tilde{\alpha}) = -1$ and $\lambda_\ell(\tilde{\alpha}), \ell = 1, 2, ...(d-1)$ with $|\lambda_\ell(\tilde{\alpha})|\neq 1$.

● Let us suppose $\xi(\alpha)$ represents the continuous curve of fixed points $\mathcal{\tilde{P}}$ near $\xi(\tilde{\alpha})$ and $\lambda_\ell(\alpha), \mathit{\text{and}}\ \ell = 1, 2, ..., (d-1)$ are the eigenvalues of the matrix $D_\xi \mathcal{\tilde{P}}|_{(\xi(\alpha))}$. Then, we have

Therefore, there exists a smooth curve of fixed points, denoted as $\xi(\alpha)$, passing through $\tilde{\xi}$ at $\tilde{\alpha}$ where the stability of the fixed point varies at $\tilde{\alpha}$ (according to varying of $\mu$).

Proof. Let us assume that the mapping or function $\Lambda(\xi, \alpha)$ is defined as:

where $\stackrel{c}{\tilde{\Sigma}}{}_1^{+}$ is an open subset of $\mathbb{R}^d$, and we have $\Lambda(\tilde{\xi}, \tilde{\alpha}) = 0$ and $\frac{\partial \Lambda}{\partial \xi}(\tilde{\xi}, \tilde{\alpha})\neq 0$ associated to the eigenvalue $\lambda_c(\tilde{\alpha}) = -1$. Given the second assumption that $D_\xi \mathcal{\tilde{P}}$ has no eigenvalue equal to 1, we can use that the implicit function theorem to find a single solution of $\Lambda(\xi, \alpha) = 0$, near $(\tilde{\xi}, \tilde{\alpha})$. This implies that there exists a curve of fixed points $\xi(\alpha)$ for $\Lambda$ near $\tilde{\alpha}$ such that $\xi(\tilde{\alpha}) = \tilde{\xi}$.

Moreover, the Lyapunov-Schmidt reduction ^[30] can be used to obtain a bifurcation function $F: \mathbb{R}\times \mathbb{R}\supset dom(F)\leftrightarrow \mathbb{R}$, whose roots correspond to the fixed points and period two points of $\mathcal{\tilde{P}}$. For further fundamental tools on the bifurcation problem, we assume that $\xi = (x, y)\in \mathbb{R}\times \mathbb{R}^{(d-1)},$ where $x$ is a right eigenvector of $D_\xi \mathcal{\tilde{P}}$ associated to the eigenvalue $\lambda_c(\alpha)$ (note that $\tilde{x}$ associated to $\lambda_c(\tilde{\alpha}) = -1)$. Our task here is to show that there exists a function $y = y(x, \alpha)$ defined on a set $V = [x_1, x_2] \times [\alpha_1, \alpha_2]$ with values in $Y\subset \mathbb{R}^{(d-1)}$ such that $y-\pi_y(\mathcal{\tilde{P}}^{(2)}(x, y)) = 0$ for all $(x, y, \alpha)\in V\times Y$, $(\mathcal{\tilde{P}}^{(2)}(x, y)$ refers to the second iterate of the map $\mathcal{\tilde{P}}(x, y))$. Then, if $D_x F|_{(\tilde{x}, \tilde{\alpha})}\neq 0$, we can use the implicit function theorem to find the solution set of the bifurcation function $F(x, \alpha) = x-\pi_x(\mathcal{\tilde{P}}^{(2)}(x, y(x, \alpha))) = 0, (x, \alpha)\in V$, and $\pi_x$ and $\pi_y$ are the coordinate projections onto the first and second factors of $\mathbb{R}\times \mathbb{R}^{(d-1)}$. In addition, the bifurcation function when $D_x F|_{(\tilde{x}, \tilde{\alpha})}$ has a non-trivial kernel, and then the situation can be simplified by applying the implicit function theorem to a new bifurcation function, that is $G(x, \alpha) = {\int}_{0}^{1}D_x F(t(x-\bar{x}(\alpha))+\bar{x}(\alpha), \alpha)dt.$ This approach allows us to evaluate the solution of the bifurcation function as a regular curve and determine the fixed points of $\mathcal{\tilde{P}}^{(2)}$, i.e., the points of period two. Figure 1 is a bifurcation diagram showing the dynamical behavior of $\mathcal{\tilde{P}}$ near the PDB point. It explains that the PDB occurs when a stable periodic orbit loses stability and undergoes bifurcation, resulting in the formation of a new periodic orbit with twice the period of the original orbit. □

3.   Perturbation of discontinuous vector fields

We consider the vector $\xi\in \mathbb{R}^n$, with its time dependence described by the vector fields $f_i:\mathbb{R}^{n}\longrightarrow\mathbb{R}^{n}$, i = 1, 2, and the switching function $\hbar(\xi; \epsilon)$ is a smooth function of both arguments, such that

where the phase space is divided into two domains separated by a hyperplane $\tilde{\Sigma} = \{\xi\in\mathbb{R}^{n} \mid\; \hbar(\xi; \epsilon) = 0\}$, $h(\xi; \epsilon)$ is a smooth function, and $\epsilon$ is a small parameter representing the size of the perturbation. In order to study (3.1), additional rules may be necessary to describe the interaction on the discontinuity surface $\tilde{\Sigma}$, so this ensures that the flow of (3.1) is uniquely defined in the forward direction of time. Let $\varrho(\xi; \epsilon) = \big(n^T(\xi; \epsilon) f_1(\xi; \epsilon)\big). \big(n^T(\xi; \epsilon) f_2(\xi; \epsilon)\big)$, where the normal vector $n(\xi; \epsilon)$ perpendicular to the manifold $\tilde{\Sigma}$ is given as $n(\xi; \epsilon) = \frac{\nabla \hbar(\xi; \epsilon)}{|\nabla \hbar(\xi; \epsilon)|_2}$, $\| n(\xi; \epsilon)\| = 1$. Consequently, the discontinuity surface $\tilde{\Sigma}$ can be partitioned as follows:

(ⅰ) direct crossing region $\tilde{\Sigma}{^{c}} = \{\xi\in \tilde{\Sigma} \mid \varrho(\xi; \epsilon) > 0\}$, more specifically $\tilde{\Sigma}^c = \tilde{\Sigma}^c_-\cup \tilde{\Sigma}^c_+$ such that $\tilde{\Sigma}{^{c}_\pm} = \{\xi\in \tilde{\Sigma}^c \mid \pm n^T(\xi; \epsilon) f_1(\xi; \epsilon)>0\}$.

(ⅱ) sliding region $\tilde{\Sigma}{^{s}} = \{\xi\in \tilde{\Sigma} \mid \varrho(\xi; \epsilon)\leq 0\}$, which in turn is divided into the attractive region $\tilde{\Sigma}{^{s}_-} = \{\xi\in \tilde{\Sigma}^s \mid n^T(\xi; \epsilon) f_1(\xi; \epsilon) > 0\}$, and escaping region $\tilde{\Sigma}{^{s}_+} = \{\xi\in \tilde{\Sigma}^s \mid n^T(\xi; \epsilon) f_2(\xi; \epsilon) > 0\}$. Further, $\tilde{\Sigma}^0_- = \{\xi \in \tilde{\Sigma} \; |\; \ \; n^T(\xi; \epsilon) f_2(\xi, \epsilon) = 0\}, \; \text{and} \; \tilde{\Sigma}^0_+ = \{\xi \in \tilde{\Sigma} \; |\; n^T(\xi; \epsilon) f_1(\xi, \epsilon) = 0\}$ are defined as the boundaries between sliding and crossing modes. Notice that $\tilde{\Sigma}$ is the disjoint union $\tilde{\Sigma}^c\cup \tilde{\Sigma}^0_\pm\cup \tilde{\Sigma}^s$.

The sliding flows on $\tilde{\Sigma}{^{s}}$ are defined by the Filippov convex ^[31] combination as:

Our objective is to examine the persistence of period doubling when (3.1) is perturbed within the category of all discontinuous vector fields with two regions separated by $\tilde{\Sigma}$. A crossing period doubling of system (3.1) forms a closed path $\Gamma$, consisting of trajectories of $f_i, \; i = 1, 2$ with matching directions, where $\Gamma\cap \tilde{\Sigma} = \tilde{\Sigma}^c$.

The following is a general system that describes the behavior of a dynamical system close to PDB. Instead of working in the full space of the dynamical system, the dynamics are described in the transversal crossing trajectories branch.

Assume that it is possible to rewrite (3.1) as:

where $\mathcal{X}$ and $\mathcal{W}$ are uniformly continuous and bounded functions $(\mathcal{X}, \mathcal{W}\in \mathcal{C}^r(\mathbb{R}^n\setminus \tilde{\Sigma}), \; r\geq1)$. Let $\Psi^{(k)}(t_i^{(k)}(\xi, 0), \xi; 0)$ be the solutions of the unperturbed system $(\epsilon = 0)$

Then, the solution of (3.3) can be expressed as:

where $\tau_i^{(k)}(\xi, 0) = t_i^{(k)}(\xi)$, $i = 1, 2, k = 1, 2$ are the number of solution intersections with $\tilde{\Sigma}$. The fundamental matrix solutions $\Phi^{(k)}_i(t_i^{(k)})$ are given by solving the following linearized equations about $(\tau_i^{(k)}, \xi; \epsilon) = (t_i^{(k)}, \xi^0, 0)$

where $\mathbb{I}$ represents the identity matrix with the same order as the number of state variables.

Hypothesis 1. The unperturbed system $(\epsilon = 0)$ is assumed to be close to a PDB, which means that the corresponding Poincaré map of (3.4) (this map is defined by the general solution (3.5) at $\epsilon = 0$) satisfies Theorem 1.

4.   Poincaré map with nonsmooth transitions

Without loss of generality, let us find a direction starting from the starting position $\xi\in \tilde{\Sigma}_{-}^{c}$. Then, the solutions is denoted by $\Psi_1^{(1)}(\tau^{(1)}_{1}(\xi, \epsilon), \xi; \epsilon)$ which is $\mathcal{C}^{r}, r\geq1$. Moreover, assume that $\Psi_1^{(1)}(\tau^{(1)}_{1}(\xi; \epsilon), \xi; \epsilon)$ is within $\tilde{\Sigma}_{+}^{c}$ of the minimum return time $\tau^{(1)}_{1}(\xi; \epsilon)$ at $\eta = \Psi_1^{(1)}(\tau^{(1)}_{1}(\xi; \epsilon), \xi; \epsilon)\in \tilde{\Sigma}_{+}^{c}$. Consequently, the flow given by $\Psi_2^{(1)}(\tau^{(1)}_{2}(\eta; \epsilon), \eta; \epsilon)$ is again within $\tilde{\Sigma}_{-}^{c}$. Thus, we define the first iterate of the Poincaré map as

where

The perturbed system (3.1) has a period- one orbit if and only if $\mathcal{\tilde{P}}(\bar{\xi}^{(1)}, \epsilon) = \bar{\xi}^{(1)}$ and $\tau^{(1)}_{i}$ exist. Since the period doubling orbit can be viewed as a fixed point of the second iteration of the Poincaré map, we have the mapping structure

where

It follows that the times of intersection exist for the first and second iteration of $\mathcal{\tilde{P}}^{(2)}(\xi; \epsilon)$ by solving the following equations:

According to Theorem 1, we remark that $\xi^{(1)}\neq\xi^{(2)}\neq \xi^{(3)}\neq\xi^{(4)}$, and that $\xi$ depends implicitly on $\epsilon$. Because the absence of the transversality conditions at the intersection did not allow the trajectory of (3.3) to be a submanifold of the crossing area, having some sort of singular point. Therefore, the following argument shows that the transversality conditions are crucial to guarantee that $\xi^{(l)}, l = 1, 2, 3, 4$ and $\tilde{\Pi}_i^{k}$ are not tangent to $\tilde{\Sigma}$. In addition, these are sufficient conditions for an intersection to be stable after a perturbation.

Corollary 1. Let $\tilde{\Pi}_i^{(k)}: \tilde{\Sigma}^c\rightarrow \tilde{\Sigma}^c$ be smooth maps, and $\tilde{\Sigma}^c$ a submanifold of $\tilde{\Sigma}$. Then, $\xi^{(l)}$ and $\tilde{\Pi}_i^{(k)}$ are within the direct crossing set if the following conditions hold:

In general, because of the non-linearity of (4.1), (4.3), and (4.5), it is not possible to obtain an explicit expression for the sub-maps $\tilde{\Pi}^{(k)}_{i}$. Thus, we are studying the characteristics of Poincaré sub-maps as they relate to the existence of period doubling orbit and the stability of the system (3.3).

Lemma 1. We assume that the transversality conditions that are given by the Corollary 1 hold. Then, for all $\bar{\xi}^{(i)}\in \tilde{\Sigma}^c_-$ and $\tilde{\Pi}_i^{(k)}(\bar{\xi}^{(l)}, \epsilon)\in \tilde{\Sigma}^c$, $i = 1, 2, k = 1, 2, l = 1, ..., 4$, the intersection function times $\tau_i^{(k)}(\xi^i; \epsilon)$ are differentiable to $\xi^{(l)}$, so we get:

$(a)$ $D_\xi \tilde{\Pi}^{(1)}_1(\tilde{\xi}^{(1)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(1)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(1)})}\big)\Phi_1^{(1)}(\tilde{t}_1^{(1)}, \tilde{\xi}^{(1)}),$

$(b)$ $D_\xi \tilde{\Pi}^{(1)}_2(\tilde{\xi}^{(2)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_2(\tilde{\xi}^{(2)})n^T}{n^T \mathcal{X}_2(\tilde{\xi}^{(2)})}\big)\Phi_2^{(1)}(\tilde{t}_2^{(1)}, \tilde{\xi}^{(2)}) \tilde{S}_1^{(1)} \times \big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(1)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(1)})}\big)^{-1}D_\xi \tilde{\Pi}^{(1)}_1(\tilde{\xi}^{(1)}; 0),$

$(c)$ $D_\xi \tilde{\Pi}^{(2)}_1(\tilde{\xi}^{(3)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(3)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(3)})}\big)\Phi_1^{(2)}(\tilde{t}_1^{(2)}, \tilde{\xi}^{(3)}) \tilde{S}_1^{(2)} \times \big(\mathbb{I}-\frac{ \mathcal{X}_2(\tilde{\xi}^{(2)})n^T}{n^T \mathcal{X}_2(\tilde{\xi}^{(2)})}\big)^{-1}D_\xi \tilde{\Pi}^{(1)}_2(\tilde{\xi}^{(2)}; 0),$

$(d)$ $D_\xi \tilde{\Pi}^{(2)}_2(\tilde{\xi}^{(4)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_2(\tilde{\xi}^{(4)})n^T}{n^T \mathcal{X}_2(\tilde{\xi}^{(4)})}\big)\Phi_2^{(2)}(\tilde{t}_2^{(2)}, \tilde{\xi}^{(4)}) \tilde{S}_2^{(2)} \times\big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(3)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(3)})}\big)^{-1}D_\xi \tilde{\Pi}^{(2)}_1(\tilde{\xi}^{(3)}; 0).$

Proof. Consider the partial derivatives of $\tilde{\Pi}^{(k)}_i(\xi; \epsilon)$ with respect to $\xi$ at the point $(\tau_i^{(k)}, \xi^{(l)}; \epsilon) = (\tilde{t}_i{(k)}, \tilde{\xi}^{(l)}, 0)$, where $D_\xi = \frac{\partial}{\partial \xi}$:

In order to find the explicit formulas of $D_\xi \tilde{\Pi}_i^{(k)}$, we have to compute the terms $D_\xi \tau_i^{(k)}$. The partial derivatives of $\tau_i^{(k)}$ are obtained by noting that $n^TD_\xi \tilde{\Pi}^{(k)}_i(\tilde{\xi}^{(l)}; 0) = 0$. Therefore, from $n^T D_\xi \tilde{\Pi}^{(1)}_1(\tilde{\xi}^{(1)}; 0) = 0$, we obtain $D_\xi\tau_1^{(1)}(\tilde{\xi}^{(1)}; 0) = -\frac{n^T \Phi_1(\tilde{t_1}^{(1)}, \tilde{\xi}^{(1)})}{n^T \mathcal{X}_1(\tilde{\xi}^{(1)})}$, and (a) holds.

At the intersection points, we have $D_{\tau_i^{(1)}}(\tilde{\xi}^{(i+1)}; 0) \Psi_{3-i}^{(i)} = -\Phi_{3-i}^{(i)} \mathcal{X}_{3-i}$. Then, the result of the equation $n^T D_\xi \tilde{\Pi}^{(1)}_2 = 0$ implies

Then, we obtain $D_\xi\tau_2^{(1)}(\tilde{\xi}^{(2)}; 0) = -\frac{n^T\Phi_2^{(1)}(\tilde{t_2}^{(1)}, \tilde{\xi}^2) \tilde{S}_1^{(1)}\Phi_1^{(1)}(\tilde{t_1}^{(1)}, \tilde{\xi}^1)}{n^T \mathcal{X}_2(\tilde{\xi}^2)}$, and the direct substitution into the second equation of (4.6) leads to statement (b).

In the same way, we are going to solve $n^T D_\xi \tilde{\Pi}^{(2)}_i(\tilde{\xi}^{(2+i)}; 0) = 0$ for $D_\xi\tau_2^{(i)}$, which gives: $D_\xi \tau_1^{(2)} = -\frac{n^T \Phi_1^{(2)}(\tilde{t}_1^{(2)}, \tilde{\xi}^3) \tilde{S}_2^{(1)} \Phi_2^{(1)}(\tilde{t}_2^{(1)}, \tilde{\xi}^2) \tilde{S}_1^{(1)}\Phi_1^{(1)}(\tilde{t}_1^{(1)}, \tilde{\xi}^1)}{n^T \mathcal{X}_1(\tilde{\xi}^3)}$, $D_\xi \tau_2^{(2)} = -\frac{n^T \Phi_2^{(2)}(\tilde{t_2}^{(2)}, \tilde{\xi}^{(4)}) \tilde{S}_1^{(2)} \Phi_1^{(2)}(\tilde{t_1}^{(2)}, \tilde{\xi}^3) \tilde{S}_2^{(1)} \Phi_2^{(1)}(\tilde{t_2}^{(1)}, \tilde{\xi}^2) \tilde{S}_1^{(1)} \Phi_1(\tilde{t_1}^{(1)}, \tilde{\xi}^1)}{n^T \mathcal{X}_2(\tilde{\xi}^4)}.$

Substituting $D_\xi \tau_i^{(2)}$ in the third and fourth formulas of (4.6), we will ultimately obtain formulas (c) and (d), respectively. □

The following theorem illustrates that, by linearizing the Poincaré map, we can carry out explicit analytical calculations.

Theorem 2. Assume that $\xi\in \tilde{\Sigma}^c_-$ and $\mathcal{\tilde{P}}^{(m)}: \tilde{\Sigma}^c_-\rightarrow \tilde{\Sigma}^c_-$, such that $\mathcal{\tilde{P}}^{(m)}(0;0) = 0$ and $\mathcal{\tilde{P}}^{(m)}(\tilde{\xi}; 0) = \tilde{\xi}$. Then, the linearization of the Poincaré map is given by

Further, $\mathcal{\tilde{P}}^{(m)}(\tilde{\xi}; 0) \mathcal{X}_1(\tilde{\xi}) = \mathcal{X}_1(\tilde{\xi})$ and attractivity of the period-m orbit is determined by the remaining (n-2) eigenvalues of $D_\xi \mathcal{\tilde{P}}^{(m)}$.

Proof. Assume that $\tilde{\xi}\in \tilde{\Sigma}{_-^c}$. The natural condition for period doubling orbit to exist in an unperturbed system is $\Psi_{1}^{(1)}(0, \tilde{\xi}; 0) = \Psi_{1}^{(1)}(T, \tilde{\xi}; 0)$, or in a compact mapping as

where $T = \sum_{k = 1}^{2} \sum_{i = 1}^{2}\tau_i^{(k)}(\tilde{\xi}; 0) = \sum_{k = 1}^{2} \sum_{i = 1}^{2}t_i^{(k)}(\xi)$. We assume that the second iteration of the Poincaré map is provided as:

Differentiating (4.7) with respect to $\xi$, we obtain

The above derivative is estimated at $(\tau, \xi, \epsilon) = (\tilde{t}, \tilde{\xi}, 0)$ as

Using the results of Lemma 1 and the value of the vector function $D_\xi \tau_2^{(2)}(\xi; 0)$, we obtain

At the second iteration of the Poincaré map, substituting $\tau_1^{(1)}$ into $\tilde{t}_2^{(2)}$,

In this situation, when the unperturbed system has period-$m$ orbits, it is simply a matter of multiplying the matrices $D_\xi \mathcal{\tilde{P}}^{(k)}_i$ in the correct order. Then, we obtain

□

Remark 1 It should be noted that, when $m = 1$, a similar result has obtained in ^[9,19].

Lemma 2. Suppose that $\Gamma\subset \tilde{\Sigma}_-^c$ is a period doubling orbit of the unperturbed system generated by (3.3), and the transversality conditions that are given by Corollary 1 hold. Then, for all $\bar{\xi}^{(i)}\in \tilde{\Sigma}^c_-$ and $\tilde{\Pi}_i^{(k)}(\bar{\xi}^{(l)}, \epsilon)\in \tilde{\Sigma}^c$, $i = 1, 2, k = 1, 2, l = 1, ..., 4$, and the intersection times functions $\tau_i^{(k)}(\xi^i; \epsilon)$ are differentiable to $\epsilon$, so we get:

$(a)$ $D_\epsilon \tilde{\Pi}^{(1)}_1(\tilde{\xi}^{(1)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(1)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(1)})}\big)D_\epsilon \Psi_1^{(1)}(\tilde{t}_1^{(1)}, \tilde{\xi}^{(1)}),$

$(b)$ $D_\epsilon \tilde{\Pi}^{(1)}_2(\tilde{\xi}^{(2)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_2(\tilde{\xi}^{(2)})n^T}{n^T \mathcal{X}_2(\tilde{\xi}^{(2)})}\big)\big(\Phi_2^{(1)}(\tilde{t}_2^{(1)}, \tilde{\xi}^{(2)}) \tilde{S}_1^{(1)} \times \big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(1)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(1)})}\big)^{-1}D_\epsilon \tilde{\Pi}^{(1)}_1(\tilde{\xi}^{(1)}; 0)+D_\epsilon \Psi_2^{(1)}(\tilde{t}_2^{(1)}, \tilde{\xi}^{(2)})\big),$

$(c)$ $D_\epsilon \tilde{\Pi}^{(2)}_1(\tilde{\xi}^{(3)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(3)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(3)})}\big)\big(\Phi_1^{(2)}(\tilde{t}_1^{(2)}, \tilde{\xi}^{(3)}) \tilde{S}_2^{(1)} \times\big(\mathbb{I}-\frac{ \mathcal{X}_2(\tilde{\xi}^{(2)})n^T}{n^T \mathcal{X}_2(\tilde{\xi}^{(2)})}\big)^{-1}D_\epsilon \tilde{\Pi}^{(1)}_2(\tilde{\xi}^{(2)}; 0)+D_\epsilon \Psi_1^{(2)}(\tilde{t}_1^{(2)}, \tilde{\xi}^{(3)})\big),$

$(d)$ $D_\epsilon \tilde{\Pi}^{(2)}_2(\tilde{\xi}^{(4)}; 0) = \big(\mathbb{I}-\frac{ \mathcal{X}_2(\tilde{\xi}^{(4)})n^T}{n^T \mathcal{X}_2(\tilde{\xi}^{(4)})}\big)\big(\Phi_2^{(2)}(\tilde{t}_2^{(2)}, \tilde{\xi}^{(4)}) \tilde{S}_2^{(2)} \times\big(\mathbb{I}-\frac{ \mathcal{X}_1(\tilde{\xi}^{(3)})n^T}{n^T \mathcal{X}_1(\tilde{\xi}^{(3)})}\big)^{-1}D_\epsilon \tilde{\Pi}^{(2)}_1(\tilde{\xi}^{(3)}; 0)+D_\epsilon \Psi_2^{(2)}(\tilde{t}_2^{(2)}, \tilde{\xi}^{(4)})\big).$

Proof. Consider the partial derivatives of $\tilde{\Pi}^{(k)}_i(\xi; \epsilon)$ with respect to $\epsilon$ at the point $(\tau_i^{(k)}, \xi^{(l)}; \epsilon) = (\tilde{t}_i{(k)}, \tilde{\xi}^{(l)}, 0)$,

To complete the computation of the above sub-maps, we need to calculate the terms $D_\epsilon \tau_i^{(k)}(\tilde{\xi}^{(l)}; 0)$ and $D_\epsilon \Psi_i^{(k)}(\tilde{t}_i^{(k)}, \tilde{\xi}^{(l)})$. The derivative of (3.3) for $\epsilon$ is given as:

Note that the recurrence Eq (4.9) describes initial value problems such as $\xi_{\epsilon}(0) = 0$. The general solution (3.5) is therefore equivalent to the following recurrence:

In addition, the partial derivatives of $D_\epsilon \tau_i^{(k)}(\tilde{\xi}^{(l)}; 0)$ are determined by solving the equations $n^T D_\epsilon \tilde{\Pi}^{(k)}_i(\tilde{\xi}^{(l)}; 0) = 0$. So, we get the following at the intersection point $\xi^{(l)}\in \tilde{\Sigma}^c$:

Through substituting (4.10) and (4.11) into (4.8), we get the equalities (a)–(d). □

The following theorem presents the derivation of the Melnikov function associated with perturbation limit cycles and PDB of a differential system (3.3).

Theorem 3. Assume that $\xi\in \tilde{\Sigma}^c_-$ and $\mathcal{\tilde{P}}^{(i)}: \tilde{\Sigma}^c_-\rightarrow \tilde{\Sigma}^c_-$, such that $\mathcal{\tilde{P}}^{(1)}(0;0) = 0$ and $\mathcal{\tilde{P}}^{(2)}(\tilde{\xi}; 0) = \tilde{\xi}$. Then,

Proof. Differentiating (4.7) with respect to $\epsilon$ at $(\tau, \xi, \epsilon) = (\tilde{t}, \tilde{\xi}, 0))$, we obtain

Using the results of Lemma 2 and the value of the vector function $D_\xi \tau_2^{(2)}(\xi; 0)$ from (4.11), the proof of the statement of this theorem is straightforward. □

It should be noted that a similar result for the expansion of the Melnikov function $D_\epsilon \mathcal{\tilde{P}}^{(1)}$ that detects the existence of limit cycles is obtained in ^[9,19]. We present the following corollary, which is based on the concepts of orthogonal projection and related Melnikov function properties discussed in ^[19].

Corollary 2. Assume that $\tilde{\xi}\in \tilde{\Sigma}^c_-$ and $e_1^T \tilde{\xi} = 0$. Then, the orthogonal projection along one-dimensional space $[ \mathcal{X}_1(\tilde{\xi})]$ is given as $\psi_{ \mathcal{X}_1} V = \frac{\langle V, \mathcal{X}_1(\tilde{\xi})\rangle \mathcal{X}_1(\tilde{\xi})}{\langle \mathcal{X}_1(\tilde{\xi}), \mathcal{X}_1(\tilde{\xi})\rangle}$, and the orthogonal projection complement is given by $\psi_{ \mathcal{X}_1^\bot} V = \mathbb{I}- \frac{\langle V, \mathcal{X}_1(\tilde{\xi})\rangle \mathcal{X}_1(\tilde{\xi})}{\langle \mathcal{X}_1(\tilde{\xi}), \mathcal{X}_1(\tilde{\xi})\rangle}$. Therefore, the linearized Poincaré map of the original system is given by

In the following section, we present a novel discontinuous systems that highlights how the nonsmooth perturbation is responsible for the period doubling orbit that emerges from an existing flat periodic orbit (or near an unperturbed one).

Because it is difficult to obtain an explicit analytical expression of the generalized Poincaré map, the theoretical results are numerically validated using path-following techniques for discontinuous dynamical systems. There are numerically sound techniques for approximating the solutions of discontinuous dynamical systems. For examples, see ^[32,33,34].

5.   Examples

Example 1. We consider a simple discontinuous system with a switching plane $\tilde{\Sigma} = \{\xi\in\mathbb{R}^{3} \mid \; e_1^T \xi = 0\}$, which is written as:

where $A_i$ are $3\times 3$ matrices defined as:

such that $\mathcal{W}_i(0) = 0$, and the nonlinear perturbations $\mathcal{W}_i(\xi) = \circ (\parallel\xi\parallel)$ are $C^k$-maps $k\geq 1$ defined on the whole phase space $\mathbb{R}^3$.

We note that the origin of coordinates is still a unique equilibrium point for both systems, and the eigenvalues of the linearized systems $\dot{\xi} = A_{i} \xi$ are given as $\lambda_1 = \mu_1, \lambda_{2, 3} = \pm \bar{i} \omega, \bar{i}^2 = -1$. Further, both linear smooth systems $(\epsilon = 0)$ possess an invariant plane $z = 0$ with constant return times $\tau_\pm = \pi$. Therefore, if $\epsilon = 0$, only one family of flat periodic orbits is generated by the eigenvector: $\xi = (0, \bar{y}, 0), \bar{y}\in \mathbb{R}$.

For $\epsilon\neq 0$, $\mu_1 = 0.15$, and $\mu_3 = -0.6$, the perturbed system (5.1) has a period doubling orbit that emerges from an existing flat periodic orbit with approximately twice the period of the flat orbit. Here, we have Figure 2 with $\omega = \mu_2 = 1$ and $\epsilon = 0.15$, Figure 3 with $\omega = \mu_2 = 1$ and $\epsilon = 0.2$, and Figure 4 with $\omega = \mu_2 = -1$ $\epsilon = 0.5$.

The following example illustrates the emergence of a perturbed period doubling orbit near an unperturbed one.

Example 2. Consider the discontinuous system

where

$a > 0$, and the phase space is divided into two domains separated by a hyperplane $\tilde{\Sigma} = \{\xi\in\mathbb{R}^{3} \mid\; \hbar(\xi; \epsilon) = x = 0\}$. Further, the perturbed system (5.2) has only the crossing mode and cannot exhibit a sliding mode due to $\tilde{\Sigma}{^{s}} = \emptyset$. By numerical computation, it is straightforward to check that the unperturbed system of (5.2) has a PDB at $a = 1.7$, see Figure 5 at $\epsilon = 0$. By setting the perturbation value to $\epsilon = 0.1$, numerical computations show that the perturbed system (5.2) has a period doubling orbit near the unperturbed period doubling orbit, see Figure 5 at $\epsilon = 0.1$.

6.   Conclusions

Practical applications in the real world motivate the development of novel mathematical techniques for analyzing the dynamics of nonsmooth systems. Existing smooth methods, for example, lack a mechanism that explains how switching manifolds generate period-doubling orbits. This paper describes a novel application of the Melnikov method for discontinuous differential systems separated by a hypersurface, with the goal of determining the effect of a nonsmooth perturbation on the period doubling bifurcation of an unperturbed system. The generalized Poincaré map for perturbed systems undergoing nonsmooth transitions has been derived, assuming transversal intersections with discontinuity-switching boundaries and the unperturbed system having a period-doubling orbit. The Melnikov function is designed for systems that have discontinuities and acts as an extension of the regular Melnikov function used in smooth scenarios. Two examples are provided to numerically validate the theoretical results, showing how non-smooth perturbations can cause a period doubling orbit to emerge from an existing flat periodic orbit, as well as how a new period doubling orbit emerges close to the original.

Author contributions

Hany A. Hosham: Conceptualization, formal analysis, investigation, software, writing-original draft; Alaa A. Alzulaibani: Data curation, funding acquisition, project administration, writing-review & editing; Tarek Sellami: Supervision, investigation, writing-review & editing; Khaled Sioud: Methodology, supervision; Thoraya N. Alharthi: Visualisation, formal analysis, methodology, writing-original draft. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number 445-9-734.

Conflict of interest

The authors declare no conflict of interest in this paper.