Analytic Matrix Method for Frequency Response Techniques Applied to Nonlinear Dynamical Systems I: Small and Medium Amplitude Oscillations

This is the first on a series of articles that deal with nonlinear dynamical systems under oscillatory input that may exhibit harmonic and non-harmonic frequencies and possibly complex behavior in the form of chaos. Frequency response techniques of nonlinear dynamical systems are usually analyzed with numerical methods because, most of the time, analytical solutions turn out to be difficult, if not impossible, since they are based on infinite series of trigonometric functions. The analytic matrix method reported here is a direct one that speeds up the solution processing compared to traditional series solution methods. In this method, we work with the invariant submanifold of the problem, and we propose a series solution that is equivalent to the harmonic balance series solution. However, the recursive relation obtained for the coefficients in our analytical method simplifies traditional approaches to obtain the solution with the harmonic balance series method. This method can be applied to nonlinear dynamic systems under oscillatory input to find the analog of a usual Bode plot where regions of small and medium amplitude oscillatory input are well described. We found that the identification of such regions requires both the amplitude as well as the frequency to be properly specified. In the second paper of the series, the method to solve problems in the field of large amplitudes will be addressed.


Introduction
Frequency response is a powerful technique to analyze the characteristics and properties of dynamical systems, describing the system response of a sinusoidal input with variable frequency. The system's output magnitude and amplitude are both functions of the frequency; the output and input are then compared. The first frequency response technique (FRT) emerged in the 1930s during the development of the feedback amplifier stability [1]. Later on, it involved transfer functions to produce Nyquist and Bode plots, which are mainly used to study feedback amplifiers. Meanwhile, its antagonist, the time response approach, involved ordinary differential equations with their associated algebraic equations to solve problems related to automatic control systems needed in mechanical, naval and chemical engineering during WW1 [2]. However, as time went by, several FRTs were developed to solve problems in chemical engineering [3] and physical chemistry [4], among other areas.
In dynamical systems, when specific changes in the input (e.g., doubling the amplitude, two inputs summed, etc.) are reflected in the output spectrum (doubled amplitude, equivalent summed outputs, etc.), the systems are linear. For the specific case of timeinvariant systems, the system response does not vary with time. Laplace transform is the method of solution of a linear system due to its simplicity by reducing ordinary differential equations to algebraic equations. For the frequency-response case, well-established stability criteria are important in the development of the frequency response methods in automatic control. To consult the classical FRT, the reader is referred to the collected works by MacFarlane [2].
Opposite to the behavior of a stable linear system under a periodic input (where an output of the same frequency is produced), a nonlinear system output oscillates with different frequencies with non-periodic harmonics that are not multiples of each other. The behavior of a nonlinear system that is not stable can be more complicated, i.e., one that is not at equilibrium nor periodic with near-periodic oscillations that evolve to chaos. Even though the system is deterministic in nature, some of these chaotic outputs exhibit randomness. The nonlinear system's output may exhibit several modes of behavior and may have more than one limit cycle if the system is unforced. Conversely, a forced system with periodic input may produce complicated steady-state outputs that exhibit harmonic or subharmonic behavior; the kind of output produced usually depends on the amplitude and frequency of the input. It is worthwhile to mention that in some cases, when the amplitude or frequency of the input is smoothly changed, the nonlinear system may even exhibit a discontinuous jump in the output mode [2].
Among the methods proposed to solve the equations that arise in frequency response techniques, the Galerkin's and Harmonic balance methods belong to the class of series solutions methods. Galerkin proposed a method to solve partial differential equations (PDE) that were too difficult to be solved analytically. This method turned out to be very useful to analyze complex PDE to find approximated solutions addressing the main concerns at the time: what residuals were to be used and how to be sure that the solutions actually converged. Galerkin was focused on the approximation to the solution by proposing a set of linearly independent test functions that are orthogonal to the residual. This method was tested first with a biharmonic equation and the results were published in 1915 [5]. Because of the flexibility of Galerkin's approach, nowadays most of the known methods fall in this category, adopting the name of the contributor who developed a particular application [6]. The method of harmonic balance is another method belonging to the category of series solutions methods. In this case, the solution is proposed as the summation of exponential terms that capture the dependence with the amplitude, frequency and harmonics or sub-harmonics. Lately, this method has been used in several mechanical applications such as a bio-inspired work [7] or the study of nonlinear vibrations of a thin, polar-orthotropic circular plate [8]. The method of harmonic balance is a well-established method for the analysis of, for example, the power-frequency response to illustrate the effects of geometrical nonlinearity [9]. However, in the literature, the application of the method of harmonic balance and other methods to analyze the behavior of nonlinear systems subjected to FRT has been mostly using numeric approaches to particular study cases (see, for instance, [3,[10][11][12][13]), and solving equations with the methods described previously can be challenging since the analytic approach may not be solvable and the numerical one does not provide an exact solution. To the best of our knowledge, there are no reported works on the general analytic solution for FRT.
In this paper, we present an analytic method for FRT with applications to nonlinear dynamical systems that is an improvement over the series solution methods available. This method provides a structured procedure to obtain information of nonlinear systems subjected to periodic inputs for the cases where the output is stable and others where the system is not at equilibrium. Hereby, we present some criteria to identify different modes of behavior such as the small amplitude oscillatory (SAOs) and the medium amplitude oscillatory (MAOs) that arise in several FRT applications. In the method reported here, we work with the invariant submanifold of the problem that is equivalent to the harmonic balance series solution. The recursive relation obtained for the coefficients in our analytical method simplifies traditional approaches to obtain the solution with the harmonic balance series method.

Preliminary
Consider a pair of cascade interconnected systems of the form (see Figure 1).
where x ∈ R n , z ∈ R m , f (0, 0) = 0, s(0) = 0, and f (x, z), s(z) are locally Lipschitz in R n × R m .ẋ andż represent the time derivative of x and z, respectively. The state variables z of system (2) can be seen as an input vector of system (1); therefore, system (2) is called the driving (or master) system of the tested (or slave) system (1), which is considered to be input-state stable. According to the input-state stability theory of dynamical systems [14], the tested system input-state is stable if for z = 0 the unforced system (1) is asymptotically stable, i.e., that x(t) → 0 as t → ∞, and for a input z(t) with finite supreme norm z(·) ∞ = sup t≥0 z(t) , it holds that where γ(·) is a class K function (a continuous function γ : [0, a) → [0, ∞) is said to belong to the class K if it is strictly increasing and γ(0) = 0 [15].) and β is a class KL function (a continuous function β is said to belong to the class KL if, for each fixed s, the mapping β(r, s) belongs to class K with respect to r and, for each fixed r, the mapping β(r, s) is decreasing with respect to s and β(r, s) → 0 as r → ∞ [15]).Then, given the properties of the class KL function β, it holds that lim t→∞ x(t) = γ( z(·) ∞ ), i.e., that the behavior of x(t) synchronizes with the behavior of z(t). This implies that the state variable, x(t), reaches an invariant submanifold, x(t) = ξ(z(t)), that depends exclusively on the behavior of z(t), i.e., any trajectory x(t) of the system (6) is the image through the mapping π(·) of a trajectory of the driving system. Note that the mapping ξ : z → x is in general an immersion, i.e., the rank is equal to the dimension of z that satisfies the partial differential equation In particular, if the driving system (2) is Poisson stable, i.e., its dynamical behavior has persistent oscillations with maximum amplitude A, then lim t→∞ x(t) = γ(|A|), and the tested system (1) might also tend to a Poisson stable behavior in the invariant submanifold, with maximum amplitude bounded by γ(|A|). This result can be applied to the analysis of FRT, as shown in the following section.

Frequency Response Analysis
Consider the input-state stable dynamical systeṁ where x ∈ X ⊂ R n is the state vector, while z ∈ R 2 is the input vector and f : X × R 2 → R n is a continuous function. In particular, in traditional FRT, z has a linear oscillatory behavior with the following dynamicsż 1 = ωz 2 (7) where ω is the frequency of oscillation. The behavior of z 1 and z 2 depends on the initial conditions and has the general form: z 1 (t) = A sin(ωt + φ) and z 2 (t) = A cos(ωt + φ), where A and φ are the amplitude and phase of oscillation. Notice that the oscillator (7) and (8) is Poisson stable and is a particular case of system (2), while system (6) is similar to system (1) assuming input-state stability; therefore, x(t) will reach an invariant submanifold, x(t) = ξ(z(t)), that satisfies a partial differential equation analogous to Equation (5), i.e., where for the oscillator (7) and (8), in particular, the time derivative ofξ iṡ For instance, consider a simple pendulum swinging under gravity and frictional viscous damping, whose dynamical behavior can be represented by the following equation where θ is the angular position, m is the mass concentrated in the tip, l is the length of the pendulum, g is the gravity acceleration, µ is the friction coefficient and T is the applied torque. For simplicity purposes, let us define the following dimensionless variables τ = g/lt, η = µ/ m gl , x 1 = θ, x 2 = dθ/dτ and u = T/(mgl), then the dynamical model becomesẋ Notice that when u = 0, the steady-state is given by x 2 = 0 and sin(x 1 ) = 0; therefore, x 1 = ±nπ, for n = 0, 1, 2, 3, . . ., and given the stability properties, the only stable equilibrium points are x 1 = ±2nπ, for n = 0, 1, 2, 3, . . .. Then, if the pendulum is submitted to an oscillatory torque: u(t) = z 1 (t), where z 1 is the first element of the state vector of the oscillator (7) and (8), for small frequencies, after a transient behavior that depends on the initial conditions, the pendulum will approach to an invariant submanifold x(t) = ξ(z(t)), which satisfies the following pair of nonlinear coupled partial differential equations: In the following section, a method is proposed to solve this class of equations.

Analytical Matrix Method
One procedure to find periodic solutions of nonlinear systems is the method of harmonic balance [15]. This method assumes that the solution can be represented by a linear combination of sinusoids (a Fourier series) of the form: where ω is the main frequency with harmonics of the form kω, for k ∈ Z. The problem reduces to find the frequency and Fourier coefficients a k . In the present work, an analogous procedure is proposed, namely, it is assumed that each element of the invariant submanifold, π(z), that satisfies Equation (9), has the following form where c k,i,j are undetermined coefficients. Then, the problem reduces to find the correct value of these coefficients, that can be found using the quadrature method, where the proposed solution (16) is substituted in Equation (9) to reduce the problem to find coefficients c k,i,j . This solution can be rewritten in matrix form as and Now it is possible to state the problem to be solved in the present work.
Problem 1. Given the tested system (6) submitted to a oscillatory regimen by the driving system (7) and (8), after the transient period, system (6) reaches a invariant submanifold, x(t) = ξ(z(t)), described by the invariance partial differential equation Then, the problem to be solved is to find an analytical series solution with the form (17) that satisfies Equation (20).
Notice that instead of using a harmonics series, such as Equation (15), in the proposed series solution (16), the driving variables, z 1 and z 2 , are elevated to integer powers and the use of the matrix form (17) allows simplifying the analysis and solution in a more systematic way. The following two definitions are instrumental for developing the proposed solution method. Definition 1. Given two square matrices Φ a and Φ b with the same dimension N × N (N could tend to infinity), the product operator is defined as where the elements of matrix P are while E is a vector of dimension N × 1, whose elements are where the elements of matrix D are while the elements of matrix P are defined in Equation (22).
To substitute the proposed solution (16) in Equation (9), it is necessary to compute the time derivative of π; given that π is the sum of coefficients of the form and considering the dynamics of z, given in Equations (7) and (8), it holds that Notice that the sum of exponents of z i 1 z j 2 is i + j and, given the linear dynamics of z, the sum of exponents of d dt z i 1 z j 2 is preserved and equal to i + j. Following this structure, the time derivative of the infinite series (16), or its matrix form (17), is presented in the following Proposition. Proposition 1. Consider the square matrix Φ a of infinite dimension with constant coefficients, associated to a function of z 1 and z 2 defined as where dΦ a is the operator described in Definition 2.
Proof. Notice that the derivative of vector Z i with respect to z i is is defined in Equation (25). Therefore, the partial derivatives of φ a with respect to z 1 and z 2 are On the other hand, notice that where is defined in Equation (22); therefore, the time derivative of φ a is Finally, with Definition 2, it holds thatφ a = ωZ T 1 (dΦ a )Z 2 , concluding the proof.
Equation (6) contains nonlinear term; therefore, it might have products of finite or infinite series of z 1 and z 2 . For instance, the product z a while the procedure to compute the product of two series is presented in the following proposition.

Proposition 2.
Consider the square matrices Φ a and Φ b with infinite dimension, associated to functions φ a and φ b of z 1 and z 2 defined as the product of these functions is which is equivalent to where the operator Φ a , Φ b N is described in Definition 1.
and considering Equation (34), the previous equation is equivalent to Notice that, using matrices E k , defined in Equation (23), Z 2 can be expressed as and finally, using Equation (34) again, it holds that as described in Equation (21), in order to rewrite Equation (43) as concluding the proof.
Finally, the following theorem presents the conditions to find the oscillatory invariant submanifolds. Theorem 1. Given the input-state stable differential Equation (6), with x ∈ X ∈ R n and z ∈ R 2 , and the oscillatory dynamics defined in Equations (7) and (8), assume that x(t) reaches an invariant submanifold that satisfies Equation (20).
Then, if function f (x(t), z(t)) can be expressed as series of integer powers of the elements of x and z, the solution of each element of vector π can be proposed to be an infinite series of integer powers of z 1 and z 2 , i.e., where matrices C k are Proof. The system of differential equations defined in Equation (20) is Assume that the elements of the invariant submanifold π(t) can be described as in Equation (45), then, according to Proposition 1, their time derivatives are where dC k = D T C k P − P T C k D. Using this time derivative, the elements of Equation (20) become for k = 1, 2, . . . , n. If functions f k (x, z) can be expressed as series of integer powers of the elements of x and z, then using identities given in Equations (37) and (42) in Proposition 2, it can be shown that functions f k (x, z) are equivalent to . . , C n )Z 2 and the differential equations can be transformed to the form where F k (C 1 , . . . , C n ) is a matrix that depends on the values of coefficient matrices C 1 , . . . , C n , defined in Equation (46). Independently of the values of vectors Z 1 and Z 2 , Equation (47) holds if ωdC k − F k (C 1 , . . . , C n ) = 0, k = 1, 2, . . . , n, then the problem reduces to finding coefficients c k,i,j , k = 1, 2, . . . , n, i, j ∈ N, by solving Equations (48) to obtain the solution given in Equation (45).

is in Medium Amplitude
Oscillations mode (MAOs) and additional terms of higher-order in solution (49) must be included to describe its oscillatory response correctly.
LAOs Finally, if the series is not convergent, i.e., c k,n ≯ A c k,n+1 as n increases, then system (9) is in Large Amplitude Oscillations mode (LAOs) and solution (49) is not longer suitable to describe the behavior of system (9).
Since the current contribution is focused on the analysis of SAOs and MAOs modes, the LAOs analysis will be left to a future contribution.

Series Solutions for a Simple Globally Asymptotically Stable Model
To illustrate the proposed method, consider the following dynamical system: where x ∈ R is the state variable and u ∈ R is the input variable. Notice that when u = 0, the steady-state is given by x = 0, and using the Lyapunov function V = x 2 , it is straightforward to verify that this equilibrium is globally asymptotically stable, be-causeV = −2x 4 = −2V 2 < 0 for V = 0. Then, if system (50) is submitted to an oscillatory input: u(t) = z 1 (t), where z 1 is the first element of the state vector of the oscillator (7)- (8), after a transient behavior that depends on the initial conditions, the system will approach an invariant submanifold, x(t) = ξ(z(t)), that satisfies the following nonlinear partial differential equation: To find the analytical solution, it is considered that π can be expressed as a series of the form of Equation (17), i.e., with unknown coefficients matrix C. The time derivative of the invariant submanifold iṡ where dC = D T CP − P T CD, while the first term of the right-hand side of Equation (51) is where Q 3 = C 1 , C 1 , C 1 ∞ ∞ , and the second term of the right-hand side of Equation (51) is z 1 = Z T 1 A 1,0 Z 2 ; therefore, Equation (51) becomes which is satisfied for any oscillatory behavior iff Equation (53) allows to compute the matrix of coefficients C, and following the procedure described in Appendix A, the solution of Equation (51) is Notice that each term in Equation (54) has the form z i 1 z j 2 , satisfying the condition that the total order n = i + j is an odd number. System (50) is in SAOs mode when the terms with order two or bigger are negligible, i.e., ξ(t) ≈ − 1 ω z 2 (t), while this system is in MAOs mode when additional terms of higher order must be included to correctly describe its oscillatory response. Given the coefficients of the analytical solution (54), SAOs occurs when A 2 ω 3 . Figure 2 shows the predictions with a numeric method solution and the analytical solution (54) truncated up to first-order, , and seventh-order (see Equation (54)).

Series Solutions for a Simple Pendulum
Now, consider the pendulum model presented in Section 3. The oscillatory invariant submanifold for this pendulum satisfies Equations (13) and (14). Now the solution is proposed to have the form of Equation (17), i.e., where Z i for i = 1, 2 is defined in Equation (18). Given that sin(π 1 ) = ∑ ∞ j=0 (−1) j π 2j+1 1 (2j+1)! , using Proposition 2, the sine function in Equation (14) is where with W 0 = C 1 , W j = W j−1 , W C , j = 1, 2, 3, . . . and W C = C 1 , C 1 . According to Proposition 1,ξ i = ωZ T 1 (dC i )Z 2 , while the term z 1 is equivalent to Z T 1 A 1,0 Z 2 ; therefore, Equations (13) and (14) can be expressed as Independently of the values of vectors Z 1 and Z 2 , Equation (13) holds if while Equation (14) holds if where dC i for i = 1, 2 are the differential operators described in Proposition 1. Equations (58) and (59) are a particular case of Equation (48). The coefficients C 1 and C 2 are obtained by solving the matrix Equations (58) and (59), as shown in Appendix B. Thus, the behavior of the pendulum angular position and angular velocity are where while the remaining coefficients can be easily obtained with Equation (A6). According to Theorem 1, the invariant submanifold will be only reached if the system is input-state stable, and given that the pendulum has two classes of equilibrium points: (1) x 1 = ±(2n + 1)π (up) and 2) x 1 = ±2nπ (down) for n = 0, 1, 2, . . .. Therefore, the family of "up" fixed points are unstable, while the "down" fixed points are stable; therefore, the correct solutions of Equations (60) and (61) are those where n is equal to zero or an even number. This analytical solution is compared with a numerical solution for two different cases in Figures 3 and 4. Figures 3a and 4a show the time evolution of the angular position, Figures 3b and 4b show the time evolution of the angular velocity, and Figures 3c and 4c show the phase plane with respect to the oscillatory driving signal.

Discussion
For the study case of the series solutions for a simple globally asymptotically stable model, Figure 2 shows that the analytical solutions truncated up to the first-order and third-order are unable to reproduce the behavior, while the analytical solution truncated up to seventh-order satisfactorily describes the system behavior. Notice that, by looking in Equation (54) at the coefficients of the series and considering that z 1 = A sin(ωt + φ) and z 2 = A cos(ωt + φ), it holds that the first-order term (n = 1) is proportional to A/ω, while the third, fifth and seventh-order terms (n = 3, 4, 5) are proportional to A 3 /ω 4 , A 5 /ω 7 , and A 7 /ω 10 , respectively. Following this sequence, the terms of order n = 2k − 1, for k = 1, 2, 3, . . ., are proportional to A 2k−1 /ω 3k−2 . Thus, the convergence of series (54) is correlated to the convergence of coefficient A 2 /ω 3 k ω 2 /A with respect to k, then, without considering the constant coefficients in series (54), an approximated limit to guarantee the convergence of series is that the relation frequency-amplitude must satisfy ω 3 > A 2 ; if this relation is satisfied then we can say that the system is in SAOs mode when A 2 /ω 3 1 or MAOs mode when ω 3 > A 2 , respectively. In Figure 2, ω = 0.5 and A = 0.2 were selected, which satisfies the aforementioned relation, while Figure 5 shows the predictions for the same frequency but with an amplitude A = 0.4, which violates the relation ω 3 > A 2 , and therefore, the analytical solution Equation (54) is no longer suitable for describing the dynamical behavior because, in this context, the system is in LAOs mode. To make the analysis in LAOs mode, it would be necessary to reconsider the assumption C 2,0 = 0 made to solve Equation (A1), which actually has an infinite number of solutions of the form C 0,2 = C 2,0 ; therefore, if C 2,0 = 0, it would appear a term of the form C 2,0 A 2 sin 2 (ωt + φ) + cos 2 (ωt + φ) that is equivalent to C 2,0 A 2 . Therefore, a new constant term would appear, modifying the equation (i, j) = (0, 0) of system (53), i.e., instead of 0 = C 3 0,0 , it would be 0 = C 3 0,0 + C 2,0 A 2 , and now C 0,0 could be different to zero and a function of the amplitude. The complete analysis of LAOs mode will be left to a future contribution. Notice that for the well-known problem of a simple pendulum that oscillates under sinusoidal potentials, the resonance characteristics are difficult to find by solving the problem analytically; hence, the numerical approach is the method of choice [16]. However, we were able to find the analytical solution for the study case of the series solutions for a simple pendulum. We found that when the system behaves almost as in the SAOs regime (Figure 3), there may be deviations from the linear behavior and additional third-order terms must be included to accurately reproduce the system's behavior. Figure 4 depicts the system in a MAOs regime with more pronounced nonlinear characteristics, where coefficients up to eleventh order must be included to correctly reproduce the system's behavior. Figure 6 shows that for ω = 0, z 1 becomes constant and the steady-state solution of Equations (11) and (12) is given by the equations x 2 = 0 and 0 = z 1 − sin(x 1 ); therefore, x 1 = ±2nπ + arcsin(z 1 ), where n may be equal to n = 0, 1, 2, . . .. However this solution is only feasible if z 1 ≤ 1, because given the domain of sine function, if z 1 > 1, this class of steady-state is not longer feasible and the pendulum reaches a permanent regime where |x 1 | keeps growing, while x 2 oscillates around a constant value, i.e., the pendulum rotates persistently. Consequently, solution (60) is not longer valid to describe this behavior, because according to Theorem 1, the system described in Equation (12) is not longer inputstate stable. An example of this behavior is shown Figure 6. In this case, the system is not longer in a SAOs or MAOs regimen. To analyze the approximated regions in which SAOs, MAOs, or LAOs hold for the angular position of the pendulum, it is possible to use Definition 3 in the coefficientŝ C 1,k , given in Equation (A6). Figure 7 sketches these regions as a function of the applied frequency and amplitude. In particular, the boundary between SAOs and MAOs regions was selected to have a maximum 1% variation with respect to the linear behavior, i.e., where A = 0.01 c 1,1 / c 1,2 and c 1,n > A c 1,n+1 for n = 1, 2, 3, . . ., while the boundary between SAOs and MAOs regions approximates the points at which solution (60) is not longer suitable to describe the oscillatory behavior because c 1,n ≯ A c 1,n+1 for some n = 1, 2, 3, . . .. As can be seen, for low frequencies and for amplitudes approximately bigger than 1, the pendulum is in LAOs mode, while for high frequencies the SAOs mode predominates regardless of the amplitude. Particularly for frequencies around 1, the boundaries between the three regions appears at smaller amplitudes, because this frequency is the natural frequency of the pendulum (given the numerical value of the damping parameter, the eigenvalues of the Jacobian around equilibrium are approximately 0.0714 ± 0.9974i), and the driving force is producing larger resonances. Another resonance region is approximately at frequencies between 0.1 and 0.4, associated with the damping constant η. Finally, notice that the diagram in Figure 7 has some similarities with the typical amplitude ratio Bode plot for a second-order damped system, which for comparison purposes was superimposed on the diagram (yellow line), suggesting that this class of diagram may be a valuable tool for frequency response techniques of nonlinear systems.

Conclusions
The analytic method for FRT with applications to nonlinear dynamical systems presented here is an improvement over the series solution methods available because it offers the possibility of obtaining the analytical solution. This method provides a structured procedure by working with the invariant submanifold of the problem that is equivalent to the harmonic balance series solution. The recursive relation obtained for the coefficients in our analytical method simplifies traditional approaches to obtain the solution with the harmonic balance series method. This proved to be useful to obtain information of nonlinear systems subjected to periodic inputs for the cases where the output is stable and others where the system is not at equilibrium. Furthermore, this method is useful to identify the modes of behavior SAOs and MAOs that arise in several FRT applications; the criteria to explore LAOs will be presented in a second paper of the series.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/math9243287/s1, the data presented in this study were generated using a Python code that is available in the supplementary material.  The following calculations were done with openware FRAMSpy module generated by the authors with Sympy package in Python 3. The codes can be consulted in the supplementary material section.

Appendix B. Detailed Solution of the Invariant Submanifold of Study Case 2
The following calculations were done with openware FRAMSpy module generated by the authors with Sympy package in Python 3. The codes can be consulted in the supplementary material section.