An Energy Conservation Algorithm for Nonlinear Dynamic Equation

1 Changan Auto Global R&D Center, State Key Laboratory of Vehicle NVH and Safety Technology, Chongqing 401120, China 2 School of Automotive Engineering, Dalian University of Technology, Dalian 116024, China 3 State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116024, China 4 College of Engineering, University of Michigan, Ann Arbor, MI 48109-2133, USA


Introduction
Numerical stability and algorithmic damping have been long recognized as two important aspects that need to be carefully handled in time integration algorithms for solving dynamic problems.Indeed, many works have been done in this area.For instance, to improve the stability, the Generalized-α method of Chung and Hulbert 1 , the HHT-α method of Hilber et al. 2 , and the WBZ-α method of Wood et al. 3 all demonstrate very good dissipation property either at low-frequency or high-frequency ranges.In references 4-8 , Fung presents a series of time-step algorithms that are based on different mathematical and mechanical principles and can be used to deal with linear dynamical problems.Recent works on

Derivation of the Integral Iteration Formula
Step-by-step time-integration algorithms are commonly used to solve dynamic equations which mostly come from actual engineering problems.By spatial discretization using the finite element method, a nonlinear system may be represented by a second order nonlinear ordinary differential equation as where M, K are n-by-n constant mass and stiffness matrices, respectively.x, ẋ, ẍ, f are vectors with rank n representing the displacement, velocity, and acceleration, respectively.t is time.The right hand term in 2.1 , f x, ẋ, t , is the force vector that includes all external forces such as the damping forces and the nonlinear forces.Using the matrix decomposition, the mass matrix M can be expressed as Substituting 2.2 into 2.1 and multiplying by L −1 , we can obtain Knowing that L, L −1 , L T , L −T are all constant matrices, a variable substitution can be executed as Substituting 2.4 into 2.3 , a new dynamic system which is equivalent to the original system can be obtained as where D is a diagonal matrix and its diagonal elements are the diagonal elements of matrix K. K 0 is a matrix whose diagonal elements are zero, and other elements are equal to those in K.Note that D ii ≥ 0 for M which is a positive definite matrix and K is a semipositive definite matrix.Let Journal of Applied Mathematics Then 2.5 can be rewritten as ÿ Dy F t − K 0 y.

2.8
The separate form of the above matrix equation can be expressed as follows: where k i 0 is the ith row of the matrix K 0 .From 2.9 , it is easy to see that the analytical solution of the displacement and the velocity can be obtained by the Duhamel integral as

2.11
In order to derive the time integral formula, let t t k τ in 2.10 and 2.11 , where τ is the integral time step, then we have

2.13
Inspecting 2.12 and 2.13 , it can be found that there are still some unknown parameters that need to be identified before the time integral formula can be carried out numerically.They are the right hand side terms consisting of the undetermined variables y t and F i t .The latter one may also be a function of y t and ẏ t .To proceed, the Taylor expansion formula is used to expand y t on the interval t k ≤ t ≤ t k τ as

2.14
In 2.14 , y t is expanded to the third order term following exactly the Taylor expansion process, while in calculating the fourth order term, a new variable vector a is Journal of Applied Mathematics 5 introduced.It should be pointed out that although the Taylor expansion is an approximation to the original variable y t , 2.14 is still an exact expression of the variable y t .This is because the last term in 2.14 which includes the newly introduced vector a can be interpreted to compensate for the difference between y t and the summation of the first three terms in 2.14 .In order to express the vector a by the variable y t , let y 1 y t t k τ , y 0 y t t k .From 2.14 , it can be obtained that

2.16
where an undetermined parameter β has been introduced in N 4 t to regulate the stability of the algorithm and will be determined by the energy conservation equation in the next section.By t t k , t t k τ in 2.14 , we can obtain

2.17
By means of multinomial interpolation, f t t k ≤ t ≤ t k τ can be written as

2.18
Here we use the third order interpolation, and generally one can choose the order of interpolation discretionarily based on solely the algorithm accuracy order that is needed.Different interpolation order will lead to different integration formulas.Now, we will derive the integration formulas first, and the discussion of the interpolation will be addressed at Section 4. Using 2.5 and 2.14 , the last two terms of the right hand of 2.11 can be expressed separately as Journal of Applied Mathematics where α i k , γ i k , k 0, 1, 2, 3 are scalar and can be obtained by follow polynomials

2.21
Knowing that y 1 y t t k τ , so we can write the above equation in a matrix form: where U 0 , U 1 , α k , γ k k 0, 1, 2, 3 are diagonal matrices and their diagonal elements are cos d i τ, sin d i τ/d i , α i k , γ i k i 1, 2, . . ., n, k 0, 1, 2, 3 , respectively.From 2.22 , the iterative solution y 1 can be expressed as Substituting 2.14 and 2.15 into 2.13 and through some mathematical manipulations, the velocity iteration formula can be obtained as

2.31
where η 0 , η 1 , η 2 , η 3 are diagonal matrices with diagonal elements η i 0 , η i 1 , η i 2 , η i 3 , respectively.c m is also a diagonal matrix and its diagonal elements are c i m m 0, 1, 2, 3 .Every term in the right hand side of 2.31 is given as bellow:

2.32
where V 0 , V 1 are diagonal matrices with diagonal elements −d i sin d i τ, cos d i τ, respectively.

Energy Conservation Equation
One reason of expressing the dynamic equation in the form of 2.1 is to establish the energy conservation equation more conveniently and more directly. where Substituting 2.30 and 2.31 into the left hand side of 3.4 , a polynomial of the undermined parameter β can be easily achieved.For the right hand side of 3.4 , the integral term can be firstly decomposed into two parts as follows: where the term f 1 is an integral of an autonomous system and can be integrated easily.The ẋ t in term f 2 can be expressed as a polynomial of time using the relationship in 2.4 and taking the derivative of 2.15 with respect to time.Two predictive methods are recommended for determining the unknown term y 1 in 2.15 .One is to let y 1 y 0 ẏ0 τ and the other is to let β 1 in 2.23 .Then through 3.5 , an algebraic equation with an undetermined parameter β can be established and β can be numerically obtained by the Newton iteration method or other algebraic methods.Finally, substituting β into 2.30 and 2.31 , a numerical result can then be achieved.

Calculations of the Interpolation
Before giving some numerical examples, choosing the proper interpolation form of 2.18 must be discussed because it will affect the accuracy and stability of the proposed algorithm.In the current study, the authors use the Hermite interpolation to approximate the r 1 , r 2 , r 3 , r 4 in 2.18 , that is, It should be noted that there are unknowns in the right hand term of 4.1 which are x i τ , ẋi τ .The prediction of the two unknowns is shown below.For example, we can let Then at every iteration of 3.4 , the parameter β can be updated.Submitting β into 2.30 and 2.31 , a new prediction of the displacement and velocity can therefore be obtained.

Numerical Examples
In this section we give some numerical examples to verify the effectiveness of the proposed algorithm, in particular, the advantage in stability of the proposed algorithm.Since 2.14 is a fourth order Taylor expansion, the energy conservation algorithm has fourth order accuracy.So we choose the Rounge-Kutta method as a numerical comparison.The numerical results show the advantages of the proposed energy conservation algorithm in terms of its integration stability and the ability to eliminate the algorithm damping inherent in the Rounge-Kutta method.

The Oscillation of a Nonlinear Simple Pendulum
The dynamic equation of a nonlinear single pendulum without damping can be written as ẍ ω 2 0 sin x 0, ω 2 0 1.0, x 0 1.57, 5.1 where x denotes the angular displacement.The numerical solutions are shown in Figure 1.
From the figure we can see that the proposed energy conservation method ECM can keep the numerical stability and have no computing damping under large-step comparing with the Rounge-Kutta RK method.The numerical result of parameter β is shown in Figure 2. Table   gives the comparison of the computing efficiency.The efficiency of the proposed method is not as good as the RK method due to the iteration of parameter β and the time needed to compute the associated matrices.Figure 3 gives the error analysis between the ECM and the RK under time step 1.0 s.

The Unforced Linear Vibration of the Cuboid Rigid Body with Two DOF
The structural diagram of the system is shown in Figure 4.The mass of the rigid body is m and the length of the hemline is a.The center of mass is collocated at the geometry center  point C .The mass moment of inertia around the center of mass is J and the stiffness of the spring is k.The deformations of the two springs are x 1 , x 2 .The displacement in vertical direction of the center of mass is x c .The angular displacement of the rigid body about the mass center is φ.Using the above parameters, the equation of motion of the system can be written as , and x 2 −1.Figures 5 and 6 compare the displacement x 1 and velocity x 2-dot results predicted by the proposed method and the RK method.It can be seen that even with a big time step 1.0 s, the proposed method still has an accurate numerical solution but the RK method does not.Table 2 gives the comparison of the computing efficiency.Again, the efficiency of the proposed method is lower than that of the RK method in calculating these two degrees of freedom problem.Furthermore, it is noticed that the RK method almost keeps the same efficiency in Sections 5.1 and 5.2.
Figure 7 gives the error analysis between the ECM and the RK under time step 1.0 s.As Figure 7 already shows that the accuracy of the ECM is almost same as the result of RK with a 0.001 time step, the comparison does not use the RK with a small time step.

The Unforced Nonlinear Oscillation of a Spring Pendulum with Two DOF
The dynamic equation of the spring pendulum can be written as  Figure 10 shows the comparison of the numerical results between ECM and RK methods under large time steps.It is obvious that the proposed method can eliminate algorithm damping better and provides better stability than the RK scheme.Parameters and initial conditions used in the calculation are as follows:    Figure 11 shows long-time response of 5.2 .Parameter and initial condition is as same as 5.4 .From this figure we can see that after long-term iteration the proposed method still keeps numerical stability under a relatively large time step.But the algorithm damping makes the RK lose accuracy.Table 3 gives the comparison of the computing efficiency between RK and the ECM algorithms.It is shown that, in solving the two degrees of freedom nonlinear problem, the efficiency of the EMC is lower than the RK method.Moreover, comparing the   2 and Table 3, it also can be seen that the computation efficiency of the ECM is worse for calculating nonlinear problems than for linear problems.
Figure 12 shows the error analysis between the ECM and the RK under time step 1.0 s.

Conclusion
1 The energy conservation algorithm has the advantage in stability and time step compared with some numerical means because the numerical solution has been corrected by the energy conservation equation.
2 All examples have shown that the energy conservation method can eliminate algorithm damping.It is also an effective means for calculating the long-term characteristics of nonlinear dynamic systems.
3 The proposed method conserves the angular momentum automatically.Although the efficiency of the energy conservation method is not as good as the RK algorithm as well as some other numerical methods discussed in the literature, the integration step is large enough to implement long-term integration with good numerical stability.4 The reason of the low efficiency of the proposed method is because the iterations need to calculate the parameter β and the time consumed in matrix computing needed by the algorithm.The efficiency of the EMC is lower in dealing with nonlinear problems compared with linear problems.

Figure 1 :
Figure 1: Angular Displacement comparison between the proposed energy conservation method ECM and the RK method.

Figure 2 :
Figure 2: Value of parameter β time step 1.0 s .

Figure 3 :
Figure 3: The log-log plot of the error between the ECM and RK.

Figures 8 Figure 6 :
Figures 8 and 9 show the numerical solution under different damping.Parameters and initial condition are given as follows:

Figure 7 :
Figure 7: The log-log plot of the error between the ECM and RK with time step 1.0 s.

c = 0 time step = 1 Figure 12 :
Figure 12:The log-log plot of the error between the ECM and RK.

Table 1 :
Comparison of computing efficiency.

Table 2 :
Comparison of computing efficiency.

Table 3 :
Comparison of computing efficiency.
elapsed time by the ECM in Table