Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 391815, 12 pages

http://dx.doi.org/10.1155/2014/391815

## An Extended Time Series Algorithm for Modal Identification from Nonstationary Ambient Response Data Only

^{1}Experimental Facility Division, National Synchrotron Radiation Research Center, Hsinchu 30076, Taiwan^{2}Department of Aeronautics and Astronautics, National Cheng Kung University, Tainan 70101, Taiwan^{3}Instrumentation Development Division, National Synchrotron Radiation Research Center, Hsinchu 30076, Taiwan

Received 10 March 2014; Revised 17 June 2014; Accepted 26 June 2014; Published 17 July 2014

Academic Editor: Jaromir Horacek

Copyright © 2014 Chang-Sheng Lin 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.

#### Abstract

Modal Identification is considered from response data of structural systems under nonstationary ambient vibration. The conventional autoregressive moving average (ARMA) algorithm is applicable to perform modal identification, however, only for stationary-process vibration. The ergodicity postulate which has been conventionally employed for stationary processes is no longer valid in the case of nonstationary analysis. The objective of this paper is therefore to develop modal-identification techniques based on the nonstationary time series for linear systems subjected to nonstationary ambient excitation. Nonstationary ARMA model with time-varying parameters is considered because of its capability of resolving general nonstationary problems. The parameters of moving averaging (MA) model in the nonstationary time-series algorithm are treated as functions of time and may be represented by a linear combination of base functions and therefore can be used to solve the identification problem of time-varying parameters. Numerical simulations confirm the validity of the proposed modal-identification method from nonstationary ambient response data.

#### 1. Introduction

Experimental identification of modal parameters of a structure is usually carried out by measuring both its input and its corresponding output. Some modal testing techniques use free or impulse responses of structures so that the input excitation need not be measured. However, there are situations where controlled excitation or initial excitation cannot be employed, such as the case of in-operation testing or in-flight measurement. Consequently, it is desirable to develop techniques for modal-parameter identification without the need of input measurement [1, 2].

Modal-parameter identification from ambient vibration data has gained considerable attention in recent years [3, 4]. The feature of experimental modal analysis lies in the emphasis on the importance of the modal-identification process of a structure from measured input/output data when compared with the classic vibration analysis. Normally, the analysis of system identification acquires the estimation of important parameters from the data measured in the experiments. Based on the historical development of time series, Yule was the first to purpose the autoregressive (AR) model of a single variable in 1927. He developed the application of AR model for investigating periodicities in disturbed series, with special reference to Wolfer’s sunspot numbers. In 1970, Box et al. [5] proposed a modeling cycle for the AR model, which assumes that future values of a time series can be expressed as a linear combination of its past values. Applications of the Box et al. methodology spread in the following decades, covering a wide range of scientific areas such as Biology, Astronomy, or Econometrics. The time-series algorithm is also extensively used in the modal identification [6, 7]. In 1980, Pandit and Wu [8] proposed a time-series method based on the mathematical theory of statistics and further developed the so-called data dependent system (DDS) using ARMA models for modal-parameter identification of structural systems [9, 10]. The main advantages of DDS are that they are accurate for extracting information from noisy signals, that estimation of the parameter covariance matrix is natural part of the estimation, and that they can be formulated as exact covariance equivalent discrete time domain models of a set of second order differential equations.

The conventional ARMA algorithm is applicable to perform modal identification, however, only for stationary-process vibration. In this paper, we will develop modal-identification techniques based on the nonstationary time series for linear systems subjected to nonstationary ambient excitation. Nonstationary ARMA model with time-varying parameters is considered because of its capability of resolving general nonstationary problems. The parameters of MA model in the nonstationary time-series algorithm are treated as functions of time and may be represented by a linear combination of base functions and therefore can be used to solve the identification problem of time-varying parameters. Numerical simulations confirm the validity of the proposed modal-identification method from nonstationary ambient response data.

#### 2. Time-Series Model of Structural Vibration Responses

The standard matrix equations of motion of a discrete linear system can be expressed as follows: In (1), , , and are, respectively, the system mass, damping, and stiffness matrices, while and are the excitation and displacement vectors, respectively. Applying Laplace transform to (1) yields where . is called the system impedance matrix, which is an inverse matrix of the transfer function matrix . According to adjoint matrix method, the transfer function can be expressed as where is the adjoint matrix of . It follows from (3) that the transfer function at the th degree of freedom (DOF) due to the input at the th DOF is where , . For a real structure vibration system with damping, the transfer function, (4), can be decomposed into where is a pole corresponding to the th mode of a system and is the residue in the form of a complex number. To obtain the time-dependent impulse response , it is, therefore, necessary to perform an inverse Laplace transformation of as follows: Define as the discrete time step, and (6) can be written as where is the time step and is the sampling interval. Let be a -transform factor; (7) can be rewritten as The transfer function of the discrete system results from the -transform of impulse response function as follows: Comparing the relationship between input and output in original transfer function, (9) can be further expressed as follows: One can, therefore, derive By inverse transforming both sides of the -domain (11) and rearranging, we obtain the following equation in the time domain: Equation (12) is a normal ARMA model. The poles of transfer function can then be obtained by solving the polynomial equation in (10) as follows: The relationship between the poles, and , of transfer function as well as natural frequencies and damping ratios of a structural system can be expressed as follows: from which the natural frequencies and damping ratios of a structural system can then be obtained to be

The preceding results show that the ARMA model is applicable to identify the modal parameters of an N-DOF structural system. However, the original ARMA model was effective only for stationary process. Remarkably, the relationship between AR and MA parts in ARMA model is approximated as a relationship between input and output. In the light of this concept, we come up with a new idea that if the nonstationary ambient vibration data can also be approximated directly as the aforementioned relationship of input and output, then the ARMA model can be extended for nonstationary process. In the following, we will propose a modification of ARMA model using the aforementioned concept to construct a nonstationary ARMA model, which can be applied to modal identification from nonstationary ambient vibration data.

#### 3. Estimation of ARMA Parameters

The ARMA algorithm is based on the assumption of a stationary input while the others assume that the input is a Gaussian white noise process. However, estimating the parameters is more involved for the ARMA model than it is for the AR model because the white noise signal in the ARMA model is also unknown, so a three-stage least-square approach or nonlinear approach has to be proposed. As mentioned in the previous section, it is important to estimate the model parameter for reconstructing the response data of a system through the ARMA algorithm. In this paper, we estimate the parameters of ARMA by using the nonlinear least-square method. However, it requires a good initial estimation as the initial parameters of the ARMA model when performing the nonlinear least-square method. We therefore propose a three-stage least-square method, which is applicable for nonstationary process, to obtain the initial estimated parameters of ARMA model.

##### 3.1. Three-Stage Least-Square Method for the Estimation of Initial Parameters of ARMA Model

Consider an ARMA model as follows: where is white noise, and are the orders of an ARMA model, and , and as well as , and are the parameters of ARMA model. By treating the linear combination of in the right-hand side of (16) as new error term , (16) can be rewritten as Equation (17) can be treated as a -order AR model, whose parameters can be determined through the least-square method as follows: where The estimations can then be obtained as follows: A new set of estimations can be composed of the obtained AR parameters, , and , as follows: Substituting the new observability data into the AR part in the original ARMA model, the following equation can be derived: Equation (22) can be treated as -order MA model. Due to the fact that the finite-order MA model can be approximated as the high-order AR model [9], (22) can be rewritten as where and is the order of high-order AR model. The parameters of high-order AR model can also be obtained through least-square method. Comparing (22) with (23), one can derive the following: from which the relationship between and can be obtained to be The above equation can be rewritten as the standard form of least-square method, and the parameters of MA model are further determined as follows: where The estimations can then be obtained as follows: In the preceding, it has been shown that the initial estimations of all parameters of ARMA model can be obtained through the three-stage least-square method without the assumption of stationary process. To obtain the relatively precise estimation of ARMA model, the nonlinear least-square method in conjunction with the Gauss-Newton algorithm will be employed, as described in detail next.

##### 3.2. Nonlinear Least-Square Method for Estimation of Parameters of ARMA Model

Consider an ARMA model, which can be expressed as where and are the orders of an ARMA model. , , and , , are constants, which are parameters of an ARMA model, as well as being white noise. It follows from (29) that can be derived as follows: By integrating the terms of parameters of ARMA, (29) can be rewritten as Let unknown be a set of the initial estimation of model parameters in conjunction with (30), whose vector form is as follows: In (32), is a function of and , where , , , , and is the number of samples in the input (or output) data. Define the prediction errors as follows: Equation (30) is necessarily to be employed to construct due to the fact that is unknown. will in general be nonlinear; that is, it is not the linear function of , and the minimization of cannot be performed. Thus, to obtain a practically applicable approach, through the Gauss-Newton algorithm, which is a method usually applied to solve nonlinear least squares problems, can be linearized using a first-order Taylor expansion at the operating point and then extract the linear part of as follows: where is an initial estimation set of ARMA parameters and can be expressed as Due to the fact that is a known observability data matrix when performing the Gauss-Newton algorithm and to emphasize that is function of , can thus be written as . Equation (34) can be therefore expressed as Substituting (36) into (34), the following equation can be derived: According to the theory of least-square method, through the minimization of , the estimated parameter can be expressed as In practice, the evaluation result of is an iterative form of running to a set of initial estimation for the beginning. Through the iteration evaluation of a set of , will be yielded better and then continue iterations until the final convergence of the results. The iteration form can be represented by and the criterion of the convergence of calculation results to stop iterations can be defined as follows: where is norm of parametric vector and is small value to be chosen and may be chosen as or .

In the preceding, it has been shown that the initial estimations of all parameters of ARMA model can be obtained through the three-stage least-square method without the assumption of stationary process. By using the nonlinear least-square method in conjunction with the Gauss-Newton algorithm, the relatively precise estimation of ARMA model can then be obtained. In the following, we will reconstruct a nonstationary ARMA algorithm and then further extend the time-series algorithm to general nonstationary ambient vibration problems.

#### 4. Nonstationary Time-Series Model for Modal Identification

The original time-series model is the well-approximately equivalent mathematical model for stationary process and is further available in practice. For vibration problems in engineering, under the assumption of stationary process, we can construct a mathematical model for the response of a system through the ARMA model. By extracting the AR parameters, the parameters of a system can then be obtained. We also treated the MA part as an equivalent excitation force for the vibration process. The aforementioned time-series algorithm was effective only for stationary process. Most ambient excitations encountered in engineering problems are, however, nonstationary in nature; that is, the characteristics of nonstationary stochastic process may change with time. The original time-series algorithm is therefore not applicable to identify the modal parameters of a structural system subjected to nonstationary ambient excitation.

In the following, we will improve the time-series algorithm to get rid of the previous restriction of stationary assumption. Without changing the AR parameters, we construct a nonstationary time-series model to effectively describe a wide variety of nonstationary behavior and then extract the modal parameters of a structural system. Furthermore, by introducing the base functions to construct a nonstationary time-series model, we extend the time-series algorithm to identify the modal parameters of a structural system subjected to nonstationary white noise.

##### 4.1. Nonstationary ARMA Algorithm Using the Temporal Root-Mean-Square Function from Structural Response Data

The conventional stationary ARMA model, whose parameters are all constants, can be expressed as In a previous study of the authors [4], it has been shown that if the excitation can be modeled as nonstationary white noise with a slowly time-varying envelope function, then the nonstationary responses of the system can also be modeled approximately as a product model with the same envelope function. In addition, if the original nonstationary response data could be represented by the product model with a slowly varying envelope function, the temporal root-mean-square functions of the data also have the same nonstationary trend as that of the original response data. Based on the above mentioned concept and assumption, in this paper, we first determine the temporal root-mean-square function and so the envelope function by using the interval average and then applying curve-fitting technique [4]. We can then acquire the approximate stationary responses by dividing the nonstationary responses of each DOF with the same , and a new set of observability data as quasistationary response data can then be obtained. The parameters of ARMA model can then be determined through the proposed algorithm in this paper. Define the nonstationary white noise in the product model by using as follows: By substituting into in (41), the nonstationary time-series model can be constructed as follows: It indicates that the right-hand side of (43) and so the MA part in the ARMA model describes a nonstationary behavior of time-varying amplitude. The parameters of AR part (in the left-hand side of (43)) are all constants and can then be employed to extract the modal parameters of a structural system using (15).

##### 4.2. Nonstationary ARMA Algorithm Using the Base Functions

The original ARMA model can be expressed as follows: where the parameters of the original ARMA model are constants due to the stationary assumption. To effectively describe a time-varying stochastic characteristic of nonstationary behavior, we propose the general nonstationary time-series model, in which the parameters of model are described as function of time as follows:

Consider a time-invariant linear dynamic system subjected to the nonstationary excitation; the MA part can be treated as an equivalent excitation force for the vibration process, and (45) can be simplified as The concept of this nonstationary time-series model describes a wide variety of nonstationary excitation by using the parametric function of MA part only, and the parameters of AR model are constants to determine the modal parameters of time-invariant system. Due to the fact that the parametric functions s in this model are all known ones, the general estimation method of model parameters cannot be employed. In this paper, we further come up with a new idea by using the known base function in conjunction with the combination of unknown coefficient instead of the original unknown parametric function in (46) as The original ARMA model can then be rewritten as Due to the fact that the base function is known, we can transfer the original problem to solve parametric function into the problem to solve the linear combination of . The choice of the base function depends on the adequate expression to nonstationary response data. In this paper, we choose the polynomial function as base function as follows: where is the number of samples in the input (or output) data and is the order of base function to be chosen. When introducing the base function to construct a nonstationary time-series model, the reconstructed ARMA model form is similar to that of the original time-series model. The major difference is that the MA parameters in the nonstationary base-function ARMA model are treated as functions of time and may be represented by a linear combination of suitable base functions. As long as by properly expanding the data matrix in the nonlinear and three-stage least-square method, respectively, the AR parameters and the coefficients of MA part corresponding to base function can be obtained. Then we construct nonstationary base-function ARMA model to describe the behavior of the nonstationary response data and further determine the modal parameters of a structural system. Note that although the structural modal parameters can be identified from the AR parameters only, the estimation of AR parameters and MA parametric function will affect each other in practice; estimation for AR parameters and MA parametric function cannot therefore be separately and solely performed.

#### 5. Numerical Simulation

To demonstrate the feasibility and effectiveness of the proposed method, we consider a linear 3-DOF chain model with viscous damping. A schematic representation of this model is shown in Figure 1. Assume that the three blocks of mass are 3, 1.5, and 1 Kg, respectively, as well as all spring constants being 4, 1, 1, and 2 N/m, respectively. The mass matrix , stiffness matrix , and the proportional damping matrix of the system are given as follows:

Note that the system has proportional damping, since the damping matrix can be expressed as a linear combination of and . Consider that the ambient vibration input can be modeled as nonstationary white noise as represented by the product model. The stationary white noise is generated using the spectrum approximation method [11] as a zero-mean band-pass noise, whose standard deviation is N^{2}·sec/rad with a frequency range from 0 to 50 Hz. The sampling interval is chosen as sec, and the sampling period is sec. The stationary white noise simulated is then multiplied by an amplitude-modulating function, as shown in Figure 2, to obtain the nonstationary white noise, which serves as the excitation input acting on the 3rd mass point of the system. The time signal of a simulated sample of the nonstationary white noise and the power spectrum of the corresponding stationary part are shown in Figures 3 and 4, respectively. The displacement responses of the system obtained through Newmark’s method [12] are shown in Figure 5. By examining the Fourier spectra associated with each of the response channel, also as shown in Figure 5, we chose the response of the 3rd channel, , which contains rich overall frequency information, as the reference response to construct a nonstationary ARMA model of the system.

According to the theory presented in the previous sections, by evaluating temporal root-mean-square function in conjunction with curve-fitting technique, as shown in Figure 6, nonstationary problem may reduce to a stationary problem if we can extract the amplitude-modulating function from the nonstationary data, as shown in Figure 7. Therefore, we can follow the same procedures as those for stationary problems, and the original time-series algorithm can thus be performed. We then consider ARMA model to perform the modal-parameter identification of the 3-DOF structural system. The initial parameters of the ARMA model can then be estimated through the three-stage least-square method, as shown in Table 1, and the relatively precise ARMA parameters can further be determined by using the nonlinear least-square methods, as listed in Table 2. We construct a nonstationary time-series model by introducing the aforementioned temporal root-mean-square function and then extracting the AR parameters from the above constructed nonstationary time-series model. The modal parameters of a structural system can then be obtained from (15), as summarized in Table 3, which show that the errors in natural frequencies are less than 1% and the maximum error in damping ratios is about 12%.

In addition, we also consider ARMA model in conjunction with a linear combination of suitable base functions to perform the modal identification of the previous 3-DOF chain-model system. The 5-order polynomial functions of time are chosen instead of the unknown MA parameter function in the original ARMA algorithm. The initial AR parameters and the coefficients of known base functions can then be estimated through the three-stage least-square method, as shown in Table 4, and the aforementioned AR parameters and coefficients of base functions can further be relatively precisely determined by using the nonlinear least-square methods, as listed in Table 5. We extract the AR parameters from the above constructed nonstationary time-series model, and the modal parameters of a structural system can then be obtained from (15), as summarized in Table 6, which show that the errors in natural frequencies are less than 1% and the maximum error in damping ratios is about 2%.

From Tables 3 and 6, it is indicated that the errors in results of natural frequency and damping ratio identification through the nonstationary ARMA algorithm with polynomial base function (as summarized in Table 3) are relatively smaller than those with curve-fitting technique (as summarized in Table 6), especially for the identified damping ratios. It should be mentioned that the errors of identified damping ratios are somewhat larger than those of identified natural frequencies. This is due to the fact that the system response generally has lower sensitivity to damping ratios than to the natural frequencies. In general, the amplitude curve-fitting technique has relatively high computational efficiency when compared with base-function method; however, it is only available under the assumption of product model for nonstationary input process. The base-function ARMA algorithm is well-performed without the assumption of product model for nonstationary input process and treatment of curve-fitting technique for extracting envelope function. In addition, the frame of mathematical model of base-function ARMA algorithm is relatively adjustable and is therefore applicable to general nonstationary ambient vibration data.

To further verify the effectiveness of the proposed method for the higher DOF structural systems, as well as for the substitution of the exact input with the data disturbed by some noise corresponding to the experience, we consider a 20-DOF chain model with nonproportional damping excited by a vibration practically recorded at Sun Moon Lake on the 21st of September, 1999, when the Chi-Chi Earthquake of a moment magnitude 7.6 occurred in central Taiwan. The sampling interval and period of this seismic record are sec and sec, respectively. A sample of the seismic record, which serves as the excitation input acting on the 20th mass of the model, is shown in Figure 8. Assume that each mass is 1 kg and all spring constants are 600 N/m. The damping matrix of the system is assumed to be In this example, we use the practical nonstationary input acting on the twentieth mass of the chain model, and then the corresponding displacement responses obtained by Newmark’s method [12] are used for modal identification. The results of identification are summarized in Table 7, which shows that the modal identification in this case is satisfactory. It is also noted that the higher modes generally are not identified so accurately as the lower ones, because the contribution of higher modes to the system response is somewhat less than that of the lower modes.

Recall that the conventional ARMA algorithm is effective only for stationary process [6, 7]. Therefore, in order to convert the original nonstationary responses into quasistationary data, we need additional treatments, such as employing the curve-fitting technique, to extract the amplitude-modulating function from the original nonstationary data in the form of a product model [4]. In contrast, in the present method, by introducing the base functions to construct a nonstationary time-series model, the MA parameters in the nonstationary time series are treated as functions of time and may be represented by a linear combination of suitable base functions. Without using input data or any additional treatment of converting the original nonstationary data into the form of quasistationary vibration, we directly perform the modal identification from nonstationary response data through the nonstationary base-function ARMA algorithm in conjunction with the proposed modified least-square method. Therefore, we avoid a distortion in the modal parameters of identification induced by the error involved in the approximate quasistationary response obtained through curve-fitting technique.

In this paper, we developed the nonstationary base-function ARMA algorithm in conjunction with the proposed modified least-square method for modal identification from nonstationary ambient response data only. We also demonstrated the validity of these methods through numerical simulations without using the practical response data. From the vibration behavior of realistic ambient excitation and the theory of the proposed method, we know that the assumptions and approximations are consistent with the time-varying nature of ambient excitation in practice. Thus, the proposed methods are generally applicable in identifying the modal parameters of a structure from the identification results obtained through numerical simulations.

#### 6. Conclusions

An extended time-series algorithm is presented in this paper for modal identification from nonstationary ambient vibration data only. Nonstationary ARMA model with time-varying parameters is considered because of its capability in resolving general nonstationary problems. Without changing the AR parameters, we construct a nonstationary time-series model to effectively describe a wide variety of nonstationary behavior and then extract the modal parameters of a structural system. Furthermore, by introducing the base functions to construct a nonstationary time-series model, the MA parameters in the nonstationary time series are treated as functions of time and may be represented by a linear combination of suitable base functions. We can then solve the identification problem of time-varying parameters and extend the time-series algorithm to identify the modal parameters of a structural system subjected to nonstationary white noise. However, identification of the mode shapes is still a challenging problem to be resolved in the proposed method from only one single sample of response signal. Through numerical simulation, applicability and effectiveness of the proposed method of modal parameter identification from nonstationary ambient vibration data are demonstrated.

#### Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This research was supported in part by the National Science Council of the Republic of China under the Grant NSC-100-2212-E-006-114 as well as the financial funding provided by the National Synchrotron Research Center (NSRRC). Chang-Sheng Lin would like to thank every member in the Optics Group of Experimental Facility Division and Precision Mechanical Engineering Group of Instrumentation Development Division, NSRRC, for their kind assistance and friendship. The authors wish to thank the anonymous reviewers for their valuable comments and suggestions in revising this paper.

#### References

- F. Shen, M. Zheng, D. Feng Shi, and F. Xu, “Using the cross-correlation technique to extract modal parameters on response-only data,”
*Journal of Sound and Vibration*, vol. 259, no. 5, pp. 1163–1179, 2003. View at Publisher · View at Google Scholar · View at Scopus - W. Ren, W. Zatar, and I. E. Harik, “Ambient vibration-based seismic evaluation of a continuous girder bridge,”
*Engineering Structures*, vol. 26, no. 5, pp. 631–640, 2004. View at Publisher · View at Google Scholar · View at Scopus - C. Lin and D. Chiang, “Modal identification from nonstationary ambient response data using extended random decrement algorithm,”
*Computers and Structures*, vol. 119, pp. 104–114, 2013. View at Publisher · View at Google Scholar · View at Scopus - D. Chiang and C. Lin, “Identification of modal parameters from nonstationary ambient vibration data using correlation technique,”
*AIAA Journal*, vol. 46, no. 11, pp. 2752–2759, 2008. View at Publisher · View at Google Scholar · View at Scopus - G. E. P. Box, G. M. Jenkins, and G. C. Reinsel,
*Time Series Analysis: Forecasting and Control*, Prentice Hall, Englewood Cliffs, NJ, USA, 3rd edition, 1994. View at MathSciNet - S. M. Moore, J. C. S. Lai, and K. Shankar, “ARMAX modal parameter identification in the presence of unmeasured excitation-I: theoretical background,”
*Mechanical Systems and Signal Processing*, vol. 21, no. 4, pp. 1601–1615, 2007. View at Publisher · View at Google Scholar · View at Scopus - S. M. Moore, J. C. S. Lai, and K. Shankar, “ARMAX modal parameter identification in the presence of unmeasured excitation-II: numerical and experimental verification,”
*Mechanical Systems and Signal Processing*, vol. 21, no. 4, pp. 1616–1641, 2007. View at Publisher · View at Google Scholar · View at Scopus - S. M. Pandit and S. M. Wu,
*Time Series and System Analysis with Applications*, John Wiley & Sons, New York, NY, USA, 1983. - S. M. Pandit and N. P. Mehta, “Data dependent systems approach to modal analysis via state space,”
*Journal of Dynamic Systems, Measurement and Control*, vol. 107, no. 2, pp. 132–138, 1985. View at Publisher · View at Google Scholar · View at Scopus - S. M. Pandit and N. P. Mehta, “Data dependent systems approach to modal analysis Part 1. Theory,”
*Journal of Sound and Vibration*, vol. 122, no. 3, pp. 413–422, 1988. View at Publisher · View at Google Scholar · View at Scopus - M. Shinozuka, “Simulation of multivariate and multidimensional random processes,”
*Journal of the Acoustical Society of America*, vol. 49, no. 1, pp. 357–367, 1971. View at Publisher · View at Google Scholar - N. M. Newmark, “A method of computation for structural dynamics,”
*Journal of the Engineering Mechanics Division*, vol. 85, no. 3, pp. 67–94, 1959. View at Google Scholar