Loading [MathJax]/jax/output/SVG/jax.js
Research article

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

  • Received: 01 March 2025 Revised: 01 May 2025 Accepted: 16 May 2025 Published: 03 June 2025
  • This paper investigates the chaotic dynamics in a fractional-order vocal fold vibration (VCV) model based on the Grünwald-Letnikov fractional derivative (GLFD). Studying the characteristics of vocal fold vibration is of great significance for revealing its vibration mechanism, the etiology of abnormal vibrations, and natural speech synthesis. Traditional vocal fold vibration models are based on integer-order systems and are unable to describe the memory effects present in real physical systems. To overcome this limitation, this paper introduces fractional derivatives and develops a high-precision numerical method to simulate the fractional-order VCV model. By incorporating nonlinear elastic and damping forces, the model can more accurately describe the complex dynamic characteristics of vocal fold vibrations, including memory effects and non-locality. The numerical simulation results reveal novel chaotic behaviors in the fractional-order VCV model, which have not been observed in integer-order models. These findings provide new insights into the possible dynamic states of vocal fold vibrations and lay the foundation for further theoretical and experimental studies on the vocal cord vibration mechanism.

    Citation: Jiaxin Zhang, Wei Zhang, Xiaoyu Li. Numerical simulation of chaotic dynamics in a fractional-order vibration model with Grünwald-Letnikov fractional derivative[J]. Networks and Heterogeneous Media, 2025, 20(2): 625-647. doi: 10.3934/nhm.2025027

    Related Papers:

    [1] Shuai Zhang, Haolu Zhang, Yulan Wang, Zhiyuan Li . Dynamic properties and numerical simulations of a fractional phytoplankton-zooplankton ecological model. Networks and Heterogeneous Media, 2025, 20(2): 648-669. doi: 10.3934/nhm.2025028
    [2] Ling Zhang, Xuewen Tan, Jia Li, Fan Yang . Dynamic analysis and optimal control of leptospirosis based on Caputo fractional derivative. Networks and Heterogeneous Media, 2024, 19(3): 1262-1285. doi: 10.3934/nhm.2024054
    [3] Yaxin Hou, Cao Wen, Yang Liu, Hong Li . A two-grid ADI finite element approximation for a nonlinear distributed-order fractional sub-diffusion equation. Networks and Heterogeneous Media, 2023, 18(2): 855-876. doi: 10.3934/nhm.2023037
    [4] Ming-Kai Wang, Cheng Wang, Jun-Feng Yin . A second-order ADI method for pricing options under fractional regime-switching models. Networks and Heterogeneous Media, 2023, 18(2): 647-663. doi: 10.3934/nhm.2023028
    [5] Yimamu Maimaiti, Zunyou Lv, Ahmadjan Muhammadhaji, Wang Zhang . Analyzing vegetation pattern formation through a time-ordered fractional vegetation-sand model: A spatiotemporal dynamic approach. Networks and Heterogeneous Media, 2024, 19(3): 1286-1308. doi: 10.3934/nhm.2024055
    [6] Jinhu Zhao . Natural convection flow and heat transfer of generalized Maxwell fluid with distributed order time fractional derivatives embedded in the porous medium. Networks and Heterogeneous Media, 2024, 19(2): 753-770. doi: 10.3934/nhm.2024034
    [7] Tong Li . Qualitative analysis of some PDE models of traffic flow. Networks and Heterogeneous Media, 2013, 8(3): 773-781. doi: 10.3934/nhm.2013.8.773
    [8] Kexin Li, Hu Chen, Shusen Xie . Error estimate of L1-ADI scheme for two-dimensional multi-term time fractional diffusion equation. Networks and Heterogeneous Media, 2023, 18(4): 1454-1470. doi: 10.3934/nhm.2023064
    [9] Leqiang Zou, Yanzi Zhang . Efficient numerical schemes for variable-order mobile-immobile advection-dispersion equation. Networks and Heterogeneous Media, 2025, 20(2): 387-405. doi: 10.3934/nhm.2025018
    [10] Wenkai Liu, Yang Liu, Hong Li, Yining Yang . Multi-output physics-informed neural network for one- and two-dimensional nonlinear time distributed-order models. Networks and Heterogeneous Media, 2023, 18(4): 1899-1918. doi: 10.3934/nhm.2023080
  • This paper investigates the chaotic dynamics in a fractional-order vocal fold vibration (VCV) model based on the Grünwald-Letnikov fractional derivative (GLFD). Studying the characteristics of vocal fold vibration is of great significance for revealing its vibration mechanism, the etiology of abnormal vibrations, and natural speech synthesis. Traditional vocal fold vibration models are based on integer-order systems and are unable to describe the memory effects present in real physical systems. To overcome this limitation, this paper introduces fractional derivatives and develops a high-precision numerical method to simulate the fractional-order VCV model. By incorporating nonlinear elastic and damping forces, the model can more accurately describe the complex dynamic characteristics of vocal fold vibrations, including memory effects and non-locality. The numerical simulation results reveal novel chaotic behaviors in the fractional-order VCV model, which have not been observed in integer-order models. These findings provide new insights into the possible dynamic states of vocal fold vibrations and lay the foundation for further theoretical and experimental studies on the vocal cord vibration mechanism.



    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.

    Figure 1.  Model of symmetrical vocal cords.

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

    {m1l¨x1l+r1l˙x1l+k1lx1l+kcl(x1lx2l)=g1,m2l¨x2l+r1l˙x2l+k2lx2l+kcl(x2lx1l)=g2,m1r¨x1r+r1r˙x1r+k1rx1r+kcr(x1rx2r)=g1,m2r¨x2r+r2r˙x2r+k2rx2r+kcr(x2rx1r)=g2, (1.1)

    where g1,g2 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.

    Table 1.  Description of parameter.
    Term Description
    m1l,m2l Mass of the left mass block
    m1r,m2r Mass of the right mass block
    k1l,k2l Stiffness of the left spring
    k1r,k2r Stiffness of the right spring
    kcl Coupling spring stiffness of on the left
    kcr Coupling spring stiffness of on the right
    r1l,r2l Damping coefficient on the left
    r1r,r2r Damping coefficient on the right

     | Show Table
    DownLoad: CSV

    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]: si(xi)=kixi(1+100x2i),i=1,2. On the other hand, for damping forces, instead of using the usual linear term rixi, nonlinear properties are used bi(xi,˙xi)=ri(1+150|xi|)˙xi,i=1,2. Therefore, in order to ensure the effectiveness of the simulation, the following fractional order VCV model is considered in this paper

    {m1lDα1l+α3ltx1l+b1l(x1l,Dα1ltx1l)+s1l(x1l)+kcl(x1lx2l)=g1,m2lDα2l+α4ltx2l+b2l(x2l,Dα2ltx2l)+s2l(x2l)+kcl(x2lx1l)=g2,m1rDα1r+α3rtx1r+b1l(x1r,Dα1rtx1r)+s1r(x1r)+kcr(x1rx2r)=g1,m2rDα2r+α4rtx2r+b1l(x2r,Dα2rtx2r)+s2r(x2r)+kcr(x2rx1r)=g2, (1.2)

    where sij(xij)=k1lxij(1+cijx2ij),i=1,2,j=l,r are the elastic recovery force. The damping forces are bij(xij,Dαijtxij)=rij(1+dij|xij|)Dαijtxij,i=1,2,j=l,r, cij,dij are coefficients. Dαtx is the GLFD. If α=[α1,α2,α3,α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.

    Definition 2.1. [37] The GLFD of α-order for the function u(t) in t[0,T] is defined as:

    Dαtu(t)=limh01hα[T/h]j=0(1)jΓ(α+1)Γ(j+1)Γ(αj+1)u(tjh). (2.1)

    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 α-order Fractional derivative for the function u(t) in t[0,T] is as follows

    Dαtu(t)=limh01hα[T/h]j=0χ(α,κ)ju(tjh), (2.2)

    where

    χ(α,κ)0=η0,χ(α,κ)k=1η0k1i=0(1i1+αk)ηiχ(α,κ)ki,k=1,2,....,κ1,χ(α,κ)k=1η0κi=0(1i1+αk)ηiχ(α,κ)ki,k=κ,κ+1,κ+2,.... (2.3)
    (η0η1η2ηκ)=(1111123κ+112232(κ+1)212p3p(κ+1)κ)1(012κ). (2.4)

    Let Dα1ltx1l=x3l,Dα2ltx2l=x4l, Dα1rtx1r=x3r,Dα2rtx2r=x4r, Eq (1.2) can be converted to Eq (2.5)

    {Dα1ltx1l=x3l,Dα2ltx2l=x4l,Dα3ltx3l=(g1b1l(x1l,x3l)s1l(x1l)kcl(x1lx2l))/m1l,Dα4ltx4l=(g2b2l(x2l,x44l)s2l(x2l)kcl(x2lx1l))/m2l,Dα1rtx1r=x3r,Dα2rtx2r=x4r,Dα3rtx3r=(g1b1r(x1r,x3r)s1r(x1r)kcr(x1rx2r))/m1r,Dα4rtx4r=(g2b2r(x2r,x4r)s2r(x2r)kcr(x2rx1r))/m2r. (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:

    {x1l(tk)=hα1lx3l(tk1)mj=1χ(α1l,κ)j(x1l(tkj)x1l(0))+x1l(0),x2l(tk)=hα2lx4l(tk1)mj=1χ(α2l,κ)j(x2l(tkj)x2l(0))+x2l(0),x3l(tk)=hα3l(g1(tk)b1l(x1l(tk1),x3l(tk))s1l(x1l(tk1))kcl(x1l(tk)x2l(tk)))/m1lmj=1χ(α3l,κ)j(x3l(tkj)x3l(0))+x3l(0),x4l(tk)=hα4l(g2(tk)b2l(x2l(tk1),x4l(tk))s2l(x2l(tk1))kcl(x2l(tk)x2l(tk)))/m2lmj=1χ(α4l,κ)j(x4l(tkj)x4l(0))+x4l(0),x1r(tk)=hα1rx3r(tk1)mj=1χ(α1r,κ)j(x1r(tkj)x1l(0))+x1r(0),x2r(tk)=hα2rx4r(tk1)mj=1χ(α2r,κ)j(x2r(tkj)x2l(0))+x2r(0),x3r(tk)=hα3r(g1(tk)b1r(x1r(tk1),x3r(tk))s1r(x1r(tk1))kcr(x1r(tk)x2r(tk)))/m1rmj=1χ(α3r,κ)j(x3r(tkj)x3r(0))+x3r(0),x4r(tk)=hα4r(g2(tk)b2r(x2r(tk1),x4r(tk))s2r(x1l(tk1))kcr(x2r(tk)x2r(tk)))/m2rmj=1χ(α4r,κ)j(x4r(tkj)x4r(0))+x4r(0). (2.6)

    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 L0 samples are retained, the summation can be approximated as:

    {mj=1χα1l,κ)jxil(tkj)min(L0,k)j=1Υ(α1,ω)jxil(tkj),i=1,2,3,4.mj=1χα1l,κ)jxir(tkj)min(L0,k)j=1Υ(α2,ω)jxir(tkj),i=1,2,3,4. (2.7)

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

    {Dα1+α3tx1+b1(x1,Dα1tx1)+s1(x1)+kc(x1x2)=Bcos(ct),Dα2+α4tx1+b2(x2,Dα2tx2)+s2(x1)+kc(x2x1)=Bcos(ct),t[0,T]. (3.1)

    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., si(xi)=kixi,bi(xi,Dαitxi)=riDαitxi,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 α[1,1,1,1], r1=0.01,k1=0.9,kc=0.5,B=2,b=2, r2=0.5,k2=0.1 are shown Figures 2 and 3. It is confirmed that the present method is effective and high precision.

    Figure 2.  Comparing time series plots of the present method with ode45 at α=[1,1,1,1], r1=0.01,k1=0.9,kc=0.5, r2=0.5,k2=0.1, B=b=2.
    Figure 3.  Comparing time series plots and phase diagrams of the present method with ode45 at α=[1,1,1,1], r1=0.01,k1=0.9,kc=0.5, r2=0.5,k2=0.1, B=b=2.

    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 α=[0.99,0.99,0.99,0.99], r1=0.01,k1=0.9,kc=0.5, r2=0.5,k2=0.1, B=b=2 are shown in Figure 4. It is confirmed that the present method is effective.

    Figure 4.  Comparing time series plots and phase diagrams of the present method with closed-form solution [39] at α=[0.99,0.99,0.99,0.99], r1=0.01,k1=0.9,kc=0.5, r2=0.5,k2=0.1, B=b=2.

    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.

    Figure 5.  Poincaré cross-section diagram at r1=1,r2=5,k1=9,k2=10,kc=0.5, α=[1,1,1,1],c=120,B=120,T=4000π/c.

    We consider the dynamic behavior of the VCVM with linear elasticity and resistance si(xi)=kixi,bi(xi,˙xi)=ri˙xi,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.

    Figure 6.  Comparing dynamical behavior of time series plots in different fractional derivative α=[1,1,1,1] and α=[0.99,1,1,1], atc=0.2,B=0.2.

    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 79. Non-stationary vocal cord oscillations, non-stationary vocal cord oscillations are typical for many speech disorders.

    Figure 7.  Comparing chaos dynamical behavior of 3D phase diagrams in integer derivatives with chaos dynamical behavior of 3D phase diagrams in fractional derivatives at c=0.2,B=0.2.
    Figure 8.  Compare 3D phase diagrams of in different parameters at α=[1,1,1,1],c=0.2,B=0.2.
    Figure 9.  Comparing chaos dynamical behavior in different parameters at α=[1,1,1,1],c=0.2,B=0.2.

    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 79. Figures 79 analyze the system's dynamic behavior using phase diagrams, which help visualize periodic, quasi-periodic, or chaotic states.

    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 (bi(xi,Dαitxi)=ri(1+ci|xi|)Dαitxi,si(xi)=kixi(1+dix2i),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 α[1,1,1,1], r1=0.01,k1=0.9,kc=0.5,B=b=0.2, r2=0.5,k2=0.1,c1=c2=0,d1=d2=100 are shown Figure 10.

    Figure 10.  Comparing numerical results of the present method with ode45 at α[1,1,1,1], r1=0.01,k1=0.9,kc=0.5,B=c=0.02, r2=0.5,k2=0.1,c1=c2=0,d1=d2=100.

    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.

    Figure 11.  C0 complexity and spectral entropy complexity for the system at α=[1,1,1,1], r1=0.01,k1=0.9,kc=0.5,b=0.8, r2=0.5,k2=0.1,c=B.

    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 12.  Lyapunov exponents at α=[1,1,1,1], r1=0.01,k1=0.9,kc=0.5,c=0.8, r2=0.5,k2=0.1.

    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.

    Figure 13.  Poincaré cross-section diagram at α=[1,1,1,1],c=B=120.

    From Figure 14, we can see that the larger the fractional derivative, the larger the amplitude. The coefficient di is greater than the coefficient ci, 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.

    Figure 14.  Comparing dynamical behavior of time series plots in different fractional derivative α[1,1,1,1] α[0.96,0.95,0.98,0.99] α[1,1.02,1.01,1], and different parameter c1,c2,d1,d2 at r1=0.01,r2=1,k1=1.5,k2=0.1,kc=1.5,c=B=0.2.

    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.

    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.

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

    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).

    The authors declare there is no conflict of interest.



    [1] Z. J. Qiao, C. L. Zhang, C. L. Zhang, X. Ma, R. H. Zhu, Z. H. Lai, et al., Stochastic resonance array for designing noise-boosted filter banks to enhance weak multi-harmonic fault characteristics of machinery, Appl. Acoust., 236 (2025), 110710. https://doi.org/10.1016/j.apacoust.2025.110710 doi: 10.1016/j.apacoust.2025.110710
    [2] Y. Yu, W. Zhou, Z. Zhang, Q. Bi, Analysis on the motion of nonlinear vibration with fractional order and time variable mass, Appl. Math. Lett., 124 (2022), 107621. https://doi.org/10.1016/j.aml.2021.107621 doi: 10.1016/j.aml.2021.107621
    [3] D. Chen, N. Wang, Z. Chen, Y. Yu, Parametrically excited vibrations in a nonlinear damped triple-well oscillator with resonant frequency, J. Vib. Eng. Technol., 10 (2022), 781–788. https://doi.org/10.1007/s42417-021-00408-5 doi: 10.1007/s42417-021-00408-5
    [4] S. Adachi, J. Yu, Two dimensional model of vocal fold vibration for sound synthesis of voice and soprano singing, J. Acoust. Soc. Am., 117 (2005), 3213–3224. https://doi.org/10.1121/1.1861592 doi: 10.1121/1.1861592
    [5] T. Wurzbacher, R. Schwarz, M. Döllinger, U. Hoppe, U. Eysholdt, J. Lohscheller, Model-based classification of nonstationary vocal fold vibrations, J. Acoust. Soc. Am., 120 (2006), 1012–1027. https://doi.org/10.1121/1.2211550 doi: 10.1121/1.2211550
    [6] D. D. Mehta, D. D. Deliyski, T. F. Quatieri, R. E. Hillman, Automated measurement of vocal fold vibratory asymmetry from high-speed videoendoscopy recordings, J. Speech Lang. Hear. Res., 54 (2011), 47–54. https://doi.org/10.1044/1092-4388(2010/10-0026) doi: 10.1044/1092-4388(2010/10-0026)
    [7] R. Cronjaeger, Die Entstehung des Primaeren Stimmklangs im Menschlichen Kehlkopf-Ein Model, Ph.D thesis, University of Braunschweig, Braunschweig, 1978.
    [8] J. C. Lucero, L. L. Koenig, Simulations of temporal patterns of oral airflow in men and women using a two-mass model of the vocal folds under dynamic control, J. Acoust. Soc. Am., 117 (2005), 1362–1372. https://doi.org/10.1121/1.1853235 doi: 10.1121/1.1853235
    [9] N. Isshiki, M. Tanabe, K. Ishizaka, D. Broad, Clinical significance of asymmetrical vocal cord tension, Ann. Otology Rhinology Laryngology, 86 (1977), 58–66. https://doi.org/10.1177/000348947708600109 doi: 10.1177/000348947708600109
    [10] I. Steinecke, H. Herzel, Bifurcations in an asymmetric vocal-fold model, J. Acoust. Soc. Am., 97 (1995), 1874–1884. https://doi.org/10.1121/1.412061 doi: 10.1121/1.412061
    [11] T. M. Onerci, Diagnosis in Otorhinolaryngology, Springer, London, 2010.
    [12] Z. J. Qiao, Y. B. He, C. R. Liao, R. H. Zhu, Noise-boosted weak signal detection in fractional nonlinear systems enhanced by increasing potential-well width and its application to mechanical fault diagnosis, Chaos, Solitons Fractals, 175 (2023), 113960. https://doi.org/10.1016/j.chaos.2023.113960 doi: 10.1016/j.chaos.2023.113960
    [13] C. Han, Y. L. Wang, Z. Y. Li, A high-precision numerical approach to solving space fractional Gray-Scott model, Appl. Math. Lett., 125 (2022), 107759. https://doi.org/10.1016/j.aml.2021.107759 doi: 10.1016/j.aml.2021.107759
    [14] X. L. Gao, H. L. Zhang, X. Y. Li, Research on pattern dynamics of a class of predator-prey model with interval biological coefficients for capture, AIMS Math., 9 (2024), 18506–18527. https://doi.org/10.3934/math.2024901 doi: 10.3934/math.2024901
    [15] H. Che, Y. L. Wang, Z. Y. Li, Novel patterns in a class of fractional reaction-diffusion models with the Riesz fractional derivative, Math. Comput. Simul., 202 (2022), 149–163. https://doi.org/10.1016/j.matcom.2022.05.037 doi: 10.1016/j.matcom.2022.05.037
    [16] Z. Li, Q. Chen, Y. Wang, X. Li, Solving two-sided fractional super-diffusive partial differential equations with variable coefficients in a class of new reproducing kernel spaces, Fractal Fract., 6 (2022), 492. https://doi.org/10.3390/fractalfract6090492 doi: 10.3390/fractalfract6090492
    [17] Q. P. Ji, J. Wang, L. X. Lu, C. F. Ge, Li-He's modified homotopy perturbation method coupled with the energy method for the dropping shock response of a tangent nonlinear packaging system, J. Low Freq. Noise Vibr. Act. Control, 40 (2021), 675–682. https://doi.org/10.1177/1461348420914457 doi: 10.1177/1461348420914457
    [18] J. H. He, Q. Yang, C. H. He, H. B. Li, E. Buhe, Pull-in stability of a fractal system and its pull-in plateau, Fractals, 30 (2022), 2250185. https://doi.org/10.1142/S0218348X22501857 doi: 10.1142/S0218348X22501857
    [19] J. H. He, C. H. He, A. A. Alsolami, A good initial guess for approximating nonlinear oscillators by the homotopy perturbation method, Facta Universitatis-Series Mech. Eng., 21 (2023), 21–29. https://doi.org/10.22190/FUME230108006H doi: 10.22190/FUME230108006H
    [20] J. H. He, Y. O. El-Dib, The enhanced homotopy perturbation method for axial vibration of strings, Facta Universitatis-Series Mech. Eng., 19 (2021), 735–750. https://doi.org/10.22190/FUME210125033H doi: 10.22190/FUME210125033H
    [21] J. H. He, Variational iteration method for autonomous ordinary differential systems, Appl. Math. Comput., 114 (2000), 115–123. https://doi.org/10.1016/S0096-3003(99)00104-6 doi: 10.1016/S0096-3003(99)00104-6
    [22] C. P. Li, Z. Q. Li, Z. Wang, Mathematical analysis and the local discontinuous Galerkin method for Caputo-Hadamard fractional partial differential equation, J. Sci. Comput., 85 (2020), 1–27. https://doi.org/10.1007/s10915-020-01353-3 doi: 10.1007/s10915-020-01353-3
    [23] C. P. Li, Z. Wang, The local discontinuous Galerkin finite element methods for Caputo-type partial differential equations: numerical analysis, Appl. Numer. Math., 140 (2019), 1–22. https://doi.org/10.1016/j.apnum.2019.01.007 doi: 10.1016/j.apnum.2019.01.007
    [24] C. P. Li, Z. Wang, The local discontinuous Galerkin finite element methods for Caputo-type partial differential equations: mathematical analysis, Appl. Numer. Math., 150 (2020), 587–606. https://doi.org/10.1016/j.apnum.2019.11.007 doi: 10.1016/j.apnum.2019.11.007
    [25] C. P. Li, Z. Wang, Non-uniform L1/discontinuous Galerkin approximation for the time-fractional convection equation with weak regular solution, Math. Comput. Simul., 182 (2021), 838–857. https://doi.org/10.1016/j.matcom.2020.12.007 doi: 10.1016/j.matcom.2020.12.007
    [26] C. P. Li, Z. Wang, Numerical methods for the time fractional convection-diffusion-reaction equation, Numer. Funct. Anal. Optim., 42 (2021), 1115–1153. https://doi.org/10.1080/01630563.2021.1936019 doi: 10.1080/01630563.2021.1936019
    [27] X. L. Gao, H. L. Zhang, Y. L. Wang, Z. Y. Li, Research on pattern dynamics behavior of a fractional vegetation-water model in arid flat environment, Fractal Fract., 8 (2024), 264. https://doi.org/10.3390/fractalfract8050264 doi: 10.3390/fractalfract8050264
    [28] X. L. Gao, Z. Y. Li, Y. L. Wang, Chaotic dynamic behavior of a fractional-order financial system with constant inelastic demand, Int. J. Bifurcation Chaos, 34 (2024), 2450111. https://doi.org/10.1142/S0218127424501116 doi: 10.1142/S0218127424501116
    [29] X. Kong, F. Yu, W. Yao, S. Cai, J. Zhang, H. Lin, Memristor-induced hyperchaos, multiscroll and extreme multistability in fractional-order HNN: Image encryption and FPGA implementation, Neural Networks, 171 (2024), 85–103. https://doi.org/10.1016/j.neunet.2023.12.008 doi: 10.1016/j.neunet.2023.12.008
    [30] F. Yu, X. Kong, W. Yao, J. Zhang, S. Cai, H. Lin, et al., Dynamics analysis, synchronization and FPGA implementation of multiscroll Hopfield neural networks with non-polynomial memristor, Chaos, Solitons Fractals, 179 (2024), 114440. https://doi.org/10.1016/j.chaos.2023.114440 doi: 10.1016/j.chaos.2023.114440
    [31] M. Chinnamuniyandia, S. Chandranb, C. J. Xu, Fractional order uncertain BAM neural networks with mixed time delays: An existence and Quasi-uniform stability analysis, J. Intell. Fuzzy Syst., 46 (2024), 4291–4313. https://doi.org/10.3233/JIFS-234744 doi: 10.3233/JIFS-234744
    [32] F. Yu, W. Zhang, X. Xiao, W. Yao, S. Cai, J. Zhang, et al., Dynamic analysis and field-programmable gate array implementation of a 5D fractional-order memristive hyperchaotic system with multiple coexisting attractors, Fractal Fract., 8 (2024), 271. https://doi.org/10.3390/fractalfract8050271 doi: 10.3390/fractalfract8050271
    [33] I. Ali, Dynamical analysis of two-dimensional fractional-order-in-time biological population model using chebyshev spectral method, Fractal Fract., 8 (2024), 325. https://doi.org/10.3390/fractalfract8060325 doi: 10.3390/fractalfract8060325
    [34] Y. W. Sha, J. Mou, S. Banerjee, Y. S. Zhang, Exploiting flexible and secure cryptographic technique for multidimensional image based on graph data structure and three-input majority gate, IEEE Trans. Ind. Inf., 20 (2024), 3835–3846. https://doi.org/10.1109/TII.2023.3281659 doi: 10.1109/TII.2023.3281659
    [35] X. H. Wang, H. L. Zhang, Y. L. Wang, Z. Y. Li, Dynamic properties and numerical simulations of the fractional Hastings-Powell model with the Grünwald-Letnikov differential derivative, Int. J. Bifurcation Chaos, 35 (2025).
    [36] K. Shah, T. Abdeljawad, On complex fractal-fractional order mathematical modeling of CO2 emanations from energy sector, Phys. Scr., 99 (2024), 015226. https://doi.org/10.1088/1402-4896/ad1286 doi: 10.1088/1402-4896/ad1286
    [37] V. Daftardar-Gejji (ed.), Fractional Calculus and Fractional Differential Equations, Springer Singapore, 2019. https://doi.org/10.1007/978-981-13-9227-6
    [38] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999.
    [39] I. Petrá˘s, Fractional-Order Nonlinear Systems, Modeling, Analysis and Simulation, Springer Science & Business Media, Heidelberg, 2011.
    [40] D. Y. Xue, Fractional Calculus and Fractional-Order Control, Science Press, Beijing, China, 2018.
    [41] D. Y. Xue, L. Bai, Numerical algorithms for Caputo fractional-order differential equations, Int. J. Control, 90 (2016), 1201–1211. https://doi.org/10.1080/00207179.2016.1158419 doi: 10.1080/00207179.2016.1158419
    [42] D. Xue, C. Zhao, Y. Chen, A modified approximation method of fractional order system, in 2006 International Conference on Mechatronics and Automation, (2006), 1043–1048. https://doi.org/10.1109/ICMA.2006.257769
    [43] Y. X. Han, J. X. Zhang, Y. L. Wang, Dynamic behavior of a two-mass nonlinear fractional-order vibration system, Front. Phys., 12 (2024), 1452138. https://doi.org/10.3389/fphy.2024.1452138 doi: 10.3389/fphy.2024.1452138
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(236) PDF downloads(12) Cited by(0)

Figures and Tables

Figures(14)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog