Multiwavelet and multiwavelet packet analysis in qualitative assessment of the chaotic states

Since chaotic phenomena are common in the behavior of nonlinear systems, there is a need to identify the ranges of parameters in the system’s mathematical model in which they occur. Chaos is present in many problems in physics, mechanics, chemistry, biology, and other fields of science [1, 2, 3, 4, 5, 6, 7, 8]. Permann and Hamilton [9] applied wavelet analysis to nonlinear dynamics problems and detected low-amplitude harmonic components in motions regarded as chaotic. Glabisz [10] used wavelet analysis to identify chaotic states in a single-degree-of-freedom system under a nonconservative load and determined parameter ranges in which the system enters critical states. Staszewski and Worden [11] analyzed the statistical measures of the coefficients in the Ueda system. Jibing et al. [12] identified various types of motions in nonlinear systems, using a wavelet harmonic transform and Poincaré maps. Billings and Coca [13] presented a method of identifying chaotic systems by means of a noisy signal wavelet transform. Using wavelet analysis, Nakao et al. [14] presented space–time relations with reference to the resolution scale adopted for the vibrations of the Ginzburg–Landan oscillator operating in the chaotic motion range. Constantine and Reinhall [15] used wavelet methods of denoising chaotic signals to estimate the fractal dimension of a noiseless system. Mastroddi and Bettoli [16] used a continuous wavelet transform to analyze a nonlinear system signal in the region of a bifurcation point. Shi et al. [17] presented an advanced synchroextracting transform based on the wavelet transform, which enables precise analysis of nonlinear and non-stationary signals. Varanis et al. [18] published a broad review of methods for analyzing chaotic signals using time-frequency tools, including discrete wavelet analysis.

This article presents a qualitative assessment of the phenomenon of chaos, based on the use of the multiwavelet transform and the packet multiwavelet transform in the analysis of the response of vibratory systems. The use of packet multiwavelet transforms for this purpose has not been applied so far. Section 2 presents a short introduction to chaotic vibrations. Section 3 describes the theoretical foundations of multiwavelet analysis and multipacket analysis. In Section 4, equations are formulated, which are then used in Section 5 to carry out several illustrative numerical analyses. Conclusions are drawn in Section 6.

Chaos frequently occurs in many nonlinear dynamical systems. Chaotic states are defined as deterministic, nonlinear system behaviors in which, after a short evolution, information about the system’s initial states is lost. Dynamical systems are highly sensitive to changes in their initial conditions. Their motion trajectories, which initially differ only infinitesimally, diverge exponentially over time, resulting in significant differences after a short period.

The qualitative methods for identifying chaos include phase portrait analysis, Poincaré cross-section analysis, bifurcation diagrams, and Fourier spectral analysis. A characteristic feature of phase portraits in chaotic systems is that the trajectory fills the phase space, passing through nearly all of its points. The Poincaré cross-sections of chaotic systems take the form of so-called strange attractors – sets of points in the phase space approached by motion trajectories with different initial conditions. These attractors exhibit self-similarity, adopting a fractal structure and maintaining the same appearance for different initial conditions. Bifurcation diagrams, which depict the value of an analyzed quantity as a function of a selected parameter, are characterized in chaotic systems by multiple splittings of the graph of the examined quantity at fixed parameter values. Fourier analysis of chaotic vibrations is distinguished by the emergence of continuous power–frequency dependencies, which indicate the presence of a continuous frequency spectrum in chaotic vibration regions.

The quantitative method for identifying chaotic states is the analysis of Lyapunov exponents, which determine the average rate of exponential divergence of initially close trajectories. A positive value of the largest Lyapunov exponent signifies the system’s sensitivity to initial conditions, indicating chaotic behavior. All the aforementioned methods, which determine whether a system is in a chaotic or nonchaotic state, are presented below alongside the approach utilizing multiwavelet and multipacket analysis for the qualitative assessment of the nature of the analyzed vibrations.

In the multiwavelet analysis, the set k 2 j k{2}^{j} of scaling functions φ j l n ( x ) {\varphi }_{jl}^{n}(x) define the space V j k \hspace{.25em}{V}_{j}^{k} and the set k 2 j k{2}^{j} of wavelet functions ψ j l n ( x ) {\psi }_{jl}^{n}(x) define the space W j k \hspace{.25em}{W}_{j}^{k} [19, 20, 21, 22]. Basis functions φ j l n ( x ) {\varphi }_{jl}^{n}(x) and ψ j l n ( x ) {\psi }_{jl}^{n}(x) are obtained by dilation ( 2 j / 2 ) ({2}^{j/2}) and translation ( l ) (l) of each from k k functions, which generate the set of fundamental multiscaling functions φ n ( x ) = { φ 0 ( x ) , φ 1 ( x ) , … , φ k − 1 ( x ) } {\varphi }^{n}(x)=\{{\varphi }^{0}(x),\hspace{.5em}{\varphi }^{1}(x),\hspace{.5em}\ldots ,\hspace{.5em}{\varphi }^{k-1}(x)\} from the space V 0 k {V}_{0}^{k} and fundamental multiwavelet functions ψ n ( x ) = { ψ 0 ( x ) , {\psi }^{n}(x)=\{{\psi }^{0}(x), ψ 1 ( x ) , … , ψ k − 1 ( x ) } {\psi }^{1}(x),\hspace{.5em}\ldots ,\hspace{.5em}{\psi }^{k-1}(x)\} from the space W 0 k {W}_{0}^{k} satisfying the following relations [23, 24]: (1) φ j l n ( x ) = 2 j / 2 φ n ( 2 j x − l ) , x ∈ [ 2 − j l , 2 − j ( l + 1 ) n s ] , {\varphi }_{jl}^{n}(x)={2}^{j/2}{\varphi }^{n}({2}^{j}x-l),\hspace{1em}x\in {[}{2}^{-j}l,\hspace{.5em}{2}^{-j}(l+1){n}_{s}]\hspace{.25em}, (2) ψ j l n ( x ) = 2 j / 2 ψ n ( 2 j x − l ) , x ∈ [ 2 − j l , 2 − j ( l + 1 ) n f ] , {\psi }_{jl}^{n}(x)={2}^{j/2}{\psi }^{n}({2}^{j}x-l),\hspace{1em}x\in {[}{2}^{-j}l,\hspace{.5em}{2}^{-j}(l+1){n}_{f}], where n s {n}_{s} is the compact support length of fundamental multiscaling functions, n f {n}_{f} is the compact support length of fundamental multiwavelet functions, n = 0 , … , k − 1 n=0,\hspace{.5em}\ldots ,\hspace{.5em}k-1\hspace{.25em} , l = 0 , … , ( 2 j − 1 ) n s , f l=0,\hspace{.5em}\ldots ,\hspace{.25em}({2}^{j}-1){n}_{s,f} , j j denotes the level of approximation, and parameter k k\hspace{.25em} determines the number of fundamental multiwavelet function sets as well as the number of fundamental multiscaling function sets. Number k k\hspace{.25em} does not depend on the assumed level j j .

The set of multiwavelet functions ψ n ( x ) = { ψ 0 ( x ) , {\psi }^{n}(x)=\{{\psi }^{0}(x), ψ 1 ( x ) , … , ψ k − 1 ( x ) } {\psi }^{1}(x),\hspace{.5em}\ldots ,\hspace{.5em}{\psi }^{k-1}(x)\} is generated on the basis of multiwavelet scaling functions φ n ( x ) = { φ 0 ( x ) , φ 1 ( x ) , … , φ k − 1 ( x ) } {\varphi }^{n}(x)=\{{\varphi }^{0}(x),\hspace{.5em}{\varphi }^{1}(x),\hspace{.5em}\ldots ,\hspace{.5em}{\varphi }^{k-1}(x)\} [19, 20, 21, 22, 23]. The equations that define the functions ψ n ( x ) {\psi }^{n}(x) and φ n ( x ) {\varphi }^{n}(x) are named the two-scale refinement equations and take the following forms: (3) ψ n ( x ) = 2 ∑ j = 0 k − 1 ( g n , j ( 0 ) φ j ( 2 x ) + g n , j ( 1 ) φ j ( 2 x − 1 ) ) + ⋯ g n , j ( r ) φ j ( 2 x − r ) ) , {\psi }^{n}(x)=\sqrt{2}\mathop{\sum }\limits_{j=0}^{k-1}({g}_{n,j}^{(0)}{\varphi }^{j}(2x)+{g}_{n,j}^{(1)}{\varphi }^{j}(2x-1))+\cdots {g}_{n,j}^{(r)}{\varphi }^{j}(2x-r)), (4) φ n ( x ) = 2 ∑ j = 0 k − 1 ( h n , j ( 0 ) φ j ( 2 x ) + h n , j ( 1 ) φ j ( 2 x − 1 ) ) + ⋯ h n , j ( r ) φ j ( 2 x − r ) ) , {\varphi }^{n}(x)=\sqrt{2}\mathop{\sum }\limits_{j=0}^{k-1}({h}_{n,j}^{(0)}{\varphi }^{j}(2x)+{h}_{n,j}^{(1)}{\varphi }^{j}(2x-1))+\cdots {h}_{n,j}^{(r)}{\varphi }^{j}(2x-r)), where n = 0 , … , k − 1 n=0,\hspace{.5em}\ldots ,\hspace{.5em}k-1\hspace{.25em} and under study in this article multiplicity r r\hspace{.25em} is equal, r = 1 r=1 .

Multiwavelet functions ψ n ( x ) {\psi }^{n}(x) and scaling functions φ n ( x ) {\varphi }^{n}(x) defined by relations (3) and (4) are presented as linear combinations of the same functions, but their differences determine, respectively, low- and high-pass filter coefficients { h n , j ( 0 ) , h n , j ( 1 ) , g n , j ( 0 ) , g n , j ( 1 ) } \{{h}_{n,j}^{(0)},\hspace{.5em}{h}_{n,j}^{(1)},\hspace{.5em}{g}_{n,j}^{(0)},\hspace{.5em}{g}_{n,j}^{(1)}\} , which in the case of multiwavelet analysis are k × k k\times k matrices.

The function space W j k \hspace{.25em}{W}_{j}^{k} created by multiwavelets ψ j l n ( x ) {\psi }_{jl}^{n}(x) is an orthogonal complementation of the function space V j k \hspace{.25em}{V}_{j}^{k} to the space of the upper level of approximation V j + 1 k \hspace{.25em}{V}_{j+1}^{k} and can be defined as (5) V j + 1 k = V j k ⊕ W j k = V j − 1 k ⊕ W j − 1 k ⊕ W j k = V j − m k ⊕ W j − m k ⊕ … ⊕ W j − 1 k ⊕ W j k . {V}_{j+1}^{k}={V}_{j}^{k}\oplus {W}_{j}^{k}={V}_{j-1}^{k}\oplus {W}_{j-1}^{k}\oplus {W}_{j}^{k}={V}_{j-m}^{k}\oplus {W}_{j-m}^{k}\oplus \ldots \oplus {W}_{j-1}^{k}\oplus {W}_{j}^{k}.

According to expression (5), any square integrable function from the j j level of approximation, when the multiscaling functions are taken from the m m level of approximation, can be expressed in multiwavelet bases in the following formula [23, 24] (6) f m j ( x ) = ∑ l = 0 2 m − 1 ∑ n = 0 k − 1 c m , l n φ m , l n ( x ) + ∑ l = 0 2 v − 1 ∑ v = 0 j − 1 ∑ n = 0 k − 1 d v , l n ψ v , l n ( x ) , {f}_{m}^{j}(x)=\mathop{\sum }\limits_{l=0}^{{2}^{m}-1}\mathop{\sum }\limits_{n=0}^{k-1}{c}_{m,l}^{n}{\varphi }_{m,l}^{n}(x)+\mathop{\sum }\limits_{l=0}^{{2}^{v}-1}\mathop{\sum }\limits_{v=0}^{j-1}\mathop{\sum }\limits_{n=0}^{k-1}{d}_{v,l}^{n}{\psi }_{v,l}^{n}(x), where the coefficients of approximation are given as follows: (7) c m , l n = ∫ 2 − m l 2 − m ( l + 1 ) f m j ( x ) φ m . l n ( x ) d x , d v , l n = ∫ 2 − m l 2 − m ( l + 1 ) f m j ( x ) ψ v . l n ( x ) d x . {c}_{m,l}^{n}=\underset{{2}^{-m}l}{\overset{{2}^{-m}(l+1)}{\int }}{f}_{m}^{j}(x){\varphi }_{m.l}^{n}(x)\text{d}x,\hspace{.5em}{d}_{v,l}^{n}=\underset{{2}^{-m}l}{\overset{{2}^{-m}(l+1)}{\int }}{f}_{m}^{j}(x){\psi }_{v.l}^{n}(x)\text{d}x.

Multiwavelet approximation coefficients j c m , l n {c}_{m,l}^{n} and d v , l n {d}_{v,l}^{n} form the basis for the calculational analyses presented in this article. The decomposition of the multiwavelet transform is shown in Figure 1.

Figure 1

Decomposition of the multiwavelet transform for three branches of basic wavelet functions.

In the frame of multiwavelet packet analysis and in contrast to multiwavelet analysis, subspaces of multiwavelet functions W j k {W}_{j}^{k} are decomposed on subsequent stages of decomposition [10, 25]. The functions that create the subspaces of decomposition W j k {W}_{j}^{k} are packet functions defined as (8) w 2 i , j n ( x ) = 2 1 / 2 j ∑ k H k w i ( 2 j x − k ) , {w}_{2i,j}^{n}(x)={2}^{1/{2}^{j}}\sum _{k}{H}_{k}{w}_{i}({2}^{j}x-k), (9) w 2 i + 1 , j n ( x ) = 2 1 / 2 j ∑ k G k w i ( 2 j x − k ) . {w}_{2i+1,j}^{n}(x)={2}^{1/{2}^{j}}\sum _{k}{G}_{k}{w}_{i}({2}^{j}x-k).

Functions w 0,0 n w{\hspace{.25em}}_{\mathrm{0,0}}^{n} and w 1,0 n w{\hspace{.25em}}_{\mathrm{1,0}}^{n} are the initial points of recursive expressions (8), (9) and enable successive functions w i , j n ( x ) w{\hspace{.25em}}_{i,j}^{n}(x) to be determined. Coefficients of matrices H k = { h n , j ( 0 ) , h n , j ( 1 ) } {H}_{k}=\{{h}_{n,j}^{(0)},\hspace{.25em}{h}_{n,j}^{(1)}\} and G k = { g n , j ( 0 ) , g n , j ( 1 ) } {G}_{k}=\{{g}_{n,j}^{(0)},\hspace{.25em}{g}_{n,j}^{(1)}\} are treated as low- and high-pass filters, respectively, and are used in classical multiwavelet analysis. The set of n = k − 1 n=k-1 multiscaling functions from space V 0 k \hspace{.25em}{V}_{0}^{k} creates the fundamental multiscaling packet functions w 0,0 n {w}_{\mathrm{0,0}}^{n} , and the set of n n multiwavelet functions from space W 0 k \hspace{.25em}{W}_{0}^{k} creates the fundamental multiwavelet packet functions w 1,0 n {w}_{\mathrm{1,0}}^{n} (10) w 0,0 n = φ n ( x ) , φ n ( x ) = { φ 0 ( x ) , φ 1 ( x ) , … , φ k − 1 ( x ) } , {w}_{\mathrm{0,0}}^{n}=\hspace{.25em}{\varphi }^{n}(x),\hspace{.5em}{\varphi }^{n}(x)=\{{\varphi }^{0}(x),\hspace{.5em}{\varphi }^{1}(x),\hspace{.5em}\ldots ,\hspace{.5em}{\varphi }^{k-1}(x)\}, (11) w 1,0 n = ψ n ( x ) , ψ n ( x ) = { ψ 0 ( x ) , ψ 1 ( x ) , … , ψ k − 1 ( x ) } . {w}_{\mathrm{1,0}}^{n}=\hspace{.25em}{\psi }^{n}(x),\hspace{.5em}{\psi }^{n}(x)=\{{\psi }^{0}(x),\hspace{.5em}{\psi }^{1}(x),\hspace{.5em}\ldots ,\hspace{.5em}{\psi }^{k-1}(x)\}.

When using formulas (3) and (4), expressions (8) and (9) take the form (12) w 2 i , j n ( x ) = 2 1 / 2 j ∑ n = 0 k ∑ j ( h n , j 0 w i n ( 2 j x ) + h n , j 1 w i n ( 2 j x − 1 ) ) , {w}_{2i,j}^{n}(x)={2}^{1/{2}^{j}}\mathop{\sum }\limits_{n=0}^{k}\sum _{j}({h}_{n,j}^{0}{w}_{i}^{n}({2}^{j}x)+{h}_{n,j}^{1}{w}_{i}^{n}({2}^{j}x-1)), (13) w 2 i + 1 , j n ( x ) = 2 1 / 2 j ∑ n = 0 k ∑ j ( g n , j 0 w i n ( 2 j x ) + g n , j 1 w i n ( 2 j x − 1 ) ) . {w}_{2i+1,j}^{n}(x)={2}^{1/{2}^{j}}\mathop{\sum }\limits_{n=0}^{k}\sum _{j}({g}_{n,j}^{0}{w}_{i}^{n}({2}^{j}x)+{g}_{n,j}^{1}{w}_{i}^{n}({2}^{j}x-1)).

The packet multiwavelet approximation coefficients are calculated similarly to expression (7), as the respective scalar products of the analyzed function f m j ( x ) {f}_{m}^{j}(x) and the functions described by relations (12) and (13).

In the presented article, the set of multiscaling functions φ n ( x ) {\varphi }^{n}(x) is defined using Legendre polynomials P n ( x ) = 1 2 n n ! ∂ n ( x 2 − 1 ) n ∂ x n {P}_{n}(x)=\frac{1}{{2}^{n}n\!}\frac{{\partial }^{n}{({x}^{2}-1)}^{n}}{\partial {x}^{n}} , which are scaled up to the range of x ∈ [ 0 , 1 ] x\in {[}0,\hspace{.25em}1] according to the expression (14) φ n ( x ) = 2 n + 1 P n ( 2 x − l ) , x ∈ [ 0 , 1 ] 0 , x ∉ ( 0 , 1 ) , n = 0 , … , k − 1 . {\varphi }^{n}(x)=\left\{\begin{array}{cc}\sqrt{2n+1}{P}_{n}(2x-l),& x\in {[}0,\hspace{.25em}1]\\ 0,& x\notin (0,\hspace{.25em}1)\end{array}\right.,\hspace{1em}n=0,\hspace{.5em}\ldots ,\hspace{.5em}k-1.

Filter coefficients, which are necessary to determine the fundamental multiwavelets and fundamental multiscaling functions, are calculated from the formula (15) h n , j ( 0 ) = 2 ∫ 0 1 2 φ n ( x ) φ j ( 2 x ) d x , g n , j ( 0 ) = 2 ∫ 0 1 2 ψ n ( x ) φ j ( 2 x ) d x , {h}_{n,j}^{(0)}=\sqrt{2}\underset{0}{\overset{\frac{1}{2}}{\int }}{\varphi }^{n}(x){\varphi }^{j}(2x)\text{d}x,\hspace{.5em}{g}_{n,j}^{(0)}=\sqrt{2}\underset{0}{\overset{\frac{1}{2}}{\int }}{\psi }^{n}(x){\varphi }^{j}(2x)\text{d}x, (16) h n , j ( 1 ) = ( − 1 ) n + j h n , j ( 0 ) , g n , j ( 1 ) = ( − 1 ) n + j + k g n , j ( 0 ) . \hspace{.25em}{h}_{n,j}^{(1)}={(-1)}^{n+j}{h}_{n,j}^{(0)},\hspace{.5em}{g}_{n,j}^{(1)}={(-1)}^{n+j+k}{g}_{n,j}^{(0)}.

The exact description of the obtained Legendre multiwavelet functions is given in previous studies [19, 20, 25].

Selected multiscaling functions and multiwavelet functions of Legendre multiwavelet and multiwavelet packets are shown in Figure 2.

Figure 2

Basic multiscaling functions of Legendre multiwavelet order k = 3 (a) and basic multiwavelet functions of Legendre multiwavelet order k = 3 (b), and second term multiwavelet packet functions of Legendre multiwavelet order k = 3 for j = 3, from final decomposition stage (c).

Two systems are considered in this article. One of them is the well-known Duffing oscillator motion described by equation (17), where the elasticity force characteristic constitutes the linear function: (17) d y ( t ) 2 d t 2 + c d y ( t ) d t − α y ( t ) + β y ( t ) 3 = F cos ( ω t ) . \frac{\text{d}y{(t)}^{2}}{\text{d}{t}^{2}}+c\frac{\text{d}y(t)}{\text{d}t}-\alpha y(t)+\beta y{(t)}^{3}=F\hspace{.25em}\cos (\omega t).

Equation (17) can be used to describe the vibration of a cantilever beam placed between two magnets and being kinematically excited [7, 26, 27]. In equation (17), the parameters F F and ω \omega are, respectively, the amplitude and frequency of the exciting force, the parameter c c describes damping, and α \alpha and β \beta are responsible for the elasticity force characteristic.

The other analyzed system [28] is shown in Figure 3. It consists of two identical massless bars connected via a hinge with mass m m placed in the hinge. The hinge is supported by an elastic spring with nonlinear characteristic S S . The system is loaded with a dynamic axial force P ( t ) P(t) [28].

Figure 3

Schematic of the analyzed system.

Taking into account the relations in Figure 3 and assuming angle α \alpha as an independent variable, one can write the system’s kinetic energy and its potential energy and the generalized force as: T = 0.5 m L 2 α ′ ( t ) T=0.5\hspace{.5em}m{L}^{2}\alpha ^{\prime} (t) , V = 0.5 k L 2 sin 2 α ( t ) + 0.25 k L 4 sin 4 α ( t ) V=0.5\hspace{.5em}k{L}^{2}{\sin }^{2}\alpha (t)+0.25\hspace{.5em}k{L}^{4}{\sin }^{4}\alpha (t) , and Q = P ( t ) 2 L sin α ( t ) Q=P(t)2L\sin \hspace{.25em}\alpha (t) . When the Lagrange equation of the second kind is used, the system vibration equation assumes the form: (18) m L 2 α ″ ( t ) + k L 2 sin α ( t ) cos α ( t ) + k 1 L 4 sin 3 α ( t ) cos α ( t ) = P ( t ) 2 L sin α ( t ) . m{L}^{2}\alpha ^{\prime\prime} (t)+k{L}^{2}\hspace{.25em}\sin \hspace{.25em}\alpha (t)\cos \hspace{.25em}\alpha (t)+{k}_{1}{L}^{4}{\sin }^{3}\alpha (t)\cos \hspace{.25em}\alpha (t)=P(t)2L\hspace{.25em}\sin \hspace{.25em}\alpha (t).

Equations (17) and (18) form the basis for the numerical analysis, and the results of which are presented in the next section.

This section presents an approach for detecting chaotic states using multiwavelet analysis and packet multiwavelet analysis. This approach constitutes an alternative to the qualitative chaos identification methods described in Section 2. The conclusions emerging from a multiwavelet analysis and a packet multiwavelet analysis in the identification of the chaotic and nonchaotic states of dynamical systems are verified by spectral analysis, Poincaré methods, bifurcation diagrams, phase portraits, and an analysis of the highest Lyapunov exponent of the investigated systems.

The range of variation of excitation amplitude F, in which chaotic vibrations are generated, is corroborated by the bifurcation diagram, constructed on the basis of the numerical integration of equation (17), together with the graph of the variation of the highest Lyapunov exponent as a function of excitation amplitude F (Figure 4). The variation of the highest Lyapunov exponent versus excitation amplitude F was determined by comparing the system response whose initial conditions differ by 10⁻⁴. When integrating equation (17), zero initial conditions y ( 0 ) = 0 , y ′ ( 0 ) = 0 y(0)=0,\hspace{.5em}y^{\prime} (0)=0 0 were imposed, and the exemplary parameters ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53 were assumed.

Figure 4

Bifurcation diagram (blue) and function of variation of max. Lyapunov exponent (red) for the Duffing oscillator described by equation (17) ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53).

An analysis of the bifurcation diagram (Figure 4) shows that if F from the interval 5.505 ≤ F ≤ 8.98 5.505\le F\le 8.98 is assumed, a system nonchaotic state can be obtained, whereas if F is taken from, e.g., the interval 15.35 ≤ F ≤ 17.0 15.35\le F\le 17.0 , a chaotic state is obtained. In the further considerations, the assumed excitation amplitude values: F = 6.5 F=6.5 and F F = 16 characterize, respectively, the nonchaotic and chaotic motion of the analyzed system.

To confirm that the assumed amplitude values are responsible for the chaotic and nonchaotic system state in Figure 5, a function of the variation of the Lyapunov exponent over time, representing the average rate of divergence of the solutions, was plotted. The average value was calculated over many initial values distributed along the motion trajectory. A chaotic state of the system occurs when the Lyapunov exponent is positive. A positive Lyapunov exponent means that trajectories that are initially close quickly diverge, leading to the unpredictability of the system.

Figure 5

Convergence of the highest Lyapunov exponent over time at F = 6.5 (a) and F = 16 (b) for the Duffing oscillator described by equation (17).

The results for qualitatively identifying the chaos phenomenon or its absence (phase trajectory, Poincaré cross sections, Fourier spectral analysis) for the Duffing oscillator described by equation (17) are depicted in Figures 6 and 7. The Poincaré cross sections in Figures 6 and 7 were plotted assuming the time distance between two consecutive points to be equal to 2 π / ω 2\pi /\omega . The nature of the results in Figures 6 and 7 also confirms that the assumed amplitude values in equation (17) are responsible for the chaotic and nonchaotic system state.

Figure 6

Phase trajectory (a) Poincaré cross sections, (b) power spectra (Fourier analysis), (c) for nonchaotic signal (F = 6.5) described by equation (17) with parameters ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53).

Figure 7

Phase trajectory (a) Poincaré cross section, (b) power spectra (Fourier analysis), (c) for chaotic signal (F = 16) described by equation (17) with parameters ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53).

In addition, using packet wavelet analysis, phase trajectory, Poincaré cross section, and Fourier spectral analysis for the Duffing oscillator described by equation (17) are obtained (Figures 8 and 9). The results presented in Figures 8 and 9 are consistent with the results in Figures 6 and 7. The input data subjected to the wavelet packet transform were a sequence of N = 1,024 numbers changing during the signal duration y(t) (for F = 6.5 and F = 16) with a constant sampling frequency of Δ t = 0.01 \text{Δ}t=0.01 . The discretized signal was subjected to packet wavelet transform using Walsh packet functions that use Haar filters. The effect of the high-pass Haar filter is the process of differentially computing the scaled derivative of the analyzed signal. Appropriate weight multipliers that scale the decomposed signal and depending on the level of decomposition and the selected package were calculated based on the formulas given in the work [29]. This made it possible to obtain the correct discrete values of the function and its derivatives. Hence, phase portraits and Poincaré sections of the discrete signals could be constructed using Walsh-wavelet packet analysis. Signals from various frequency levels obtained as a result of the packet wavelet transform of the input signal y(t) were subjected to Fourier analysis. This made it possible to select a signal representation in which it is easier to detect more subtle changes in the distribution of the signal power spectrum, indicating more critical changes in the dynamics of the system (Figures 8c and 9c). Presented in this way, the packet-based Walsh wavelet analysis seems to be a particularly valuable tool, among others for qualitative identification of critical states, because it can be used for any discrete measurement data and can be a valuable source of information about the current state of measured real objects (designation: y(t)(n,k), where: “n” – nth decomposition level of the packet wavelet analysis and “j” – jth packet number from the nth decomposition level).

Figure 8

Phase trajectory (a), Poincaré cross sections (b), and power spectra (Fourier analysis) (c) obtained from wavelet packet transform of nonchaotic signal (F = 6.5) described by equation (17).

Figure 9

Phase trajectory (a), Poincaré cross section (b), and power spectra (Fourier analysis) (c) obtained from wavelet packet transform of chaotic signal (F = 16) described by equation (17).

A multiwavelet analysis is based on a system response in the form of a discrete set of displacements y ( t ) y(t) with number N = 2 9 = 512 N={2}^{9}=512 of equidistant readings y i \hspace{.25em}{y}_{i} from the first 200 s of the response duration of the system described by equation (17). Legendre wavelets of order 3 (L_k3) were used in the analysis of the signal. Since the multiwavelet functions from space V 0 {V}_{0} are determined over interval t ∈ [ 0 , 1 ] t\in {[}0,1] while the signal was determined over any interval t ∈ [ t 0 , t k ] t\in {[}{t}_{0},\hspace{.5em}{t}_{k}] , the analyzed signal was rescaled by multiplying it by parameter τ = ( t − t 0 ) / ( t k − t 0 ) \tau =(t-{t}_{0})/({t}_{k}-{t}_{0}) in order to avoid problems with boundary conditions. As a result, the considered problem comes down to a problem determined over interval t ∈ [ 0 , 1 ] t\in {[}0,1] .

Figure 10a and b depicts the time distribution of the wavelet expansion coefficients of the response of system (17) in, respectively, nonchaotic state ( F = 6.5 ) (F=6.5) and chaotic state ( F = 16 ) (F=16) at the particular signal perception levels. The presented resolution levels j j correspond to the analyzed system’s frequencies in a range of 0.712–1.282 Hz. The higher the level j j , the higher the vibration frequency.

Figure 10

Selected resolution levels j of multiwavelet analysis coefficients of the Duffing oscillator described by equation (17) for nonchaotic signal F = 6.5 (left part of figure) and chaotic signal F = 16 (right part of figure).

As one can see, there are differences in the distribution of multiwavelet analysis coefficients between the chaotic signal and the nonchaotic signal, especially at the lower resolution levels, whereby the two states can be easily distinguished from each other. In the chaotic motion range, the coefficients are less ordered than for nonchaotic motions, which is corroborated by Figures 11 and 12. Figure 11 applies to the Duffing oscillator at the equation parameters ( ω \omega = 3.3, c c = 0.8, α \alpha = 12, β \beta = 100).

Figure 11

Multiwavelet expansion coefficients obtained at resolution j = 7 of nonchaotic signal F = 6.5 (a) and chaotic signal F = 16 (b) for the Duffing system described by equation (17) ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53).

Figure 12

Multiwavelet expansion coefficients obtained at resolution level j = 7 of nonchaotic signal F = 2.5 F=2.5\hspace{.25em} (a) and chaotic signal F = 1.329 (b) for the Duffing oscillator described by equation (17) ( ω \omega = 3.3, c c = 0.8, α \alpha = 12, β \beta = 100).

When plotting a diagram of the energy of the analyzed signals, calculated as the sum square of the wavelet expansion coefficients as a function of the number of the coefficients (ordered from the lowest to the highest absolute value) taken into account, one can notice that a larger number of expansion coefficients than in the case of nonchaotic states is needed for a chaotic state at a given energy level (Figure 13).

Figure 13

Nonchaotic signal energy (F = 6.5, F = 2.5 – dashed line) and chaotic signal energy (F = 16, F = 1.329 – solid line) versus the number of multiwavelet expansion coefficients for the Duffing oscillator: ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53) (a) and ( ω \omega = 3.3, c c = 0.8, α \alpha = 12, β \beta = 100) (b).

Significant differences between the nonchaotic signal and the chaotic signal observed in the multiwavelet analysis are even more apparent in a packet multiwavelet analysis. Figures 14 and 15 show selected levels of the wavelet decomposition of, respectively, a nonchaotic signal (F = 6.5) and a chaotic signal (F = 16). The resolution levels j j correspond to signal frequencies from a range of 0.176 to 0.318 Hz. The differences between the chaotic state and the nonchaotic state are similar to the case of the multiwavelet analysis, which confirms the effectiveness of multiwavelet packet analysis in the qualitative identification of the chaos phenomenon.

Figure 14

Selected resolution levels j of multiwavelet packet analysis coefficients for nonchaotic signal F = 6.5 of the Duffing oscillator described by equation (17) using packets of Legender’s multiwavelets k3_2.

Figure 15

Selected resolution levels j of multiwavelet packet analysis coefficients for chaotic signal F = 16 of the Duffing oscillator described by equation (17) using packets of Legender’s multiwavelets k3_2.

A classical wavelet analysis and packet wavelet analysis of a chaotic signal and a nonchaotic signal carried out using a Daubechies filter of order 12 and taking into account periodic boundary conditions lead to the same conclusions in the qualitative identification of the chaos phenomenon as does a multiwavelet analysis. The expansion coefficients of the nonchaotic signal in the base of wavelet or packet wavelet functions are more regular, and the number of nonzero expansion coefficients is, as a rule, higher than at chaotic states, particularly at high frequency levels of the investigated signals (Figures 16 and 17).

Figure 16

Selected resolution levels j = {4, 5, 6, 7} of wavelet analysis coefficients for nonchaotic signal F = 6.5 (a) and chaotic signal F = 16 (b) of Duffing oscillator described by equation (17) ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53).

Figure 17

Expansion coefficients of wavelet analysis (a) and (b) and packet wavelet analysis (c) and (d), obtained for nonchaotic signal F = 6.5 (a), (c) and chaotic signal F = 16 (b), (d) of Duffing oscillator described by equation (17) ( ω \omega = 0.7, c c = 0.1, α \alpha = 0.2, β \beta = 0.53).

If the system shown in Figure 3, described by equation (18), is analyzed for the exemplary parameters L = 0.5 m , L=0.5\hspace{.5em}\text{m}, m = 1 kg , k = 100 N m , k 1 = 1,000 N m 3 , ω = 10 rad/s m=1\hspace{.5em}\text{kg},\hspace{.5em}k=\frac{100\hspace{.5em}\text{N}}{m},\hspace{.25em}{k}_{1}=\frac{\mathrm{1,000}\hspace{.5em}\text{N}}{{m}^{3}},\hspace{.25em}\omega =10\hspace{.5em}\text{rad/s} , it can be demonstrated that stability is lost at a force P cr ≥ 28.49 N {P}_{\text{cr}}\ge 28.49\hspace{.5em}\text{N} . The value of the force at which stability loss occurs also depends on the frequency with which force P P changes, i.e., on ω \omega . For example, for ω = 15 rad/s \omega =15\hspace{.5em}\text{rad/s} , the stability of the system is lost at P cr ≥ 53.5 N {P}_{\text{cr}}\ge 53.5\hspace{.5em}\text{N} . The critical load refers to the state in which the motion of the system becomes unbounded. One indicator of stability loss is a sudden jump in the value of the Lyapunov coefficient at the point where the critical force occurs. This phenomenon is shown for different values of ω \omega in Figures 18 and 19. The bifurcation diagrams (Figures 18a and 19a) indicate that for the assumed parameters of equation (18), the system responses in both cases ( ω = 10 rad/s \omega =10\hspace{.5em}\text{rad/s} , ω = 15 rad/s \omega =15\hspace{.5em}\text{rad/s} ) are within the area of chaotic solutions.

Figure 18

Bifurcation diagram (a) and variation function of max. Lyapunov exponent (b) for the system described by equation (18) L = 0.5 m , m = 1 kg , k = 100 N m , k 1 = 1,000 N m 3 , ω = 10 rad/s \left(\phantom{\rule[-0.75em]{}{0ex}},L=0.5\hspace{.5em}\text{m},\hspace{.5em}m=1\hspace{.5em}\text{kg},\hspace{.5em}k=\frac{100\hspace{.5em}\text{N}}{m},\hspace{.5em}{k}_{1}=\frac{\mathrm{1,000}\hspace{.5em}\text{N}}{{m}^{3}},\hspace{.5em}\omega =10\hspace{.5em}\text{rad/s}\right) .

Figure 19

Bifurcation diagram (a)and variation function of max. Lyapunov exponent (b) for the system described by equation (18) L = 0.5 m , m = 1 kg , k = 100 N m , k 1 = 1,000 N m 3 , ω = 15 rad/s \left(\phantom{\rule[-0.75em]{}{0ex}},L=0.5\hspace{.5em}\text{m},\hspace{.5em}m=1\hspace{.5em}\text{kg},\hspace{.5em}k=\frac{100\hspace{.5em}\text{N}}{m},\hspace{.5em}{k}_{1}=\frac{\mathrm{1,000}\hspace{.5em}\text{N}}{{m}^{3}},\hspace{.5em}\omega =15\hspace{.5em}\text{rad/s}\right) .

The distribution of the expansion coefficients of the signal multiwavelet analysis for the case when the system is in the pre-critical phase and the post-critical phase (Figure 20) was checked for the second system described by equation (18). Due to the fact that in the two states the analyzed system is sensitive to a change in the initial conditions, and so shows the signs of the chaos phenomenon, the distribution of the multiwavelet coefficients in chaotic states is irregular and unpredictable, similarly to the previous problem. One can notice that a larger number of nonzero wavelet expansion coefficients occur when the system is unstable (Figures 21 and 22).

Figure 20

Selected resolution levels j of multiwavelet signal analysis coefficients in pre-critical state P = 25 N (a) and in post-critical state P = 30 N (b) for the system described by equation (18) L = 0.5 m , m = 1 kg , k = 100 N m , k 1 = 1,000 N m 3 , ω = 10 rad/s \left(\phantom{\rule[-0.75em]{}{0ex}},L=0.5\hspace{.5em}\text{m},\hspace{.5em}m=1\hspace{.5em}\text{kg},\hspace{.5em}k=\frac{100\hspace{.5em}\text{N}}{m},\hspace{.5em}{k}_{1}=\frac{\mathrm{1,000}\hspace{.5em}\text{N}}{{m}^{3}},\hspace{.25em}\omega =10\hspace{.5em}\text{rad/s}\right) .

Figure 21

Multiwavelet expansion coefficients obtained at resolution level j = 6 in pre-critical state P = 25 N (a) and in post-critical state P = 30 N (b) for the system described by equation (18) L = 0.5 m , m = 1 kg , k = 100 N m , k 1 = 1,000 N m 3 , ω = 10 rad/s \left(\phantom{\rule[-0.75em]{}{0ex}},L=0.5\hspace{.5em}\text{m},\hspace{.5em}m=1\hspace{.5em}\text{kg},\hspace{.5em}k=\frac{100\hspace{.5em}\text{N}}{m},\hspace{.5em}{k}_{1}=\frac{\mathrm{1,000}\hspace{.5em}\text{N}}{{m}^{3}},\hspace{.5em}\omega =10\hspace{.5em}\text{rad/s}\right) .

Figure 22

Multiwavelet expansion coefficients obtained at resolution level j = 6 in pre-critical state P = 40 N (a) and in post-critical state P = 55 N (b) for the system described by equation (18) L = 0.5 m , m = 1 kg , k = 100 N m , k 1 = 1,000 N m 3 , ω = 15 rad/s \left(\phantom{\rule[-0.75em]{}{0ex}},L=0.5\hspace{.5em}\text{m},\hspace{.5em}m=1\hspace{.5em}\text{kg},\hspace{.5em}k=\frac{100\hspace{.5em}\text{N}}{m},\hspace{.5em}{k}_{1}=\frac{\mathrm{1,000}\hspace{.5em}\text{N}}{{m}^{3}},\hspace{.5em}\omega =15\hspace{.5em}\text{rad/s}\right) .

This article demonstrates the effectiveness of multiwavelet analysis and packet multiwavelet analysis in the qualitative evaluation of chaotic behavior in the responses of nonlinear oscillatory systems. Multiwavelet analysis enables an accurate decomposition of the signal into components of different frequencies, allowing one to capture subtle changes in the signal structure that make it easier to indicate a chaotic state. Employing packet multiwavelet analysis further decomposes the signal into smaller subspaces, thereby increasing the sensitivity for detecting such changes. The system response, subjected to multiwavelet and packet multiwavelet analysis, is evaluated by comparing the distribution of multiwavelet coefficients of the signal expansions in nonchaotic and chaotic states. In addition, the cumulative energy of the multiwavelet approximation is calculated, which serves to distinguish between these two states. The results obtained via multiwavelet methods are verified using traditional qualitative approaches, such as Fourier analysis, Poincaré sections, bifurcation diagrams, as well as a quantitative method like Lyapunov exponent analysis. Furthermore, in the qualitative identification of critical states, Walsh wavelet packet analysis was applied as an innovative tool particularly useful for detecting critical states having at its disposal only a set of discrete measurement data.

The numerical analyses conducted have led to the following detailed conclusions:

Chaotic states observed through multiwavelet analysis are characterized by a relatively small number of coefficients with comparable absolute values, resulting in an uneven distribution of these coefficients compared to nonchaotic states.

It should be emphasized that the analytical answers obtained using qualitative methods in the assessment of chaotic states are not always clear, unlike quantitative methods. Therefore, they should be confirmed by several methods. Hence, the methods presented in the article based on multiwavelet and packet multiwavelet analysis constitute an additional qualitative tool for identifying chaotic states of the analyzed systems and make a significant contribution to the development of fast qualitative methods for assessing critical states in nonlinear systems.

Author states no funding involved.

The author Kamila Jarczewska confirms the sole responsibility for the conception of the study, presented results and manuscript preparation. Kamila Jarczewska confirms the sole responsibility for Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Visualization, Writing – original draft, Writing – review & editing, Supervision.

The author declares that there are no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Not applicable: the work involves no human participants, animals or proprietary datasets.

All data generated or analysed during this study are contained within the article; further details (e.g., numerical code) are available from the author upon reasonable request.