Research Article  Open Access
Bowen Li, Jieyu Ding, Yanan Li, "LStable Method for DifferentialAlgebraic Equations of Multibody System Dynamics", Mathematical Problems in Engineering, vol. 2019, Article ID 9283461, 11 pages, 2019. https://doi.org/10.1155/2019/9283461
LStable Method for DifferentialAlgebraic Equations of Multibody System Dynamics
Abstract
An Lstable method over time intervals for differentialalgebraic 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 Lstable solution formula are obtained by comparison with Pade approximation. Newton iteration method is used during the solution process. Taking the planar twolink manipulator system as an example, the results of Lstable 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 RungeKutta method, Astable 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 longterm 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 Astable, that is, , , and all . Subsequently, Axellson [2] determines the of the by the quadrature coefficients at the zeros of or , constructing a class of Astable 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 Astable concept to nonlinear problems. In 1975, Dahlquist [5] proposed the Gstable concept of singlemethod linear multistep method. Butcher [6, 7] studied the nonlinear stability of RungeKutta method by using nonlinear equations as model equations, proposed Bstable and algebraic stability of the RungeKutta 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 Astable is not enough. Ehle [9] proposed Lstable. Hairer [10] pointed out that the Lstable numerical method will introduce damping effect, but it is very effective for suppressing highfrequency oscillation, and some scholars have proposed implicit singlestep and multistep methods based on Lstable for solving rigid equations, but these methods are based on ordinary differential equations. Lstable 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 longterm 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 Lstable 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 twolink manipulator system verified and compared.
2. DAEs of Multibody System Dynamics
The DAEs of index3 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 index2 and the index1 can be obtained by the combination of (1).When the initial values are known, the RungeKutta method can be used to directly solve the index1 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 index3 equation (1) and the index2 equation (2) cannot be solved directly using the RungeKutta method. The Lstable method is constructed below and can be used to solve the multibody system dynamics DAEs of index1, 2, and 3.
3. Construction of the DAEs LStable Format
3.1. General Form of LStable Method
In 1972, H. A. Watts and L. F. Shampine [13] proposed a block implicit onestep 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 Astable. 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 singlestep interval, and in the iterative solution process, a singlestep 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 index1 equation (3) as an example, after the order reduction, the following firstorder 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 rorder unit matrix, and . For the stability function, using the Pade approximation, we can construct Lstabilized solution formula.
Theorem 1 ( theorem and conjecture). For approximation, if , then this Pade approximation is Astable; if , then this Pade approximation is Lstable [14].
Perform a Taylor expansion on at , let the coefficient of the term be a pending constant, and omit the following highorder 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. LStable 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 Lstable solution formats for other r values.
3.3. LStable 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 Lstability solution format for nonequidistant nodes of other r values can be constructed using the same method.
4. Numerical Example
Planar twolink 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 twolink 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, speedlevel constraints, and accelerationlevel constraints are shown in Figure 5.
It can be seen from Figures 2–5 that the Lstable 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, speedlevel constraints, and accelerationlevel 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, speedlevel constraints, and accelerationlevel 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, speedlevel constraints, and accelerationlevel 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 Lstable method when r=4 with fourthorder Explicit RungeKutta method (ERK4) in Table 3. In the Lstable method (L) for solving differentialalgebraic equations, the maximum energy error, the displacement constraints, and the speedlevel constraints are obviously smaller than those in the fourthorder Explicit RungeKutta method, although the accelerationlevel constraints are not as good as those in the fourthorder Explicit RungeKutta method; the total energy of the Lstable method keeps stable with time, and the total energy of the fourthorder Explicit RungeKutta 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 differentialalgebraic equations of index1, index2, and index3 are solved in Lstable method in Table 4. The time taken for the differentialalgebraic equations of the Lstability block format to solve index1, index2, and index3 is roughly the same, and the maximum energy errors for index1 and index2 are significantly smaller than those for index3. The displacement constraints for index3 are the smallest, the speedlevel constraints for index2 are the smallest, and the acceleration level constraints for index1 are the smallest. By comprehensive comparison, the overall accuracy of index1 is higher.

The step size is h=0.01; under the same step size, the Lstable method of index3 selects equidistant nodes to be compared with nonequidistant nodes (Chebyshev nodes, Legendre nodes) in Table 5. By comparison, it is found that when the Lstable 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 differentialalgebraic equation of index3 is used to solve the planar twolink model in time interval , taking the result of fourthorder Explicit RungeKutta method as its approximate exact solution, and calculate the logarithmic relationship between the overall error of and the step size of the Lstable 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 Lstable method overall error can guarantee the fourthorder accuracy. When r=3, the Lstable method overall error can guarantee the secondorder accuracy; i.e., the accuracy of r is greater.
6. Stability Analysis
According to Theorem 1, one can construct Astable method and get matrix and .
Taking the planar twolink manipulator as an example, with the simulation time being t=10s and the step size being h=0.01, the differentialalgebraic equations for the index1, 2, and 3 are obtained in Table 6. It can be seen that the results of solving the index1, 2 are basically the same as those obtained by the Lstable method, but when solving the index3, there is a singularity that cannot be solved. It can be seen that the Lstable method is also helpful for solving highindex problems.

The Lstable method in this paper can be written into the corresponding RungeKutta method [15]. Taking the approximation as an example, the corresponding fourthorder RungeKutta coefficient table is as follows.
The Radau IA, Radau IIA, and Lobatto IIIC methods in the Implicit RungeKutta method are also Lstabilized [16]. The above method is used to calculate the planar twolink manipulator dynamics model. In order to compare the results with the results of this paper, the fourthorder 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, speedlevel constraints, and accelerationlevel constraints of the Lstable 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 RungeKutta, although these methods are Lstable.
7. Conclusion
This paper constructs Lstable method for solving differentialalgebraic 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 RungeKutta method, the Lstable method better maintains the displacement constraints and speedlevel level constraints, the maximum energy error is significantly smaller than the Explicit RungeKutta method, and the total energy retention effect is better under longterm simulation. When compared with Astable method, Lstable method has no singularity in calculating highindex systems of equations, and each index keeps good results. When compared with Radau IA, Radau IIA, and Lobatto IIIC method, Lstable 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.
References
 G. G. Dahlquist, “A special stability problem for linear multistep methods,” BIT Numerical Mathematics, vol. 3, pp. 27–43, 1963. View at: Publisher Site  Google Scholar  MathSciNet
 O. Axelsson, “A class of Astable methods,” BIT Journal, vol. 9, no. 3, pp. 185–199, 1969. View at: Publisher Site  Google Scholar
 O. B. Widlund, “A note on unconditionally stable linear multistep methods,” BIT Numerical Mathematics, vol. 7, pp. 65–70, 1967. View at: Publisher Site  Google Scholar  MathSciNet
 C. W. Gear, The Automatic Integration of Stiff Ordinary Differential Equations, North Holland publishing Company, Amsterdam, Netherlands, 1963. View at: MathSciNet
 G. Dahlquist, Error Analysis for a Class of Methods for Stiff Nonlinear Initial Value Problems, SpringerVerlag, Berlin, Germany, 1975. View at: MathSciNet
 J. C. Butcher, “A stability property of implicit RungeKutta methods,” BIT Journal, vol. 15, no. 4, pp. 358–361, 1975. View at: Publisher Site  Google Scholar
 K. Burrage and J. C. Butcher, “Stability criteria for implicit RungeKutta methods,” SIAM Journal on Numerical Analysis, vol. 16, no. 1, pp. 46–57, 1979. View at: Publisher Site  Google Scholar  MathSciNet
 S. F. Li, “Nonlinear stability of general linear methods,” Journal of Computational Mathematics, vol. 9, no. 2, pp. 97–104, 1991. View at: Google Scholar  MathSciNet
 B. L. Ehle, “Astable methods and Pade Approximations to the exponential,” SIAM Journal on Mathematical Analysis, vol. 4, no. 5, pp. 671–680, 1973. View at: Publisher Site  Google Scholar  MathSciNet
 E. Hairer and G. Wanner, Solving Ordinary Differential Equations II stiff and differentialAlgebraic problems, Science Press, Beijing, China, 2nd edition, 2006. View at: Publisher Site  MathSciNet
 J. Y. Ding and Z. K. Pan, “Generalizedα projection method for differentialalgebraic equations of multibody dynamics,” Engineering Mechanics, vol. 30, no. 4, pp. 380–384, 2013. View at: Google Scholar
 L. P. Wen, C. H. Yang, and H. Y. Wen, “Stability and asymptotic stability of multistep rungekutta methods for nonlinear functionalintegrodifferential equations,” Natural Science Journal of Xiangtan University, vol. 40, no. 1, pp. 1–5, 2018. View at: Google Scholar
 H. A. Watts and L. F. Shampine, “Astable block implicit onestep methods,” Bit Numerical Mathematics, vol. 12, no. 2, pp. 252–266, 1972. View at: Publisher Site  Google Scholar  MathSciNet
 Z. D. Yuan, G. Fei, and D. G. Liu, Numerical Solution of Initial Value Problems for Stiff Ordinary DifferentialEquations, Science Press, Beijing, China, 2016.
 Z. Q. Wu, “LStable methods for numerical solution of structural dynamics equations and multibody dynamics equations,” in Changsha, Hunan: National University of Defense Technology, pp. 1–132, 2009. View at: Google Scholar
 F. S. Zhu, “A Class of RungeKutta methods with Lstability and Bstability,” Journal of Mathematics, vol. 21, no. 2, pp. 183–188, 2001. View at: Google Scholar  MathSciNet
Copyright
Copyright © 2019 Bowen Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.