Journal of Applied Mathematics and Physics
Vol.06 No.08(2018), Article ID:86549,11 pages
10.4236/jamp.2018.68138

Improving Wilson-θ and Newmark-β Methods for Quasi-Periodic Solutions of Nonlinear Dynamical Systems

G. Liu, Z. R. Lv, Y. M. Chen*

Department of Applied Mechanics and Engineering, School of Aeronautics and Astronautics, Sun Yat-sen University, Guangzhou, China

Received: May 24, 2018; Accepted: August 6, 2018; Published: August 9, 2018

ABSTRACT

Quasi-periodic responses can appear in a wide variety of nonlinear dynamical systems. To the best of our knowledge, it has been a tough job for years to solve quasi-periodic solutions, even by numerical algorithms. Here in this paper, we will present effective and accurate algorithms for quasi-periodic solutions by improving Wilson-θ and Newmark-β methods, respectively. In both the two methods, routinely, the considered equations are re-arranged in the form of incremental equilibrium equations with the coefficient matrixes being updated in each time step. In this study, the two methods are improved via a predictor-corrector algorithm without updating the coefficient matrixes, in which the predicted solution at one time point can be corrected to the true one at the next. Numerical examples show that, both the improved Wilson-θ and Newmark-β methods can provide much more accurate quasi-periodic solutions with a smaller amount of computational resources. With a simple way to adjust the convergence of the iterations, the improved methods can even solve some quasi-periodic systems effectively, for which the original methods cease to be valid.

Keywords:

Wilson-θ Method, Newmark-β Method, Quasi-Periodic Solution, Predictor-Corrector Algorithm

1. Introduction

Both Wilson-θ and Newmark-β methods are generally used approaches for solving numerical solutions of dynamical systems [1]-[6]. As we know, a fundamental assumption of the Wilson-θ method lies in that, the acceleration changes linearly in a single time step [7]. For this reason, the Wilson-θ is considered to be one of linear acceleration methods. One of the best merits of this method is that its result is unconditionally stable when the internal parameter θ is chosen to be larger than 1.37 [8]. The Newmark-β method is in essence an extension of linear acceleration method [9]. When solving linear problems, both the two methods can provide accurate results by adjusting the internal parameters.

For some nonlinear problems, sometimes the two methods would provide false results, especially when there are quasi-periodic or chaotic responses [10]. Even worse, the methods may even cease to be valid for some cases, for instance leading to numerical divergence as shown later. Fu et al. [11] proposed an approximate method to linearize a quasi-periodic nonlinear equation, and suggested a gradual integration process based on the Wilson-θ method. To the best of our knowledge, it is tough to make such a linearization to complex quasi-periodic systems especially when strongly nonlinearities are taken into account. Even though the linearization is done successfully, it will pose restriction to the computational accuracy [12] [13].

In this paper, both the Wilson-θ and Newmark-β methods will be improved with the help of a predictor-corrector algorithm. This algorithm is based on the incremental process which makes Wilson-θ and Newmark-β methods free of repeated linearizations of the nonlinear terms. The initial solution is predicted at the previous time point and corrected recursively to be the true one at the next point. Numerical solutions have been successfully obtained for quasi-periodic responses of both autonomous and non-autonomous nonlinear dynamical systems. Interestingly, the results show that the iterative correction of the predicted solution can converge within two iterations, which ensures the improved methods with high computational efficiency. Most importantly, the improved methods can directly deal with a variety of quasi-periodic dynamical systems with even strong nonlinearities, as there is no requirement of linearizing the considered equations.

2. Classical Wilson-θ and Newmark-β Methods

2.1. Linearization

The nonlinear dynamical system is expressed as the following form

P [ x ¨ ( t ) , x ˙ ( t ) , x ( t ) ] = Q ( t ) (1)

where x ¨ ( t ) , x ˙ ( t ) , x ( t ) represent the generalized acceleration, velocity and displacement, respectively. After a step length Δ t , the following equation will be obtained as

P [ x ¨ ( t + Δ t ) , x ˙ ( t + Δ t ) , x ( t + Δ t ) ] = Q ( t + Δ t ) (2)

The incremental balance equation can be obtained by subtracting Equations ((1) from (2)). In practice, the linearization equation can also be obtained by expanding the incremental equilibrium equation as Taylor series at time t with the higher order terms being neglected

P ( x ¨ , x ˙ , x ) x ¨ | t Δ x ¨ + P ( x ¨ , x ˙ , x ) x ˙ | t Δ x ˙ + P ( x ¨ , x ˙ , x ) x | t Δ x = Q ( t + Δ t ) Q ( t ) (3)

which can be expressed in matrix form as

M ˜ Δ x ¨ + C ˜ Δ x ˙ + K ˜ Δ x = Δ Q (4)

with coefficient matrixes as

M ˜ = P ( x ¨ , x ˙ , x ) x ¨ | t , C ˜ = P ( x ¨ , x ˙ , x ) x ˙ | t , K ˜ = P ( x ¨ , x ˙ , x ) x | t , Δ P = Q ( t + Δ t ) Q (t)

At the beginning of each time step, one should calculate the coefficient matrixes. One possible way is to substitute into Equation (4) the velocity, acceleration and displacement of the previous time step, i.e., x ¨ , x ˙ ( t ) , x ( t ) . It is worth noting that, however, this treatment will bring serious numerical errors as the higher-order terms are all removed from the Taylor expansion. On the other hand, Equation (1) will be unbalanced if we substitute x ¨ ( t + Δ t ) , x ˙ ( t + Δ t ) , x ( t + Δ t ) into the coefficient matrixes. Denote the unbalanced term as

R = P [ x ¨ ( t + Δ t ) , x ˙ ( t + Δ t ) , x ( t + Δ t ) ] Q ( t + Δ t ) (5)

where R is also the truncation error for each calculation. The precision of the computing result can be controlled by adjusting R to a relatively low magnitude [14].

2.2. Solving Incremental Equation

A basic assumption of the Wilson-θ method is that the acceleration varies linearly between time [ t , t + θ Δ t ] , so that the velocity and displacement can be expressed as

{ x ˙ ( t + Δ t ) = x ˙ ( t ) + θ Δ t 2 [ x ¨ ( t + Δ t ) + x ¨ ( t ) ] x ( t + Δ t ) = x ( t ) + θ Δ t x ˙ ( t ) + ( θ Δ t ) 2 6 [ x ¨ ( t + θ Δ t ) + 2 x ¨ ( t ) ] (6)

Rewrite Equation (6) equivalently in the incremental form and substitute it into (4), we can deduce the incremental displacement equation at [ t , t + θ Δ t ]

K ¯ Δ x θ t = Δ Q ¯ (7)

in which

K ¯ = K ˜ + 3 θ Δ t C ˜ + 6 ( θ Δ t ) 2 M ˜ , Δ Q ¯ = Δ Q ˜ + [ θ Δ t 2 x ¨ ( t ) + 3 x ˙ ( t ) ] C ˜ + [ 6 θ Δ t x ˙ ( t ) + 3 x ¨ ( t ) ] M ˜ (8)

With Δ x θ t being obtained, we can update the velocity, acceleration, and displacement as

Δ x ¨ = a 0 Δ x θ t + a 1 x ˙ + a 2 x ¨ Δ x ˙ = a 3 Δ x θ t + a 4 x ˙ + a 5 x ¨ Δ x = a 6 Δ x θ t + a 7 x ˙ + a 8 x ¨ (9)

where the coefficients a i , i = 0 , 1 , , 8 are all listed in Appendix. And, the system response at time t + Δ t could be finally obtained as

x ¨ ( t + Δ t ) = x ¨ ( t ) + Δ x ¨ x ˙ ( t + Δ t ) = x ˙ ( t ) + Δ x ˙ x ( t + Δ t ) = x ( t ) + Δ x (10)

Similar to Equation (6), when the Newmark-β method is employed we can obtain

{ x ¨ ( t + Δ t ) = 1 α Δ t 2 [ x ( t + Δ t ) x ( t ) ] 1 α Δ t x ˙ ( t ) + ( 1 1 2 α x ¨ ( t ) ) x ˙ ( t + Δ t ) = x ˙ ( t ) + ( 1 δ ) Δ t x ¨ ( t ) + δ Δ t x ¨ ( t + Δ t ) (11)

with α and δ are priorly given parameters. Also, rewrite Equation (11) in the incremental form and substitute it into (4), then we can obtain

Δ x ¨ = b 0 Δ x + b 1 x ˙ + b 2 x ¨ Δ x ˙ = b 3 Δ x + b 4 x ˙ + b 5 x ¨ (12)

with the expressions of b i , i = 0 , 1 , , 5 in Appendix. Substituting Equation (12) into (10) yields the response at time t + Δ t .

3. Improving Wilson-θ and Newmark-β Methods

To elucidate the improved methods, we consider the following dynamic equation

M x ¨ ( t ) + C x ˙ ( t ) + K ( x ) x ( t ) = F ( x , x ˙ , t ) (13)

If Equation (13) is linear such that F ( x , x ˙ , t ) is independent upon either displacement x or velocity x ˙ , the Wilson-θ method can also be straightforwardly implemented by carrying out the following procedures

Step 1a: Given any initial conditions ( x ¨ 0 , x ˙ 0 and x 0 ), time step Δ t as well as θ 1.37 , calculate the integral constants c 0 ~ c 8 shown in Appendix.

Step 2a: Compute the payload, F ˜ t + θ Δ t , at time t + Δ t

F ˜ t + θ Δ t = F t + θ ( F t + Δ t F t ) + M ( c 0 x t + c 2 x ˙ t + 2 x ¨ t ) + C ( c 1 x t + 2 x ˙ t + c 3 x ¨ t ) (14)

Step 3a: Solve the displacement at time t + θ Δ t according to the following equation

K * x t + θ Δ t = F ˜ t + θ Δ t (15)

where K * = K + c 0 M + c 1 C .

Step 4a: Finally, calculate the displacement, velocity and acceleration at time t + Δ t according to the rules

{ x ¨ ( t + Δ t ) = c 4 [ x ( t + θ Δ t ) x ( t ) ] + c 5 x ˙ ( t ) + c 6 x ¨ ( t ) x ˙ ( t + Δ t ) = x ˙ ( t ) + c 7 [ x ¨ ( t + Δ t ) + x ¨ ( t ) ] x ( t + Δ t ) = x ( t ) + Δ t x ˙ ( t ) + c 8 [ x ¨ ( t + Δ t ) + 2 x ¨ ( t ) ] (16)

Then go to Step 1a to calculate the responses at the next time step.

Also, the Newmark-β method is capable of solving Equation (13) as long as it is linear, by employing the following procedures

Step 1b: Given any initial conditions ( x ¨ 0 , x ˙ 0 and x 0 ), time step Δ t as well as δ 0.50 , α 0.25 ( 0.5 + δ ) 2 , compute the constants d 0 ~ d 7 as shown in Appendix.

Step 2b: Compute the payload, F ˜ t + θ Δ t , at time t + Δ t

F ˜ t + Δ t = F t + Δ t + M ( d 0 x t + d 2 x ˙ t + d 3 x ¨ t ) + C ( d 1 x t + d 4 x ˙ t + d 5 x ¨ t ) (17)

Step 3b: Solve the displacement at time t + Δ t according to the following equation

K * x t + Δ t = F ˜ t + Δ t (18)

where K * = K + d 0 M + d 1 C .

Step 4b: Finally, calculate the displacement, velocity and acceleration at time t + Δ t according to the rules

{ x ¨ ( t + Δ t ) = d 0 [ x ( t + Δ t ) x ( t ) ] d 2 x ˙ ( t ) d 3 x ¨ ( t ) x ˙ ( t + Δ t ) = x ˙ ( t ) + d 6 x ¨ ( t ) + d 7 x ¨ ( t + Δ t ) x ( t + Δ t ) = F ˜ t + Δ t / K * (19)

Then go to Step 1b to calculate the responses at the next time step.

Note that, the above two methods are designed to solve linear systems. For nonlinear systems such as a cubic stiffness being included [15], the payload cannot be directly determined using either Step 2a or Step 2b, because the term F t + Δ t in both Equations (14) and (17) is a function of the responses (i.e., x ¨ t + Δ t , x ˙ t + Δ t , x t + Δ t ), which are to be determined at the next time step. For this reason, Fu et al. [11] suggested an approximate method by transforming the nonlinear problem into linear one during one time step. As will be shown in the numerical examples, this method requires a large amount of computational resources. Even worse, it ceases to be valid in some cases possibly because of large errors accumulation in the procedure of linearization.

The key procedure of the improved methods is to calculate the responses by a predictor-corrector algorithm, without transforming the nonlinearities into linear ones.

Denote the state responses of the system at time t as

{ u } t = { x ¨ ( t ) , x ˙ ( t ) , x ( t ) } (20)

After the time Δ t , the state responses are rewritten as

{ u } t + Δ t = { x ¨ ( t + Δ t ) , x ˙ ( t + Δ t ) , x ( t + Δ t ) } (21)

In both the improved Wilson-θ and Newmark-β methods, specifically, the second steps (Step 2a and Step 2b) which aim at calculating the payload are replaced by the following predictor-corrector algorithm

Step A: Predict the solution { u } t + Δ t at t + Δ t according to priorly given responses { u } t with a small increment Δ u .

Step B: Compute the payload F ˜ t + θ Δ t by applying the predicted solution { u } t + Δ t to Equation (14) or (17).

Step C: Obtain a new solution { u } t + Δ t at t + Δ t by applying F ˜ t + θ Δ t to Eqs. (15-16) or (18-19);

Step D: Correct { u } t + Δ t by proceeding Steps B-C-D if necessary, until the residual Re = | M x ¨ ( t ) + C x ˙ ( t ) + K ( x ) x ( t ) F ( x , x ˙ , t ) | is smaller than a given tolerance T o l .

For more details, Step 2a and Step 2b will be realized by the procedures displayed in Figure 1.

In real practice, there are two key variables to be elaborately chosen such as the initial increment Δ u and the tolerance T o l . In principle, a small value of Δ u is helpful to the convergence of the predictor-corrector algorithm. In our

Figure 1. A predictor-corrector algorithm for realizing the second steps in the improved Wilson-θ and Newmark-β methods, respectively.

numerical study, Δ u is chosen randomly chosen to be at the order of 10-6. Differently, the value of T o l should be chosen to be relatively large to guarantee the convergence. In this paper, T o l is chosen as T o l = O ( Δ t 2 ) because both the Wilson-θ and Newmark-β methods generally have a precision of the second order.

4. Numerical Examples

In this section, a nonlinear dynamical system will be presented as numerical examples to validate the presented methods for solving quasi-periodic solutions.

It should be pointed out that, both the classical Wilson-θ and Newmark-β methods are incapable of solving nonlinear dynamical systems. As mentioned above, Fu et al. [11] suggested an approximate approach to enable Wilson-θ be capable of doing a nonlinear problem. The Wilson-θ modified by Fu et al. [11] will also be employed in all the following numerical studies, for the purpose of comparing the effectiveness and precision of the presented improved methods.

Considering the forced system is subjected to external excitations varying over time domain such as [16]

M x ¨ + C x ˙ + ( K 1 + K 3 ( x ) ) x = F (22)

where x = [ x 1 x 2 ] , M = [ 1 0 0 1 ] , C = [ c 0 0 c ] , K 1 = [ 1 δ δ 1 ] , K ( 3 ) ( x ) = [ a 1 x 1 2 + a 2 x 1 x 2 , a 3 x 1 x 2 + a 4 x 2 2 b 1 x 1 2 + b 2 x 1 x 2 , b 3 x 1 x 2 + b 4 x 2 2 ] , F = [ f 1 cos ω t f 2 sin ω t ] with the constants ω = 2.3112 , c = 0.12 , δ = 0.1299 , f 1 = f 2 = 0.1587 , a 2 = a 4 = 0 , b 1 = b 3 = 0 , a 1 = a 3 = b 2 = b 4 = 1 .

Also, the time step is chosen as Δ t = 0.01 in this example. Figure 2 shows the comparisons of the displacement obtained by RK and the Wilson-θ methods, respectively. There is a slight discrepancy even at the early stage of the solution domain. As t increases further, the errors accumulate gradually resulting into noticeable differences. For quasi-periodic solutions, the coefficient matrix of the incremental balance equation should be corrected at every iteration step during calculation [17], which is possibly responsible for the rapid growing errors. The quasi-periodic solutions obtained by the improved methods are shown in Figure 3, also compared with the RK results. Nice agreement between the obtained soltuions domenstrate that both the improved Wilson-θ and Newmark-β methods can provide much more accurate quasi-periodic solutions.

To further check the computaion accuracy, the error curves provided by the obtained results are shown in Figure 4. The error of the classical Wilson-θ is nearly of the same magnitue of the quasi-periodic responses, whereas the errors of the improved methods are at a much lower level.

Figure 5 shows the CPU running times for different methods, versus the number of time steps. The improved methods are more computational efficient than the classical Wilson-θ method. The reason consists in that the iteration of

Figure 2. Displacement of Equation (23) obtained by RK and classical Wilson-θ methods [11], respectively.

Figure 3. Displacement of Equation (23) obtained by RK, improved Wilson-θ and improved Newmark-β methods, respectively.

classical method requires more iterations during Step 2a for convergence, that is designed to approximate the preload associated with nonlinearities. Interestingly, in most cases only two iterations are needed for the convergence of the second steps of the improved methods.

Figure 4. Error curves of the obtained results for Equation (23) with RK results as the references.

Figure 5. CPU running time versus the number of time steps for the classical Wilson-θ, improved Wilson-θ and improved Newmark-β methods, respectively.

5. Conclusions

Two improved methods have been proposed based on Wilson-θ and Newmark-β methods, respectively, for solving quasi-periodic responses of nonlinear dynamical systems. In order to avoid the same linearization procedure of Wilson-θ and Newmark-β methods, an improved fast algorithm is proposed in this paper, which enables both the two methods be capable of solving the quasi-periodic responses of nonlinear dynamical systems. In this algorithm, the solution at the next time step is first predicted according to the known one at the present time. And the predicted solution is corrected by iterative manner, which together provides us with a predictor-corrector algorithm.

Numerical examples have showed that, compared with the classical method, the improved methods can provide much more accurate quasi-periodic responses with less computational resources. Moreover, the improved methods are sometimes valid for systems for which the classical method is not. In addition, the convergence criterion is simple and effective for the presented methods, which guarantees both the computation accuracy and efficiency.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (11572356, 11672337).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Liu, G., Lv, Z.R. and Chen, Y.M. (2018) Improving Wilson-θ and Newmark-β Methods for Quasi-Periodic Solutions of Nonlinear Dynamical Systems. Journal of Applied Mathematics and Physics, 6, 1625-1635. https://doi.org/10.4236/jamp.2018.68138

References

  1. 1. Fang, D.P. and Wang, Q.F. (2010) Accuracy Analysis of Wilson-θ Method with Modified Acceleration. J. Vib. Shock, 29, 216-218.

  2. 2. Dai, H.L., Dai, T. and Yan, X. (2015) Thermoelastic Analysis for Rotating Circular HSLA Steel Plates with Variable Thickness. Appl. Math. Comput., 268, 1095-1109. https://doi.org/10.1016/j.amc.2015.07.017

  3. 3. Zupan, E., Saje, M. and Zupan, D. (2013) Dynamics of Spatial Beams in Quaternion Description Based on the Newmark Integration Scheme. Comput. Mech., 51, 47-64. https://doi.org/10.1007/s00466-012-0703-0

  4. 4. Dabaghi, F., Petrov, A., Pousin, J., et al. (2016) A Robust Finite Element Redistribution Approach for Elastodynamic Contact Problems. Appl. Numer. Math., 103, 48-71. https://doi.org/10.1016/j.apnum.2015.12.004

  5. 5. Krause, R. and Walloth, M. (2012) Presentation and Comparison of Selected Algorithms for Dynamic Contact Based on the Newmark Scheme. Appl. Numer. Math., 62, 1393-1410. https://doi.org/10.1016/j.apnum.2012.06.014

  6. 6. An, C. and Su, J. (2011) Dynamic Response of Clamped Axially Moving Beams: Integral Transform Solution. Appl. Math. Comput., 218, 249-259. https://doi.org/10.1016/j.amc.2011.05.035

  7. 7. Wilson, E.A. and Parsons, B. (1970) Finite Element Analysis of Elastic Contact Problems Using Differential Displacements. Int. J. Numer. Method Eng., 2, 387-395. https://doi.org/10.1002/nme.1620020307

  8. 8. Liu, J.B. and Du, X.L. (2005) Structural Dynamics. China Machine Press, Beijing.

  9. 9. Cho, I.H. and Porter, K. (2015) Three-Stage Multiscale Nonlinear Dynamic Analysis Platform for Building-Level Loss Estimation. Earthquake Spectra, 31, 1021-1042. https://doi.org/10.1193/092712EQS293M

  10. 10. Janin, O. and Lamarque, C.H. (2001) Comparison of Several Numerical Methods for Mechanical Systems with Impacts. Int. J. Numer. Method Eng., 51, 1101-1132. https://doi.org/10.1002/nme.206

  11. 11. Fu, Z.G. and Ren, F.C. (2001) A Linearization and an Integral Method Step by Step for the General Strong Non-Linear Vibration Differential Equation. J. Nonlinear Dyn. Sci. Tech., 8, 359-364.

  12. 12. Mirzaee, A., Shayanfar, M. and Abbasnia, R. (2014) Damage Detection of Bridges Using Vibration Data by Adjoint Variable Method. Shock Vib., 2014, 1-17. https://doi.org/10.1155/2014/698658

  13. 13. Grammont, L., Vasconcelos, P.B. and Ahues, M. (2016) A Modified Iterated Projection Method Adapted to a Nonlinear Integral Equation. Appl. Math. Comput., 276, 432-441. https://doi.org/10.1016/j.amc.2015.12.019

  14. 14. Ferziger, J.H. and Peric, M. (2012) Computational Methods for Fluid Dynamics. Springer Science & Business Media.

  15. 15. Chen, Y.M., Liu, Q.X. and Liu, J.K. (2014) Limit Cycle Analysis of Nonsmooth Aeroelastic System of an Airfoil by Extended Variational Iteration Method, Acta Mech, 225, 1-9. https://doi.org/10.1007/s00707-013-1026-8

  16. 16. Kirrou, I. and Belhaq, M. (2016) On the Quasi-Periodic Response in the Delayed Forced Duffing Oscillator, Nonlinear Dyn, 84, 2069-2078. https://doi.org/10.1007/s11071-016-2629-0

  17. 17. Ren, F.C. (1997) Research on the Shafting Bending Torsional Coupling Vibration in Big Turbogenerator Unit. Ph.D. Dissertation, North China Electric Power University.

Appendix

Coefficients for Wilson-θ method based on the incremental equation

a 0 = 6 θ 3 Δ t 2 , a 1 = 6 θ 2 Δ t , a 2 = 3 θ , a 3 = 3 θ 2 Δ t , a 4 = 3 θ 2 , a 5 = ( 1 3 2 θ ) Δ t , a 6 = 1 θ 3 , a 7 = ( 1 3 θ 2 ) Δ t , a 8 = 1 2 ( 1 1 θ ) Δ t 2

Coefficients for Newmark-β method based on the incremental equation

b 0 = 1 α Δ t 2 , b 1 = 1 α Δ t , b 2 = 1 2 α , b 3 = δ α Δ t , b 4 = δ α , b 5 = 1 + δ 2 α

Integral constants in classical Wilson-θ method

c 0 = 6 ( θ Δ t ) 2 , c 1 = 3 θ Δ t , c 2 = 2 c 1 , c 3 = Δ t 2 θ , c 4 = c 0 θ , c 5 = c 2 θ , c 6 = 1 3 θ , c 7 = Δ t 2 , c 8 = Δ t 2 6

Integral constants in classical Newmark-β method

d 0 = 1 α Δ t 2 , d 1 = δ α Δ t , d 2 = 1 α Δ t , d 3 = 1 2 α 1 , d 4 = δ α 1 , d 5 = Δ t 2 ( δ α 2 ) , d 6 = Δ t ( 1 δ ) , d 7 = Δ t δ