Numerical simulation of chaotic dynamics in a fractional-order vibration model with Grünwald-Letnikov fractional derivative

1.   Introduction

Nonlinear vibration remains a hot research topic due to its complex dynamics, including chaos, bifurcations, and frequency-energy dependence ^[1,2,3]. Qiao et al. ^[1] proposed a stochastic resonance array to enhance weak multi-harmonic fault features via noise-boosted nonlinear filtering. Yu et al. ^[2] investigated the dynamics of a fractional-order nonlinear oscillator with time-varying mass, analyzing vibration under parameter variations. Chen et al. ^[3] examined parametrically excited vibrations in a damped triple-well oscillator, exploring resonant responses in a multi-stable system, and so on. Vocal cord vibration is a kind of high-speed, complex and subtle three-dimensional nonlinear vibration, which makes it impossible to apply various advanced medical testing and imaging techniques in vocal research. So far, the mechanism of vocal cord vibration is still in the theoretical stage. Studying the characteristics of vocal cord vibration has far-reaching significance to reveal the mechanism of vocal cord vibration, the etiology of abnormal vibration and natural speech synthesis. Research on VCV model ^[4,5] can reveal the characteristics of vocal cord vibrations and is instrumental in understanding the mechanics of vocal cord vibrations. This has theoretical guidance significance and practical application value for medical fields such as Laryngopathology, artificial organs, and voice therapy.

Researchers have also been attempting to mimic human voice production, developing various mechanical and mathematical models to describe the human organs associated with voice generation, as well as studying the process of phonation. Essentially, sound is produced by the movement of the two vocal cords located on opposite sides of the larynx, which is caused by airflow generated by the lungs and transmitted through the trachea. The complexity of the vocal cords, including their histology, shape, position, etc., necessitates the analysis of specific issues in different ways. There are many arguments for simulating the voice production process, and many questions remain to be studied, with some research results already published.

The vocal cord vibration (VCV) can be simplified into various physical models, among which the mass-spring-damper system is a typical representation of the vocal cord. These VCV models simulate the vibration of vocal cords and polyps using different numbers of mass-spring blocks, such as the single-mass block model ^[6]. Early research on vocal cord vibration models (VCVM) was primarily based on this single-mass block model ^[6]. One of the simplest models is the single-degree-of-freedom single-mass block model, which was first proposed by Flanagan ^[7] and later extensively studied by numerous researchers in related fields. Flanagan ^[7] described the dynamics of the vocal cords using electroacoustic analysis and introduced a vocal cord model based on a single-mass vibration system. The self-excited vibration of this model is driven by airflow that varies with pressure. Through experimental observations and theoretical analysis, the self-excited force of this VCV model can be derived. Although this model is linear and effectively explains the kinematic vibration characteristics of the vocal cords, it falls short in elucidating the underlying principles of the self-excited force.

Lucero and Koenig ^[8] incorporated the nonlinear characteristics of vocal cord tissue into the fundamental two-mass model of the vocal cords, proposing an enhanced symmetric nonlinear two-mass model. Isshiki et al. ^[9] explored the clinical implications of asymmetric vocal cord tension. Such asymmetry may indicate pathological conditions of the vocal cords, the mechanical properties of lesions, and the instability of asymmetric vocal cords, which can qualitatively result in the production of different sounds by the larynx. Dynamic models often include typical effects, such as the impact arising from the tension imbalance between the left and right vocal cords. Building on this foundation, Steinecke and Herzel ^[10] developed an asymmetric VCV model based on the two-mass framework, as illustrated in Figure 1. As previously mentioned, the differential equations governing the motion of the upper and lower masses remain unchanged. However, the left and right vocal cords are asymmetrical, allowing for various types of motion.

According to the motion characteristics of the mass spring, the motion equation of the asymmetric VCV model is

where $g_1, g_2$ respectively represent the air pressure on the low side and high side of the mass block, which is determined by the glottic air pressure equation.

The traditional VCV model is described as an integer-order system, which means it does not account for the memory effects that are present in real-world physical systems. To address the limitations of the traditional model, the paper introduces fractional derivatives into the VCV model. In the VCVM, memory refers to the ability of the system to retain information about past events and use it to influence future behavior. In this model (1.1), Lucero and Koenig introduce the cubic feature of biological tissue elasticity ^[11]: $s_i(x_i) = k_ix_i(1+100x_i^2), i = 1, 2$. On the other hand, for damping forces, instead of using the usual linear term $r_ix_i$, nonlinear properties are used $b_i(x_i, \dot{x_i}) = r_i(1+150|x_i|)\dot{x_i}, i = 1, 2$. Therefore, in order to ensure the effectiveness of the simulation, the following fractional order VCV model is considered in this paper

where $s_{ij}(x_{ij}) = k_{1l}x_{ij}(1+c_{ij}x_{ij}^2), i = 1, 2, j = l, r$ are the elastic recovery force. The damping forces are $b_{ij}(x_{ij}, D_t^{\alpha_ij}{x_{ij}}) = r_{ij}(1+d_{ij}|x_{ij}|)D_t^{\alpha_{ij}}{x_{ij}}, i = 1, 2, j = l, r$, $c_{ij}, d_{ij}$ are coefficients. $D_t^\alpha{x}$ is the GLFD. If $\alpha = [\alpha_1, \alpha_2, \alpha_3, \alpha_4] = [1, 1, 1, 1]$, the model (1.2) is the model (1.1).

Fractional calculus represents a research domain that is both ancient and modern. Its origins can be traced back to the late 17th century, lending it a historical depth. However, it is also considered novel because, for nearly three centuries after its conceptualization, research in this area was predominantly confined to the realm of pure mathematics. In recent decades, a wealth of studies and experiments have demonstrated that fractional calculus is particularly adept at characterizing substances and processes that exhibit memory and hereditary properties. This capability stands as its primary advantage over classical calculus. As fractional calculus has found applications across various scientific and engineering disciplines, the proliferation of fractional differential equation models has surged. Consequently, this field has garnered significant attention from both theoretical and applied researchers globally.

Some scholars have applied fractional micro-product to memory modeling damping vibration systems. Qiao et al. ^[12] demonstrated that fractional-order derivatives enhance weak signal detection in nonlinear systems. Their simulations and experiments showed that adjusting fractional-order parameters optimizes system dynamics, while increasing potential-well width improves noise-boosted detection, leading to an effective weak fault diagnosis method for machinery. The authors ^[13] gave a high-precision numerical approach to solving space fractional Gray-Scott model. The authors ^[14] investigated how uncertainty in biological parameters, represented as intervals, affects the dynamics of predator and prey populations. This could involve analyzing the stability of equilibria, bifurcation points, and the overall behavior of the system. In ^[15], authors gave some novel patterns for a class of fractional reaction-diffusion models with the Riesz fractional derivative. The study in ^[16] addressed the solution of PDEs that describe super-diffusive processes with variable coefficients. These equations are solved in the context of a new type of reproducing kernel space, which might offer novel insights into the behavior of super-diffusive systems. In ^[17], authors applied Li-He's modified homotopy perturbation method for doubly-clamped electrically actuated micro beams-based micro electromechanical system and the dropping shock response of a tangent nonlinear packaging system. He ^{[18,19,20,21]} is particularly prominent in his research on the numerical solution of fractional differential equations, and has proposed a variety of mathematical methods, such as variational iteration method, homotopy perturbation method, Exp-function method, and variational semi-inverse method. These methods are widely used to solve nonlinear problems and numerical solutions of fractional differential equations. In ^[18], authors studied pull-in stability of a fractal system and its pull-in plateau. In ^[19], authors studied approximating nonlinear ossciliators by using the homotopy perturbation method. In ^[20], authors gave the enhanced homotopy perturbation method for axial vibration of strings. In ^[21], authors studied autonomous ordinary differential systems. Li ^[22,23,24] has conducted in-depth research on the theory of fractional partial differential equations, especially in the asymptotic and regularity of solutions and the reliability of numerical algorithms, and proposed a variety of numerical algorithms to solve fractional partial differential equations ^[25,26], which have achieved significant improvements in accuracy and efficiency. In ^[22], authors used the local discontinuous Galerkin method for Caputo-Hadamard fractional partial diferential equation. In ^[23,24], authors used the local discontinuous Galerkin finite element methods for Caputo-type partial diferential equations. The research in ^[27] focused on the interaction between vegetation and water resources in arid regions, using a fractional-order model to capture the complex dynamics that may not be well represented by traditional integer-order models. The work in ^[28] explored the behavior of financial systems that are modeled using fractional calculus, particularly when there is a constant inelastic demand. This could involve studying the emergence of chaotic behavior under certain conditions, which is important for understanding financial market fluctuations. and so on fractional-order equation have investigated, such as fractional-order HNN ^[29], fractional-order neural networks ^[30,31], fractional memristive hyperchaotic system ^[32], fractional-order biological population model ^[33], etc. ^[34,35,36].

The innovations of this article are mainly reflected in the following aspects:

● The article introduces the Grünwald-Letnikov fractional-order derivative into the vocal cord vibration (VCV) model for the first time, overcoming the limitations of traditional integer-order models that fail to describe memory effects present in real-world physical systems. By incorporating fractional-order derivatives, the model can more accurately describe the complex dynamic characteristics of vocal cord vibrations, including memory effects and non-locality.

● The article develops a high-precision numerical method to simulate the fractional-order vocal cord vibration model. This method effectively captures the complex dynamic behaviors arising from memory effects, providing a new tool for understanding the intricate mechanisms of vocal cord vibrations.

● Through numerical simulations, the article reveals novel chaotic behaviors in the fractional-order vocal cord vibration model, which were not observed in integer-order models. These findings provide new insights into the possible dynamic states of vocal cord vibrations and lay the foundation for further theoretical and experimental studies.

● The article introduces nonlinear elastic and damping forces into the model to more accurately describe the complex dynamic characteristics of vocal cord vibrations. By incorporating nonlinear terms, the model can better simulate the actual physical processes of vocal cord vibrations, particularly in regulating vibration amplitude as the glottal width increases.

2.   Brief description of the methodology

Definition 2.1. ^[37] The GLFD of $\alpha$-order for the function $u\left(t \right)$ in $t\in[0, T]$ is defined as:

The Grünwald-Letnikov (GL) fractional derivative has numerous numerical approximation methods ^[37,38,39]. such as The standard GL definition forms ^[37,39], short memory principle ^[38], Lubich's higher-order approximations ^[38], centered difference approximation ^[38], Riemann-Liouville equivalence form ^[40,42], Closed-form solution method ^[39], and so on ^[40,41,42].

The standard GL definition forms ^[37,39], expressing the derivative as a limit of a weighted sum involving generalized binomial coefficients and function values at backward-shifted points; while conceptually simple, its first-order accuracy and slow convergence limit its utility in high-precision applications. To address the infinite memory requirement, The short memory principle ^[38] was developed, truncating the summation to a fixed number of recent terms, thereby reducing storage and operations while introducing a controllable truncation error that scales, making it suitable for long-time simulations where early history has diminishing influence. For improved accuracy, Lubich's higher-order approximations ^[38] employ carefully designed convolution weights derived from generating functions, achieving high-precision ^[40,41,42] by matching Taylor expansions up to order, though this comes at the cost of increased complexity in weight computation and potential stability issues for non-smooth functions. The centered difference approximation offers second-order accuracy by symmetrically distributing the fractional derivative stencil around the target point, reducing phase errors in oscillatory problems, but requires specialized boundary treatments and may exhibit instability for fractional derivative near 1. Recognizing the oscillatory behavior inherent in GL weights, damped weight modifications introduce exponential or polynomial decay factors to the binomial coefficients, significantly improving stability for larger fractional derivative while maintaining reasonable accuracy, particularly useful in viscoelastic modeling and anomalous diffusion problems. The Riemann-Liouville equivalence form ^[40,42] incorporates an initial value correction term to align the GL derivative with the Riemann-Liouville definition for smooth functions, crucial for consistent initialization in fractional differential equations. Closed-form solution method ^[39] is obtained by combining the standard Grünwald-Letnikov definition approximation scheme with the short memory principle.

Theorem 2.2. ^[40,41,42] The GLFD of with $p$-order generating function and $\alpha$-order Fractional derivative for the function $u\left(t \right)$ in $t\in[0, T]$ is as follows

where

Let $D_t^{\alpha_{1l}} x_{1l} = x_{3l}, D_t^{\alpha_{2l}} x_{2l} = x_{4l}$, $D_t^{\alpha_{1r}} x_{1r} = x_{3r}, D_t^{\alpha_{2r}} x_{2r} = x_{4r}$, Eq (1.2) can be converted to Eq (2.5)

Xue ^[40,42] has put forward a set of high-precision numerical calculation methods of fractional calculus, which are not only innovative in theory, but also very effective in practical application. Xue ^[40,42] has also developed a general simulation framework based on block diagrams, which provides a new way to solve linear and nonlinear fractional differential equations. His research also includes general solutions to fractional implicit differential equations, delayed differential equations, and boundary value problems of fractional differential equations. In ^[43], authors gave numerical calculation scheme of two equation. In this paper, we give the following numerical calculation scheme of Eq (2.5) by using Theorem 2.2:

If the number of simulation points is large, the short-memory effect can be utilized to approximate the summation terms. Assuming only the most recent $L_0$ samples are retained, the summation can be approximated as:

3.   Simulating analyzed

Since the left and right models have the same structure, the following models are simulated in this article for convenience

3.1.   Linear elastic force and damping force

First, we compare the present method with ode45 and closed-form solution ^[39] to confirm the effectiveness of the present method. We consider linear elastic force and damping force, i.e., $s_i(x_i) = k_ix_i, b_i(x_i, D_t^{\alpha_i}{x_i}) = r_iD_t^{\alpha_i}{x_i}, \; i = 1, 2$. Comparing numerical results of the present method (time step $h = 0.01$) with four-five Runge-Kutta method (time step is $0.000001$) at $\alpha\in[1, 1, 1, 1]$, $r_1 = -\; 0.01, \; k_1 = -\; 0.9, \; kc = 0.5, \; B = 2, \; b = 2$, $r_2 = -\; 0.5, \; k_2 = -\; 0.1$ are shown Figures 2 and 3. It is confirmed that the present method is effective and high precision.

It can be seen from the Figures 2 and 3 that the results of this paper are almost consistent with those of the Runge-Kutta method. Since Runge-kutta method is a powerful numerical method for solving initial value problems, it is a classical numerical solution of ordinary differential equations, which can reach the fourth order accuracy, and the method has good numerical stability, which means that the stability and reliability of the numerical solution can be maintained even under large steps. It is shown that this method is of high precision and good stability. Comparing result for time series plots and phase diagrams of the present method with closed-form solution ^[39] at $\alpha = [0.99, 0.99, 0.99, 0.99]$, $r_1 = -\; 0.01, k_1 = -\; 0.9, kc = 0.5$, $r_2 = -\; 0.5, k_2 = -\; 0.1$, $B = b = 2$ are shown in Figure 4. It is confirmed that the present method is effective.

The Poincaré cross-section diagram is widely used in nonlinear dynamics and chaos theory, and it can help researchers intuitively understand the dynamic behavior of complex systems. Next, we give the Poincaré cross-section diagram 5 of the model. If there is only one fixed point and a few discrete points on the Poincare surface of section, the motion can be judged to be periodic. If there is a closed curve in the section, the motion is quasi-periodic. If there are dense points in the section, and there is a hierarchy, the motion can be judged to be in a chaotic state.

We consider the dynamic behavior of the VCVM with linear elasticity and resistance $s_i(x_i) = k_ix_i, \; b_i(x_i, \dot{x_i}) = r_i\dot{x_i}, \; i = 1, 2$. We choose different fractional derivatives and different parameters to numerical simulate. Figure 6 shows time series plots at $c = 0.2, \; B = 0.2.$ From the perspective of the vibration waveform of the mass block Figure 6(a), the amplitude of the motion of the high-side mass block (red line) is always greater than that of the low-side mass block (green line). From the view of the vibration waveform of the mass block Figure 6(b), the amplitude of the motion of the high fractional derivative mass block (red line) is always greater than that of the low fractional derivative mass block (green line), and both the high and low orders conform to the physiological phenomenon that the motion of the upper edge of the vocal cord always lags behind the lower edge. It is the phase difference between the lower and upper edges of the vocal cords that converts the kinetic energy of the air flow into the kinetic energy of the VCV. The vibration frequency of the two mass blocks is almost the same. This shows that the introduction of fractional derivatives into VCV model is effective.

The introduction of fractional derivatives adds additional degrees of freedom to the model, allowing the system to describe more complex dynamic properties such as memory effects and non-locality. Next, we show some novel non-stationary vocal cord oscillation in Figures 7–9. Non-stationary vocal cord oscillations, non-stationary vocal cord oscillations are typical for many speech disorders.

Figure 7 compares 3D phase diagrams for different fractional derivatives. The phase diagrams help visualize the system's dynamic behavior, including periodic, quasi-periodic, or chaotic states. Different fractional derivatives lead to different dynamic behaviors. The results show that the fractional derivatives significantly influence the system's dynamics, leading to novel chaotic behaviors that were not observed in integer-order models.

Figure 8 compares 3D phase diagrams for different parameters. The phase diagrams demonstrate how varying parameters (e.g., damping coefficients, stiffness) affect the system's dynamics. Some parameter sets lead to chaotic behavior. The results highlight the sensitivity of the system's dynamics to parameter changes, emphasizing the need for careful parameter selection in the model.

Figure 9 compares 2D phase diagrams for different parameters. Similar to Figure 8, the 2D phase diagrams help visualize the system's dynamic behavior under different parameter sets. The results further confirm the influence of parameters on the system's dynamics, with some parameter sets leading to chaotic behavior. Some novel chaos dynamical behavior are shown in Figures 7–9. Figures 7–9 analyze the system's dynamic behavior using phase diagrams, which help visualize periodic, quasi-periodic, or chaotic states.

3.2.   Nonlinear elastic force and damping force

Next, for the elastic force and damping force, the usual linear terms are not used, but the nonlinear terms of the elastic force and damping force ($b_i(x_i, D_t^{\alpha_i}{x_i}) = r_i(1+c_i|x_i|)D_t^{\alpha_i}{x_i}, s_i(x_i) = k_ix_i(1+d_ix_i^2), \; i = 1, 2$) are considered. The nonlinear term is used because when the width of the glottis increases, it is necessary to limit the amplitude of vocal cord vibration to ensure the effectiveness of the simulation.

Comparing numerical results of the present method (time step $h = 0.01$) with Runge-Kutta method (ode45, time step is $0.000001$) at $\alpha\in[1, 1, 1, 1]$, $r_1 = -\; 0.01, \; k_1 = -\; 0.9, \; kc = 0.5, \; B = b = 0.2$, $r_2 = -\; 0.5, \; k_2 = -\; 0.1, \; c_1 = c_2 = 0, \; d_1 = d_2 = 100$ are shown Figure 10.

Figure 10 compares the numerical results of the proposed high-precision method with the well-established Runge-Kutta method (ode45) for the case of linear elastic and damping forces. The close agreement between the two methods validates the accuracy and effectiveness of the proposed numerical scheme. The results show that the proposed method can achieve high precision even with a larger time step (0.01) compared to the Runge-Kutta method (0.000001). The simulation confirms that the proposed method is reliable for solving the fractional-order VCV model, even when nonlinear terms are introduced. It is confirmed that the present method is effective and high precision. $C0$ Complexity and spectral entropy complexity for the system at $c = 1.5, b = 0.2$ are shown Figure 11.

$C0$ complexity reflects the randomness of a time series. In Figure 11, the higher the value, the system is the more complex (such as chaotic), high spectral entropy indicates a wide frequency distribution (such as chaos), while low spectral entropy indicates energy concentrated in a few frequencies (such as periodic signals).

Figure 12 shows the Lyapunov exponents for the nonlinear elastic and damping forces case. Lyapunov exponents are used to detect chaotic behavior. Positive Lyapunov exponents indicate chaotic dynamics.

Figure 13 shows the Poincaré cross-section for the nonlinear elastic and damping forces case with specific parameters. The dense points in the Poincaré cross-section indicate chaotic behavior, further validating the presence of chaos in the system under nonlinear conditions. The results demonstrate that the system's chaotic behavior is sensitive to parameter changes, highlighting the importance of parameter selection in the model.

From Figure 14, we can see that the larger the fractional derivative, the larger the amplitude. The coefficient $d_i$ is greater than the coefficient $c_i$, and the amplitude is increasing with time. The chaotic phenomenon disappears gradually with the increase of nonlinear terms. The results show that the system's dynamic behavior is influenced by both fractional derivatives and nonlinear terms, with larger fractional derivatives leading to increased amplitude and reduced chaotic behavior.

According to the above research results, it can be concluded that in order to improve the model, we propose to introduce fractional derivatives, so that the fractional model can get more accurate and richer results. The mathematical model of vocal cord vibration gives a very good qualitative description of vocal cord pronunciation. Although it can be applied to modeling, the obtained research results are quite different from the real situation, so it is still necessary to improve the accuracy of the model. For all the mass blocks of simulated vocal cords in the system, the motion should be divided into three directions: horizontal, vertical, and vertical. If there is such a perfect fit model, it will give a better explanation of the VCV.

The proposed method is validated as accurate and effective for both linear and nonlinear cases. The fractional-order VCV model exhibits novel chaotic behaviors, which are influenced by fractional derivatives and nonlinear terms. The system's dynamics are highly sensitive to parameter changes, with some parameter sets leading to chaotic behavior. The introduction of fractional derivatives adds memory effects and non-locality to the model, allowing it to capture more complex dynamic properties of vocal cord vibrations.

4.   Conclusions

This paper investigates the chaotic dynamics in a fractional-order vocal fold vibration (VCV) model based on the Grünwald-Letnikov fractional derivative. By introducing fractional derivatives, the model accounts for memory effects present in real-world physical systems. This advancement allows the model to capture more complex dynamic properties of vocal cord vibrations, such as memory effects and non-locality, which are not addressed by traditional integer-order models. A high-precision numerical scheme was developed specifically for simulating the fractional-order VCV model. This scheme accurately captures the complex dynamics arising from the inclusion of memory effects, which are crucial for understanding the intricate behavior of vocal cord vibrations. The simulation results reveal novel chaotic behaviors in the fractional-order VCV model that were not observed in integer-order models. These findings contribute to a deeper understanding of the possible dynamic states of vocal cord vibrations and provide new insights for further theoretical and experimental studies. The study also explored the impact of nonlinear elastic and damping forces on the system's dynamics. The results show that the system's dynamic behavior is highly sensitive to parameter changes, with some parameter sets leading to chaotic behavior. The introduction of nonlinear terms further enriches the model's ability to describe real-world vocal cord vibrations. The proposed numerical method was validated against the well-established Runge-Kutta method, demonstrating its accuracy and effectiveness for both linear and nonlinear cases. The method achieves high precision even with larger time steps, making it a reliable tool for solving the fractional-order VCV model. The incorporation of nonlinear terms and the GLFD is motivated by their ability to regulate the amplitude and chaotic behavior of VCV vibrations as the glottis widens. This regulation is crucial for producing distinct vocal sounds, thereby enhancing the effectiveness of the simulation. In future work, we will continue to combine the physiological structure of the vocal cords with the vocalization process to further clarify the potential role and significance of these chaotic behaviors in sound production.

Use of AI tools declaration

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

Acknowledgments

This paper is supported by Doctoral research start-up fund of Inner Mongolia University of Technology (DC2400003130) and the Natural Science Foundation of Inner Mongolia (2024LHMS06025).

Conflicts of interest

The authors declare there is no conflict of interest.