# L-Stable Method for Differential-Algebraic Equations of Multibody System Dynamics

Authors: Bowen Li; Jieyu Ding; Yanan Li
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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 y′=μy as a model equation and proposed A-stable, that is, yr→0, r→∞, and all h>0. Subsequently, Axellson [2] determines the aik of the yi,k=yn,r-1+h∑k=1naikfyk,r by the quadrature coefficients at the zeros of pn-pn-1 or pn-pn-2, constructing a class of A-stable methods, where pn is the Legendre polynomial. Widlund [3] proposed A(α)-stable, that is, if all solutions of αkxn+k+…+α0xn=hβkfn+k+…+β0fn tend to zero as n→∞ when the method is applied, with a fixed h, to any differential equation of the form y′=μy,where μ is a complex constant which lies in the set Sα=z;arg⁡(-z)<α,z≠0. 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 y′=f(y) 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 are(1)Mqq¨+ΦqTq,tλ-Fq˙,q,t=0Φq,t=0where q=[q1(t)q2(t)…qn(t)]T are generalized coordinates of the system, q˙=dq/dt are generalized speed of the system, q¨=d2q/dt2 are generalized acceleration of the system, λ=λ1,λ2,…,λmT are the Lagrange multiplier vectors, M∈Rn×n is a symmetric positive definite mass matrix, F∈Rn are generalized force vectors, and Φ=[Φ1,Φ2,…,Φm]T are m complete displacement constraints.For the constraint equationΦ, the first derivative and the second derivative are, respectively, obtained for the timet, and the DAEs form of the index-2 and the index-1 can be obtained by the combination of (1).(2)Mqq¨+ΦqTq,tλ-Fq˙,q,t=0Φqq,tq˙+Φtq,t=0(3)Mqq¨+ΦqTq,tλ-Fq˙,q,t=0Φqq,tq¨+Φqq,tq˙qq˙+2Φqtq,tq˙+Φttq,t=0When the initial values q0=q(0),q˙0=q˙(0) 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. Therdiscrete moments are made into a block vector, the output values ofrmoments 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 [tk,tk+1] and a time step h=tk+1-tk are set, selecting tk<c1<c2<…<cr=tk+1 and constructing DAEs solution format, where ci(i=1,2,…,r) 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. (4)q˙=zMqz˙+ΦqTq,tλ-Fz,q,t=0Φqq,tz˙+Φqq,tzqz+2Φqtq,tz+Φttq,t=0Let y=qTzTλTT∈R2n+m, yi=y(ci),y˙i=y˙(ci),i=1,2,…,r, Y=y1Ty2T⋯yrTT,Y˙=y˙1Ty˙2T⋯y˙rTT; the construction solution formula is as follows(5)Y=e⊗yk+hd⊗y˙k+hB⊗Y˙where e=11⋯1T∈Rr, d∈Rr, B∈Rr×r, yk=y(tk),y˙k=y˙(tk), and ⊗ is Kronecker.Bring (5) into the test equation(6)y′t=αytyt0=y0,Reα<0and get stability function(7)Rλ=erTIr-λB-1e+λdwhere λ=αh, Ir is ther-order unit matrix, and er=0⋯01∈Rr. For the stability function, using thePade approximation, we can construct L-stabilized solution formula.Theorem 1 (Ehle theorem and conjecture).
Fork,j-Pade approximation, if k≤j≤k+2, then thisPadeapproximation is A-stable; if k<j≤k+2, then thisPade approximation is L-stable [14].
Perform a Taylor expansion onyi,i=1,2,…,r at tk, let the coefficient of the yk(r+1)=y(r+1)(tk) term be a pending constant, and omit the following high-order terms; one can get(8)yi≈yk+hiy˙k+hi22!y¨k+⋯+hirr!ykr+μir+1!ykr+1,i=1,2,⋯,rwhere μi,i=1,2,⋯,r, is pending constant and hi=ci-tk,i=1,2,⋯,r. Bringing (8) into (5), one can get the following.(9)Bd=h1h12⋯h1rμ1h2h22⋯h2rμ2⋮hrhr2⋯hrrμr·12h13h12⋯r+1h1r12h23h22⋯r+1h2r⋮12hr3hr2⋯r+1hrr100⋯0-1From (9), B,d, which are expressed using the undetermined constant, can be obtained bringing these into the stability function (7). The value of the undetermined constant μi,i=1,2,⋯,r, can be obtained by comparison with the correspondingPade approximation form, so that the value of B,dcan be obtained bringing into formula (5) a specific solution formula.


### 3.2. L-Stable Format Based on Equidistant Nodes

For the equidistant nodeci(i=1,2,...,r), there are ci=tk+i/rh,i=1,2,…,r. For example, when r=3,h=1, the equidistant node is ci=tk+i/3,i=1,2,3; at this time hi=i/3,i=1,2,3, and (9) becomes as follows.(10)Bd=1319127μ12349827μ2111μ3·1231342714343322712341000-1Calculate B,d, substitute stability function (7), and compare it with 2,3-Pade approximation formula; one can solve undetermined constants μ1=1/45,μ2=8/45,μ3=1, and bringing them into the expression B,d, one can obtain(11)B=13107120-371203401715815-115989838,d=13411202538Bringing them into (5) can be used to obtain a solution formula. The same method can be used to construct L-stable solution formats for otherr values.


### 3.3. L-Stable Format Based on Nonequidistant Nodes

When the nodeci(i=1,2,…,r) of the interval is the zero point of the Chebyshev polynomial, within the interval [tk,tk+1], one gets the following.(12)ci=tk+h21+cos⁡2i-12r-1π,i=1,2,…,r-1Take r=3, h=1 as an example, get c1=tk+1/2-2/4, c2=tk+1/2+2/4,c3=1, and bring it into (9); one can get(13)Bd=h1h12h13μ1h2h22h23μ2111μ3·12h13h124h1312h23h224h2312341000-1where h1=1/2-2/4, h2=1/2+2/4.BringB,d into the stability function (7) and make it equal to the 2,3-Pade approximation to get μ1=10-52/32, μ2=10+52/32, and μ3=1 and get the corresponding B,d:(14)B=11-22245-822472-1485+822411+2224-72+1482323-16,d=2-748-2+748-16T,and bringing these into (5), one can get the formula.When the nodeci(i=1,2,…,r) is the zero point of the Legendre polynomial, take r=3,h=1 as an example, get c1=1/2-3/3+1, c2=1/23/3+1, c3=1, get μ1=3-3/12,μ2=3+3/12,μ3=1, and get the corresponding B,d:(15)B=27+3729-173723+53729+1737227-3723-537212120d=-3+3723-3720T,and bringing these into (5) can be used to obtain a solution formula. The L-stability solution format for nonequidistant nodes of otherrvalues can be constructed using the same method.


## 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. Therdiscrete moments are made into a block vector, the output values ofrmoments 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 [tk,tk+1] and a time step h=tk+1-tk are set, selecting tk<c1<c2<…<cr=tk+1 and constructing DAEs solution format, where ci(i=1,2,…,r) 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. (4)q˙=zMqz˙+ΦqTq,tλ-Fz,q,t=0Φqq,tz˙+Φqq,tzqz+2Φqtq,tz+Φttq,t=0Let y=qTzTλTT∈R2n+m, yi=y(ci),y˙i=y˙(ci),i=1,2,…,r, Y=y1Ty2T⋯yrTT,Y˙=y˙1Ty˙2T⋯y˙rTT; the construction solution formula is as follows(5)Y=e⊗yk+hd⊗y˙k+hB⊗Y˙where e=11⋯1T∈Rr, d∈Rr, B∈Rr×r, yk=y(tk),y˙k=y˙(tk), and ⊗ is Kronecker.Bring (5) into the test equation(6)y′t=αytyt0=y0,Reα<0and get stability function(7)Rλ=erTIr-λB-1e+λdwhere λ=αh, Ir is ther-order unit matrix, and er=0⋯01∈Rr. For the stability function, using thePade approximation, we can construct L-stabilized solution formula.Theorem 1 (Ehle theorem and conjecture).
Fork,j-Pade approximation, if k≤j≤k+2, then thisPadeapproximation is A-stable; if k<j≤k+2, then thisPade approximation is L-stable [14].
Perform a Taylor expansion onyi,i=1,2,…,r at tk, let the coefficient of the yk(r+1)=y(r+1)(tk) term be a pending constant, and omit the following high-order terms; one can get(8)yi≈yk+hiy˙k+hi22!y¨k+⋯+hirr!ykr+μir+1!ykr+1,i=1,2,⋯,rwhere μi,i=1,2,⋯,r, is pending constant and hi=ci-tk,i=1,2,⋯,r. Bringing (8) into (5), one can get the following.(9)Bd=h1h12⋯h1rμ1h2h22⋯h2rμ2⋮hrhr2⋯hrrμr·12h13h12⋯r+1h1r12h23h22⋯r+1h2r⋮12hr3hr2⋯r+1hrr100⋯0-1From (9), B,d, which are expressed using the undetermined constant, can be obtained bringing these into the stability function (7). The value of the undetermined constant μi,i=1,2,⋯,r, can be obtained by comparison with the correspondingPade approximation form, so that the value of B,dcan be obtained bringing into formula (5) a specific solution formula.


## 3.2. L-Stable Format Based on Equidistant Nodes

For the equidistant nodeci(i=1,2,...,r), there are ci=tk+i/rh,i=1,2,…,r. For example, when r=3,h=1, the equidistant node is ci=tk+i/3,i=1,2,3; at this time hi=i/3,i=1,2,3, and (9) becomes as follows.(10)Bd=1319127μ12349827μ2111μ3·1231342714343322712341000-1Calculate B,d, substitute stability function (7), and compare it with 2,3-Pade approximation formula; one can solve undetermined constants μ1=1/45,μ2=8/45,μ3=1, and bringing them into the expression B,d, one can obtain(11)B=13107120-371203401715815-115989838,d=13411202538Bringing them into (5) can be used to obtain a solution formula. The same method can be used to construct L-stable solution formats for otherr values.


## 3.3. L-Stable Format Based on Nonequidistant Nodes

When the nodeci(i=1,2,…,r) of the interval is the zero point of the Chebyshev polynomial, within the interval [tk,tk+1], one gets the following.(12)ci=tk+h21+cos⁡2i-12r-1π,i=1,2,…,r-1Take r=3, h=1 as an example, get c1=tk+1/2-2/4, c2=tk+1/2+2/4,c3=1, and bring it into (9); one can get(13)Bd=h1h12h13μ1h2h22h23μ2111μ3·12h13h124h1312h23h224h2312341000-1where h1=1/2-2/4, h2=1/2+2/4.BringB,d into the stability function (7) and make it equal to the 2,3-Pade approximation to get μ1=10-52/32, μ2=10+52/32, and μ3=1 and get the corresponding B,d:(14)B=11-22245-822472-1485+822411+2224-72+1482323-16,d=2-748-2+748-16T,and bringing these into (5), one can get the formula.When the nodeci(i=1,2,…,r) is the zero point of the Legendre polynomial, take r=3,h=1 as an example, get c1=1/2-3/3+1, c2=1/23/3+1, c3=1, get μ1=3-3/12,μ2=3+3/12,μ3=1, and get the corresponding B,d:(15)B=27+3729-173723+53729+1737227-3723-537212120d=-3+3723-3720T,and bringing these into (5) can be used to obtain a solution formula. The L-stability solution format for nonequidistant nodes of otherrvalues can be constructed using the same method.


## 4. Numerical Example

Planar two-link manipulator is shown in Figure1; the length of the link is L1,L2, where L1=1m,L2=3m, and the quality is m1,m2, where m1=1kg,m2=2kg. The corners are denoted by θ1 and θ2, and at the initial time, θ1=π/3,θ2=-π/6. x1,y1,x2,y2 represent the coordinate position of the connecting rod centroid.Figure 1
Planar two-link manipulator.Based on the Lagrange function, planar two-link dynamics model is established, whereq=x1,y1,θ1,x2,y2,θ2T, M=diagm1,m1,I1,m2,m2,I2 is a generalized mass matrix, I1=1/12m1L12 and I2=1/12m2L22 are the moments of inertia of the connecting rod, and F=[0,-m1g,0,0,-m2g,0]T is a generalized force vector.Here using the2,4-Pade, we taker=4 as an example to divide the cells into four equal parts for the solution, and the step sizeh=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 10-3. 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.Figure 2
Displacement trajectory at the end of the connecting rod.Figure 3
Displacement time history of the end of the rod.Figure 4
Total energy, kinetic energy, and potential energy time history of the system.Figure 5
Time history of displacement and speed-level and acceleration-level constraints.It can be seen from Figures2–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.( 1 ) The simulation time ist=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 (max⁡εH) 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.Table 1
Results of L-stable method under different steps.
Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints h=0.01 5.8906 7.0007e-05 6.8459e-07 3.7480e-07 5.6595e-12 h=0.005 10.8750 8.8078e-07 1.9154e-09 6.6686e-09 1.9846e-12 h=0.001 42.6406 5.5074e-11 7.6766e-13 4.1800e-13 1.4211e-13( 2 ) The simulation time ist=10s, and the step size ish=0.01, selecting the number of different equidistant nodes in each time interval for comparison (r=3,4) in Table 2. Whenr=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 thanr=3; therefore, an increase in the number of nodes can improve the accuracy.Table 2
Results of the L-stable method under different nodes.
Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints r=3 1.6875 0.0021 5.8212e-05 3.9267e-04 6.2629e-12 r=4 5.8906 7.0007e-05 6.8459e-07 3.7480e-07 5.6595e-12( 3 ) The simulation time ist=10s, and the step size ish=0.01. We compare L-stable method whenr=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.Table 3
Comparison of L-stable method and Explicit Runge-Kutta method at the same step size,h=0.01.
Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints L 5.8906 7.0007e-05 6.8459e-07 3.7480e-07 5.6595e-12 ERK4 0.3906 0.2628 0.0087 0.0021 9.9476e-14Figure 6
L-stable method and fourth-order Explicit Runge-Kutta method for total system time history (10s).Figure 7
L-stable method and fourth-order Explicit Runge-Kutta method for total system time history (100s).( 4 ) The simulation time ist=10s, and the step size ish=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.Table 4
Comparison of differential-algebraic equations for L-stable method solving indices-1,2,3.
Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints index -1 5.8906 7.0007e-05 6.8459e-07 3.7480e-07 5.6595e-12 index -2 5.4219 3.2163e-05 4.3949e-09 1.7764e-14 0.0704 index -3 5.6094 0.0011 2.7756e-15 5.9405e-04 0.9854( 5 ) The step size ish=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.Table 5
Comparison of the results of the L-stable method taking uniform and nonuniform nodes.
Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints equidistant nodes 5.6094 0.0011 2.7756e-15 5.9405e-04 0.9854 Chebyshev nodes 5.2188 4.3245e-05 2.9976e-15 1.8094e-04 0.6880 Legendre nodes 5.1094 2.5327e-04 2.8866e-15 2.9047e-04 0.7560


## 5. Accuracy Experiments and Analysis

The differential-algebraic equation of index-3 is used to solve the planar two-link model in time intervalt∈0,1, taking the result of fourth-order Explicit Runge-Kutta method h=0.000001 as its approximate exact solution, and calculate the logarithmic relationship between the overall error of q and the step size of the L-stable method att=1s,r=4, as shown in Figure 8; the logarithm of the overall error and step size for each component of q is shown in Figures 9 and 10, whenr=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 log⁡e/log⁡h is used to calculate the accuracy. Whenr=4, the L-stable method overall error can guarantee the fourth-order accuracy. Whenr=3, the L-stable method overall error can guarantee the second-order accuracy; i.e., the accuracy ofris greater.Figure 8
The overall error ofq and the logarithm ofh.Figure 9
The overall error ofx1,y1,θ1and the logarithm ofh.Figure 10
The overall error ofx2,y2,θ2 and the logarithm ofh.Figure 11
The overall error ofq and the logarithm ofh.Figure 12
The overall error ofx1,y1,θ1 and the logarithm ofh.Figure 13
The overall error ofx2,y2,θ2and the logarithm ofh.


## 6. Stability Analysis

According to Theorem1, one can construct A-stable method and get matrix B and d.(16)B=1420812520-1094201912520-4350406245415245-190377280111140167280-31560644581564451445,d=141847504029901795601445TTaking the planar two-link manipulator as an example, with the simulation time beingt=10s and the step size beingh=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.Table 6
Comparison of DAEs for L-stable and A-stable method solving indices-1,2,3.
index-1 Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints A-stable 6.2188 3.8045e-05 2.3863e-07 2.0615e-07 5.6648e-12 L-stable 5.8906 7.0007e-05 6.8459e-07 3.7480e-07 5.6595e-12  index-2 Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints  A-stable 5.4688 1.1145e-05 1.6059e-09 1.7729e-14 0.0214 L-stable 5.4219 3.2163e-05 4.3949e-09 1.7764e-14 0.0704  index-3 Runtime/s max ⁡ ε H displacement constraints speed-level constraints acceleration-level constraints  A-stable Singular L-stable 5.6094 0.0011 2.7756e-15 5.9405e-04 0.9854The L-stable method in this paper can be written into the corresponding Runge-Kutta method [15]. Taking the 1,3-Pade approximation as an example, the corresponding fourth-order Runge-Kutta coefficient table is as follows.(17)
0
0
0
0
0
1
3
1
72
43
72
-
29
72
1
8
2
3
1
12
19
36
1
36
1
36
1
1
8
3
8
3
8
1
8
1
8
3
8
3
8
1
8The 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 ish=0.01, and the simulation time ist=10s. The comparison results are shown in Table 7.Table 7
Comparison of L-stable method (r=3) with Radau IA, Radau IIA, and Lobatto IIIC methods.
Runtime/s max ⁡ | ε H | displacement constraints speed-level constraints acceleration-level constraints L-stable 1.6875 0.0021 5.8212e-05 3.9267e-04 6.2629e-12 Radau IA 14.9063 1.4075 0.1803 0.1087 9.8100 Radau IIA 13.5625 0.1035 2.1396 8.2918e-04 9.8100 Lobatto IIIC 13.9375 0.1038 2.1396 0.0010 9.8100By 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.



---

*Source: 9283461-2019-05-29.xml*