Both Wilson-θ and Newmark-β methods are generally used approaches for solving numerical solutions of dynamical systems -. As we know, a fundamental assumption of the Wilson-θ method lies in that, the acceleration changes linearly in a single time step . 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 . The Newmark-β method is in essence an extension of linear acceleration method . 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 . 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.  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  .
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
The nonlinear dynamical system is expressed as the following form
where represent the generalized acceleration, velocity and displacement, respectively. After a step length , the following equation will be obtained as
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
which can be expressed in matrix form as
with coefficient matrixes as
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., . 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 into the coefficient matrixes. Denote the unbalanced term as
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 .
2.2. Solving Incremental Equation
A basic assumption of the Wilson-θ method is that the acceleration varies linearly between time , so that the velocity and displacement can be expressed as
Rewrite Equation (6) equivalently in the incremental form and substitute it into (4), we can deduce the incremental displacement equation at
With being obtained, we can update the velocity, acceleration, and displacement as
where the coefficients are all listed in Appendix. And, the system response at time could be finally obtained as
Similar to Equation (6), when the Newmark-β method is employed we can obtain
with and are priorly given parameters. Also, rewrite Equation (11) in the incremental form and substitute it into (4), then we can obtain
with the expressions of in Appendix. Substituting Equation (12) into (10) yields the response at time .
3. Improving Wilson-θ and Newmark-β Methods
To elucidate the improved methods, we consider the following dynamic equation
If Equation (13) is linear such that is independent upon either displacement x or velocity , the Wilson-θ method can also be straightforwardly implemented by carrying out the following procedures
Step 1a: Given any initial conditions ( and ), time step as well as , calculate the integral constants shown in Appendix.
Step 2a: Compute the payload, , at time
Step 3a: Solve the displacement at time according to the following equation
Step 4a: Finally, calculate the displacement, velocity and acceleration at time according to the rules
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 ( and ), time step as well as , , compute the constants as shown in Appendix.
Step 2b: Compute the payload, , at time
Step 3b: Solve the displacement at time according to the following equation
Step 4b: Finally, calculate the displacement, velocity and acceleration at time according to the rules
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 , the payload cannot be directly determined using either Step 2a or Step 2b, because the term in both Equations (14) and (17) is a function of the responses (i.e., ), which are to be determined at the next time step. For this reason, Fu et al.  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
After the time , the state responses are rewritten as
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 at according to priorly given responses with a small increment .
Step B: Compute the payload by applying the predicted solution to Equation (14) or (17).
Step C: Obtain a new solution at by applying to Eqs. (15-16) or (18-19);
Step D: Correct by proceeding Steps B-C-D if necessary, until the residual is smaller than a given tolerance .
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 and the tolerance . In principle, a small value of 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, is chosen randomly chosen to be at the order of 10-6. Differently, the value of should be chosen to be relatively large to guarantee the convergence. In this paper, is chosen as 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.  suggested an approximate approach to enable Wilson-θ be capable of doing a nonlinear problem. The Wilson-θ modified by Fu et al.  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 
where , , , , , with the constants , , , , , , .
Also, the time step is chosen as 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 , 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 , 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 5. CPU running time versus the number of time steps for the classical Wilson-θ, improved Wilson-θ and improved Newmark-β methods, respectively.
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.
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.
Coefficients for Wilson-θ method based on the incremental equation
, , , , , , , ,
Coefficients for Newmark-β method based on the incremental equation
, , , , ,
Integral constants in classical Wilson-θ method
, , , , , , , ,
Integral constants in classical Newmark-β method
, , , , , , ,
 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
 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
 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
 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
 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
 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
 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
 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