Third order complex-time-step methods for transient analysis

Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

www.elsevier.com/locate/cma

Third order complex-time-step methods for transient analysis T.C. Fung * School of Civil and Structural Engineering, Nanyang Technological University, Blk N1, No. 1A-8, Nanyang Avenue, Singapore 639 798, Singapore Received 9 December 1998; received in revised form 15 February 2000

Abstract Using complex-time-steps and the h-method for the ®rst order differential equations, third order accurate A-stable step-by-step time integration algorithms with controllable numerical dissipation are constructed for transient analysis. The results at two sub-step locations are evaluated and combined linearly with the corresponding weighting factors to give higher order accurate results at the end of a time-step. The sub-step locations and the weighting factors are complex numbers. The ampli®cation factors at in®nity are controlled by another algorithmic parameter l. A fourth order accurate algorithm is obtained when l ˆ 1. The algorithmic characteristics and the computational procedure for the present third order algorithms are discussed. Equivalent algorithms and comparisons with other algorithms are given. The present algorithms are found to be competitive with the other third and fourth order accurate algorithms. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction It is common to employ numerical techniques to evaluate approximate solutions of a given linear timedependent problem. The ®nite element method is widely used for the spatial discretization. For parabolic problems (e.g. transient heat conduction, di€usion, ground water ¯ow, etc.), the resultant semi-discretized equations are ®rst order di€erential equations in time and can be written as _ ‰CŠfu…t†g ‡ ‰KŠfu…t†g ˆ fF…t†g

…1†

with initial condition fu…0†g ˆ fu0 g. For hyperbolic problems (e.g. wave transmission in ¯uids, dynamic vibration of structures, etc.), the resultant semi-discretized equations are second order di€erential equations in time and can be written as _ ‰MŠf u…t†g ‡ ‰CŠfu…t†g ‡ ‰KŠfu…t†g ˆ fF…t†g _ with initial condition fu…0†g ˆ fu0 g and fu…0†g ˆ fv0 g: Eq. (2) can be reduced to an equivalent ®rst order system as        d u…t† 0 I u…t† 0 ‡ ˆ ; 1 1 1 M K M C v…t† M F…t† dt v…t† _ where fv…t†g ˆ fu…t†g: *

Tel.: +65-790-5305; fax: +65-791-0676. E-mail address: [email protected] (T.C. Fung).

0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 0 ) 0 0 2 7 0 - X

…2†

…3†

2790

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

Eqs. (1) and (3) can be written in the following standard form: _ fy…t†g ‡ ‰BŠfy…t†g ˆ ff…t†g

…4†

with appropriate initial condition fy…0†g ˆ fy0 g: The numerical solutions of Eq. (4) can be obtained by using step-by-step time integration methods. The time integration algorithms with desirable characteristics can be constructed via ®nite di€erence methods, weighted residual methods, Taylor series collocation methods, HamiltonÕs principle, HamiltonÕs law or least-squares methods (see, for example, books by Hairer et al. [1], Hairer and Wanner [2], Hughes [3], Wood [4], and Zienkiewicz and Taylor [5]). 1.1. h-Method [4] The h-method is perhaps the most widely used time integration algorithm for the numerical solutions of ®rst order systems of equations. It is also called the generalized trapezoidal rule and is an SS11 algorithm [5]. The recurrence relation for the h-method relating the approximate solution fyn‡1 g at t ˆ tn‡1 to fyn g at t ˆ tn for Eq. (4) can be written as …‰IŠ ‡ hDt‰BŠ†fyn‡1 g ˆ …‰IŠ

…1

h†Dt‰BŠ†fyn g ‡ Dt…hff n‡1 g ‡ …1

h†ff n g†;

…5†

where h is an algorithmic parameter, Dt ˆ tn‡1 tn is the time-step, [I] is an identity matrix, [B] is de®ned in Eq. (4) and ff n g ˆ ff…tn †g: When h ˆ 0; 1; 2=3 and 1/2, the h-method reduces to the well-known Euler method (or the forward difference method), the backward difference method, the Galerkin method and the Crank±Nicolson method, respectively. The h-method is second order accurate only if h is set to 1/2 and ®rst order accurate otherwise. 1.2. Higher order algorithms Higher order accurate algorithms give very accurate numerical results and are good for long-term prediction of system responses and preservation of system invariant (such as energy and momentum). Fast and accurate algorithms are useful for dynamical design, analysis and control of mechanical and structural systems. Higher order algorithms based on the Pade approximations to the exponential function have been widely used in solving di€erential equations [6,7]. However, powers of the matrices appear in the higher order formulations. This makes the computational e€ort very high. It is also well known that the diagonal and ®rst sub-diagonal Pade approximations are A-stable and L-stable, respectively. 1.3. Outline of the paper In this paper, third order accurate A-stable step-by-step time integration algorithms with controllable numerical dissipation are constructed for transient analysis using the h-method and sub-steps in complex form. As shown in Section 3, the results at two sub-step locations are evaluated and combined linearly with the corresponding weighting factors to give higher order accurate results at the end of a time-step. The substep locations and the weighting factors are complex numbers. The ampli®cation factors at in®nity are controlled by another algorithmic parameter l. A fourth order accurate algorithm is obtained when l ˆ 1. The algorithmic characteristics and the computational procedure for the present third order algorithms are discussed in Sections 4 and 5. Equivalent algorithms and comparisons with other algorithms are given in Sections 6 and 7. The present algorithms are found to be competitive with the other third and fourth order accurate algorithms.

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

2791

2. Single-degree-of-freedom systems The semi-discretized matrix di€erential equations in Eq. (4) are coupled in general. It would be cumbersome and dicult to study the characteristics of an algorithm by applying it to Eq. (4) directly. It has been rigorously established that the integration of the uncoupled equations is equivalent to the integration of the original system [3,5]. It is therefore more convenient and sucient for the purpose of investigating the characteristics of a proposed algorithm by considering the equation of motion of a single-degree-offreedom system in the form y_ …t† ‡ ky…t† ˆ f …t†;

…6†

where k may be complex. The conclusions in stability and accuracy are valid for the original multi-degreeof-freedom system as well. The exact solution of Eq. (6) can be written as Z t y…t† ˆ y0 exp … kt† ‡ exp … k…t s††f …s† ds: …7† 0

The ®rst expression in y(t) is related to the given initial condition y0 at t ˆ 0. It is known as the homogeneous solution. All the numerical methods are trying to approximate the homogeneous solution in various ways. Poor approximations would give unstable numerical results. 3. Formulation of third order algorithms When the h-method is applied to Eq. (6), the recurrence relation becomes yn‡1 ˆ

1 ‡ …1 h†q Dt……1 h†fn ‡ hfn‡1 † yn ‡ ˆ A…h; q†yn ‡ L…h; q; fn ; fn‡1 †; 1 hq 1 hq

…8†

where q ˆ k Dt, A and L are the ampli®cation factor and load operator respectively. It can be seen that A is an approximation to the homogeneous solution exp …q†: The Taylor series expansion of A and exp…q† can be written as: A…h; q† ˆ 1 ‡ q ‡ hq2 ‡ h2 q3 ‡


‡ hk 1 qk ‡


ˆ 1 ‡

1 X

hj 1 qj ;

…9a†

jˆ1

exp …q† ˆ 1 ‡ q ‡

1 X 1 2 1 3 1 q ‡ q ‡


‡ qk ‡


ˆ qj =j!: 2! 3! k! jˆ0

…9b†

It can be seen that if h 6ˆ 1=2, the truncation error is O(Dt2 ) and hence the algorithm is ®rst order accurate only. To reduce the truncation error, an sth order ampli®cation factor As …h; q† with truncation error O…Dts‡1 † is obtained by linearly combining the ampli®cation factors in Eq. (9a) at different sub-step location bj Dt with weighting factor aj , i.e. X aj A…h; bj q† As …h; q† ˆ j

ˆ

X

aj ‡

X

aj bj q ‡

X

aj b2j hq2 ‡

X

aj b3j h2 q3 ‡


‡

X

aj bkj hk 1 qk ‡


;

…10†

where aj and bj are algorithmic parameters to be determined. A similar procedure has been adopted by Fung [8±10] to construct unconditionally stable higher order accurate algorithms for second order equations. The algorithmic parameters are chosen so that the Taylor series expansions of Eqs. (10) and (9b) match for the ®rst (s + 1) terms, i.e. from Dt0 and Dts . Comparing the Taylor series expansions, it can be shown that the s + 1 conditions are

2792

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

X

ai ˆ 1

and hk

1

X

…11a†

aj bkj ˆ 1=k! for k ˆ 1; . . . ; s:

…11b†

3.1. Ampli®cation factor at in®nity It can be seen that for b0 ˆ 0, A…h; 0† ˆ 1 and no evaluation is required. Therefore, it is at no cost to assume one of the sub-step locations to be at the beginning of the time interval, i.e. b0 ˆ 0. The corresponding weighting factor a0 can be related to the asymptotic value of the ampli®cation factor at in®nity (i.e. Dt ! 1). Using Eq. (11a), the asymptotic value is given by  X  1 1 …1 a0 †; aj 1 …12† lim As ˆ a0 ‡ ˆ1 Dt!1 h h jˆ1 which is real and independent of k. For stable algorithms, the magnitude of this value must be less than unity. If the value in Eq. (12) is de®ned as l, a0 can be expressed in terms of l as a0 ˆ 1

h…1

l†:

…13†

3.2. Algorithmic parameters Since b0 ˆ 0 and a0 is given by Eq. (13), four more parameters (a1 , b1 , a2 and b2 ) are required to satisfy four equations (Eqs. (11a), and (11b) with s ˆ 1; 2 and 3) in order to construct third order accurate algorithms. In the following, h ˆ 1=2 is assumed. It can be veri®ed that the algorithmic parameters are given explicitly as: 1 a0 ˆ …1 ‡ l†; 2

a1 ˆ

1

l 4

4 ‡ 7l ‡ l2 ‡ i p ; 4 2 ‡ 2l l2

a2 ˆ

1

l 4

4 ‡ 7l ‡ l2 i p ; 4 2 ‡ 2l l2

…14a†

p p 2 ‡ 2l l2 2 ‡ 2l l2 2‡l 2‡l i ‡i ; b2 ˆ ; …14b† b0 ˆ 0; b1 ˆ 3…1 ‡ l† 3…1 ‡ l† 3…1 ‡ l† 3…1 ‡ l† p where i ˆ 1. It is interesting to note that b1 and b2 are in fact the roots of the quadratic equation 3…1 ‡ l†x2 …4 ‡ 2l†x ‡ 2 ˆ 0. It can be seen that the algorithmic parameters a1 and a2 , and b1 and b2 are complex conjugate pairs when 0 6 l 6 1. The resultant ampli®cation factor in Eq. (10) is given by A3 …q† ˆ

lq2 ‡ 2…2l ‡ 1†q ‡ 6…l ‡ 1† : q2 2…l ‡ 2†q ‡ 6…l ‡ 1†

…15†

Furthermore, it can be shown that the ampli®cation factor A3 would be the same for other values of h (6ˆ 0) [11]. 4. Algorithmic characteristics 4.1. Parabolic equations Consider the transient heat conduction problems. k is real and positive. The ampli®cation factor A3 in Eq. (15) is real as well. The ampli®cation factors for various l values against k Dt are plotted in Fig. 1. It can be seen that for p 1 6 l 6 1, the magnitude of the ampli®cation factors would not exceed unity. Besides, A3 …q† P 0 if l P … 3 1†=2. In this case, the sign of the amplitude of the high-frequency response would

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

2793

Fig. 1. Comparison of ampli®cation factor with third order complex-time-step algorithms.

not change from step to step. When l ˆ 0, A3 ! 0 as k Dt ! 1. The oscillation limit is k Dt ˆ 3, i.e. the solution could produce oscillations if k Dt P 3. When l ˆ 1, A3 ! 1 as k Dt ! 1. For other values of l, as expected, A3 ! l as k Dt !p1.  The minimum values of A3 for l ˆ 0 and 1 are )0.0981 and 0.0718, respectively. When l ˆ 7 4 3, the minimum value would be equal to )l so that the numerical dissipation for the high frequency range is maximized. In Fig. 1, the ampli®cation factors for the Euler method, Crank± Nicolson method, Galerkin method and backward difference method are shown for comparison. The present methods are obviously of higher order accuracy and agree better with the exact solution exp … k Dt†: 4.2. Hyperbolic equations Consider the undamped dynamic vibration problems. k ˆ ix is imaginary. It can be veri®ed that A3 in Eq. (15) is complex in general and is given by A3 …ix Dt† ˆ

lx4 Dt4

2…7 ‡ 16l ‡ 7l2 †x2 Dt2 ‡ 36…1 ‡ l†

2

x4 Dt4 ‡ 4…1 ‡ l ‡ l2 †x2 Dt2 ‡ 36…1 ‡ l†2 ‡i

2…1 ‡ 4l ‡ l2 †x3 Dt3 ‡ 36…1 ‡ l†2 x Dt x4 Dt4 ‡ 4…1 ‡ l ‡ l2 †x2 Dt2 ‡ 36…1 ‡ l†

2

:

…16†

The magnitude of the ampli®cation factors is the spectral radius of the ampli®cation matrix for the dynamic vibration problems. It can be veri®ed that 2

jA3 …ix Dt†j ˆ ‰A3 …ix Dt†  A3 … ix Dt†Š ˆ ˆ1

…1

l2 x4 Dt4 ‡ 4…1 ‡ l ‡ l2 †x2 Dt2 ‡ 36…1 ‡ l† x4 Dt4 ‡ 4…1 ‡ l ‡ l2 †x2 Dt2 ‡ 36…1 ‡ l†

l2 †x4 Dt4

x4 Dt4 ‡ 4…1 ‡ l ‡ l2 †x2 Dt2 ‡ 36…1 ‡ l†

2

:

2

2

…17†

The spectral radii of A3 for various l values against Dt/T are shown in Fig. 2, where T ˆ 2p=x. It can be seen that the spectral radii stays close to the unity level for a while and then decrease monotonically to the ultimate spectral radius l. In Fig. 2, the spectral radii for the unconditionally stable second order accurate

2794

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

Fig. 2. Comparison of spectral radii with third order complex-time-step algorithms (n ˆ 0).

Hilber±Hughes±Taylor (HHT) method, the Houbolt method and the Park method are shown for comparison. The present third order complex-time-step algorithms show good numerical dissipation property as they stay close to the unity level longer. Any desirable value of ultimate spectral radius can be obtained by setting l directly. For dynamic vibration problems, the numerical errors related to dissipation and dispersion are usually measured in terms of the algorithmic damping ratio and relative period error [3]. If A3 …ix Dt† is expressed in the form P ‡ iQ, where P and Q are real as in Eq. (16), it can be shown that the algorithmic damping ratio and relative period error are given by: Algorithmic damping ratio ˆ

Relative period error ˆ

1 ln…P 2 ‡ Q2 † 1 l x3 Dt3 ‡ O…Dt5 †; ˆ 2 tan 1 …Q=P † 72…1 ‡ l†

x Dt tan 1 …Q=P †

1ˆ

2

l ‡ 2l2

540…1 ‡ l†

2

x4 Dt4 ‡ O…Dt6 †:

…18a†

…18b†

The algorithmic damping ratio and relative period error for various l values against Dt/T are shown in Figs. 3 and 4, respectively. It can be seen that the performance error is very different for various values of Dt and obviously deteriorates at large values of it. Such values in a real multi-variable problem correspond of course to the `high-frequency' responses which are less important. For smooth solutions, algorithms with numerical dissipation in the high-frequency range are preferred. From these ®gures, the third order characteristics of the present complex-time-step algorithms are observed. In Fig. 4, it can be seen that the relative period error for the present algorithms is better than the second order Newmark method (or the trapezoidal rule), which has the smallest relative period error among all unconditionally stable second order accurate linear multi-step algorithms.

5. Computational procedure From Eq. (15), an implementation similar to the direct application of Pade approximations can be constructed as

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

2795

Fig. 3. Comparison of algorithmic damping ratios with third order complex-time-step algorithms (n ˆ 0).

Fig. 4. Comparison of relative period errors with third order complex-time-step algorithms (n ˆ 0).

2

2

…‰BŠ Dt2 ‡ 2…l ‡ 2†‰BŠDt ‡ 6…l ‡ 1†‰IŠ†fyn‡1 g ˆ …l‰BŠ Dt2

2…2l ‡ 1†‰BŠDt ‡ 6…l ‡ 1†‰IŠ†fyn g:

…19†

In the formulation, the computation of ‰BŠ2 ˆ ‰CŠ 1 ‰KŠ‰CŠ 1 ‰KŠ is required. This is very computationally intensive. Furthermore, the computation would destroy the band and symmetry properties of the original [C] and [K] matrices. Alternatively, for the present third order complex-time-step algorithms, the solution fX ‡ iYg at substep location b1 can be obtained from ‰Kr ‡ iKi ŠfX ‡ iYg ˆ fFr ‡ iFi g;

…20†

2796

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

where Kr , Ki , Fr and Fi are real and for parabolic problems, they are given as ‰Kr Š ˆ ‰CŠ ‡

2‡l ‰KŠDt; 6…1 ‡ l†

…21a†

p 2 ‡ 2l l2 ‰KŠDt; 6…1 ‡ l†

‰Ki Š ˆ

 fFr g ˆ

‰CŠ

…21b†

 2‡l ‰KŠDt fyn g; 6…1 ‡ l†

…21c†

p 2 ‡ 2l l2 ‰KŠDtfyn g: fFi g ˆ 6…1 ‡ l†

…21d† 1

It can be seen that the matrix ‰BŠ ˆ ‰CŠ ‰KŠ is not required. The solution corresponding to b2 is simply fX iYg and therefore no additional computation is required. The result at t ˆ tn‡1 is then given by fyn‡1 g ˆ a0 fyn g ‡ a1 fX ‡ iYg ‡ a2 fX

iYg ˆ

1‡l 1 l fyn g ‡ fXg 2 2

4 ‡ 7l ‡ l2 p fYg: 2 2 ‡ 2l l2

…22†

The computational procedure in Eq. (20) is very similar to the h-method except complex numbers are used. Eq. (20) can also be solved in terms of real number calculations by rewriting it as      Kr Ki X Fr ˆ : …23† Ki Kr Fi Y fyn‡1 g can then be obtained from Eq. (22) as well. Eq. (23) needs more storage and does not possess symmetry property. Besides, nowadays many computer languages already support the complex number data type implementation. Programming using complex number has no additional diculties. 6. Equivalence to other algorithms Eq. (15) treats many existing algorithms as special cases. 6.1. Pad e approximations The …1; 2† and …2; 2† Pade approximations to exp …q† are given by: R12 ˆ

R22 ˆ

6 ‡ 2q ; 4q ‡ q2

…24a†

12 ‡ 6q ‡ q2 ; 12 6q ‡ q2

…24b†

6

respectively. It can be veri®ed that Eqs. (24a) and (24b) can be obtained by substituting l ˆ 0 and l ˆ 1 into Eq. (15). Hence, the Pade approximations are reproduced by the present complex-time-step methods. 6.2. Lambert approximation Lambert [6] gave an expression approximating exp …q† in terms of two parameters a and b as RLambert ˆ

1 ‡ 12 …1 a†q ‡ 14 …b a†q2 : 1 12 …1 ‡ a†q ‡ 14 …b ‡ a†q2

…25†

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

2797

The order of accuracy of the approximation is 2 in general, 3 if b ˆ 1=3, a 6ˆ 0 and 4 if a ˆ 0 and b ˆ 1=3. The Lambert third order approximation with b ˆ 1=3 can be obtained from A3 …q† in Eq. (15) if l is set to (1 ) 3a)/(1 + 3a). 6.3. Calahan approximation p When l ˆ 1 3, it can be shown that A3 (q) in Eq. (15) is equivalent to the Calahan [12] approximation p or the Nùrsett [13] approximation with (in his notation) p ˆ 2 and ap ˆ …3 ‡ 3†=6. However, this approximation p cannot be constructed from the complex-time-step method as a1 and a2 are unde®ned when lˆ1 3 in Eq. (14a) even though A3 is well de®ned in Eq. (15). 6.4. Argyris±Vaz±Willam approximation Argyris et al. [14] proposed a series of approximations in terms of the values of the function and its derivatives at the end of the time interval. The ®rst member is equivalent to the h-method. The second member can be written as RArgyris ˆ

6

2…f 2†q ‡ …1 f†q2 ; 6 2…1 ‡ f†q ‡ fq2

…26†

where f represents the collocation position. It is order 3 in general and order 4 when f ˆ 1=2. This approximation can be obtained from A3 (q) in Eq. (15) by setting l ˆ …1 f†=f: 6.5. Weighted residual method Unconditionally stable higher order accurate time-step integration algorithms can be constructed from the weighted residual method [15]. The third and fourth order accurate algorithms for Eq. (1) can be written as      Y1 C W1 K Dt Kyn Dt ; …27† ˆ K Dt 2C ‡ W2 K Dt 0 Y2 where W1 ˆ

…1=3†…1=…1 ‡ l††, W2 ˆ …2=3†……2 ‡ l†=…1 ‡ l†† and fyn‡1 g at t ˆ tn‡1 is given by

fyn‡1 g ˆ fyn g ‡ fY1 g ‡ fY2 g:

…28†

Eq. (27) is also equivalent to many existing algorithms, for example, the …1; 2† and …2; 2† Pade approximations [16±18], the fully implicit Runge±Kutta methods [1,2], time-discontinuous Galerkin method [19] or bi-discontinuous time Galerkin method [20]. It can be shown that the numerical results given by Eqs. (28) and (22) are equivalent and corresponding to the generalized Pade approximations [15]. Besides, the computational e€ort to solve Eq. (27) is comparable to Eqs. (23) and (20) for third order accurate algorithms. 7. Comparison with other algorithms 7.1. Tarnow and Simo's method Tarnow and Simo [21] proposed a procedure to construct a fourth order accurate algorithm from an underlying second order accurate algorithm by stepping forward a Dt, backward (1 ) 2a)Dt and then forward a Dt again to advance a time-step Dt. The value for a is found to be 2 1 1 a ˆ ‡ 22=3 ‡ 21=3 ˆ 1:351207 . . . 3 6 3 The algorithm is suitable for Hamiltonian (non-dissipative) systems.

…29†

2798

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

Fig. 5. Ampli®cation factors for Tarnow and Simo's and Nùrsett methods.

When the procedure employs the Crank±Nicolson method as the underlying second order accurate algorithm, the ampli®cation factor for the resultant fourth order accurate algorithm is given by   2 …2 ‡ …1 2a†q†…2 ‡ aq† 1 1 2=3 5 1=3 5 2 ˆ exp …q† …30† ATS ˆ ‡ 2 ‡ q ‡ O…q6 †: 45 72 288 …2 …1 2a†q†…2 aq†2 It can be shown that ATS is only I-stable and not A-stable (hence only suitable for non-dissipative Hamiltonian systems). In particular, it is not even A0 -stable as the magnitude of the ampli®cation factor is unde®ned when k Dt ˆ 2=…2a 1† and the magnitude of ATS would exceed unity for real k when pp 2 p 8 ‡ 3  22=3 ‡ 4  21=3 …31† 2 3 4 ‡ 22=3 ‡ 2  21=3 6 k Dt 6 p 5 or 1:13442 6 k Dt 6 1:20063: Fig. 5 shows the comparison between ATS and exp()k Dt). It can be seen that the two functions does not agree very well when k Dt P 0:5: Hence, this fourth order algorithm would not be very useful for parabolic problems. When the procedure is applied to an undamped vibration problem, the resultant algorithm is unconditionally stable and non-dissipative (i.e. the spectral radius is unity and the algorithmic damping ratio is zero for all Dt/T). However, as shown in Fig. 4, the relative period error is found to be bigger than those for the present third order methods with various l values. 7.2. Nrsett method Nùrsett [13] proposed to construct rational approximations to the exponential function with denomip nators having a single linear factor, i.e. in the form …1 ‡ ap k Dt† : In other words, the numerical results for Eq. (4) at the end of a time-step can be obtained by evaluating {Yj }, for j from 1 to p using ‰I ‡ ap DtBŠfYj g ˆ ap Dt‰BŠfYj 1 g;

…32†

with fY0 g ˆ fyn g; and fyn‡1 g is given by fyn‡1 g ˆ fyn g ‡

p X cj fYj g; jˆ1

…33†

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

2799

where ap and cj are algorithmic parameters. ap and cj are real numbers and are related to the Laguerre polynomials [4]. It can be shown that the maximum order of accuracy achievable by the Nùrsett method is p + 1 only [2,4]. Besides, A-stable algorithms exist only for p ˆ 1; 2; 3 and 5. The computational procedure is found to be very similar to the present complex-time-step methods. However, the present methods achieve higher order accuracy and remain A-stable by combining the results at various complex sub-step locations which can be evaluated independently. It can be veri®ed that the algorithmic parameters for the third and fourth order accurate unconditionally stable Nùrsett algorithms are given by p p p p ˆ 2 : a2 ˆ …3 ‡ 3†=6; c1 ˆ 3 ‡ 3; c2 ˆ 3 2 3 …34a† and pˆ3:

a3 ˆ 1:068579021;

c1 ˆ

0:9358222279;

c2 ˆ

0:4979406069;

c3 ˆ

0:1966521039; …34b†

respectively. The corresponding ampli®cation factors are p p   36 12 3q 6…1 ‡ 3†q2 1 1 p 4 p 2 N2 ˆ 3 q ‡ O…q5 †; ˆ exp…q† ‡ 24 36 …6 …3 ‡ 3†q† N3 ˆ

1

…35a†

2:205737063q ‡ 0:7198463093q2 ‡ 0:7692126656q3 …1

1:068579021q†3

ˆ exp…q† ‡ 0:1643928967q5 ‡ O…q6 †:

…36a†

The computational e€orts for the third and fourth order algorithms are about only half and threequarter of the present complex-time-step method. However, it can be seen from Fig. 5 that the accuracy of N2 and N3 is not as high when compared to the present third order algorithms. The algorithmic damping ratios and relative period errors for N2 and N3 for hyperbolic problems are shown in Figs. 3 and 4, respectively as well for comparison. It can be seen that the accuracy is not as good as those for the present third order algorithms. Hence, the present third order algorithms are more competitive than the Nùrsett method.

8. Numerical examples 8.1. Example 1: transient heat conduction problem Consider a bar of unit length at uniform initial temperature subjected suddenly to zero temperature applied at one end. The governing equation is ou…x; t† ot

o2 u…x; t† ˆ0 ox2

for 0 < x < 1;

…37†

with boundary conditions u…0; t† ˆ 0; …ou=ot†…1; t† ˆ 0 and initial condition u…x; 0† ˆ 1. The exact solution is given by  2 !    1 X 2 1 1 u…x; t† ˆ exp k p2 t sin k px : …38† …k …1=2††p 2 2 kˆ1 Using 16 linear heat conduction elements, the maximum and minimum k are 3050 and 2.469, respectively. If the Euler method (which is only conditionally stable) is used, the maximum time-step allowed is

2800

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

Fig. 6. Numerical results for Example 1 using 16 elements.

6.557 ´ 10 4 to maintain numerical stability. The numerical results at x ˆ 1 for various methods with different time-steps are shown in Fig. 6. From Fig. 6, it can be seen that for the Crank±Nicolson method with Dt ˆ 0:5, the solutions are oscillating and even giving negative temperature. It is well known that the oscillation can be reduced by using a larger value of h or reducing the time step. Using the present higher order methods, good results can be obtained by just using Dt ˆ 0:5. In Fig. 6, it can be seen that Tarnow and SimoÕs fourth order algorithm also gives oscillation results when Dt ˆ 0:5. Reasonable results are obtained when Dt ˆ 0:25. This can be explained by studying the corresponding value of kmin Dt. When Dt ˆ 0:5, kmin Dt ˆ 1:2345 and the corresponding ampli®cation factor from Fig. 5 is negative and hence the oscillating results. Furthermore, it can be veri®ed that when Dt ˆ 0:47, unstable numerical results are obtained since kmin Dt ˆ 1:1604 is within the unstable range. As a result, it can be seen that the use of Tarnow and SimoÕs method should be restricted to non-dissipative Hamiltonian systems. 8.2. Example 2: free vibration of a single-degree-of-freedom system Consider a single-degree-of-freedom system governed by  u…t† ‡ u…t† ˆ 0;

…39†

with initial conditions u0 ˆ 1 and v0 ˆ 0. The period of the system is 2p. A time-step of about one-tenth of the period of the system is usually recommended for second order accurate algorithms, i.e. Dt ˆ 0:6 is usually suggested. Fig. 7 shows the calculated results by the Newmark method (or the trapezoidal rule) with Dt ˆ 0:5. Results are not very accurate. Using the fourth order algorithm, accurate results can still be obtained with Dt ˆ 1. The results are better than the trapezoidal rule with Dt ˆ 0:25! The accuracy of the present higher order accurate algorithms is therefore demonstrated.

9. Conclusions Using complex-time-steps and the h-method for the ®rst order differential equations, third order accurate A-stable step-by-step time integration algorithms with controllable numerical dissipation are con-

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

2801

Fig. 7. Numerical results for Example 2.

structed for transient analysis. The results at two sub-step locations are evaluated and combined linearly with the corresponding weighting factors to give higher order accurate results at the end of a time-step. The sub-step locations and the weighting factors are complex numbers. The algorithmic characteristics and computational procedure for the third order algorithms are discussed. Equivalent algorithms and comparison with other higher order algorithms are given. It is found that many existing algorithms can be recovered by the present algorithms. The algorithmic characteristics for the present algorithms are shown to be better than the fourth order accurate Tarnow and SimoÕs method and the Nùrsett method. References [1] E. Hairer, S.P. Nùrsett, G. Wanner, Solving Ordinary Di€erential Equations I, Springer, Berlin, 1987. [2] E. Hairer, G. Wanner, Solving Ordinary Di€erential Equations II, Springer, Berlin, 1991. [3] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cli€s, NJ, 1987. [4] W.L. Wood, Practical Time-Stepping Schemes, Clarendon Press, Oxford, 1990. [5] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, 4th edn., McGraw-Hill, New York, 1991. [6] J.D. Lambert, Computational Methods in Ordinary Di€erential Equations, Wiley, New York, 1973. [7] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall, Englewood Cli€s, NJ, 1962. [8] T.C. Fung, Unconditionally stable higher-order Newmark methods by sub-stepping procedure, Computational Methods in Applied Mechanics and Engineering 147 (1997) 61±84. [9] T.C. Fung, Complex-time-step Newmark methods with controllable numerical dissipation, International Journal for Numerical Methods in Engineering 41 (1998) 65±93. [10] T.C. Fung, Higher order time-step integration methods with complex time steps, Journal of Sound and Vibration 201 (1) (1998) 69±89. [11] T.C. Fung, Complex-time-step methods for transient analysis, International Journal for Numerical Methods in Engineering 46 (1999) 1253±1271. [12] D.A. Calahan, A stable accurate method of numerical integration for non-linear systems, Proceedings of IEEE 56 (1968) 744. [13] S.P. Nùrsett, One-step methods of Hermite type for numerical solution of sti€ systems, BIT 14 (1974) 63±77. [14] J.H. Argyris, L.E. Vaz, K.J. Willam, Higher order methods for transient di€usion analysis, Computational Methods in Applied Mechanics and Engineering 12 (1977) 243±278. [15] T.C. Fung, Weighting parameters for unconditionally stable higher-order accurate time step integration algorithms: Part 1. First order equations, International Journal for Numerical Methods in Engineering 45 (1999) 941±970. [16] P.W. M oller, High-order hierarchical A- and L-stable integration methods, International Journal for Numerical Methods in Engineering 36 (1993) 2607±2624.

2802

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 190 (2001) 2789±2802

[17] M. Gellert, A new algorithm for integration of dynamic systems, Computer and Structures 9 (1978) 401±408. [18] T.C. Fung, Unconditionally stable higher order accurate Hermitian time ®nite elements, International Journal for Numerical Methods in Engineering 39 (1996) 3475±3495. [19] G.M. Hulbert, A uni®ed set of single-step asymptotic annihilation algorithms for structural dynamics, Computer Methods in Applied Mechanics and Engineering 113 (1994) 1±9. [20] M. Borri, C. Bottasso, A general framework for interpreting time ®nite element formulations, Computational Mechanics 13 (1993) 133±142. [21] N. Tarnow, J.C. Simo, How to render second order accurate time-stepping algorithms fourth order accurate while retaining the stability and conservation properties, Computer Methods in Applied Mechanics and Engineering 115 (1994) 233±252.