Efficient, universal and easy-to-use solution method for periodic response of high-dimensional complex nonlinear differential equation

Abstract

The invention relates to a high-dimensional complex nonlinear differential equation periodic response efficient universal easy-to-use solving method, which is characterized in that an automatic differential technology and a harmonic balance method are fused, and the method comprises the following steps: step 1, solution setting: representing the periodic response of a system by using discrete Fourier series; step 2, harmonic balance: substituting periodic response in a Fourier series form of the system back to a system kinetic equation, and converting a differential equation into an algebraic equation through a harmonic balance program; 3, calculating a Jacobian matrix, namely, quickly obtaining an accurate Jacobian matrix of the algebraic equation after harmonic balance by utilizing an automatic differential technology; and 4, solving an algebraic equation, namely calculating by using a Newton-Raphson method to obtain a system periodic response. The method solves the problems that a high-dimensional complex nonlinear differential equation response solving method is low in efficiency, not universal and not easy to use, can be widely applied to calculation of high-dimensional complex nonlinear system dynamic response, and has the advantages of being efficient, universal, capable of being used immediately after being opened and modularized.

Description

Technical Field

[0001] This invention belongs to the field of dynamic simulation algorithms. It is an efficient, general and easy-to-use solution method for the periodic response of high-dimensional complex nonlinear differential equations. It integrates automatic differentiation and harmonic balance methods and utilizes the parallel computing and CUDA acceleration capabilities of deep learning frameworks to construct an efficient, general and out-of-the-box solution method suitable for high-dimensional complex nonlinear systems. It is applicable to the rapid dynamic simulation of complex systems in engineering, such as aero-engines, gas turbines, large spacecraft structures, and power grid systems. Background Technology

[0002] The evolution of dynamic systems is usually described by second-order ordinary differential equations. In order to accurately characterize the dynamic characteristics of the system, the actual model not only needs high degrees of freedom, but also needs to introduce nonlinear factors, such as nonlinear stiffness, nonlinear damping and parametric excitation. Therefore, the dynamic model has the characteristics of high dimensionality and nonlinear complexity, making it difficult to obtain its steady-state periodic response.

[0003] Methods for solving nonlinear ordinary differential equations mainly fall into three categories: numerical methods, analytical methods, and semi-analytical methods. Numerical methods yield discrete and isolated solutions, failing to capture the full picture of the solution space by providing unstable periodic responses. Furthermore, numerical methods exhibit poor stability, making them unsuitable for rigid equations. Implicit numerical methods suffer from low computational efficiency, hindering large-scale parameter influence analysis. Analytical methods lack versatility, proving difficult to apply to high-dimensional systems or handle complex nonlinear factors (such as fractional exponential functions). Semi-analytical methods, such as the Newton-Raphson iteration, rely on the Jacobian matrix, which remains within the realm of symbolic or numerical differentiation. Symbolic differentiation suffers from sign inflation and requires extensive manual derivation, making it unsuitable for high-dimensional and complex nonlinear problems. Numerical differentiation suffers from truncation and rounding errors, making accurate calculation of bifurcation and instability issues in complex nonlinear systems difficult.

[0004] Currently, there is a lack of efficient, universal, and out-of-the-box solution methods for high-dimensional nonlinear differential equations of complex dynamic systems. This limits the analysis of the dynamic characteristics of complex systems to linearization, numerical calculations, or low-dimensional simplified models, thus restricting further optimization and improvement of the dynamic performance of complex systems. Therefore, developing efficient, universal, and out-of-the-box solution methods for the periodic responses of high-dimensional complex nonlinear systems has significant theoretical value and broad engineering application prospects. Summary of the Invention

[0005] This invention provides an efficient, general, and easy-to-use solution method for the periodic response of high-dimensional complex nonlinear differential equations. It constructs a basic solution framework through harmonic balance-Newton iteration, uses automatic differentiation to obtain the Jacobian matrix required for Newton iteration, and relies on the parallel computing and CUDA acceleration capabilities of deep learning frameworks to completely eliminate the manual derivation of complex formulas in the solution process, thus constructing an efficient, general, and out-of-the-box program architecture.

[0006] This invention provides an efficient, general, and easy-to-use method for solving periodic responses of high-dimensional complex nonlinear differential equations, comprising the following steps:

[0007] Step 1: Solution—Express the system's periodic response using a discrete Fourier series. Specifically, for a general nonlinear ordinary differential equation, express the system's periodic response X using a discrete Fourier series, i.e.: X = AH, where H is the Fourier series harmonic basis matrix, and A is the corresponding Fourier coefficient matrix, as shown below:

[0008] (4)

[0009] (5)

[0010] Where s represents the order of the Fourier series, the value of which depends on the nature of the external excitations and nonlinear terms of the system.

[0011] Step 2: Harmonic Balancing—Substitute the system's periodic response in Fourier series form back into the system's dynamic equations, and transform the differential equations into algebraic equations using a harmonic balancing procedure. Specifically, this includes: substituting the system's periodic response in Fourier series form back into the system's dynamic equations, performing a fast Fourier transform on the dynamic equations, discretizing the residuals of the dynamic equations into Fourier series, extracting the coefficients of each harmonic pair, obtaining the algebraic residuals of the dynamic equations, and transforming the system's motion differential equations into algebraic equations. Simultaneously, using system parameters as extension parameters (e.g., rotor speed as an extension parameter when calculating the amplitude-frequency response), add corresponding algebraic constraint equations, and assemble these constraint equations into the residuals of the system's algebraic equations, forming a global residual square algebraic equation with arc length extension.

[0012] Step 3: Calculate the Jacobian matrix—Use automatic differentiation technology to quickly obtain the accurate Jacobian matrix of the harmonic-balanced algebraic equations. Specifically, this includes: using an automatic differentiation system to calculate the accurate Jacobian matrix of the residual algebraic equations, completely eliminating the manual derivation process of complex formula differentiation, overcoming the problem of cumbersome and even difficult-to-calculate Jacobian matrix derivation for complex nonlinear systems, and breaking the inherent defects of symbolic differentiation (expression inflation) and numerical differentiation (truncation and rounding errors), thus improving the method's versatility and ease of use. Furthermore, based on advanced deep learning frameworks and CUDA support, a better program architecture is designed to directly call the inverse automatic differentiation function (jacrev) to efficiently calculate the Jacobian matrix, improving computational efficiency and achieving out-of-the-box usability.

[0013] Step 4: Solving the Algebraic Equation—Iterative calculations are performed using the Newton-Raphson iterative method until convergence, obtaining the Fourier series coefficients of the periodic response. These coefficients are then multiplied by the Fourier series trigonometric harmonic basis in the solution to output the system's periodic response. Specifically, this includes iterative calculations using the Newton-Raphson iterative method to update the Fourier coefficients of the system's periodic response, i.e.:

[0014] (6)

[0015] When the norm of residual B is less than the given error That is, satisfying When the iteration converges, the periodic response of the system is obtained as X = AH. For the evolution law of the system response in the parameter space, the parameters can be specified as extension parameters. For example, when calculating the amplitude frequency response, the rotor speed is used as the extension parameter. The system parameters are updated by using the arc length extension method. By repeating steps (1)-(4), the response of the system in the entire parameter space can be obtained.

[0016] Beneficial effects

[0017] (1) High efficiency: This invention benefits from the deep optimization of parallel computing and CUDA by advanced deep learning frameworks such as PyTorch and JAX. It directly uses automatic differentiation functions, such as jacrev, to achieve efficient computation of the Jacobian matrix. In this process, vmap is used to vectorize the function, so the computation of the Jacobian matrix no longer needs to be looped, which greatly improves the computation efficiency. Furthermore, it can achieve efficient parallel computing on GPUs, showing significant advantages when dealing with large-scale problems. Especially in engineering practice, when large-scale parameter influence analysis is required, its parallel architecture can significantly improve computation efficiency and greatly shorten the analysis cycle.

[0018] (2) Universality: The automatic differentiation in this invention takes the chain rule as the mathematical core and the computation graph as the engineering framework. It provides two complementary modes, forward and backward, to adapt to different computation scenarios. It can accurately and efficiently handle derivative calculation problems from simple functions to complex neural networks. It is applicable to all nonlinear problems and realizes the universality of this invention.

[0019] (3) Ease of use: This invention utilizes automatic differentiation to automatically calculate the accurate Jacobian matrix, completely eliminating the manual derivation process of complex formulas, overcoming the problem that the derivation of Jacobian matrix for complex nonlinear systems is complicated or even difficult to calculate, and relying on the deep learning framework, it can directly call the automatic differentiation function to calculate the Jacobian matrix, achieving out-of-the-box use;

[0020] (4) Global Analyticalness: This invention is based on the framework of the harmonic balance method, which is essentially an analytical method. Automatic differentiation can obtain the accurate Jacobian matrix required for Newton iteration. Combined with the arc length extension method, it can obtain all periodic response solutions of the system, including unstable periodic solutions, laying the foundation for in-depth analysis of nonlinear dynamic mechanisms. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the process of the present invention.

[0022] Figure 2 This is a schematic diagram of the squeeze oil film damper-rotor system model.

[0023] Figure 3 This is a comparison chart of the solution results of the method proposed in this invention and the fourth-order Runge-Kutta (RK4). Detailed Implementation

[0024] The embodiments of this invention will be described in detail below with reference to the accompanying drawings and specific examples. The technical solutions will be clearly and thoroughly described. Obviously, the described embodiments are only individual embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art without making innovative achievements are within the scope of protection of this invention.

[0025] This invention is an efficient, general, and easy-to-use solution method for the periodic response of high-dimensional complex nonlinear differential equations. Its implementation process is as follows: Figure 1 As shown. Now, as follows Figure 2 Taking the squeeze film damper-rotor system shown as an example, the implementation steps of the present invention are discussed in detail. Before solving, the motion differential equations of the squeeze film damper-rotor system are first given as follows:

[0026] (7)

[0027] Where m represents the rotor mass, c represents the damping coefficient at the support, k represents the support stiffness, and J d J represents the rotor's diameter and moment of inertia. p Let l represent the polar moment of inertia of the rotor, l1 represent the length of the shaft segment from the left support to the turntable, l2 represent the length of the shaft segment from the right support to the turntable, and e represent the rotor eccentricity. The value represents the rotor speed, and x and y represent the rotor's horizontal and vertical displacement responses, respectively. , This indicates the rotor's speed response in the horizontal and vertical directions. and This represents the rotor's acceleration response in the horizontal and vertical directions. and This indicates the rotor's rotation angle in the x and y directions. and This represents the rotor's angular velocity response in the x and y directions. and This represents the rotor's angular acceleration response in the x and y directions. and These represent the components of the oil film force in the x and y directions of the extrusion oil film damper, respectively. The calculation method for the nonlinear restoring force of the extrusion oil film damper is as follows:

[0028] (8)

[0029] in, Let ψ represent the dimensionless radial displacement of the journal center, ψ be the precession angle of the journal, and I be the Sommerfeld integral constant. Furthermore, we have:

[0030] (9)

[0031] (10)

[0032] (11)

[0033] Where θ1 and θ2 represent the start and end points of the positive pressure zone of the oil film, and we have:

[0034] (12)

[0035] Transforming the oil film force shown in formula (5) into a rectangular coordinate system, we can obtain the expressions for the oil film force in the x and y directions as follows:

[0036] (13)

[0037] As can be seen from equations (8)–(13), the oil film force of the extrusion oil film damper exhibits complex nonlinear characteristics; it is (x, y, , The nonlinear functions of the system and their derivatives, including fractional exponential functions, negative exponential functions, and inverse trigonometric functions, as well as integral functions, are coupled together. When using implicit numerical or analytical methods based on Newton iteration to solve the periodic response of the system, the derivation of the Jacobian matrix required for Newton iteration is not only cumbersome but also error-prone. This makes it almost impossible to use these methods to solve the dynamic equations of the squeeze film damper-rotor system to obtain the periodic response. The method proposed in this invention is efficient, universal, and ready to use, capable of solving such complex nonlinear differential equations. The specific implementation process is described below:

[0038] Step 1: Solution—Express the periodic response of the differential equation shown in formula (1) using discrete Fourier series, as follows:

[0039] Let X = [x, y, , ] T Let X = AH, where H is the Fourier series harmonic basis matrix, A is the corresponding Fourier coefficient matrix, and we have:

[0040] (14)

[0041] (15)

[0042] Step 2: Harmonic Balancing—Substitute the system's periodic response in Fourier series form back into the system's dynamic equations, and use the harmonic balancing procedure to transform the differential equations into algebraic equations, as follows:

[0043] According to X = AH, we have , Substitute back into the system dynamic equation (4) and perform a fast Fourier transform to discretize the residuals of the dynamic equations into a Fourier series. Extract the coefficients of each harmonic pair to obtain the algebraic form of the residuals of the dynamic equations. Transform the system's differential equations of motion into algebraic equations. Meanwhile, assume the rotor speed... To extend the parameters, corresponding algebraic constraint equations are added, and these constraint equations are assembled into the residuals of the system algebraic equations to form the global residual square algebraic equations with arc length extension, expressed as follows:

[0044] B(A) = 0(16)

[0045] Where B is the Fourier coefficient of the residuals of the system dynamics equations.

[0046] Step 3: Calculate the Jacobian matrix—Use automatic differentiation techniques to quickly obtain the exact Jacobian matrix of the algebraic equations after harmonic balance, as detailed below:

[0047] The completely open-source inverse automatic differentiation function jacrev provided by PyTorch is directly called. The algebraic residual B(A) of the system's dynamic equation is input, and inverse automatic differentiation is performed to obtain the Jacobian matrix of the residual relative to the Fourier series A of the system's periodic response.

[0048] Step 4: Solving the algebraic equation—Use the Newton-Raphson iterative method to perform iterative calculations until the iteration converges, obtain the Fourier series coefficients of the periodic response, and take the inner product with the Fourier series trigonometric function harmonic basis in the solution to output the periodic response of the system, as follows:

[0049] The Newton-Raphson iterative method is used to iteratively calculate and update the Fourier coefficients of the system's periodic response, i.e.:

[0050] (17)

[0051] When the norm of residual B is less than the given error That is, satisfying When the iteration converges, the periodic response of the system is obtained as X = AH. For the amplitude-frequency response of the system, the rotor speed is taken as... To obtain the system's amplitude-frequency response, the rotor speed is updated using the arc length extension method, and steps (1)–(4) are repeated.

[0052] The periodic response results of the squeeze film damper-rotor system obtained according to the solution example of this invention are plotted in... Figure 3 The solution results are consistent with the fourth-order Runge-Kutta (RK4) method. This invention can achieve accurate solution of complex nonlinear dynamic equations of the extrusion oil film damper-rotor system.

[0053] The above description of the present invention is merely one embodiment of the present invention and is only used to illustrate the implementation steps of the present invention. It is not intended to limit the implementation of the present invention. Those skilled in the art can easily make corresponding modifications or alterations based on the main concept and spirit of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of protection claimed in the claims.

Claims

1. A highly efficient, general, and easy-to-use solution method for the periodic response of high-dimensional complex nonlinear differential equations, characterized in that... The method integrates automatic differentiation technology with harmonic balance, including the following steps: Step 1: Solution—Express the periodic response of the system using discrete Fourier series; Step 2: Harmonic balance—Substitute the system's periodic response in Fourier series form back into the system's dynamic equations, and use the harmonic balance procedure to transform the differential equations into algebraic equations; Step 3: Calculate the Jacobian matrix—Use automatic differentiation techniques to quickly obtain the exact Jacobian matrix of the algebraic equations after harmonic balance; Step 4: Solving the algebraic equation—Use the Newton-Raphson iterative method to perform iterative calculations until the iteration converges, obtain the Fourier series coefficients of the periodic response, and perform an inner product with the Fourier series trigonometric function harmonic basis in the solution to output the periodic response of the system.

2. The efficient, general, and easy-to-use solution method for periodic responses of high-dimensional complex nonlinear differential equations according to claim 1, characterized in that: For a general nonlinear ordinary differential equation, the periodic response X of the system can be represented by a discrete Fourier series, i.e.: X = AH, where H is the Fourier series harmonic basis matrix and A is the corresponding Fourier coefficient matrix, as shown below: (1) (2) Where s represents the order of the Fourier series, the value of which depends on the nature of the external excitations and nonlinear terms of the system.

3. The efficient, general, and easy-to-use solution method for periodic responses of high-dimensional complex nonlinear differential equations according to claim 1, characterized in that: Substituting the system's periodic response in Fourier series form back into the system's dynamic equations, and performing a fast Fourier transform on the dynamic equations, the residuals of the dynamic equations are discretized into Fourier series. The coefficients of each harmonic pair are extracted to obtain the algebraic form of the residuals of the dynamic equations, thus transforming the system's differential equations of motion into algebraic equations. Simultaneously, using system parameters as extension parameters (e.g., rotor speed as an extension parameter when calculating the amplitude-frequency response), corresponding algebraic constraint equations are added, and these constraint equations are assembled into the residuals of the system's algebraic equations, forming a global residual square algebraic equation with arc length extension.

4. The efficient, general, and easy-to-use solution method for periodic responses of high-dimensional complex nonlinear differential equations according to claim 1, characterized in that: By utilizing the accurate Jacobian matrix of residual algebraic equations obtained through an automatic differential calculation system, the manual derivation process for complex formulas is completely eliminated. This overcomes the problem of cumbersome and even computationally difficult derivation of Jacobian matrices for complex nonlinear systems, and breaks through the inherent defects of symbolic differentiation (expression inflation) and numerical differentiation (truncation and rounding errors), improving the method's versatility and ease of use. Furthermore, based on advanced deep learning frameworks and CUDA support, a better program architecture is designed to directly call the inverse automatic differential function (jacrev) for efficient Jacobian matrix calculation, improving computational efficiency and achieving out-of-the-box usability.

5. The efficient, general, and easy-to-use solution method for periodic responses of high-dimensional complex nonlinear differential equations according to claim 1, characterized in that: The Newton-Raphson iterative method is used to iteratively calculate and update the Fourier coefficients of the system's periodic response, i.e.: (3) When the norm of residual B is less than the given error That is, satisfying When the iteration converges, the periodic response of the system is obtained as X = AH. For the evolution law of the system response in the parameter space, the parameters can be specified as extension parameters. For example, when calculating the amplitude frequency response, the rotor speed is used as the extension parameter. The system parameters are updated by using the arc length extension method. By repeating steps (1)-(4), the response of the system in the entire parameter space can be obtained.

Images