INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 39, 3985-3998 (1996) A TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS USING GAUSS POINT COLLOCATION B . W. GOLLEY Department of Civil Engineering, University College, University of NS W, Australian Defmce Force Academy, Canberra, ACT 2600, Australia SUMMARY When a cubic function is interpolated between the prescribed initial displacement and velocity and the exact displacement and velocity at the end of a time step for a single degree of freedom system, the error, or residual, in the governing equation is zero at a number of times. It is shown that for a general undamped system, in the limit as the time step approaches zero these times correspond to Gauss points. This observation is verified by considering a general collocation procedure in which the displacement in any time step is approximated as a cubic function of time, with two coefficients chosen to satisfy the displacement and velocity at the beginning of the time step with the other two coefficients being chosen to satisfy the governing differential equation at any two times. It is shown that optimum accuracy is obtained if these points are the Gauss points. Detailed expressions are then presented for this particular case, and stability of the algorithm is investigated showing that the procedure is conditionally stable. For time steps which are a small proportion of the least period of vibration of the structure, the algorithm is considered to be the most accurate possible procedure based on cubic approximation of the displacement. KEY WORDS: structural dynamics; time stepping; collocation; Gauss points 1. INTRODUCTION In structural dynamics, the equations of motion for a spatially discretized system are of the form MX+ CX+KX - P = 0 (1) where M, C and K are the mass, damping and stiffness matrices, x, x and x are acceleration, velocity and displacement vectors, P is a prescribed force vector and an overdot denotes differentiation with respect to time t . If the elements of M,C and K are constant and the damping is of a restricted type, the equations may be uncoupled into a series of equations, each being of the same form as the equation for a single degree of freedom system, namely mx + c x + kx - p = 0 in which m is the mass, c is the damping coefficient, k is the stiffness, x is the displacement and p is the prescribed force. In this paper, the solution of equation (2) is addressed, but extension to the solution of equation (1 ) is obvious. Many algorithms are available for the numerical integration of equation (2). A number of these are discussed by Hughes and Belytschko’ and Dokainish and S~bbaraj.’.~In single-step algorithms, the displacement and velocity are specified at time t =O and are estimated during some period 0 <t <At where At is the time step. Values of displacement and velocity are calculated at CCC 0029-598 1/96/233985-14 0 1996 by John Wiley & Sons, Ltd. Received 7 February I995 Revised 5 January 1996 3986 B. W.GOLLEY t = At, and become the initial values for the next time step. The order of errors in velocity and displacement during the time step are generally not important, the important errors for accuracy when large numbers of time steps are involved are the values at the end of the time step. Many, but by no means all, single-step algorithms approximate displacements within the time step by functions which are cubic in time. These include the general Newmark method! presented in cubic form by Burnett? the Wilson-8 method6 and A number of cubic displacement algorithms use cubic Hermite polynomials as shape f~nctions,~-'~ with the final recurrence relationships being determined using weighted residuals, Hamilton's principle or Hamilton's law of varying action. Zienkiewicz9 gives a detailed description of the weighted residual method, and considers a range of weight functions to rederive a number of well-known algorithms using this method. In this paper, a highly accurate single-step procedure using cubic displacement functions is presented. In the first section, a cubic function is interpolated between the known displacement and velocity at the beginning of the time step and the exact displacement and velocity at the end of the time step, obtained by solving equation (2). This cubic function satisfies the governing differential equation at discrete dimensionless times T = tfAt only, and an investigation shows 0 these dimensionless times correspond to Gauss points in the general that in the limit as At damped case. In a following section, an algorithm is developed in which the coefficients of a cubic displacement function are chosen to satisfy the initial displacement and velocity and to satisfy the governing differential equation at two general collocation points. The locations of these points are then selected to maximize accuracy at the end of the time step. It is shown that fifth order accuracy is achieved for both displacement and velocity at the end of the time step if the collocation points are selected as Gauss points. For the specific case of Gauss point collocation, formulas required to implement the method are presented. Continuous weight functions leading to identical results are also presented. The spectral radius obtained from the amplification matrix indicates that the algorithm is conditionally stable. While it is recognized that this algorithm does not satisfy three of the five attributes suggested by Hilber and Hughes13 for a time-stepping procedure to be competitive, it is presented as a method which is simple to understand and implement, and which is felt to be useful when high accuracy is warranted. 2. CUBIC INTERPOLATION BETWEEN EXACT END VALUES In this section, a cubic function is interpolated between the displacement xo and velocity i o which are prescribed at time t = 0 and the exact displacement X; and velocity x; at time t = At obtained by solving equation (2). Equation (2) may be written in the alternative form in which o = step O<t<At, = cf2mw and 3 = pfm. In the following, it is assumed that in the time m, may be approximated by a cubic $=TP (4) where T = [ l t t2 t 3 ] , P = [PO p1 p2 p31T and the coefficients PO, ..., p3 are determined from the prescribed force p and the mass m. In the following procedure, the coefficients PO,..., p 3 TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS 3987 are assumed to be independent of the time step At. If @ is inadequately approximated by a cubic function, At should be decreased. It should be noted that in most time-stepping algorithms, the force jj is approximated by a linear function within each step. The exact solution for the displacement, xe, may be written as where the functions X ( ( , wt),i = 1, ...,6 are readily determined. The exact velocity xe is obtained from xe by differentiation. The exact displacement and velocity at t = 0 and t = At may be written in matrix form where 6 = oat. The cubic x = Tb" where be = [b: b: be2 b;lT is now interpolated between the exact values of displacement and velocity at t = 0 and t = A t given by equation (6). Substitution of this cubic function x and its derivatives x and x into the left hand side of equation (3) gives a residual Re which may be expressed in the form Re = T[Rbe - P] where R= -02 25w 2 0 0 w2 4(0 6 0 0 w2 6(0 -o o o w 2 3988 B. W.GOLLEY or after some manipulation where z = tfAt and the coefficients A1,A2,...,f46 of matrix F are functions of 5 and oat. In general, Re is a cubic function of 7 , although there are cases where Re degenerates to a quadratic function of z. For given values of C, w, At, xo, . . . , p 3 , values of 7 for which Re = 0 may be determined by solving the general cubic equation. Depending on the coefficients of the cubic function, there are one, two or three real roots. If we now consider a converse approach, in which the displacement is again approximated by a cubic function, with coefficients chosen to satisfy the initial conditions xo and i o , and the governing differential equation at two times when Re = 0, if there are two or three such times, the exact displacement and velocity should be reproduced at 7 = 1. Such a procedure is of course impractical, as it relies on determining the exact solution to obtain a solution which is approximate except at the end points of the time step. However, the times z when Re = 0 are considered in some detail, as they reveal optimum collocation points used in the method discussed later. If there is only one time when Re = 0, the exact displacement and velocity at z = 1 cannot be reproduced using a cubic approximation and two-point collocation. Expanding the coefficients J i l , f i 2 , . ..,f4 as Taylor series in ant, Re may be written in the form xo ~ , x o-, , w PO p1 p2 -, -, - (1 - 67 + 6 ~ ~ ) 0 ( w A+t 0) ~( ~ A t ) ~5, # 0 0 2 w3 w4 (11) and (1 - 62 +w2r3 (4 w ’ w3’ E) 0 5 (2 - +6 ~ ~ ) O ( o A t ) ~ 97 + 1 0 ~ ~ ) 0 ( w A+t )o~( ~ A t ) ~5, = 0 (12) where w Hence as At 7 -+ 0, two cases emerge. (a) Residual is O(At)2. This occurs in the general damped case, i.e. when for the residual reduces to a quadratic, with zeros obtained by solving 1 - 67 + 6r2 = 0 5 # 0. The cubic in (16) TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS 3989 which has roots 71 = (1 - &/3)/2 x 0.21 1325, 72 = (1 + &/3)/2 0.788675 (17) which are Gauss points. (b) Residual is O(At)’. This occurs in undamped systems, i.e. when [ = 0, and only when xo = po = p2 = 0. Hence the case only arises under unusual conditions during the time-stepping procedure, but is considered for completeness. The times when the residual is zero are obtained by solving the cubic 2 - 99 + lor3 = 0 (18) Roots of this equation cannot be obtained in simple closed form, but can be expressed in trigonometric form as where L o s - 1 6 (20) Numerical values of these roots are r3 = -1.044720, 74 = 0-237016, r5 = 04307703 (21 1 where it should be noted that two of the roots are within the range of the time step and are reasonably close to the Gauss points. Graphs showing the times 7 when Re = 0 are shown in Figure l(a)-l(e) taking one variable at a time to be non-zero. Two values of the damping coefficient are considered, the undamped case (( = 0) and a damped case with [ = 0.1. In addition, horizontal dashed lines are shown on the graphs to indicate the Gauss points. Values are plotted for O G o A t G 8 . Referring to those graphs, it may be seen that limiting cases as At + 0, times 7 for which Re = 0 are as given by equations (17) or (21) depending on the conditions described earlier. It is also noted that for O<o At G 1 ‘0,the residual is zero at times very close to the Gauss points, and hence if a cubic interpolation function is adopted with coefficients selected to satisfy the initial conditions and the governing equation at the Gauss points, a very accurate solution would be expected for displacement and velocity at 7 = 1. The graphs also indicate that for certain values of o At there is only one time 7 when Re = 0, and hence there are time steps for which it is not possible to obtain exact displacements and velocities at the end of the time step with cubic interpolation and two-point collocation. 3. TIME STEPPING USING POINT COLLOCATION In this section, a procedure is investigated in which a cubic displacement function is assumed, with coefficients chosen to satisfy the initial conditions and the governing equation at two general collocation points. The locations of these points are selected to give maximum accuracy for displacement and velocity at the end of the time step. 3990 T B. W.GOLLEY Z L 0.0 -0.5- -05 -1.07 -1 0 - . -Gauss 0 =405 A = ot - - -5 5 6 i = 0.1 - . -Gauss pts -1.5 3 pts .1.5 (b) 1.5 1.0- 0.5- z 0.0 I I I I I 1 2 3 4 0 =@At 5 I 6 I 7 - - - 5 = 0.1 - - -Gauss pts - - -Gauss -1.0- pts -1.5 -1.5 (C) 1.5 1.0- 0.5- z 0.0 1 2 3 4 5 6 7 D = mAt -5-0 -1.0 - - -Gauss pts -1.5 (el Figure 1. Times when residual equals zero: (a) xo # 0 or po # 0; (b) .to # 0; (c) p~ # 0; (d) p2 # 0; ( e ) p3 # 0 It is assumed that x may be approximated by the cubic function x=Tb (22) where b = [bo bl bz b3IT and the coefficients bo, ..., b3 are to be determined. The velocity is determined by differentiating equation (22) giving X= Tb (23) TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS 399 1 The initial displacement xo and velocity x o and displacement XI and velocity x l at time 7 = 1 are obtained in terms of bo, ..., b3 by substituting t = 0 and t = At into equations (22) and (23) as appropriate. We thus have four equations containing six unknowns, bo, bl, 62, b3, XI and .?I. Two additional equations are obtained by setting the residual to zero at two collocation points tcl and tc2. The six equations may be collected in matrix form as [ I,: I: [ [;][%; = K l: + [PI K3I K32 where Kzl = 1 At 1 [o [ ] KZ2= (At)2 (At)3 2At 3(At)2 ] Solving for the unknowns B1, B2 and XI gives BI = Xo B2 = A'Xo + L* P and x1= A X ~ + L P where (34) (35) 3992 B. W.GOLLEY and L = K22KG1D3 Matrix A is the amplification matrix. 4. ERROR ORDER IN GENERAL CASE - i, may be expressed in The error in displacement, E = xe - x , and the error in velocity i = ie where the coefficients E I , ..., 66 are functions of [, o,t , At, tcl and tc2. The exact solution xe given by equation ( 5 ) is a transcendental function of t while the approximate solution, x , is a cubic function of t , and a rational function of At. Expanding xe as a Taylor series in t and x as a Taylor series in At, coefficients ~ i ,ii, i = 1, ..., 6 may be expressed as polynomials in t and At. Hence at the end of the time step, values of 61, ..., i6 may be determined in polynomial form. In Appendix I it is shown that in the general undamped case if the collocation points are chosen as the Gauss points zCl = tc,/At= (1 - h / 3 ) / 2 , tc2 = tc2/At= (1 + h/3)/2 (43 ) values of el, . . ., i 6 are O(At)’, giving accuracies for both displacement and velocity of that order. In the special case where ( = xo = po = p2 = 0, it is also shown in the appendix that choosing collocation points ~~1 and zc2 as any pair of solutions of the cubic 2 - 9 ~ + 10z3 = O (44) which has been considered previously gives displacements with errors 0(At)7and velocities with errors of O(At)6.As this is a special and unusual case, collocation at the Gauss points is appropriate for consistent accuracy. 3993 TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS 5 . GAUSS POINT COLLOCATION In this section, terms required in the application of the procedure with collocation points taken as the Gauss points are presented, and stability of the algorithm is investigated. Hence, collocation points given by equations (43) are considered. It is of interest to note that identical equations to the equations obtained by setting the residual, R ( t ) , to zero at the Gauss points are obtained using weight functions W l ( z ) and W2(z), where 5Js W,(t) = 1 - -( 1 9 - 30t + 84t2 - 56t3), W2(t) = 1 + $(1 - 30t + 8 4 -~5 6~t 3 ) (45) and setting The weight function Wl(t)is shown in Figure 2 for reference. The location of the Gauss point leading to the same weighted residual equation is shown on the same figure as a vertical line. The function W,(t)is a mirror image of W l ( z ) . Formulas for the elements of matrices A, A*, L and L* given in equations (37)-(40) determined using Gauss point collocation are given in Appendix 11. Expressions for the eigenvalues of the amplification matrix A from which the spectral radius, p, may be determined are also given in that appendix. A plot of the region where p< 1, indicating numerical stability of the algorithm, is shown in Figure 3. The procedure is conditionally stable, being stable in the undamped case if O < C3<3 or WG6. If 120.077, the algorithm is stable for OG6. To satisfy the condition that C3<3, the time step should be less than 3/2x x 0.477 times the undamped period. m< 6. EXAMPLE A simple system with w = 1 and p = 0 is considered to demonstrate the accuracy of the method. The displacement, XT, and velocity, XT, at time T = 12, which is about 1.9 cycles, are determined. The time is divided into N time steps. Using the amplification matrix A given in Appendix 11, XT 3.0 2.0 1 .o 0.o -1 .o Figure 2. Equivalent weight function 3994 9. W.GOLLEY 0.20 0.18 0.16 0.14 < 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0 5 9 12 1 5 20 25- 30 36 4 0 -z 0 45 ’ 0 Figure 3. Stability regions for algorithm Table I. Ratios of terms in amplification matrix: x 6.283 4 6 8 10 20 30 40 50 60 80 3.00 2.00 1.50 1.20 0.60 0.40 0.30 0.24 0.20 0.15 1.185039 0.981219 0.992653 0.996752 0-999777 0.999955 0.999986 0.999994 0.999997 0.999999 2.484906 1.079273 1.025865 1.010863 1.000706 1.000141 1.000044 1.OOOO18 1.000009 1.oooO03 5 = 0.0, Period 0~000000 1.01 1818 1.O10082 1.0051 12 1 .OW393 1~000080 1.000025 1.oooo 10 1~000005 1-000002 1.185039 0.981219 0.992653 0.996752 0.999778 0.999955 0.999986 0.999994 0.999997 0.999999 and i~ are obtained as where the superscripts ‘a’ refer to the solution being approximate. An exact amplification matrix relating displacements and velocities at t = 0 and t = 2‘ can be determined from equation ( 5 ) . The ratios rij = A$/A$,i,j = 1,2, where the superscript ‘e’ refers to the exact solution, are presented in Tables 1 and 2 for = 0 and ( = 0.1, respectively, for various values of N. Examination of the tables indicates that for small values of At, halving the time step reduces the error by a factor of 16 for all the ratios rij. As halving the time step doubles the number of time steps, this is consistent with an error of 0 ( A t ) 5for a single time step for displacement and velocity as previously discussed. 3995 TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS Table 11. Ratios of terms in amplification matrix. C = 0.1,Period = 6.252 N 4 6 8 10 20 30 40 50 60 80 At rI I 3.00 2.00 1.50 1.20 0.60 0.40 0.30 0.24 0.20 0.15 1.511500 1.020964 1.005829 1402311 1 .OOOI 41 1.000028 1400009 1 ~ooooo4 1~000002 1 ~000001 ~ I Z r2I 2.615383 1.091205 19029229 1 -012186 1 .OW786 1 .0001 56 1.000050 1 ~000020 1 .oooo10 1.000003 0.000000 1.023004 1.013394 1-006435 1 400474 1.000096 1~000031 1.oooo1 2 1.000006 1 ~000002 1-22 1.306505 1.021241 1406855 1-002870 1.000186 1 *000037 1400012 1*000005 1~000002 1~00000 1 7. CONCLUSIONS Time stepping using displacement functions cubic in time has been investigated. It has been shown that for general damped systems, maximum accuracy is achieved by satisfying the governing equation at two Gauss points. An implicit conditionally stable algorithm based on this observation has been presented, and has been shown to lead to the accurate solution of structural dynamics problems. APPENDIX I Optimum collocation points The following terms up to O(At)5 for [ # 0 are obtained by expanding the error coefficients .,&6 as polynomials . ~ 1 , .. + 0(At)5 EI = &w4(l - 4(2)41(At)4 El = &w4(l - 4c2)$2(Af)3 - &w5(($3 - 4(2$4)(At)4 ~2 = ao4[(l- 2c2)41(At)4 &2 = +4((1 + 0(At)5 - 2c2)$2(At)3 + & 0 5 ( 4 5 (48) + O(At)’ (49) (50) - 4c2& + 16c444)(At)4+ O(At)’ (51 1 E3 = -El (52) E3 = -i, (53) + E4 = - i ~ ~ ( 4 2 ( A t+) &w5(-95 ~ + 41244)(At)4+ O(At)5 ~5 = h 0 1 ~ 4 1 ( A t + ) ~0(At)5 = - i ~ ~ c 4 4 ( A+ t 0(At)5 )~ (57) &6 = 0(At )5 (58) ~4 &j = - k ~ ~ [ 4 1 ( A t ) ~O(At)5 + = f ~ ’ 4 5 ( A t ) ~O(At)5 (54) (55) (56) (59) 3996 B. W.GOLLEY where 41= 1 - ~ ( G I + tc2) + 6 ~ 1 ~ 2 42 = 2 - ~ ( G4I G2) + 6~17tc2 5(?%1-k rf2) - 2'k1zc2 -k 10~cl'k2(~cl -k '49) = 1 - 3(& & ) 6 L ~ ' k 2 ( %f~k2) = 1 - 2 ( 4 -k & ) - ~ T C I ~ ~C ' ~~ I ~ c ~ (Tc2) ZCI = 3 - g ( ~ f 1 42) - ~ Z C I T C ~1 6 t c 1 ~ (-k~t c12 ) 43 =2 44 45 46 + + + + + (60) (61) (62) (63) (64) (65) and z,1 = ?,,/At and z,2 = tc2/At. The collocation points may now be selected to achieve maximum accuracy in the limiting case as At + 0. Selecting ~~1 and ~~2 such that 41 and 42 are zero would be expected to lead to displacements being accurate to 0 ( A t ) 5 and velocities being accurate to 0 ( A t ) 4 at the end of the time step. However, setting 41 = 42 = 0 and solving equations (60) and (61) gives t , ~= (1 - ~ / 3 ) / 2and z,2 = (1 6 / 3 ) / 2 and substitution of these values in equations (62)-(65) gives $3 = 4 4 = 45 = 4 6 = 0. Hence collocation at these times also gives accuracies for velocity of 0 ( A t ) 5 . Higher accuracy is attainable in certain cases, which are considered for completeness. Thus we consider in more detail the undamped case (C = 0) associated with f o , p1 and p3. In these cases, we find + + i 2 = & ~ J 4 5 ( A t+ )~ 0(At)6 E2 = &jw547(At)' 0 ( A t ) 7 E4 = -E2 i4 = -62 =&~~47(A -I-to)( A ~ t)7 &,j = ~ 0 ' 4 5 ( A t 4) ~0 ( A t ) 6 &j where 47 = 3 - lo(?:, + &) - 102C15C2 + 30?,1?,2(~Cl + zc2) (72) In these special cases, collocation points can be selected such that 4s = 4 7 = 0. The non-linear equations can be solved giving solutions for ~~1 and t c 2 which are any pair of the three solutions of the cubic 2-9t+10t3=0 (73) which has been considered previously. Selection of pairs of collocation points according to this solution then gives displacements with errors of O(At)7and velocities with errors of O(At)6. APPENDIX I1 Matrix terms and eigenvalues The terms of the matrices required for time stepping when the collocation points are the Gauss points, together with the eigenvalues from which the spectral radius may be determined are given 3997 TIME-STEPPING PROCEDURE FOR STRUCTURAL DYNAMICS below. d A=-[ + 6W2(a2- 12cw - 36) (is2- 36)(is2 - 12)At 48G2(is2 - 9)/At -6is2(502 A* = d 3 6 d + 6is( G2 - 36)( 0 + 41) + 36(0 + 36)/(At)2 ~is~+ (41)/(At)3 + 16(is2 + 72c2is + 72()/At 6G2(G2+ 1216- 12(1 -4[2))/(At)2 -6w(is3 36ga 1 + 2)(at)3 -6( (is2 + 2 4 ( ~ 5 + 3 6 ) ( A t ) ~ - (is2+ 18(~5+24)/(At)~ -6(5G2 + 3 6 ( 0 + 3 6 ) 6is(is+65)At 1 L* = - - -6((0 ] (75) a2- 12(a - 36)(At)’ 2 ( a 2 +36(6+72)(Ar)’ 3(is2+20(a+36)(At)4 -36is(a++()/At - -6(is2 - 12) + 3)At + ~ ) ( A z ) ~ ( G2 + 2450 + 60)(At)2- ( is2 - 36)(At)2 (74) -48( is2- 9)At --6(is2 - 12(6- 36)(At)2 L=- 1 24((6 where + a = 7G4 - 48(4 - 3C2)G2 432 and /.?= 48G2(a2- 9)( a4- 48( 1 - C2)is2 + 432( 1 - C2)) The spectral radius, p. is the maximum absolute value of 21 and 12. (77) 3998 B. W. GOLLEY REFERENCES 1. T. J. R. Hughes and T. Belytschko, ‘A precis of developments in computational methods for transient analysis’, J. Appl. Mech. ASME, 50, 1033-1041 (1983). 2. M. A. Dokainish and K. Subbaraj, ‘A survey of direct time-integration methods in computational structural dynamics-I. Explicit methods’, Comp. Slrucr., 32, 1371-1386 (1989). 3. K. Subbaraj and M. A. Dokainish, ‘A survey of direct time-integration methods in computational structural dynamics11. Implicit methods’, Comp. Struct., 32, 1387-1401 (1989). 4. N. M. Newmark, ‘A method of computation for structural dynamics’, J. Eng. Mech. Diu. ASCE, 85, (EM3) 67-94 (1959). 5. D. S . Bumett, Finite Element Analysis from Concepts to Applications, Addison-Wesley, Reading, MA, 1987. 6. K. J. Bathe and E. L. Wilson, ‘Stability and accuracy analysis of direct integration methods’, Earthquake eng. struct. dyn., 1, 283-291 (1973). 7. C . Hoff and P. J. Pahl, ‘Development of an implicit method with numerical dissipation from a generalised single-step algorithm for structural dynamics’, Comp. Methods Appl. Mech. Eng., 67, 367-385 (1988). 8. D. Karamanlidis, ‘On the equivalence of alternative finite element schemes in the time domain’, Eng. Anal., 3, 4 5 4 0 (1986). 9. 0.C. Zienkiewicz, The Finite Element Method, 3rd edn, McGraw-Hill, London, 1977. 10. G. F. Howard and J. E. T. Penny, ‘The accuracy and stability of time domain finite element solutions’, J. Sound Vib., 61, 585-595 (1978). 11. R. Riff and M. Baruch, ‘Time finite element discretization of Hamilton’s law of varying action’, AIAA J., 22, 13101318 (1983). 12. T. E. Simkins, ‘Finite elements for initial value problems in dynamics’, AZAA J., 19, 1357-1362 (1981). 13. H. M. Hilber and T. J. R. Hughes, ‘Collocation, dissipation and ‘overshoot’ for time integration schemes in structural dynamics’, Eurthquake eng. struct. dyn., 6, 99-1 17 (1978).