# Robustness analysis of externally driven damped solitons in the presence of uncertainties and disorders

Arthur Barbosa, Najib Kacem, Noureddine Bouhaddi

► To cite this version:

Arthur Barbosa, Najib Kacem, Noureddine Bouhaddi. Robustness analysis of externally driven damped solitons in the presence of uncertainties and disorders. Physica D: Nonlinear Phenomena, 2025, 476, pp.134612. 10.1016/j.physd.2025.134612. hal-05221435

HAL Id: hal-05221435

https://hal.science/hal-05221435v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Highlights

## Robustness analysis of externally driven damped solitons in the presence of uncertainties and disorders

A. Barbosa, N. Kacem, N. Bouhaddi

* Robustness of stationary solitons under frequency and parameter uncertainties

* Natural frequency uncertainties as the main factor affecting soliton stability

* Lyapunov exponents for analyzing soliton nucleation sensitivity to initial conditions
Robustness analysis of externally driven damped
solitons in the presence of uncertainties and disorders

A. Barbosa<sup>a,*</sup>, N. Kacem<sup>a</sup>, N. Bouhaddi<sup>a</sup>

<sup>a</sup>University of Franche-Comté, CNRS, FEMTO-ST Institute, Department of Applied Mechanics, 26 Rue de l'Épitaphe, Besançon, 25000, Franche-Comté, France

Abstract

This paper investigates the sensitivity of localized vibrations phenomena in externally driven Duffing oscillator chains. Such investigation is conducted by generalizing the Nonlinear Schrödinger Equation (NLSE) to accommodate disorder functions in all physical parameters, beyond impurities commonly found in the literature, limited to the natural frequency of the components. Given the absence of closed-form solutions for externally driven damped systems, we employ a numerical method, followed by statistical analysis, to elucidate the effects of parameter uncertainties across the lattice on solitons behavior. Our findings highlight the diverse effects of the physical nature of uncertainties within the mechanical structure, offering insights into possible experimental investigations. Additionally, we illustrate how specific impurities along the chain, capable to nucleate oscillations, mitigate resonant chaotic behaviors, reinforcing soliton stability. The results affirm the feasibility of generating standing waves in nonlinear lattices, emphasizing

<footnote>*Corresponding author
Email addresses: arthur.barbosa@femto-st.fr (A. Barbosa), najib.kacem@femto-st.fr (N. Kacem), noureddine.bouhaddi@femto-st.fr (N. Bouhaddi)</footnote>
their relevance beyond traditional periodic assumptions, where uncertainties
in physical parameters are commonly disregarded.

*Keywords:* Disordered Nonlinear Schrödinger Equation, Duffing lattices, Soliton sensitivity, Largest Lyapunov Exponent, Externally driven damped solitons.

# 1. Introduction

The Nonlinear Schrödinger Equation (NLSE) is a topic of investigation across diverse domains, spanning analyses pertinent to micro-mechanical systems [1], plasma physics [2], biology [3], optics [4], photonics [5], and metamaterials [6], thermal mediums [7], among others. Concurrently, research into the dynamics of arrays of nonlinear oscillators has demonstrated significant relevance in numerous applications [8, 9, 10, 11, 12, 13], and akin to other fields of study, the NLSE assumes a pivotal role in elucidating such systems.

One main distinction between nonlinear oscillator chains and their linear counterparts lies in their capability to manifest Intrinsic Localized Modes (ILMs) of vibration, even while adhering to the periodicity assumption. However, in linear systems, this phenomenon is commonly observable upon the introduction of impurities, as evidenced in [14, 15, 16]. Indeed, the association between ILMs and standing solitons (which remain stationary in space) is established in the literature [1, 17, 18, 19, 20, 21, 22, 23], contributing to the understanding of energy localization mechanisms in periodic structures [24]. Despite the descriptions of breathing oscillations, interest in the dynamics of oscillations in nonlinear lattices continues to grow.
[25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36].

Even though the non-disorder hypothesis (absence of impurities and uncertainties in the system) holding particular significance in the context of soliton stabilization [20, 37] and propagation [29, 38], the introduction of perturbations in the spatial domain of the NLSE presents a diverse spectrum of opportunities for novel research endeavors. For instance, Villarroel et al. [39] examined the propagation of optical beams in fibers containing random impurities, while Ezzoug et al. [40] utilized impurities to enhance the stability of soliton solutions. Chen et al. [41] investigated soliton collisions within a fiber Bragg grating, elucidating their behavior under such conditions. Additionally, Brazhnyi et al. [42] explored interactions between moving discrete solitons and Gaussian defects, shedding light on their effects during propagation. Parasuraman [43] reported on the instability of localized modes in a chain of Heisenberg ferromagnetic spins with impurities, contributing to the understanding of impurity-induced phenomena. Ribama et al. [44] documented interactions between breathers and impurities, as well as modulational instability within a two-dimensional Klein-Gordon quantum chain. In the realm of biology, Gninzanlong et al. [45] investigated the dynamics of DNA subjected to damping and periodic excitation.

Regardless of the parallels among various physical domains where NLSE appears, mechanical solitons add an additional layer to the problem by imposing energy dissipations that constrain the emergence of localized oscillations [46]. Despite the presence of damping, analytical solutions and well-defined stability diagrams are known for stationary nonlinear waves in parametrically excited scenarios [1, 47, 48, 49]. As highlighted in a previous study
[50], the nature of excitation significantly influences soliton dynamics, and to date, despite the existence of numerical solutions and stability diagrams [51], no analytical solutions have been reported for stationary damped solitons subjected to external excitation. Regarding discrete NLSE, besides periodic lattices, quasiperiodic systems incorporating localized impurities have been demonstrated to effectively generate multi-solitons, reduce the basin of attraction for non-soliton solutions, and broaden the stability diagram for standing solitons, evidenced in [52]. In fact, the impurities reported in the literature are typically related to the natural frequency of oscillators. Although useful for uncertainty analysis, for a more realistic assessment of what could be encountered experimentally, other physical parameters should also be perturbed. In [53], the authors introduced random impurities to non-linearity throughout the NLSE spatial domain. While the article does not specify the physical origin of these perturbations, the work lays the computational foundations for studying potential experimental uncertainties in mechanical field.

Given the context of impurities’ significance in non-dispersive wave formation, this study seeks to deepen our understanding of standing solitons. It does so by statistically analyzing how disorder affects the generation of damped standing solitons in nonlinear mechanical lattices, focusing solely on external excitation. As a result of this study, we present a hierarchization of the physical nature of periodicity disruption in oscillator chains, demonstrating the sensitivity of the phenomenon to the physical parameters. It is also demonstrated how breaks in the periodicity of the natural frequency of the components can expand the basin of attraction of the standing waves,
confirming the indifference of the phenomenon to initial conditions of the integration. We illustrate the robustness of standing solitons under a level of uncertainty in physical parameters, providing designers a guide for applications of the NLSE. The rest of the paper is organized as follows: Section 2 introduces the mathematical foundations, wherein the derivation of the Disordered Damped Nonlinear Schrödinger Equation is conducted. Section 3 presents a sensitivity analysis and a statistical study on the uncertainty of physical parameters. In Section 4, a method based on Lyapunov Exponents is used to draw conclusions about potential chaotic dynamics. Lastly, Section 5 offers the conclusion and outlines possible future directions. Regardless of the context being within a mechanical scope (coupled oscillator chains), due to the general nature of the NLSE, the conclusions presented can also be extrapolated to other physical domains [39, 43, 54].

## 2. Ac-driven NLSE in Non-periodic Oscillator Chains

The mathematical treatment employed in this article is extensively documented in the literature. Among the pioneering works that investigated soliton equations in the presence of impurities on mechanical contexts, we can mention [55], where the authors formulated mathematical expressions that identified impurities in mass, damping, parametric excitation, and coupling stiffness, thereby breaking the assumption of periodicity.

A chain of Duffing oscillators can be mathematically characterized by the interplay of forces among the *N* elements constituting the system (see Figure 1). Initially, a heterogeneity in physical parameters across each oscillator is presumed, acknowledging the potential non-periodic behavior of the system,
<figure><img src="image_2.png" /><figcaption>Figure 1: Equivalent mechanical system describing coupled Duffing Oscillators.</figcaption></figure>

from which the following relationship is deduced:

$$m_n \ddot{\eta}_n + b_n \dot{\eta}_n + k_{ln} \eta_n - k_{nl_n} \eta_n^3 + k_{c_{n+1}} (\eta_n - \eta_{n+1}) + k_{c_n} (\eta_n - \eta_{n-1}) \\ + m_n A(0.5e^{i\omega t} + 0.5e^{-i\omega t}) = 0, \quad (1)$$

where $m_n, b_n, k_{ln}, k_{nl_n}, k_{c_n}, A, \omega$, and $\eta$ represent the mass, damping coefficient, linear stiffness, nonlinear stiffness, coupling stiffness, base acceleration amplitude, excitation frequency, and displacement of each oscillator with index $n$.

Without loss of generality, the physical parameters can be expressed as the sum of a mean value and a deviation from this mean. From the perspective of uncertainty analysis, the fluctuation around the mean should be accompanied by a function that incorporates a probability distribution describing the occurrence of these uncertainties. It is worth noting that, besides a probabilistic approach, variations around the mean could serve as parameters selected by the designer. The notation adopted in this paper will be
formulated according to the following relationships:

$$k_{l_n} = \check{k}_l + \Delta k_l p_{k_l}(n), \quad b_n = \check{b} + \Delta b p_b(n), \quad k_{c_n} = \check{k}_c + \Delta k_c p_{k_c}(n),$$

$$A = \check{A} + \Delta A p_A, \quad m_n = \check{m} + \Delta m p_m(n), \quad k_{nl_n} = \check{k}_{nl} + \Delta k_{nl} p_{k_{nl}}(n), \quad (2)$$

where, for any arbitrary variable $\rho$, a mean value $\check{\rho}$ and a deviation $\Delta\rho$ are associated. As mentioned, a probabilistic distribution $p_\rho$ follows these variations, thus generalizing the mathematical modeling to the case where uncertainties are not controlled.

Taking into account the relationships established by (2), the dynamic equation of the nonlinear chain can be reformulated as follow:

$$\left(1 + \frac{\Delta m p_m(n)}{\check{m}}\right) \ddot{\eta}_n + 2\check{\omega}_0 (\check{\gamma} + \Delta\gamma p_b(n)) \dot{\eta}_n + (\check{\omega}_0^2 \pm \Phi p_{k_l}(n)) \eta_n \\ + \left(-\frac{\check{k}_{nl}}{\check{m}} - \frac{\Delta k_{nl} p_{k_{nl}}(n)}{\check{m}}\right) \eta_n^3 + \left(\frac{\check{k}_c}{\check{m}} + \frac{\Delta k_c p_{k_c}(n)}{\check{m}}\right) (2\eta_n - \eta_{n+1} - \eta_{n-1}) \\ \left(\frac{m_n}{\check{m}}\right) (\check{A} + p_A \Delta A) (0.5e^{itw} + 0.5e^{-itw}) = 0, \quad (3)$$

where:

$$\check{\omega}_0^2 = \frac{\check{k}_l}{\check{m}}, \quad 2\check{\gamma}\check{\omega}_0 = \frac{\check{b}}{\check{m}}, \quad 2\check{\omega}_0\Delta\gamma = \frac{\Delta b}{\check{m}}, \quad \Phi = \frac{\Delta k_l}{\check{m}}, \quad \gamma_n = \check{\gamma} + \Delta\gamma p_\gamma(n). \quad (4)$$

It is important to note that while changes in mass can impact the natural frequency of the oscillators, the symbol $\Phi$ is specifically used to denote variations in linear stiffness. This distinction is made because Equation (1) has been normalized based on the average/nominal mass of the components. For
algebraic simplicity, variations in the nominal mass of the oscillators are neglected ($m_n = \tilde{m}, \forall n$), which could be incorporated into variations of other physical parameters by normalization via $m_n$.

For solving Equation (3), the method of multiple scales is employed, introducing a perturbation parameter $\epsilon$. Considering an arbitrary variable $\rho$, the notation $\rho_\epsilon$ denotes the rescaling of this variable in accordance with $\epsilon$ by the relation: $\rho = \epsilon\rho_\epsilon$. Consequently, Equation (3) can be reformulated alternatively as:

$$ 
\begin{aligned}
\ddot{\eta}_n + 2\dot{\omega}_0 (\dot{\gamma}_\epsilon \epsilon + \Delta\gamma_\epsilon \epsilon p_b(n)) \dot{\eta}_n + (\dot{\omega}_0^2 + \Phi_\epsilon \epsilon p_{kl}(n)) \eta_n \\
+ \left( -\frac{\tilde{k}_{nl}\epsilon}{\tilde{m}} - \frac{\Delta k_{nl}\epsilon p_{knl}(n)}{\tilde{m}} \right) \eta_n^3 \\
+ \left( \frac{\tilde{k}_{c\epsilon}\epsilon}{\tilde{m}} + \frac{\Delta k_{c\epsilon}\epsilon p_{k\epsilon}(n)}{\tilde{m}} \right) (2\eta_n - \eta_{n+1} - \eta_{n-1}) \\
+ (\tilde{A}_\epsilon \epsilon + p_A \Delta A_\epsilon \epsilon) (0.5e^{itw} + 0.5e^{-itw}) = 0.
\end{aligned}
\quad (5) $$

The solution of Equation (5) can be described as the addition of terms with magnitudes that are sorted according to powers of $\epsilon$:

$$ \eta_n(t) = \left( -\psi_n(t\epsilon)e^{-i\omega} - e^{i\omega} \overline{\psi_n(t\epsilon)} \right) + O(\epsilon). \quad (6) $$

While unconventional, the function $\psi_n$ is accompanied by a negative sign due to mathematical simplifications applied in the algebraic approach.

Assuming an excitation frequency close to the expected natural frequency of each oscillator, $w = \dot{\omega}_0(1 - \epsilon)$, the algebraic manipulation of Equation
(5), based on the solution given in expression (6), results in the emergence of secular terms that must be canceled to ensure convergence. Such cancellation can be achieved through the following Equation:

$$ 
\begin{aligned}
& \left( \frac{\check{k}_{c\epsilon}}{\check{k}_l} + \frac{\Delta k_{c\epsilon} p_{k_c}(n)}{\check{k}_l} \right) (2\psi_n - \psi_{n+1} - \psi_{n-1}) \\
& + \left( \frac{3\check{k}_{nl\epsilon}}{2\check{k}_l} + \frac{3\Delta k_{nl\epsilon} p_{k_{nl}}(n)}{2\check{k}_l} \right) \psi_n^2(T) \overline{\psi_n(T)} \\
& + \left( \check{\gamma}_\epsilon i + \Delta\gamma_\epsilon ip_b(n) - \frac{\Phi_\epsilon p_{k_l}(n)}{2\check{\omega}_0^2} - 1 \right) \psi_n(T) + \frac{0.25\check{m}}{\check{k}_l} (\check{A}_\epsilon + p_A \Delta A_\epsilon) \\
& = -\frac{i}{\check{\omega}_0} \frac{d}{dT} \psi_n(T),
\end{aligned}
\quad (7) $$

where the time variable $T$ arises from the time scale $T = \epsilon t$.

Equation (7) can be identified as a representation of the discretized non-linear Schrödinger equation with cubic nonlinearity. In order to broaden the scope of the analysis, an adimensionalization is performed using the following functions:

$$ \Psi_n(\tau) = \frac{\sqrt{3}\sqrt{\check{k}_{nl\epsilon}}}{2\sqrt{\check{k}_l}} \psi_n(T), \quad \tau = T\check{\omega}_0. \quad (8) $$

Thus Equation (7) can be rewritten in the form:
$$ \frac{\frac{\Delta s^2}{2\check{k}_l} (\check{k}_{c\epsilon} + \Delta k_{c\epsilon} p_{k_c}(n)) (-2\Psi_n(\tau) + \Psi_{n+1}(\tau) + \Psi_{n-1}(\tau))}{\Delta s^2} \\ + \left(2 + \frac{2\Delta k_{nl\epsilon} p_{k_{nl}}(n)}{\check{k}_{nl\epsilon}}\right) \Psi_n(\tau) |\Psi_n(\tau)|^2 \\ + \left(\check{\gamma}_\epsilon i + \Delta\gamma_\epsilon ip_b(n) - \frac{\Phi\epsilon p_{k_l}(n)}{2\check{\omega}_0^2} - 1\right) \Psi_n(\tau) \\ + \frac{0.125\sqrt{3}\check{m}\sqrt{\check{k}_{nl\epsilon}} (\check{A}_\epsilon + p_A\Delta A_\epsilon)}{\check{k}_l^{\frac{3}{2}}} = -i \frac{d}{d\tau} \Psi_n(\tau), \quad (9) $$

where, in parallel with the process of nondimensionalization, a spatial variable, $\Delta s$, is introduced. This variable is related to the length, $L$, of the chain through the expression $\Delta s = L/(N - 1)$, representing the distance between components.

## 2.1. Considerations on the Continuity Hypothesis

As previously discussed, the propagation of soliton waves in chains of nonlinear oscillators is extensively documented in the literature, and it is typically described through the Nonlinear Schrödinger Equation. Modeling this phenomenon can be approached in two distinct ways: either by approximating the system as a continuous entity ($\Psi_n(\tau) \to \Psi(s, \tau)$), or by investigating the dynamic behavior in a discrete domain. The choice between these two strategies depends on the evaluation of the subsequent approximation:

$$ \frac{(-2\Psi_n + \Psi_{n+1} + \Psi_{n-1})}{\Delta s^2} \to \frac{\partial^2}{\partial s^2} \Psi(s, \tau). \quad (10) $$

From a numerical computing perspective, if the second-order spatial
derivative of $\Psi(s, \tau)$ exists, we also could assume that these derivatives could be estimated using the finite difference method, leading to the approximation presented in Equation (10) and, through this technique, the continuity hypothesis could be used in perfectly periodic lattices.

For singular impurities, the assessment of the continuity assumption has been conducted from both a mechanical standpoint [52, 23, 8] and other physical domains [55, 54, 43]. Typically, an approximation is made for the first derivative, both to the left and right of the impurity position, but with a discontinuity at the point where the impurity occurs. To exemplify this issue, let's explore the following Schrödinger Equation, incorporating an impurity situated at $s = 0$ by the Dirac Delta function $\delta_d$:

$$i \frac{\partial}{\partial t} \vartheta + \frac{\partial^2}{\partial s^2} \vartheta + 2|\vartheta|^2 \vartheta - \vartheta + 2q\delta_d(s)\vartheta = -i\gamma\vartheta. \quad (11)$$

In accordance with results reported in the literature [56, 57], it is coherent to assume that:

$$2q\vartheta(0) + \vartheta_X|_{-0}^{+0} = 0, \quad (12)$$

providing the mathematical support to the application of the continuous modeling strategy while the periodicity is not respected.

The temporal nondimensionalization was performed based on the mean natural frequency of the oscillators (Equations (8)). However, for the spatial domain, a new variable needs to be formulated, as the allocation of the structure's length in the dimensionless domain of the NLSE is also open to an infinite range of possibilities. With this aim in consideration, the variable $\delta$ is built, where the following expression is valid for the dimensionless length
$\hat{L}$:

$$ \hat{L} = \frac{L}{\delta}. \qquad (13) $$

In accordance with the findings in [51], the NLSE exhibits energy localization near $\hat{s} = 0$, where $\hat{s}$ denotes the dimensionless spatial coordinate ($\hat{s} = \frac{s}{\delta}$). This knowledge is crucial for determining the appropriate value for $\delta$ because as $\hat{L}$ values increase, there are fewer oscillators in the surroundings of the region with high oscillation amplitudes ($\hat{s} = 0$). For a better understanding of the choice of $\delta$, a detailed study is given in [58].

Without compromising the generality of the model, and considering spatial normalization and continuity notation, Equation (9) is reformulated in the following format:

$$ \begin{aligned} & \zeta(\hat{s}) \frac{\partial^2}{\partial \hat{s}^2} \Psi(\hat{s}, \tau) + i \frac{\partial}{\partial \tau} \Psi(\hat{s}, \tau) \\ & + hQ(\hat{s}) + \Psi(\hat{s}, \tau)F(\hat{s}) + 2\Psi(\hat{s}, \tau)P(\hat{s})|\Psi(\hat{s}, \tau)|^2 = 0, \end{aligned} \qquad (14) $$

where:

$$ \begin{aligned} Q(\hat{s}) &= 1 + \frac{p_A \Delta A_{\epsilon}}{\check{A}_{\epsilon}}, \quad \zeta(\hat{s}) = 1 + \frac{\Delta k_{c\epsilon} p_{k_c}(\hat{s})}{\check{k}_{c\epsilon}}, \\ P(\hat{s}) &= 1 + \frac{\Delta k_{nl\epsilon} p_{k_{nl}}(\hat{s})}{\check{k}_{nl\epsilon}}, \end{aligned} \qquad (15) $$

$$ F(\hat{s}) = -\frac{\Phi_{\epsilon} p_{k_l}(\hat{s})}{2\check{\omega}_0^2} - 1 + i\zeta \left( 1 + \frac{\Delta\gamma_{\epsilon} p_b(\hat{s})}{\check{\gamma}_{\epsilon}} \right), $$

and
$$ \zeta = \check{\gamma}_{\epsilon}, \quad \Delta s = \Delta \hat{s} \delta, \quad \check{k}_{c\epsilon} = \frac{2\check{k}_l \delta^2}{\Delta s^2}, \quad h = \frac{0.125\sqrt{3}\check{m}\check{A}_{\epsilon}\sqrt{\check{k}_{nl\epsilon}}}{\check{k}_l^{\frac{3}{2}}}. \quad (16) $$

We should note that an extrapolation of the continuity notation is also adopted for the $p(n)$ functions:

$$ p(n) \rightarrow p(s) \rightarrow p(\hat{s}). \quad (17) $$

The term “continuity notation” instead of “continuity assumption” stems from the fact that the extrapolation of the relationship (10) extends beyond the scope of the performed investigations, and all the analyses presented are grounded in numerical strategies. For the sake of notation economy, the disorder functions will be rewritten according to the following convention:

$$ Q(s) = 1 + p_{q_1}(s)q_1, \quad \zeta(s) = 1 + p_{q_2}(s)q_2, $$

$$ P(s) = 1 + p_{q_3}(s)q_3, \quad (18) $$

$$ F(s) = p_{q_4}(s)q_4 - 1 + i\zeta(1 + p_{q_5}(s)q_5), $$

where the functions $p$ determine the location and probability of occurrence, and the variables $q$ represent the magnitude of periodicity disruption. It should be noted that the values of $q_i$ for $i = 1, 2, 3$, and $5$ represent the proportional variations of their associated physical parameters. The same reasoning cannot be directly applied to $q_4$, since, despite also being a function of the periodicity breaking of $k_l$, the perturbation factor $\epsilon$ must also be taken into account. The physical interpretation of $q_4$ is related to percentage variations of the natural frequency of the oscillators but in light of the excitation frequency: $q_4 = \frac{\Delta k_l}{2(1 - \frac{w}{\check{\omega}_0})\check{k}_l}$
Equation (14) represents the generalization of the damped and externally driven NLSE and was first reported in [59], in an experimental study on energy nucleation in a chain of oscillators. The nomenclature employed the term disorder, thus it is referred to as the Disordered Damped Nonlinear Schrödinger Equation (DDNLSE). The specific case without the presence of impurities, perfectly periodic, can be obtained by examining the functions accompanying $\Psi$. The mathematical treatment remains unaffected by the physical architecture of the oscillators. This consideration indicates the broad applicability of Equation (14), as it can be employed in pendulum [20], linear [23], or circular [21] architectural configurations.

### 3. Study of NLSE sensitivity to uncertainties in physical parameters

#### 3.1. Numerical simulation strategy

In pursuit of achieving stationary waves, localized in space for any value of $\tau$, it is immediate to assume that $\frac{\partial}{\partial\tau}\Psi(\hat{s}, \tau) = 0$. The particular case without energy dissipation ($\gamma = 0$), as initially proposed by [51], accommodates the following solutions:

$$ \Psi_{\pm}(\hat{s}) = \Psi_0 \left( 1 + \frac{2 \sinh^2 \alpha}{1 \pm \cosh \alpha \cosh(C\hat{s})} \right), \quad (19) $$

where

$$ h = \frac{\sqrt{2} \cosh^2 \alpha}{(1 + 2 \cosh^2 \alpha)^{3/2}}, \quad \Psi_0 = \frac{1}{\sqrt{2(1 + 2 \cosh^2 \alpha)}}, \quad C = 2\Psi_0 \sinh \alpha. \quad (20) $$
Although closed-form solutions for damped situations are not reported, for
absence of lattice disorder, soliton stability is known (see Fig. 2). As eluci-
dated in [51], only Dark Solitons Ψ_ are temporally stable when damping is
imposed, therefore these will be the objective of study in this article.

<figure><img src="image_3.png" /><figcaption>Figure 2: Stability diagram for Equation (14) in the absence of disorders, based on [51].

Curve I, in blue ($0.3 \le \zeta \le 0.5$): $h - \left( \frac{1}{3} \left( \frac{1}{3} \left( -\zeta^2 + \frac{1}{3} \right)^3 \right)^{0.5} + \frac{1}{3} \left( \zeta^2 + \frac{1}{9} \right) \right)^{0.5} = 0$.

Curve II, in orange ($0.5 \le \zeta \le 0.66$): $h - \left( \zeta^3 - \zeta^2 + \frac{\zeta}{2} \right)^{0.5} = 0$. Curve III (Hopf bifurcation), in purple ($0 \le \zeta \le 0.3$): $h - (1.8728\zeta^3 + 1.6822\zeta^2 + 0.0997\zeta + 0.0777) = 0$.

Curve IV, in black ($0 \le \zeta \le 0.66$): $h - \frac{2\zeta}{\pi} = 0$. Curve III (Hopf bifurcation) is responsible for separating the region of stable parameter values (S) from the unstable parameter region (U).</figcaption></figure>

Strategies have been proposed over the years for the numerical investiga-
tion of Equation (14). Some approaches have proven effective in addressing
this issue, as the one reported by [51], employing the Newton’s method and
the numerical continuation to seek potential parameters that would ensure
stationary solutions. Another possible approach was employed by [60], where
the use of the pseudo-spectral method, for low damping values also proved
to be efficient. In [53], the authors employed a method based on the Crank-Nicolson algorithm to address the analyses of undamped yet chaotic cases. Recently, in [61], the use of a hybrid variational method has also proven to be effective as a computational strategy.

With the aim of optimizing computational time, a neural network was created in our previous work, significantly reducing the computational time required for obtaining numerical solutions [50]. The data used can be freely accessed in the provided reference and the solution obtained by the metamodel operates with a discretization of 600 points. In this work, the discretization will comprise of 601 points. This adjustment is accomplished by adding an extra discretization at the origin ($\hat{s} = 0$) and setting its amplitude as the arithmetic mean of the amplitudes of the adjacent points. This discretization ensures that only 1 point will have the highest value for $|\Psi(\hat{s}, \tau)|$.

The numerical process follows the subsequent sequence of steps:

1. Choice of pairs $h \times \zeta$, defining an undisturbed equation ($q_i = 0$, for $i = 1, 2, 3, 4, 5$);

2. Obtaining a solution by the metamodel;

3. Calibration of the metamodel solution error through the Newton method [51] and addition of one extra point in the discretization;

4. Numerical integration of Equation (14), starting from the previously calibrated steady-state solution with no disorder;

5. Identification of the integrated solution under steady-state conditions;

6. Categorization of the solution with respect to localization.

The computational parameters used are presented in Table 1.
The classification of the steady-state dynamic behavior of Equation (14) could be approached in more than one way. In the proposed analysis, this classification is performed by comparing the integrated solution with the metamodel-calibrated one, thus verifying the temporal stability of the standing wave. To ensure that solutions do not disperse or diverge to infinity, the oscillation amplitude will be classified as localized at $\hat{s} = \hat{s}_i$ if $0.25|\Psi(\hat{s} = 0, \tau = 0)| \le |\Psi(\hat{s}_i, \tau = 30)| \le 4|\Psi(\hat{s} = 0, \tau = 0)|$. In parallel with the criterion of intensity of localization, a spatial criterion is also adopted. Once the peaks of $|\Psi(\hat{s}, \tau)|$ are identified (points where the function $|\Psi(\hat{s}, \tau)|$ transitions from increasing to decreasing behavior), the classification is considered a soliton if the ratio between the value of the second-highest peak and the highest peak is less than or equal to 0.5; otherwise, the solution is categorized as either possessing a multi-soliton or a dispersive behavior.

<table><caption>Table 1: Simulation parameters used in Numerical Integration.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Newton method calibration error</td><td>10<sup>-5</sup></td></tr><tr><td>Integration domain</td><td>-50 ≤ <i>s</i> ≤ 50</td></tr><tr><td>Number of points</td><td>601</td></tr><tr><td>Spacial domain discretization Δ<i>s</i></td><td>0.17</td></tr><tr><td>Integration method</td><td>DOP853</td></tr><tr><td>Temporal integration domain <i>τ</i></td><td>30</td></tr><tr><td>Temporal domain discretization</td><td>0.05</td></tr><tr><td>Integration tolerance</td><td>10<sup>-4</sup></td></tr></tbody></table>

The core idea behind the classification algorithm revolves around the
<figure><img src="image_4.png" /><figcaption>Figure 3: Graphical comparison between $|\Psi(s, 30)|$ in the presence of impurities described by disorder functions. Simulation parameters: $h = 0.25$ and $\zeta = 0.375$. Assumed uncertainties: $q_i = 0.05$ and $q_j = 1$, where $i = 1, 2, 3, 5$ and $j = 4$. Associated probability distribution: $p_j(s) = p_i(s) = U(-1, 1)$ for all points in $\hat{s}$.</figcaption></figure>

postulate that, despite the presence of disturbances, the soliton behavior should not significantly differ from what would be obtained if no disorder is considered. Fig. 3 provides a visual perception of this premise, drawing a parallel between the ideal, non disordered solution, and the one accounting for potential uncertainties. The spatial domain under consideration extends from $-50 \le \hat{s} \le 50$, so that the boundary conditions ($|\Psi(-50, \tau)| = |\Psi(50, \tau)| = 0$) had minimal impact on the numerical integration.

## 3.2. Uncertainties unrelated to the natural frequency of oscillators

The physical origin of uncertainty related to nonlinear and coupling stiffness typically arises from the experimental difficulty in handling them, whether in controlling the distance of magnetic interactions [62] or in managing the geometry of the structure [63]. In certain applications, nonlinearity can be artificially induced outside of the oscillator, as reported by Alfahmi et al. [64], where a synthetic impedance-based Duffing shunt with hardening
and softening nonlinearities, was studied. On the other hand, the physical
origin of damping uncertainties may emerge due to difficulties in measuring
or modeling it, as is often the case in viscoelastic systems [65, 66].

The study of this of type disorder is conducted by the numerical scanning
of Equation (14), while disregarding other ones ($q_4 = 0$). A statistical anal-
ysis is implemented, and a Uniform probability distribution is adopted for
the functions $p, U[-1, 1]$. The choice of the Uniform distribution, although
arbitrary, aims to generalize $p$. In parallel with the choice of the probability
density function, the number of samples is also crucial for the robustness<fnref n="1" />
of the analysis. In this work, 1000 randomly selected configurations for dis-
orders are applied on the numerical integration at randomly points near the
boundaries of the regions delimiting stable parameters for the periodic NLSE
(Figures 4, 5 and 6).

Intuitively, we can hypothesize that increasing of disorders would gradu-
ally reduce the U area displayed in Fig. 2. To investigate this premise, three
progressively decreasing areas are examined, and parameters on the borders
of these new areas are explored while considering uncertainties throughout
the lattice. The reduction of the areas is determined by an arbitrary param-
eter $r$, where the new curves are adjusted using the relationships:

<footnote n="1">¹We corroborate the notion that the level of robustness is linked to the level of uncertainties and disorders that the model is able to tolerate without changing the decision in term of design: the more the model can tolerate uncertainties and desorders, the more it will be considered robust.</footnote>
$$ (\text{new } h) = (\text{old } h) \pm r, \qquad (21) \\ (\text{new } \zeta) = (\text{old } \zeta) \pm r. $$

The sign associated with $r$ is determined based on the curves shown in Fig. 2 with the objective of reduce the area enclosed by these curves.

Table 2 displays the selected points and the resulting outcomes of the numerical investigation. According to the simulations, at a 5% uncertainty level, as we move away from the boundaries, we can generate stationary solitons with greater probability. As anticipated in the case without disorder, high damping values notably impede the achievement of stable solutions and, for a choice at the midle of the diagram ($h = 0.25$ and $\zeta = 0.3$), 100% of the simulations exhibited soliton behavior.

<table><caption>Table 2: Description of simulation parameters and results obtained in the numerical investigation of disorders not related to the natural frequency of the oscillators.</caption><thead><tr><th>Area</th><th>r</th><th>Results</th><th>q<sub>i</sub> (i = 1, 2, 3, 5)</th></tr></thead><tbody><tr><td>Fig. 4</td><td>8.5 × 10<sup>-3</sup></td><td>Table 3</td><td>5%</td></tr><tr><td>Fig. 5</td><td>10 × 10<sup>-3</sup></td><td>Table 4</td><td>5%</td></tr><tr><td>Fig. 6</td><td>20 × 10<sup>-3</sup></td><td>Table 5</td><td>5%</td></tr></tbody></table>

### 3.3. Uncertainties in all physical parameters

The breakdown of periodicity induced by changes in the natural frequency of oscillators, unlike other types of disorder, has been widely studied in the literature. In linear systems, such breakdowns can be manifested through Anderson Localization, whereas in nonlinear systems, such disruptions can
<figure><img src="image_5.png" /><figcaption>Figure 4: Stability boundaries and parameter pairs used in statistical analyses. Boundaries reduction performed from $r = 8.5 \times 10^{-3}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Stability boundaries and parameter pairs used in statistical analyses. Boundaries reduction performed from $r = 10 \times 10^{-3}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Stability boundaries and parameter pairs used in statistical analyses. Boundaries reduction performed from $r = 20 \times 10^{-3}$.</figcaption></figure>
<table><caption>Table 3: Percentage of soliton results in numerical investigation in Table 2. Associated boundaries: Fig. 4.</caption><thead><tr><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th></tr></thead><tbody><tr><td>p1</td><td>100%</td><td>p7</td><td>100%</td><td>p13</td><td>100%</td><td>p19</td><td>92.3%</td></tr><tr><td>p2</td><td>100%</td><td>p8</td><td>100%</td><td>p14</td><td>100%</td><td>p20</td><td>90.8%</td></tr><tr><td>p3</td><td>100%</td><td>p9</td><td>100%</td><td>p15</td><td>100%</td><td>p21</td><td>91.7%</td></tr><tr><td>p4</td><td>100%</td><td>p10</td><td>99.9%</td><td>p16</td><td>100%</td><td>p22</td><td>90.9%</td></tr><tr><td>p5</td><td>100%</td><td>p11</td><td>99.2%</td><td>p17</td><td>100%</td><td>p23</td><td>88.4%</td></tr><tr><td>p6</td><td>100%</td><td>p12</td><td>100%</td><td>p18</td><td>91.3%</td><td></td><td></td></tr></tbody></table>

<table><caption>Table 4: Percentage of soliton results in numerical investigation in Table 2. Associated boundaries: Fig. 5.</caption><thead><tr><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th></tr></thead><tbody><tr><td>p1</td><td>100%</td><td>p6</td><td>100%</td><td>p11</td><td>100%</td><td>p16</td><td>96.6%</td></tr><tr><td>p2</td><td>100%</td><td>p7</td><td>100%</td><td>p12</td><td>100%</td><td>p17</td><td>98.3%</td></tr><tr><td>p3</td><td>100%</td><td>p8</td><td>100%</td><td>p13</td><td>100%</td><td>p18</td><td>98.6%</td></tr><tr><td>p4</td><td>100%</td><td>p9</td><td>99.9%</td><td>p14</td><td>100%</td><td>p19</td><td>97.5%</td></tr><tr><td>p5</td><td>100%</td><td>p10</td><td>100%</td><td>p15</td><td>100%</td><td>p20</td><td>95.4%</td></tr></tbody></table>

be analyzed from a soliton perspective. In [67], it was demonstrated that the center of mass of a soliton may be affected by Anderson localization, a finding supported by [68]. Studies examining the degree of disorder that Anderson localization can withstand in the presence of nonlinearities are presented in [69] and [70].

Similar to what was shown in [71] (but in another physical domain),
<table><caption>Table 5: Percentage of soliton results in numerical investigation in Table 2. Associated boundaries: Fig. 6.</caption><thead><tr><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th></tr></thead><tbody><tr><td>p1</td><td>100%</td><td>p6</td><td>100%</td><td>p11</td><td>100%</td><td>p16</td><td>100%</td></tr><tr><td>p2</td><td>100%</td><td>p7</td><td>100%</td><td>p12</td><td>100%</td><td>p17</td><td>100%</td></tr><tr><td>p3</td><td>100%</td><td>p8</td><td>100%</td><td>p13</td><td>100%</td><td>p18</td><td>100%</td></tr><tr><td>p4</td><td>100%</td><td>p9</td><td>100%</td><td>p14</td><td>100%</td><td>p19</td><td>100%</td></tr><tr><td>p5</td><td>100%</td><td>p10</td><td>100%</td><td>p15</td><td>100%</td><td>p20</td><td>100%</td></tr></tbody></table>

in [59], we reported simulations showing that the allocation of impurities can tune mechanical lattice solitons, expanding the stability diagram and promoting spontaneous nucleation. The experimental results indicated the effectiveness of the new stability diagram, supporting previously reported findings in the literature regarding the imposition of nonlinearities on Anderson localization. The outcomes allow for the exploration of disorder as a strategy to control nonlinear energy localization in damped and non-periodic systems. In the present study, the analysis does not focus on impurities capable of promoting spontaneous nucleation, but rather in the opposite direction: examining the influence of impurities that do not localize the oscillation energy, thereby in conjunction with disorders of other physical natures, given the existence of works studying isolated influences of $q_4$.

The assumption of the absence of one type of uncertainty in favor of the other, despite being important for understanding the role of each type of disorder, is not feasible in many contexts. Therefore, separating them is not a realistic hypothesis to adopt in a general perspective. In this context, a
<table><caption>Table 6: Description of simulation parameters and results obtained in the numerical investigation of both types of disorders considered simultaneously.</caption><thead><tr><th>Area</th><th>r</th><th>Results</th><th>q<sub>i</sub> (i = 1, 2, 3, 5)</th><th>q<sub>4</sub></th></tr></thead><tbody><tr><td>Fig. 4</td><td>8.5 × 10<sup>-3</sup></td><td>Table 7</td><td>8%</td><td>1% (ε = 0.01)</td></tr><tr><td>Fig. 6</td><td>20 × 10<sup>-3</sup></td><td>Table 8</td><td>5%</td><td>1% (ε = 0.01)</td></tr></tbody></table>

numerical examination regarding the influence of all uncertainties, simultaneously present, is performed. Table 6 shows the simulation results under this condition.

As expected, an increase in the level of disorders significantly decreases the probability of obtaining stationary solitons, particularly with the increase in damping. In fact, besides observing a greater sensitivity of the phenomenon concerning $q_4$, we can also infer a correlation between it and other physical disorders. It is interesting to note that the difference between $q_4$ and $q_{\neq 4}$ can be up to 8 times, clearly indicating a hierarchy among the levels of uncertainty. This is due to the perturbation factor $\epsilon$, which, being small by definition, can considerably increase the value of $q_4$ (see Equation (16)).

The physical implication of the results presented in Tables 7 and 8 provides a design guideline for experimental work. Given that the robustness of the phenomenon is more sensitive to periodicity disruptions in the natural frequency, controlling nonlinearity, damping, and lattice acceleration allows a greater margin for error, whether inherent to the model or due to measurement inaccuracies. A physical perspective on this statement can be obtained through the frequency response of Duffing oscillators shown in Figure
<table><caption>Table 7: Percentage of soliton results in the numerical investigation in Table 6. Associated boundaries: Fig. 4.</caption><thead><tr><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th></tr></thead><tbody><tr><td>p1</td><td>97.8%</td><td>p7</td><td>96.4%</td><td>p13</td><td>95.5%</td><td>p19</td><td>6.6%</td></tr><tr><td>p2</td><td>98%</td><td>p8</td><td>95.7%</td><td>p14</td><td>99.2%</td><td>p20</td><td>6.2%</td></tr><tr><td>p3</td><td>99.4%</td><td>p9</td><td>93.3%</td><td>p15</td><td>100%</td><td>p21</td><td>5.3%</td></tr><tr><td>p4</td><td>99.5%</td><td>p10</td><td>87.7%</td><td>p16</td><td>99.8%</td><td>p22</td><td>6.2%</td></tr><tr><td>p5</td><td>99.4%</td><td>p11</td><td>65.3%</td><td>p17</td><td>79.2%</td><td>p23</td><td>5.3%</td></tr><tr><td>p6</td><td>97.9%</td><td>p12</td><td>76.9%</td><td>p18</td><td>5.1%</td><td></td><td></td></tr></tbody></table>

<table><caption>Table 8: Percentage of soliton results in the numerical investigation in Table 6. Associated boundaries: Fig. 5.</caption><thead><tr><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th><th>Point</th><th>Solitons</th></tr></thead><tbody><tr><td>p1</td><td>94.8%</td><td>p6</td><td>100%</td><td>p11</td><td>83.9%</td><td>p16</td><td>96.5%</td></tr><tr><td>p2</td><td>97.3%</td><td>p7</td><td>100%</td><td>p12</td><td>97%</td><td>p17</td><td>27.3%</td></tr><tr><td>p3</td><td>99.7%</td><td>p8</td><td>100%</td><td>p13</td><td>99.5%</td><td>p18</td><td>33%</td></tr><tr><td>p4</td><td>100%</td><td>p9</td><td>99.7%</td><td>p14</td><td>100%</td><td>p19</td><td>28.8%</td></tr><tr><td>p5</td><td>100%</td><td>p10</td><td>93.3%</td><td>p15</td><td>100%</td><td>p20</td><td>33.4%</td></tr></tbody></table>

7 [72, 73]. Variations in $k_{nl}$, A, or $\zeta$ affect the slope of the frequency response backbone, while changes related to $\omega_0$ additionally shift it. As shown in [58], once the design parameters are fixed, the coupling stiffness imposes a relationship between the physical spatial domain $s$ and the dimensionless domain of the Schrödinger Equation, $\hat{s}$, enabling the design of infinite spatial configurations. This fact provides a physical perspective on why the system is less sensitive to coupling variations compared to its sensitivity to the natural
<figure><img src="image_8.png" /><figcaption>Figure 7: Frequency response curves of Duffing oscillators described by Equation (1). Blue curves are obtained using the parameters $m = 1$, $A = 1$, $\omega_0 = 1$, $\zeta = 0.01$, and $k_{nl} = 10^{-5}$. Parameter changes for generating the other curves are detailed in the individual legends. Units omitted for simplicity of notation.</figcaption></figure>

frequency of the oscillators.

4. Sensitivity of solitons to initial conditions

As pointed out by [1] and [58], obtaining soliton dynamics from initially stationary oscillators is not a trivial task. Both studies suggest an alternative to this challenge by driving the system to an unstable parameter configuration $^2$, and after a certain number of oscillation periods, imposing soliton stability on the system. As expected, this methodology suggests a sensitivity of nonlinear chains to initial conditions.

Alexeeva et al. [52] revealed the role that $q$ plays in parametrically driven chains, raising questions about solitons externally excited, which were answered in [59]. As a synthesis of the two works, we can extract two important conclusions about such impurities: nucleation capacity and expansion of the stability diagram. These findings suggest an enhancement in obtaining sta-

<footnote>$^2$Parameter configuration $h \times \zeta$ above the upper threshold set by the stability diagram.</footnote>
tionary solitons by the allocation of $q_4$ impurities, as the prerequisite initial instability suggested by [1] and [58] would no longer be necessary.

For an understanding of the phenomenon, we should rewrite Equation (14) from the perspective of a single impurity localized at position $s = s_i$:

$$ (i\zeta + q\delta(-\hat{s}_i + \hat{s}) + 2|\Psi(\hat{s}, \tau)|^2 - 1) \Psi(\hat{s}, \tau) + h + 2\Psi(\hat{s}, \tau)|\Psi(\hat{s}, \tau)|^2 + \frac{\partial^2}{\partial \hat{s}^2} \Psi(\hat{s}, \tau) + i \frac{\partial}{\partial \tau} \Psi(\hat{s}, \tau) = 0, \quad (22) $$

where

$$ p_{q_4}(s_i) = \delta(-s_i + s) \quad \text{and} \quad q = q_4. \quad (23) $$

The function $p_{q_4}$ is now being interpreted as being deterministically defined, thus ensuring the allocation of $q = q_4$ at the desired position. Despite the well-known properties of $q_4$, to the best of the authors' knowledge, there is no proof that localized oscillation would occur regardless of initial conditions, with other basins of attraction.

The investigation of this problem can be conducted through the analysis of the Lyapunov Exponents (or Lyapunov Exponents Spectrum for systems with multiple degrees of freedom), where sensitivity to initial conditions is associated with positive values of these quantities. Obtaining the Largest Lyapunov Exponent (LLE) can be computationally expensive and several methods have been proposed over the years [74, 75, 76, 77, 78, 79] for various types of problems [80, 81, 82, 83, 84, 85, 86, 87, 88, 89].

The approach proposed in [90] is employed, a practical way to obtain the LLE is established without the need for calculating the entire spectrum. For the application of the methodology, let us consider Equation (1) with
constant physical parameters and a variation $Q_n$ in the linear stiffness of the nth component (for simplicity of notation, index omitted):

$$ \ddot{\eta} + 2\gamma\dot{\eta} + \left(1 + \frac{Q_n}{k_l}\right) \eta - \eta^3 + \frac{k_c}{k_l}(2\eta - \eta_{n+1} - \eta_{n-1}) + \sqrt{\frac{knl}{kl}} \frac{A}{\omega_0^2} \cos\left(\frac{\omega}{\omega_0}t\right) = 0, \quad (24) $$

where, with the aim of generalizing the expression, the displacement was reformulated through the normalization $\eta(t) \rightarrow \sqrt{\left(\frac{k_{nl}}{k_l}\right)}\eta(\omega_0 t)$.

For the application of the method, we must write the system in a matrix format, in order to generate the following expression:

$$ \frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x}, t), \quad (25) $$

where the vector $\mathbf{x}$ is defined as:

$$ \begin{bmatrix} x_1(t) \\ x_2(t) \\ x_3(t) \\ x_4(t) \\ x_5(t) \\ x_6(t) \\ \vdots \end{bmatrix} = \begin{bmatrix} \eta_1(t) \\ \frac{d}{dt}\eta_1(t) \\ \eta_2(t) \\ \frac{d}{dt}\eta_2(t) \\ \eta_3(t) \\ \frac{d}{dt}\eta_3(t) \\ \vdots \end{bmatrix}, \quad (26) $$

and $\mathbf{F}(\mathbf{x}, t)$ defined from Equation (24). Alongside Equation (25), a vector $\mathbf{z}(t)$ (related to the difference between undisturbed and the disturbed solution) is defined as:
$$ \frac{d\mathbf{z}}{dt} = U(\mathbf{x}(t))\mathbf{z}, \qquad (27) $$

where $U(\mathbf{x}(t))$ is the Jacobian matrix of $F(\mathbf{x}, t)$ calculated in $\mathbf{x}(t)$. The LLE is reached asymptotically through numerical simulation, with convergence being described by the following expression:

$$ \frac{\mathbf{z} \circ \frac{d\mathbf{z}}{dt}}{|\mathbf{z}|^2} = \lambda^*, \qquad (28) $$

where the average value $\hat{\lambda}^*$ of $\lambda^*$ is the parameter allowing to determine stability of the system in the direction of the disturbance $\mathbf{z}$:

$$ \hat{\lambda}^* \rightarrow \text{LLE} \text{ when } t \rightarrow \infty. \qquad (29) $$

Figures 8 and 9 illustrate the computed LLE.

The numerical integration of Equation (25), when carried out directly from $\eta_n(0) = 0$ and $\dot{\eta}_n(0) = 0$ for $q_4 = 0$, would lead to flat solutions, even when imposing stable parameters for stationary solitons. Indeed, the existence of flat solutions is predicted for Equation (14) in the absence of disorder [91]. The attainment of an LLE< 0 (see Figure 8) or a LLE> 0 (see Figure 9 confirms the findings in the literature concerning the reduction of the basin of attraction of non-soliton solutions at appropriate levels of impurity.

## 5. Concluding remarks

This paper investigates the robustness of standing solitons in chains of Duffing oscillators in presence of uncertainties and disorder. The study is
<figure><img src="image_9.png" /><figcaption>Figure 8: Numerical integration of the Equation 24. Stable parameters are depicted, and a negative Largest Lyapunov Exponent (LLE) is converged. Simulation parameters are presented in Table 9</figcaption></figure>

conducted around the Damped Disordered Nonlinear Schrödinger Equation, derived using the method of multiple scales without the assumption of periodicity. The analysis methodology categorized the disorders into two groups: one unrelated and one related group to the natural. This categorization is
<figure><img src="image_10.png" /><figcaption>Figure 9: Numerical integration of Equation 24. Unstable parameters are used through the increase of acceleration (and consequent increase of $h$), where a positive Largest Lyapunov Exponent (LLE) is observed. Simulation parameters are presented in Table 9</figcaption></figure>

based on the properties reported in the literature regarding the nucleation capability that impurities in the stiffness of the oscillators possess. The results presented allow us to conclude that both types of disorder exhibit correlation; precisely, disorders in the natural frequencies decrease the probability
<table><caption>Table 9: Physical parameters used in the integration of Equation (24).</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>ε</td><td>0.01</td></tr><tr><td>N</td><td>45</td></tr><tr><td>m</td><td>1 [kg]</td></tr><tr><td>ω<sub>0</sub></td><td>100 2π [rad/s]</td></tr><tr><td>Q</td><td>-0.2k<sub>l</sub></td></tr><tr><td>γ</td><td>0.3%</td></tr><tr><td>k<sub>c</sub></td><td>0.7k<sub>l</sub></td></tr><tr><td>k<sub>nl</sub></td><td>10<sup>4</sup> k<sub>l</sub>/η<sup>2</sup></td></tr><tr><td>ω</td><td>ω<sub>0</sub>(1 - ε)</td></tr><tr><td>A</td><td>Equations (16) [m/s<sup>2</sup>]</td></tr><tr><td>Time step</td><td>0.02 [s]</td></tr><tr><td>Total integration time</td><td>2 10<sup>3</sup> [s]</td></tr><tr><td>method</td><td>DOP853</td></tr><tr><td>Tolerance</td><td>1 10<sup>-12</sup></td></tr><tr><td>Initial z vector</td><td>U[-0.2, 0.2]</td></tr></tbody></table>

of obtaining stationary solitons when compared to disorders related to other physical parameters when investigated alone. This is due to the fact that the expression for $q_4$ (stiffness disorder) also considers the perturbation factor used in the method of multiple scales.

In addition to the sensitivity to the lattice design parameters, the chaotic dynamics is investigated through the Largest Lyapunov Exponents. As shown by the integration of the physical equation, the nucleation of solitons
by $q_4$ also eliminates the basins of attraction of other solutions. Although we cannot assert that a negative LLE implies soliton dynamics, we can infer that soliton dynamics obtained through stationary oscillators (in disturbed lattices) are not sensitive to the initial conditions of integration.

This study serves as a design guide for potential applications of the NLSE, such as those mentioned in the introduction. The synthesis of this guide can be summarized as follows: the robustness of stationary solitons is ensured at uncertainty levels below 5% for physical parameters unrelated to the natural frequency and 0.5% for the natural frequency of the oscillators, for an excitation at 99% of the natural frequency (damping of 0.3% and $\varepsilon = 0.01$). To ensure greater robustness in soliton dynamics, it is recommended to impose $h$ and $\zeta$ from central zones of the periodic stability chart, where stability remains higher than near the boundaries.

**Declaration of competing interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgments**

We gratefully acknowledge support from EUR EIPHI program, Europe (Contract No. ANR 17-EURE-0002).
References

[1] E. Kenig, B. A. Malomed, M. Cross, R. Lifshitz, Intrinsic localized modes in parametrically driven arrays of nonlinear resonators, Physical Review E 80 (4) (2009) 046202. doi:10.1103/PhysRevE.80.046202.

[2] S. El-Tantawy, A. H. Salas, H. A. Alyousef, M. Alharthi, Novel approximations to a nonplanar nonlinear schrödinger equation and modeling nonplanar rogue waves/breathers in a complex plasma, Chaos, Solitons & Fractals 163 (2022) 112612. doi:10.1016/j.chaos.2022.112612.

[3] C. Boopathy, L. Kavitha, R. Ravichandran, Nonlinear modelling of soliton collision dynamics for blood flow in a stenotic artery, Materials Today: Proceedings 51 (2022) A1–A13. doi:10.1016/j.matpr.2022.03.213.

[4] K. Shehzad, A. R. Seadawy, J. Wang, M. Arshad, Multi peak solitons and btreather types wave solutions of unstable nlses with stability and applications in optics, Optical and Quantum Electronics 55 (1) (2023) 7. doi:10.1007/s11082-022-04252-z.

[5] S. M. Hernández, J. Bonetti, N. Linale, D. F. Grosz, P. I. Fierens, Soliton solutions and self-steepening in the photon-conserving nonlinear schrödinger equation, Waves in Random and Complex Media 32 (5) (2022) 2533–2549. doi:10.1080/17455030.2020.1856970.

[6] A. S. Ali, A. Govindarajan, M. Lakshmanan, Filamentation and stabilization of vortex solitons in nonlinear metamaterial waveguides, Physics Letters A 451 (2022) 128416. doi:10.1016/j.physleta.2022.128416.
[7] F. Guo, W. Dai, Model and numerical method for soliton propagation through thermal medium based on nonlinear schrödinger and heat transfer equations, Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107790. doi:10.1016/j.cnsns.2023.107790.

[8] S. H. Strogatz, R. E. Mirollo, P. C. Matthews, Coupled nonlinear oscillators below the synchronization threshold: relaxation by generalized landau damping, Physical review letters 68 (18) (1992) 2730. doi:10.1103/PhysRevLett.68.2730.

[9] A. M. Morgante, M. Johansson, G. Kopidakis, S. Aubry, Standing wave instabilities in a chain of nonlinear coupled oscillators, Physica D: Nonlinear Phenomena 162 (1-2) (2002) 53-94. doi:10.1016/S0167-2789(01)00378-5.

[10] V. Chandrasekar, J. H. Sheeba, R. G. Pradeep, R. Divyasree, M. Lakshmanan, A class of solvable coupled nonlinear oscillators with amplitude independent frequencies, Physics Letters A 376 (32) (2012) 2188-2194. doi:10.1016/j.physleta.2012.04.058.

[11] M. Antonelli, L. De Micco, H. A. Larrondo, Measuring the jitter of ring oscillators by means of information theory quantifiers, Communications in Nonlinear Science and Numerical Simulation 43 (2017) 139-150. doi:10.1016/j.cnsns.2016.05.002.

[12] A. D. Pano-Azucena, B. Ovilla-Martinez, E. Tlelo-Cuautle, J. M. Muñoz-Pacheco, L. G. de la Fraga, Fpga-based implementation of different families of fractional-order chaotic oscillators applying grünwald-
letnikov method, Communications in Nonlinear Science and Numerical Simulation 72 (2019) 516-527. doi:10.1016/j.cnsns.2019.01.014.

[13] K. Sone, Y. Ashida, T. Sagawa, Topological synchronization of coupled nonlinear oscillators, Physical Review Research 4 (2) (2022) 023211. doi:10.1103/PhysRevResearch.4.023211.

[14] P. W. Anderson, Absence of diffusion in certain random lattices, Physical review 109 (5) (1958) 1492. doi:10.1103/PhysRev.109.1492.

[15] Z. Zergoune, N. Kacem, N. Bouhaddi, On the energy localization in weakly coupled oscillators for electromagnetic vibration energy harvesting, Smart Materials and Structures 28 (7) (2019) 07LT02. doi:10.1088/1361-665X/ab05f8.

[16] S. Dowlati, N. Kacem, N. Bouhaddi, Optimal design for vibration energy harvesters based on quasi-periodic structures, Physica Scripta 97 (8) (2022) 085212. doi:10.1088/1402-4896/ac7fc5.

[17] A. Grolet, N. Hoffmann, F. Thouverez, C. Schwingshackl, Traveling and standing envelope solitons in discrete non-linear cyclic structures, Mechanical systems and signal processing 81 (2016) 75-87. doi:10.1016/j.ymssp.2016.02.062.

[18] F. Fontanela, A. Grolet, L. Salles, A. Chabchoub, N. Hoffmann, Dark solitons, modulation instability and breathers in a chain of weakly non-linear oscillators with cyclic symmetry, Journal of Sound and Vibration 413 (2018) 467-481. doi:10.1016/j.jsv.2017.08.004.
[19] A. T. Savadkoohi, C.-H. Lamarque, J. D. Goutte, On the nonlinear interactions and existence of breathers in a chain of coupled pendulums, International Journal of Non-Linear Mechanics 94 (2017) 375-379. doi:10.1016/j.ijnonlinmec.2017.01.018.

[20] A. Jallouli, N. Kacem, N. Bouhaddi, Stabilization of solitons in coupled nonlinear pendulums with simultaneous external and parametric excitations, Communications in Nonlinear Science and Numerical Simulation 42 (2017) 1-11. doi:10.1016/j.cnsns.2016.05.012.

[21] F. Fontanela, A. Grolet, L. Salles, A. Chabchoub, A. Champneys, S. Pat-sias, N. Hoffmann, Dissipative solitons in forced cyclic and symmetric structures, Mechanical Systems and Signal Processing 117 (2019) 280-292. doi:10.1016/j.ymssp.2018.08.011.

[22] Z. I. Djoufack, E. Tala-Tebue, J. Nguenang, Quantum breathers and intrinsic localized excitation associated with the modula-tional instability in 1d bose-hubbard chain, Communications in Nonlinear Science and Numerical Simulation 69 (2019) 134-147. doi:10.1016/j.cnsns.2018.07.018.

[23] A. D. Adile, F. Kenmogne, A. K. S. Tewa, H. Simo, A. M. Tahir, S. Kumar, Dynamics of a mechanical network consisting of discontinuous coupled system oscillators with strong irrational nonlinearities: resonant states and bursting waves, International Journal of Non-Linear Mechanics 137 (2021) 103812. doi:10.1016/j.ijnonlinmec.2021.103812.

[24] C. Wang, S. Tawfick, A. F. Vakakis, Irreversible energy transfer, lo-
calization and non-reciprocity in weakly coupled, nonlinear lattices
with asymmetry, Physica D: Nonlinear Phenomena 402 (2020) 132229.
doi:10.1016/j.physd.2019.132229.

[25] I. Grinberg, O. V. Gendelman, Localization in finite vibroimpact chains: Discrete breathers and multibreathers, Physical Review E 94 (3) (2016) 032204. doi:10.1103/PhysRevE.94.032204.

[26] A. J. Martínez, H. Yasuda, E. Kim, P. Kevrekidis, M. A. Porter, J. Yang, Scattering of waves by impurities in precompressed granular chains, Physical Review E 93 (5) (2016) 052224. doi:10.1103/PhysRevE.93.052224.

[27] B. Yousefzadeh, C. Daraio, Complete delocalization in a defective periodic structure, Physical Review E 96 (4) (2017) 042219. doi:10.1103/PhysRevE.96.042219.

[28] M. Hwang, A. F. Arrieta, Solitary waves in bistable lattices with stiffness grading: Augmenting propagation control, Physical Review E 98 (4) (2018) 042205. doi:10.1103/PhysRevE.98.042205.

[29] H. Yamaguchi, S. Houri, Generation and propagation of topological solitons in a chain of coupled parametric-micromechanical-resonator arrays, Physical Review Applied 15 (3) (2021) 034091. doi:10.1103/PhysRevApplied.15.034091.

[30] Y. Miyazawa, C. Chong, P. G. Kevrekidis, J. Yang, Rogue and solitary waves in coupled phononic crystals, Physical Review E 105 (3) (2022) 034202. doi:10.1103/PhysRevE.105.034202.
[31] J. R. Tempelman, K. H. Matlack, A. F. Vakakis, Topological protection in a strongly nonlinear interface lattice, Physical Review B 104 (17) (2021) 174306. doi:10.1103/PhysRevB.104.174306.

[32] H. Gao, H. Xue, Z. Gu, L. Li, W. Zhu, Z. Su, J. Zhu, B. Zhang, Y. Chong, Anomalous floquet non-hermitian skin effect in a ring resonator lattice, Physical Review B 106 (13) (2022) 134112. doi:10.1103/PhysRevB.106.134112.

[33] M. A. Mohammed, P. Grover, Phase space analysis of nonlinear wave propagation in a bistable mechanical metamaterial with a defect, Physical Review E 106 (5) (2022) 054204. doi:10.1103/PhysRevE.106.054204.

[34] A. Upadhyaya, M. N. Semenova, A. A. Kudreyko, S. V. Dmitriev, Chaotic discrete breathers and their effect on macroscopic properties of triangular lattice, Communications in Nonlinear Science and Numerical Simulation 112 (2022) 106541. doi:10.1016/j.cnsns.2022.106541.

[35] Y. V. Bebikhov, E. Naumov, M. Semenova, S. Dmitriev, Discrete breathers in a β-fput square lattice from in-band external driving, Communications in Nonlinear Science and Numerical Simulation 132 (2024) 107897. doi:10.1016/j.cnsns.2024.107897.

[36] S. Shcherbinin, Y. V. Bebikhov, D. Abdullina, A. Kudreyko, S. Dmitriev, Delocalized nonlinear vibrational modes and discrete breathers in a body centered cubic lattice, Communications in Nonlinear Science and Numerical Simulation 135 (2024) 108033. doi:10.1016/j.cnsns.2024.108033.
[37] L. Zhang, H. Han, N-soliton solutions of coupled schrödinger-boussinesq equation with variable coefficients, Communications in Nonlinear Science and Numerical Simulation 138 (2024) 108185. doi:10.1016/j.cnsns.2024.108185.

[38] M. N. Islam, O. A. İlhan, M. A. Akbar, F. B. Benli, D. Soybaş, Wave propagation behavior in nonlinear media and resonant nonlinear interactions, Communications in Nonlinear Science and Numerical Simulation 108 (2022) 106242. doi:10.1016/j.cnsns.2021.106242.

[39] J. Villarroel, M. Montero, On the effect of random inhomogeneities in kerr media modelled by a nonlinear schrödinger equation, Journal of Physics B: Atomic, Molecular and Optical Physics 43 (13) (2010) 135404. doi:10.1088/0953-4075/43/13/135404.

[40] E. Ezzoug, W. Kechiche, E. Zahrouni, Finite dimensional global attractor for a semi-discrete nonlinear schrödinger equation with a point defect, Applied mathematics and computation 217 (19) (2011) 7818-7830. doi:10.1016/j.amc.2011.02.096.

[41] P. Y. Chen, B. A. Malomed, P. L. Chu, Trapping bragg solitons by a pair of defects, Physical Review E 71 (6) (2005) 066601. doi:10.1103/PhysRevE.71.066601.

[42] V. A. Brazhnyi, C. P. Jisha, A. Rodrigues, Interaction of discrete nonlinear schrödinger solitons with a linear lattice impurity, Physical Review A 87 (1) (2013) 013609. doi:10.1103/PhysRevA.87.013609.
[43] E. Parasuraman, Effect of impurity on modulational instability of localized modes in a discrete quantum ferromagnetic spin chain, Journal of Magnetism and Magnetic Materials 539 (2021) 168348. doi:10.1016/j.jmmm.2021.168348.

[44] R. A. A. Ribama, Z. Djoufack, J. Nguenang, Breather-impurity interactions and modulational instability in a quantum 2d klein-gordon chain, The European Physical Journal B 95 (5) (2022) 86. doi:10.1140/epjb/s10051-022-00337-6.

[45] C. L. Gninzanlong, F. T. Ndjomatchoua, C. Tchawoua, Taming intrinsic localized modes in a dna lattice with damping, external force, and inhomogeneity, Physical Review E 99 (5) (2019) 052210. doi:10.1103/PhysRevE.99.052210.

[46] W. Hu, Z. Deng, T. Yin, Almost structure-preserving analysis for weakly linear damping nonlinear schrödinger equation with periodic perturbation, Communications in Nonlinear Science and Numerical Simulation 42 (2017) 298-312. doi:10.1016/j.cnsns.2016.05.024.

[47] E. Kenig, R. Lifshitz, M. Cross, Pattern selection in parametrically driven arrays of nonlinear resonators, Physical Review E 79 (2) (2009) 026203. doi:10.1103/PhysRevE.79.026203.

[48] M. Bondila, I. V. Barashenkov, M. M. Bogdan, Topography of attractors of the parametrically driven nonlinear schrödinger equation, Physica D: Nonlinear Phenomena 87 (1-4) (1995) 314-320. doi:10.1016/0167-2789(95)00126-O.
[49] I. Barashenkov, M. Bogdan, V. Korobov, Stability diagram of the phase-locked solitons in the parametrically driven, damped nonlinear schrödinger equation, Europhysics Letters 15 (2) (1991) 113. doi:10.1209/0295-5075/15/2/001.

[50] A. Barbosa, J. P. Sena, N. Kacem, N. Bouhaddi, An artificial intelligence approach to design periodic nonlinear oscillator chains under external excitation with stable damped solitons, Mechanical Systems and Signal Processing 205 (2023) 110879. doi:10.1016/j.ymssp.2023.110879.

[51] I. Barashenkov, Y. S. Smirnov, Existence and stability chart for the ac-driven, damped nonlinear schrödinger solitons, Physical Review E 54 (5) (1996) 5707. doi:10.1103/PhysRevE.54.5707.

[52] N. Alexeeva, I. Barashenkov, G. Tsironis, Impurity-induced stabilization of solitons in arrays of parametrically driven nonlinear oscillators, Physical review letters 84 (14) (2000) 3053. doi:10.1103/PhysRevLett.84.3053.

[53] W. Cardoso, S. Leão, A. Avelar, D. Bazeia, M. S. Hussein, Non-linear schrödinger equation with chaotic, random, and nonperiodic nonlinearity, Physics Letters A 374 (45) (2010) 4594-4598. doi:10.1016/j.physleta.2010.09.037.

[54] R. Kamburova, S. Varbev, M. Primatarowa, Interaction of solitons with impurities in anisotropic ferromagnetic chains, Physics Letters A 383 (5) (2019) 471-476. doi:10.1016/j.physleta.2018.11.009.

[55] O. M. Braun, Y. S. Kivshar, Nonlinear dynamics of the frenkel-
kontorova model with impurities, Physical Review B 43 (1) (1991) 1060.
doi:10.1103/PhysRevB.43.1060.

[56] N. Alexeeva, I. Barashenkov, G. Tsironis, Taming spatiotemporal chaos by impurities in the parametrically driven damped nonlinear schrödinger equation, Journal of Nonlinear Mathematical Physics 8 (sup1) (2001) 5–12. doi:10.2991/jnmp.2001.8.s.2.

[57] D. Witthaut, S. Mossmann, H. Korsch, Bound and resonance states of the nonlinear schrödinger equation in simple model systems, Journal of Physics A: Mathematical and General 38 (8) (2005) 1777. doi:10.1088/0305-4470/38/8/013.

[58] A. Barbosa, N. Kacem, N. Bouhaddi, A design methodology for non-linear oscillator chains enabling energy localization tuning and soliton stability enhancement with optimal damping, Mechanical Systems and Signal Processing 213 (2024) 111358. doi:10.1016/j.ymssp.2024.111358.

[59] A. Barbosa, N. Kacem, N. Bouhaddi, Standing solitons tuned by impurities in damped nonlinear lattices under external excitation, Phys. Rev. Lett. 133 (2024) 187201. doi:10.1103/PhysRevLett.133.187201.
URL https://link.aps.org/doi/10.1103/PhysRevLett.133.187201

[60] P. P. Galuzio, S. Benkadda, S. R. Lopes, Characterization of intermittency at the onset of turbulence in the forced and damped nonlinear schrödinger equation, Communications in Nonlinear Science and Numerical Simulation 42 (2017) 404–419. doi:10.1016/j.cnsns.2016.05.031.
[61] T. K. Vi, N. D. Q. Huy, T. C. Quy, B. D. Tinh, D. Q. Tri, M. Trippenbach, N. V. Hung, et al., A hybrid variational method for beam propagation and interaction in a graded-index nonlinear waveguide, Communications in Nonlinear Science and Numerical Simulation (2024) 108124doi:10.1016/j.cnsns.2024.108124.

[62] A. O. Oyelade, Experiment study on nonlinear oscillator containing magnetic spring with negative stiffness, International Journal of Non-Linear Mechanics 120 (2020) 103396. doi:10.1016/j.ijnonlinmec.2019.103396.

[63] S. Leadenham, A. Erturk, M-shaped asymmetric nonlinear oscillator for broadband vibration energy harvesting: Harmonic balance analysis and experimental validation, Journal of Sound and Vibration 333 (23) (2014) 6209-6223. doi:10.1016/j.jsp.2014.06.046.

[64] O. Alfahmi, A. Erturk, Programmable hardening and softening cubic inductive shunts for piezoelectric structures: Harmonic balance analysis and experiments, Journal of Sound and Vibration 571 (2024) 118029. doi:10.1016/j.jsp.2023.118029.

[65] G. Failla, A. Pirrotta, On the stochastic response of a fractionally-damped duffing oscillator, Communications in Nonlinear Science and Numerical Simulation 17 (12) (2012) 5131-5142. doi:10.1016/j.cnsns.2012.03.033.

[66] P. Balasubramanian, G. Ferrari, M. Amabili, Identification of the viscoelastic response and nonlinear damping of a rubber plate in nonlinear
vibration regime, Mechanical Systems and Signal Processing 111 (2018)
376-398. doi:10.1016/j.ymssp.2018.03.061.

[67] M. Płodzień, K. Sacha, Breakdown of anderson localization of interact-
ing quantum bright solitons in a disorder potential, Physical Review A
86 (3) (2012) 033617. doi:10.1103/PhysRevA.86.033617.

[68] K. Sacha, C. A. Müller, D. Delande, J. Zakrzewski, Anderson localization of solitons, Physical review letters 103 (21) (2009) 210402. doi:10.1103/PhysRevLett.103.210402.

[69] D. Basko, Local nature and scaling of chaos in weakly nonlin-
ear disordered chains, Physical Review E 86 (3) (2012) 036202.
doi:10.1103/PhysRevE.86.036202.

[70] B. Senyange, C. Skokos, Identifying localized and spreading chaos in nonlinear disordered lattices by the generalized alignment index (gali) method, Physica D: Nonlinear Phenomena 432 (2022) 133154. doi:10.1016/j.physd.2022.133154.

[71] A. Boudjemâa, Spontaneous formation of bright solitons in self-localized impurities in bose-einstein condensates, Communications in Nonlinear Science and Numerical Simulation 33 (2016) 85-90. doi:10.1016/j.cnsns.2015.09.003.

[72] I. S. Sawkmie, D. Kharkongor, Theoretical and numerical study of vibrational resonance in a damped softening duffing oscillator, International Journal of Non-Linear Mechanics 144 (2022) 104055. doi:10.1016/j.ijnonlinmec.2022.104055.
[73] H. Xiao, H. Qie, C. R. Bowen, Modelling of the circular edge-clamped interface of a hydraulic pressure energy harvester to determine power, efficiency and bandwidth, *Mechanical Systems and Signal Processing* 146 (2021) 107013. doi:10.1016/j.ymssp.2020.107013.

[74] M. T. Rosenstein, J. J. Collins, C. J. De Luca, A practical method for calculating largest lyapunov exponents from small data sets, *Physica D: Nonlinear Phenomena* 65 (1-2) (1993) 117-134. doi:10.1016/0167-2789(93)90009-P.

[75] H. Kantz, A robust method to estimate the maximal lyapunov exponent of a time series, *Physics letters A* 185 (1) (1994) 77-87. doi:10.1016/0375-9601(94)90991-1.

[76] B. J. Kim, G. H. Choe, High precision numerical estimation of the largest lyapunov exponent, *Communications in Nonlinear Science and Numerical Simulation* 15 (5) (2010) 1378-1384. doi:10.1016/j.cnsns.2009.05.064.

[77] A. Stefanski, Estimation of the largest lyapunov exponent in systems with impacts, *Chaos, Solitons & Fractals* 11 (15) (2000) 2443-2451. doi:10.1016/S0960-0779(00)00029-1.

[78] A. Stefański, T. Kapitaniak, Estimation of the dominant lyapunov exponent of non-smooth systems on the basis of maps synchronization, *Chaos, Solitons & Fractals* 15 (2) (2003) 233-244. doi:10.1016/S0960-0779(02)00095-4.

[79] A. Stefanski, A. Dabrowski, T. Kapitaniak, Evaluation of the largest lyapunov exponent of a time series, *Physical Review E* 85 (2017) 021003. doi:10.1103/PhysRevE.85.021003.
punov exponent in dynamical systems with time delay, Chaos, Solitons & Fractals 23 (5) (2005) 1651–1659. doi:10.1016/j.chaos.2004.06.051.

[80] R. Gharavi, V. Anantharam, An upper bound for the largest lyapunov exponent of a markovian product of nonnegative matrices, Theoretical computer science 332 (1-3) (2005) 543–557. doi:10.1016/j.tcs.2004.12.025.

[81] H. Rong, G. Meng, X. Wang, W. Xu, T. Fang, Largest lyapunov exponent for second-order linear systems under combined harmonic and random parametric excitations, Journal of Sound and Vibration 283 (3-5) (2005) 1250–1256. doi:10.1016/j.jsv.2004.07.025.

[82] R. O. Vallejos, C. Anteneodo, Largest lyapunov exponent of long-range xy systems, Physica A: Statistical Mechanics and its Applications 340 (1-3) (2004) 178–186. doi:10.1016/j.physa.2004.04.005.

[83] A. Petry, D. A. C. Barone, Preliminary experiments in speaker verification using time-dependent largest lyapunov exponents, Computer Speech & Language 17 (4) (2003) 403–413. doi:10.1016/S0885-2308(03)00029-9.

[84] D. Yılmaz, N. F. Güler, Analysis of the doppler signals using largest lyapunov exponent and correlation dimension in healthy and stenosed internal carotid artery patients, Digital Signal Processing 20 (2) (2010) 401–409. doi:10.1016/j.dsp.2009.06.016.

[85] K. E. Chlouverakis, M. J. Adams, Stability maps of injection-locked
laser diodes using the largest lyapunov exponent, Optics communica-
tions 216 (4-6) (2003) 405-412. doi:10.1016/S0030-4018(02)02357-X.

[86] A. N. Pavlov, N. B. Janson, V. S. Anishchenko, V. I. Gridnev, P. Y. Dovgalevsky, Diagnostic of cardio-vascular disease with help of largest lyapunov exponent of rr-sequences, Chaos, Solitons & Fractals 11 (5) (2000) 807-814. doi:10.1016/S0960-0779(98)00212-4.

[87] A. Giovanni, M. Ouaknine, J.-M. Triglia, Determination of largest lyapunov exponents of vocal signal: application to unilateral laryngeal paralysis, Journal of Voice 13 (3) (1999) 341-354. doi:10.1016/S0892-1997(99)80040-X.

[88] Z. Feng, H. Hu, Largest lyapunov exponent and almost certain stabil-
ity analysis of slender beams under a large linear motion of basement
subject to narrowband parametric excitation, Journal of Sound and Vi-
bration 257 (4) (2002) 733-752. doi:10.1006/jsvi.2002.5041.

[89] I. Bashkirtseva, L. Ryashko, J. M. Seoane, M. A. Sanjuán, Chaotic transitions in a tumor-immune model under chemotherapy treatment, Communications in Nonlinear Science and Numerical Simulation 132 (2024) 107946. doi:10.1016/j.cnsns.2024.107946.

[90] A. Dabrowski, Estimation of the largest lyapunov exponent from the perturbation vector and its derivative dot product, Nonlinear Dynamics 67 (2012) 283-291. doi:10.1007/s11071-011-9977-6.

[91] I. Barashenkov, E. Zemlyanaya, Travelling solitons in the ex-
ternally driven nonlinear schrödinger equation, Journal of
Physics A: Mathematical and Theoretical 44 (46) (2011) 465211.
doi:10.48550/arXiv.1107.0483.
<figure><img src="image_11.png" /></figure>

A cos(ωt)