Abstract
An L-stable method over time intervals for differential-algebraic equations (DAEs) of multibody system dynamics is presented in this paper. The solution format is established based on equidistant nodes and nonequidistant nodes such as Chebyshev nodes and Legendre nodes. Based on Ehle’s theorem and conjecture, the unknown matrix and vector in the L-stable solution formula are obtained by comparison with Pade approximation. Newton iteration method is used during the solution process. Taking the planar two-link manipulator system as an example, the results of L-stable method presented are compared for different number of nodes in the time interval and the step size in the simulation, and also compared with the classic Runge-Kutta method, A-stable method, Radau IA, Radau IIA, and Lobatto IIIC methods. The results show that the method has the advantages of good stability and high precision and is suitable for multibody system dynamics simulation under long-term conditions.
1. Introduction
The stability of the numerical solution of differential equations is one of the important indicators to judge the quality of a numerical method. In the 1960s, after the innovation and development of Dahlquist, Widlund, Gear, and other scholars, the stability theory was deeply developed. Dahlquist [1] proposed a linear equation as a model equation and proposed A-stable, that is, , , and all . Subsequently, Axellson [2] determines the of the by the quadrature coefficients at the zeros of or , constructing a class of A-stable methods, where is the Legendre polynomial. Widlund [3] proposed A(α)-stable, that is, if all solutions of tend to zero as when the method is applied, with a fixed , to any differential equation of the form ,where is a complex constant which lies in the set . Gear [4] proposed stiffness stable, in order to generalize the A-stable concept to nonlinear problems. In 1975, Dahlquist [5] proposed the G-stable concept of single-method linear multistep method. Butcher [6, 7] studied the nonlinear stability of Runge-Kutta method by using nonlinear equations as model equations, proposed B-stable and algebraic stability of the Runge-Kutta method, and established the algebraic stability criteria for general linear methods. Li Shoufu [8] established the (k,p,q) stability theory of the general linear method and the weak algebraic stability criterion of (k,p,q). For a very rigid differential/algebraic equations, it is hoped that the rigid component can be attenuated as soon as possible, and the stable domain contains the left half complex plane; only A-stable is not enough. Ehle [9] proposed L-stable. Hairer [10] pointed out that the L-stable numerical method will introduce damping effect, but it is very effective for suppressing high-frequency oscillation, and some scholars have proposed implicit single-step and multistep methods based on L-stable for solving rigid equations, but these methods are based on ordinary differential equations. L-stable methods for DAEs are relatively rare.
The typical form of multibody system dynamic equations is DAEs, which consist of ordinary differential equations and algebraic constraint equations. The numerical solution is a matter of concern in the study of multibody system dynamics. At present, most of the numerical solutions of the DAEs mainly include Direct integration method and Newmark method. Although the calculation process is relatively simple, the stability is poor. In the long-term simulation process, the error accumulation is more serious, and it needs to be corrected by combining constrained projection [11]. In this paper, the general form of L-stable method for multibody system dynamics DAEs is established. Based on the equidistant nodes, Chebyshev nodes, and Legendre nodes, the concrete construction process is given and the planar two-link manipulator system verified and compared.
2. DAEs of Multibody System Dynamics
The DAEs of index-3 in multibody systems dynamics arewhere are generalized coordinates of the system, are generalized speed of the system, are generalized acceleration of the system, are the Lagrange multiplier vectors, is a symmetric positive definite mass matrix, are generalized force vectors, and are complete displacement constraints.
For the constraint equation , the first derivative and the second derivative are, respectively, obtained for the time t, and the DAEs form of the index-2 and the index-1 can be obtained by the combination of (1).When the initial values are known, the Runge-Kutta method can be used to directly solve the index-1 equation (3); however, if the selected step size is large, the accuracy will be reduced and the total energy will increase with time [12], smaller step will lead to lower computation efficiency and poor stability, and the index-3 equation (1) and the index-2 equation (2) cannot be solved directly using the Runge-Kutta method. The L-stable method is constructed below and can be used to solve the multibody system dynamics DAEs of index-1, -2, and -3.
3. Construction of the DAEs L-Stable Format
3.1. General Form of L-Stable Method
In 1972, H. A. Watts and L. F. Shampine [13] proposed a block implicit one-step method for solving ordinary differential equations. The r discrete moments are made into a block vector, the output values of r moments are obtained from the input values of one moment, and it is proved that the method is A-stable. In this method, the time step in each block vector is equidistant. This article improves on this method; the block vector is modified to correspond to the point in the single-step interval, and in the iterative solution process, a single-step interval and a time step are set, selecting and constructing DAEs solution format, where can be an equidistant node or a nonequidistant node.
Taking the multibody system dynamics DAEs of index-1 equation (3) as an example, after the order reduction, the following first-order differential equations can be obtained. Let , , ; the construction solution formula is as followswhere , , , , and is Kronecker.
Bring (5) into the test equationand get stability functionwhere , is the r-order unit matrix, and . For the stability function, using the Pade approximation, we can construct L-stabilized solution formula.
Theorem 1 ( theorem and conjecture). For approximation, if , then this Pade approximation is A-stable; if , then this Pade approximation is L-stable [14].
Perform a Taylor expansion on at , let the coefficient of the term be a pending constant, and omit the following high-order terms; one can getwhere , is pending constant and . Bringing (8) into (5), one can get the following.From (9), , which are expressed using the undetermined constant, can be obtained bringing these into the stability function (7). The value of the undetermined constant , can be obtained by comparison with the corresponding Pade approximation form, so that the value of can be obtained bringing into formula (5) a specific solution formula.
3.2. L-Stable Format Based on Equidistant Nodes
For the equidistant node , there are . For example, when , the equidistant node is ; at this time , and (9) becomes as follows.Calculate , substitute stability function (7), and compare it with approximation formula; one can solve undetermined constants , and bringing them into the expression , one can obtainBringing them into (5) can be used to obtain a solution formula. The same method can be used to construct L-stable solution formats for other r values.
3.3. L-Stable Format Based on Nonequidistant Nodes
When the node of the interval is the zero point of the Chebyshev polynomial, within the interval , one gets the following.Take , as an example, get , ,, and bring it into (9); one can getwhere , .
Bring into the stability function (7) and make it equal to the approximation to get , , and and get the corresponding :and bringing these into (5), one can get the formula.
When the node is the zero point of the Legendre polynomial, take as an example, get , , , get , and get the corresponding :and bringing these into (5) can be used to obtain a solution formula. The L-stability solution format for nonequidistant nodes of other r values can be constructed using the same method.
4. Numerical Example
Planar two-link manipulator is shown in Figure 1; the length of the link is , where , and the quality is , where . The corners are denoted by and , and at the initial time, . , represent the coordinate position of the connecting rod centroid.
Based on the Lagrange function, planar two-link dynamics model is established, where , is a generalized mass matrix, and are the moments of inertia of the connecting rod, and is a generalized force vector.
Here using the , we take r=4 as an example to divide the cells into four equal parts for the solution, and the step size h=0.01 and the time history 10s are selected and programmed using MATLAB. The Newton iteration method is used for calculation and the maximum allowable error is . The displacement trajectory at the end of the connecting rod can be obtained as shown in Figure 2; the starting position of the connecting rod end is indicated by a circle symbol, and the ending position is represented by a star. The displacement time course of the connecting rod end is shown in Figure 3. The total energy, kinetic energy, and potential energy of the system are shown in Figure 4. Displacement constraints, speed-level constraints, and acceleration-level constraints are shown in Figure 5.
It can be seen from Figures 2–5 that the L-stable method given in this paper can maintain good stability during the simulation process, and the total energy error and the constraint error of each level are small, which is suitable for multibody system dynamics simulation. The numerical verification and comparison of the method are carried out in different situations.
The simulation time is t=10s, the number of identical equidistant nodes (r=4) is taken in each time interval, and the results of different step sizes are compared in Table 1. When the step size is reduced, the time spent on program operation is increasing, but the maximum energy error is decreasing; at the same time, displacement constraints, speed-level constraints, and acceleration-level constraints are better when the step size is smaller. Therefore, in the case of a short simulation time, choosing a smaller step size will improve the accuracy of the solution, and accuracy better maintains displacement constraints, speed-level constraints, and acceleration-level constraints.
The simulation time is t=10s, and the step size is h=0.01, selecting the number of different equidistant nodes in each time interval for comparison (r=3,4) in Table 2. When r=4, although the time spent on program operation is increasing, the maximum energy error, displacement constraints, speed-level constraints, and acceleration-level constraints are less than r=3; therefore, an increase in the number of nodes can improve the accuracy.
The simulation time is t=10s, and the step size is h=0.01. We compare L-stable method when r=4 with fourth-order Explicit Runge-Kutta method (ERK4) in Table 3. In the L-stable method (L) for solving differential-algebraic equations, the maximum energy error, the displacement constraints, and the speed-level constraints are obviously smaller than those in the fourth-order Explicit Runge-Kutta method, although the acceleration-level constraints are not as good as those in the fourth-order Explicit Runge-Kutta method; the total energy of the L-stable method keeps stable with time, and the total energy of the fourth-order Explicit Runge-Kutta method has been increasing rapidly with time, as in Figures 6 and 7.
The simulation time is t=10s, and the step size is h=0.01. Under the same number of nodes (r=4), the differential-algebraic equations of index-1, index-2, and index-3 are solved in L-stable method in Table 4. The time taken for the differential-algebraic equations of the L-stability block format to solve index-1, index-2, and index-3 is roughly the same, and the maximum energy errors for index-1 and index-2 are significantly smaller than those for index-3. The displacement constraints for index-3 are the smallest, the speed-level constraints for index-2 are the smallest, and the acceleration- level constraints for index-1 are the smallest. By comprehensive comparison, the overall accuracy of index-1 is higher.
The step size is h=0.01; under the same step size, the L-stable method of index-3 selects equidistant nodes to be compared with nonequidistant nodes (Chebyshev nodes, Legendre nodes) in Table 5. By comparison, it is found that when the L-stable method selects nonequidistant nodes, the maximum energy error is smaller than equidistant nodes, and Chebyshev nodes have higher overall accuracy.
5. Accuracy Experiments and Analysis
The differential-algebraic equation of index-3 is used to solve the planar two-link model in time interval , taking the result of fourth-order Explicit Runge-Kutta method as its approximate exact solution, and calculate the logarithmic relationship between the overall error of and the step size of the L-stable method at t=1s, r=4, as shown in Figure 8; the logarithm of the overall error and step size for each component of is shown in Figures 9 and 10, when r=3; the logarithmic relation curve of the overall error is shown in Figure 11; the logarithm of the overall error and step size for each component is shown in Figures 12 and 13. The formula is used to calculate the accuracy. When r=4, the L-stable method overall error can guarantee the fourth-order accuracy. When r=3, the L-stable method overall error can guarantee the second-order accuracy; i.e., the accuracy of r is greater.
6. Stability Analysis
According to Theorem 1, one can construct A-stable method and get matrix and .
Taking the planar two-link manipulator as an example, with the simulation time being t=10s and the step size being h=0.01, the differential-algebraic equations for the index-1, -2, and -3 are obtained in Table 6. It can be seen that the results of solving the index-1, -2 are basically the same as those obtained by the L-stable method, but when solving the index-3, there is a singularity that cannot be solved. It can be seen that the L-stable method is also helpful for solving high-index problems.
The L-stable method in this paper can be written into the corresponding Runge-Kutta method [15]. Taking the approximation as an example, the corresponding fourth-order Runge-Kutta coefficient table is as follows.
The Radau IA, Radau IIA, and Lobatto IIIC methods in the Implicit Runge-Kutta method are also L-stabilized [16]. The above method is used to calculate the planar two-link manipulator dynamics model. In order to compare the results with the results of this paper, the fourth-order method is selected for calculation. The step size is h=0.01, and the simulation time is t=10s. The comparison results are shown in Table 7.
By comparison, it is found that running time, maximum energy error, displacement constraints, speed-level constraints, and acceleration-level constraints of the L-stable method constructed in this paper for solving the multibody system dynamics model are better than those in the Radau IA, Radau IIA, and Lobatto IIIC methods of the Implicit Runge-Kutta, although these methods are L-stable.
7. Conclusion
This paper constructs L-stable method for solving differential-algebraic equations of multibody system dynamics. The numerical results show that, under the equidistant node, when the same number of nodes is selected, the smaller the step size, the higher the accuracy, and when the step size is fixed, the more the nodes, the higher the accuracy. When compared with the Explicit Runge-Kutta method, the L-stable method better maintains the displacement constraints and speed-level level constraints, the maximum energy error is significantly smaller than the Explicit Runge-Kutta method, and the total energy retention effect is better under long-term simulation. When compared with A-stable method, L-stable method has no singularity in calculating high-index systems of equations, and each index keeps good results. When compared with Radau IA, Radau IIA, and Lobatto IIIC method, L-stable method is more efficient and has better performance than the above methods. The numerical results show that the nonequidistant nodes (Chebyshev nodes and Legendre nodes) result accuracy is higher than equidistant node. The method is simple in mathematics, easy to use, and easy to program. It is easy to implement in computers and worthy of research and promotion. How to improve the calculation efficiency while ensuring the numerical accuracy is the content that needs further study in the future.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study has been supported by the National Natural Science Foundations of China under Grant no. 11472143, 11772166.