Stability and accuracy of differential quadrature method in solving dynamic problems

Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

www.elsevier.com/locate/cma

Stability and accuracy of di€erential quadrature method in solving dynamic problems T.C. Fung

*

School of Civil and Structural Engineering, Nanyang Technological University, Blk N1, No. 1C-88, Nanyang Avenue, Singapore 639798, Singapore Received 1 May 2001

Abstract In this paper, the di€erential quadrature method is used to solve dynamic problems governed by second-order ordinary di€erential equations in time. The Legendre, Radau, Chebyshev, Chebyshev±Gauss±Lobatto and uniformly spaced sampling grid points are considered. Besides, two approaches using the conventional and modi®ed di€erential quadrature rules to impose the initial conditions are also investigated. The stability and accuracy properties are studied by evaluating the spectral radii and truncation errors of the resultant numerical ampli®cation matrices. It is found that higher-order accurate solutions can be obtained at the end of a time step if the Gauss and Radau sampling grid points are used. However, the conventional approach to impose the initial conditions in general only gives conditionally stable time step integration algorithms. Unconditionally stable algorithms can be obtained if the modi®ed di€erential quadrature rule is used. Unfortunately, the commonly used Chebyshev±Gauss±Lobatto sampling grid points would not generate unconditionally stable algorithms. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Higher order algorithms; Unconditionally stable time step integration algorithms; Stability of time step integration algorithms; In¯uence of sampling grid points; Modi®ed di€erential quadrature rules

1. Introduction Recently, there is an increasing interest in applying the di€erential quadrature method to solve dynamic problems. Since the governing di€erential equations are second order in time with two initial conditions given at the initial time, the application of the di€erential quadrature method is not as straight-forward as the ®rst-order ordinary di€erential equations. The main concerns are (1) how to incorporate the initial/ boundary conditions and (2) how to select the sampling grid points so that the resultant algorithms would be accurate and stable. 1.1. Approaches to impose the initial/boundary conditions It is well known that the imposition of the initial/boundary conditions for the di€erential quadrature method can be very tricky [1]. Conventionally, the boundary conditions are expressed as di€erential quadrature analog equations at the sampling grid points on or near the boundaries. These analog equations are then used to replace the di€erential quadrature analog equations of the governing di€erential equations at these points in order to solve the problems. This procedure becomes very tedious when higher-order derivatives or multi-boundary conditions are involved along the boundaries. *

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

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

1312

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Bert et al. [2] and Jang et al. [3] proposed a d-technique to impose the ®rst derivative boundary conditions. In this approach, one boundary condition is exactly discretized while the other is approximately discretized by applying it at a small distance d away from the boundary edge. Hence, d must be small ( 10 5 in dimensionless value) for accurate solutions. However, this procedure may also create ill-conditioned weighting coecient matrices and unexpected oscillation behavior of the solutions [1]. On the other hand, the higher-order derivatives initial/boundary conditions can also be imposed exactly by modifying the trial functions to incorporate the degrees of freedom of the higher-order derivatives at the boundary [4] or by using the di€erential quadrature element method directly [5]. For initial value problems, Wu and Liu [6,7] have employed the Hermite interpolation functions to incorporate the initial conditions into the approximate solutions. It has been shown [8] that this approach is in fact equivalent to the conventional approach by using the di€erential quadrature analog equations of the initial conditions at the initial time. Besides, a modi®ed di€erential quadrature rule has also been presented in [8] to impose the initial conditions. The ®rst derivatives at the sampling grid points are related to the function values at all the sampling grid points and the given function value at the initial time, as in the conventional di€erential quadrature method. For the second derivatives, however, they are related to the values of the ®rst derivatives at the sampling grid points and the ®rst derivative given at the initial time. This approach is di€erent from the conventional approach where the second derivatives are obtained by di€erentiating the trial functions twice. Alternatively, Wang and Bert [9], Wang et al. [10] and Malik and Bert [11] had proposed several methods to incorporate the boundary conditions by modifying the elements in the weighting coecient matrices. However, as reported by Shu and Du [12], these techniques have some major limitations and cannot be used to tackle general boundary conditions. For initial value problems, Tanaka and Chen [13] have recently used this approach to incorporate the given initial conditions for transient responses for elastodynamic problems. In this paper, it is shown that this approach is in fact equivalent to the modi®ed di€erential quadrature rule presented in [8]. Furthermore, it is well known that the second-order equations can also be expressed as an equivalent system of ®rst-order equations in time and then solved by the di€erential quadrature method. This is called an indirect method. It is shown in [8] that the numerical results given by the modi®ed di€erential quadrature rule would be equivalent to the indirect method but without doubling the number of unknowns. Hence, the accuracy and stability properties of the numerical solutions obtained by using the modi®ed di€erential quadrature rule are related to those for the ®rst-order equations in [14]. 1.2. Choices of sampling grid points Obviously, the accuracy, stability and rate of convergence of the numerical solutions depend on the choices of the sampling grid points. It is well known that the equally spaced sampling grid points are not desirable (e.g. [1]). It has been suggested that non-uniformly spaced sampling grid points could give better results. The zeros of some orthogonal polynomials are commonly adopted as the sampling grid points. It is not surprising since the solutions of di€erential equations are closely related to quadrature where the particular signi®cance of orthogonal polynomials is well known [15]. In fact, Bellman et al. [16] had proposed to use the zeros of the Legendre polynomials as the sampling grid points in one of his introductory papers. Quan and Chang [17,18] had also studied the use of the zeros of orthogonal polynomials on many typical chemical engineering problems. The performances of using the zeros of the Chebyshev polynomials were consistently better than the other choices of orthogonal polynomials considered in their work. On the other hand, it was also reported that the zeros of Legendre polynomials could deliver more accurate results than the zeros of Chebyshev polynomials for interaction problems of a plate and an elastic half space [1]. However, these Legendre and Chebyshev sampling grid points do not include the end points of the domain under consideration. This makes the imposition of initial/boundary conditions quite dicult. As a result, the Chebyshev±Gauss±Lobatto sampling grid points are commonly used. It was reported [1] that the Chebyshev±Gauss±Lobatto sampling grid points performed consistently better than the equally spaced,

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1313

Legendre and Chebyshev sampling grid points in a variety of problems. They have been demonstrated to enhance the accuracy of the quadrature solutions and can tackle the problem of deterioration of the quadrature solutions [17±19] as well. Furthermore, the Chebyshev±Gauss±Lobatto sampling grid points had been used in [6,7,13] to solve dynamic problems. For initial value problems, if the whole time domain of interest is discretized simultaneously (for example, as in [1, Example 6] and [7,20]), many unknowns have to be solved simultaneously. As a result, it is more natural to adopt a step-by-step time integration approach to advance the solutions progressively over the time domain of interest [8]. Hence, it is important to ensure that the resultant time step integration algorithms are numerically stable, i.e. small initial errors would not grow with the time steps. Besides, since the results at the end of a time step will be used as the initial conditions for the next time step, it is desirable to have the results at the end of the time step as accurate as possible. Recently, Fung [8] has found that higher-order accurate and unconditionally stable algorithms can be obtained if the Legendre points (obtained by setting l ˆ 1 in [8]) or the Radau points (obtained by setting l ˆ 0 in [8]) are used with the modi®ed di€erential quadrature rule to impose the initial conditions. However, up to now, there is no systematic study on how the choices of the sampling grid points and the ways to impose the initial conditions would a€ect the accuracy and stability of the numerical solutions for dynamic problems. In this paper, special attentions will be focused on the spectral radii of the numerical ampli®cation matrices and the order of accuracy of the solutions at the end of a time step. Five types of sampling grid points and the two approaches using the conventional and modi®ed di€erential quadrature rules to impose the initial conditions are considered. If s0 ; s1 ; . . . ; sn are the sampling grid points considered in this paper, the ®ve types of sample grid points are: 1. Legendre sampling grid points s0 ˆ 0;

s1 ; . . . ; sn are the zeros of

dn n …s …s dsn

1†n † ˆ 0:

…1†

2. Radau sampling grid points s0 ˆ 0;

s1 ; . . . ; sn are the zeros of

dn 1 n 1 …s …s dsn 1

3. Chebyshev sampling grid points    1 2k 1 1 cos p s0 ˆ 0; sk ˆ 2 2n

for 1 6 k 6 n:

4. Chebyshev±Gauss±Lobatto sampling grid points    1 k 1 cos p for 0 6 k 6 n: sk ˆ 2 n

n

1† † ˆ 0:

…2†

…3†

…4†

5. Equally spaced sampling grid points k for 0 6 k 6 n: …5† n Note that s0 ˆ 0 for all the ®ve types of sampling grid points and sn ˆ 1 for the Radau, Chebyshev±Gauss± Lobatto and equally spaced sampling grid points. Furthermore, the di€erential quadrature analog equations of the di€erential equations will be established at s1 ; . . . ; sn only. The manuscript is organized as follows. In Section 2, the conventional approach to impose the initial conditions is reviewed. In Section 3, the modi®ed di€erential quadrature rule to impose the initial conditions is discussed. The approach to incorporate the initial conditions by modifying the weighting coecient matrices is shown to be equivalent to this approach. In Sections 4 and 5, the accuracy and stability properties of the di€erential quadrature method using various types of sampling grid points and approaches to impose the initial conditions are discussed. Conclusions are then given in Section 6. sk ˆ

1314

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

2. Conventional approaches to impose the initial conditions In this section, the conventional approaches to impose the given initial conditions are presented. Without loss of generality, one time step of length Dt is considered. 2.1. Imposing initial conditions by choosing suitable test/trial functions The initial conditions can be incorporated directly by employing suitable test/trial functions. If the initial conditions are given by du u…t ˆ 0† ˆ u0 and ˆ v0 ; …6† dt tˆ0 the approximate solution can be written as u…s† ˆ u0 Jun …s† ‡ v0 DtJvn …s† ‡ u1 J1n …s† ‡


‡ un Jnn …s†;

…7†

where s ˆ t=Dt is the non-dimensional time, Dt is the time step, u1 ; . . . ; un are unknowns to be determined. The functions Jun …s†, Jvn …s†, J1n …s†; . . . ; Jnn …s† can be chosen such that J_ un …s0 † ˆ 0;

Jun …s0 † ˆ 1;

Jun …si † ˆ Jvn …si † ˆ 0

Jvn …s0 † ˆ 0;

J_ vn …s0 † ˆ 1;

…8a†

for i ˆ 1; . . . ; n;

…8b†

for i ˆ 0; 1; . . . ; n; k ˆ 1; . . . ; n and J_ kn …s0 † ˆ 0;

Jkn …si † ˆ dik

k ˆ 1; . . . ; n;

…8c†

where dik ˆ 1 if i ˆ k and 0 otherwise, s0 ˆ 0 and s1 ; . . . ; sn are the sampling grid points within a time step and an over dot denotes the di€erentiation with respect to s. In this case, u1 ; . . . ; un can be interpreted as the approximate values of u…s† at the sampling grid points s1 ; . . . ; sn . If the functions Jun …s†, Jvn …s† and Jkn …s† are constructed from polynomials of degree n ‡ 1, they can be expressed as Jun …s† ˆ …1

L_ n0 …s0 †…s

Jvn …s† ˆ …s

s0 †Ln0 …s†;

…9b†

s0 n L …s† for k ˆ 1; . . . ; n; s0 k

…9c†

Jkn …s† ˆ

s sk

s0 ††Ln0 …s†;

…9a†

where Lnk …s† is the nth order Lagrange polynomial and is given by Lnk …s† ˆ

n Y s s jˆ0 k

sj : sj

…10†

j6ˆk

The di€erential quadrature rules for the ®rst and Eq. (7) directly and can be written as 8 …1† 9 2 …1† 8 9 8 …1† 9 > > > u_ 1 > > B B11 B > > > > > 1u 1v = = < < = < 6 . .. . . 6 u0 ‡ v0 Dt ‡ 4 .. ˆ .. .. . > > > > > > ; > > > : > ; ; : : …1† …1† …1† u_ n Bnv Bn1 Bnu 8 …2† 9 2 …2† 8 9 8 …2† 9 > > >  u1 > > B B11 B > > > > > 1u 1v = = < < = < 6 . .. . . 6 u0 ‡ v0 Dt ‡ 4 .. ˆ .. .. . > > > > > > ; > > > : > ; ; : : …2† …2† …2†  un Bnv Bn1 Bnu

second derivatives can be obtained by di€erentiating













3 9 …1† 8 u1 > B1n > = < .. 7 7 .. ; . . 5> ; : > …1† un Bnn 3 9 …2† 8 u1 > B1n > = < .. 7 7 .. ; . . 5> ; : > …2† un Bnn

…11a†

…11b†

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

where …r†

Biu ˆ

dr n J …s† ; dsr u sˆsi

…r†

Biv ˆ

dr n J …s† ; dsr v sˆsi

…r†

Bij ˆ

1315

dr n J …s† : dsr j sˆsi

…12† …r†

…r†

…r†

A more detailed discussion on how to compute the weighting coecients Biu , Biv and Bij explicitly can be found in [8]. 2.2. Imposing initial conditions by di€erential quadrature rules Another conventional approach is to impose the initial conditions by establishing the appropriate differential quadrature analog equations for the initial conditions by using the di€erential quadrature rules. If the approximate solution is assumed to be a polynomial of degree n ‡ 1 as in Eqs. (7) and (9a)±(9c), then n‡1 n‡1 n‡1 u…s† ˆ u0 Ln‡1 0 …s† ‡ u1 L1 …s† ‡


‡ un Ln …s† ‡ un‡1 Ln‡1 …s†;

…13†

where Lkn‡1 …s† ˆ

n‡1 Y s s jˆ0 k

sj sj

for 0 6 k 6 n ‡ 1

…14†

j6ˆk

are the Lagrange polynomials with sampling grid points s0 …ˆ 0†, s1 ; . . . ; sn ; sn‡1 . The di€erential quadrature rules for the ®rst and second derivatives can be obtained by di€erentiating Eq. (13) and can be written as 3 9 2 …1† 8 …1† …1† A00 A01


A0;n‡1 8 u0 9 u_ 0 > > > > > > > > 7> …1† …1† …1† = 6 < u_ 1 > 6 A10 A11


A1;n‡1 7< u1 = 6 7 ˆ …15a† .. 7> .. >; .. .. . > 6 > > 5> . . > 4 .. ; > . > : . > ; : …1† …1† …1† u_ n‡1 un‡1 An‡1;0 An‡1;1


An‡1;n‡1 3 9 2 …2† 8 …2† …2† A00 A01


A0;n‡1 8 u0 9  u0 > > > > > > 6 …2† > > > 7> …2† …2† <  u1 = 6 A10 A11


A1;n‡1 7< u1 = 7 6 ˆ6 . …15b† .. 7> .. >; .. .. > . > . > > > 5> 4 .. . . > ; > : ; : …2† …2† …2†  un‡1 un‡1 An‡1;0 An‡1;1


An‡1;n‡1 where …r† Aij

dr n‡1 ˆ r Lj …s† : ds sˆsi

…16†

In Eq. (15a), u0 is given and u_ 0 ˆ v0 Dt. Hence, the ®rst equation in Eq. (15a) is v0 Dt ˆ u_ 0 ˆ

n‡1 X jˆ0

…1†

A0j uj :

From Eq. (17), un‡1 can be expressed in terms of u0 , v0 Dt, u1 ; . . . ; un as ! n X 1 …1† …1† A00 u0 ‡ v0 Dt A0j uj : un‡1 ˆ …1† A0;n‡1 jˆ1

…17†

…18†

If the analog equations of the governing di€erential equations are to be established at s1 ; . . . ; sn only, u_ n‡1 and  un‡1 are not required in setting up the equations. Using Eq. (18), u_ 1 . . . ; u_ n in Eq. (15a) and u1 . . . ; un in Eq. (15b) can be expressed in terms of u0 ; v0 Dt; u1 ; . . . ; un after eliminating un‡1 . It has been shown [8] that the relations are independent of the actual value of sn‡1 and are in fact given by Eqs. (11a) and (11b) as well.

1316

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Once u1 ; . . . ; un are solved, un‡1 can be obtained from Eq. (18). u_ 0 …ˆ v0 Dt†, u_ 1 , u_ 2 . . . ; u_ n‡1 can then be obtained from Eq. (15a). The interpolated solutions for u…s† and v…s† are then given by u…s† ˆ

n‡1 X kˆ0

v…s† ˆ

Lkn‡1 …s†uk ;

…19a†

n‡1 n‡1 1 1 X 1 X _ L_ kn‡1 …s†uk ˆ Ln‡1 …s†u_ k : u…s† ˆ Dt Dt kˆ0 Dt kˆ0 k

…19b†

In other words, u…s† and v…s† are polynomials of degree n ‡ 1 and n, respectively, and are related by differentiation. 3. Modi®ed di€erential quadrature rules As proposed in [8], u…s† and v…s† can be both approximated by polynomials of degree n with sampling grid points at s0 …ˆ 0†, s1 ; . . . ; sn as u…s† ˆ u0 Ln0 …s† ‡ u1 Ln1 …s† ‡


‡ un Lnn …s†; v…s†Dt ˆ

v0 DtLn0 …s†

‡

u_ 1 Ln1 …s†

‡


‡

u_ n Lnn …s†;

…20a† …20b†

where Lnk …s† ˆ

n Y s s jˆ0 k

sj : sj

…21†

j6ˆk

The di€erential quadrature rules for the ®rst and second derivatives at the n sampling grid points can be related to the function values at the n sampling grid points and the initial conditions at s0 by 8 9 8 …1† 9 2 …1† 38 9 …1† > A11


A1n > = > = = < u_ 1 > < A10 > < u1 > 6 . 7 .. .. . . _ ˆ fG0 gu0 ‡ ‰GŠfug; ˆ u or fug ‡ …22a† . . . 4 5 0 . : . > : > > > . ; > .…1† ; > . ; : …1† …1† u_ n un An1


Ann An0 8 9 8 9 8 9  u > u_ > u > > > > > = = = < 1> < 1> < 1> .. .. . ˆ fG0 gv0 Dt ‡ ‰GŠ . ˆ fG0 gv0 Dt ‡ ‰GŠfG0 gu0 ‡ ‰GŠ‰GŠ .. or . > > > > > > …22b† > > ; > > : : ; : ;  un un u_ n 2

f ug ˆ fG0 gv0 Dt ‡ ‰GŠfG0 gu0 ‡ ‰GŠ fug; where …1† Aij

d n L …s† ˆ : ds j sˆsi

…23†

As before, if the n unknowns u1 ; u2 ; . . . ; un in fug are solved, u_ 1 ; u_ 2 ; . . . ; u_ n can be obtained from Eq. (22a). The interpolated solutions are then given by u…s† ˆ

n X kˆ0

v…s† ˆ

Lnk …s†uk ;

n 1 X Ln …s†u_ k ; Dt kˆ0 k

where u_ 0 ˆ v0 Dt is used in Eq. (24b) and is not obtained from the di€erential quadrature rule.

…24a† …24b†

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1317

Note that in general, v…s†Dt 6ˆ du…s†=ds as both u…s† and v…s† in Eqs. (24a) and (24b) are polynomials of _ ˆ 0† 6ˆ v0 Dt in general. However, it can be veri®ed that the degree n. Moreover, while v…s ˆ 0† ˆ v0 ; u…s order of accuracy at the initial time is comparable to the order of accuracy at the other locations with the time interval [8]. 3.1. Imposing initial conditions by modifying the weighting coecient matrices Recently, Tanaka and Chen [13] have proposed to incorporate the initial conditions by modifying the weighting coecient matrices. Instead of using the original variable u…s† with non-homogenous initial conditions, a new variable z…s† de®ned as z…s† ˆ u…s†

u0

v0 Dts

…25†

z_ …0† ˆ 0:

…26†

is introduced so that z…0† ˆ 0 and

If z…s† is assumed to be a polynomial of degree n as z…s† ˆ z0 Ln0 …s† ‡ z1 Ln1 …s† ‡


‡ zn Lnn …s†

…27a†

the di€erential quadrature rule relating z_ 0 , z_ 1 ; . . . ; z_ n to z0 ; z1 ; . . . ; zn is given by 3 9 8 9 2 …1† …1† …1† 8 A00 A01


A0n > z0 > z_ 0 > > > > > > > 6 …1† > > = …1† …1† 7> < z_ 1 6 A10 A11


A1n 7< z1 = 7 6 ˆ .. 6 .. .. 7> .. >: .. > . > . > > > 4 . . 5> . > > > ; : ; : > …1† …1† z_ n zn A A


A…1† n0

n1

…28†

nn

In the conventional approach, 8 9 2 …1† A00 z0 > > > > > > < z1 = 6 6 A…1† 10 ˆ6 .. 6 .. > > . > 4 . > > > : ; …1† zn An0

…1†

A01 …1† A11 .. . …1† An1












3 9 …1† 2 8 A0n > z0 > > > > …1† 7 > A1n 7 < z1 = 7 . .. 7 > . >: . > . 5> > ; : > …1† zn A

…29†

nn

However, since z0 ˆ 0, following the idea in [9], it was proposed to replace all the elements in the ®rst column by zero, i.e. 8 9 2 0 z_ 0 > > > > > 6 > > = 60 < z_ 1 > .. > ˆ 6 6. > > . > 4 .. > > > > ; : z_ n 0

A01

…1†






A11 .. . …1† An1

…1†









38 9 z > > > 0> > …1† 7> < z A1n 7 1 = 7 . : .. 7> . > . > . 5> > ; : > …1† zn A …1†

A0n

…30†

nn

Similarly, the di€erential quadrature rule relating z0 , z1 ; . . . ; zn to z0 ; z1 ; . . . ; zn is proposed as 8 9 2 0 z0 > > > > > > 6 > > < z1 = 6 0 6 .. > ˆ 6 > 6. > . > > 4 .. > > ; : > zn 0

A01

…1†






A11 .. . …1† An1

…1†









38 9 2 0 > z_ 0 > > > > 7> 6 > …1† 7> A1n 7< z_ 1 = 6 60 ˆ 7 . .. 7> . > 6 6. > 4 .. . 5> > . > > > ; : z_ n 0 A…1† A0n

…1†

A01

…1†






…1†






nn

A11 .. . …1† An1






32 8 9 > z0 > > > > 7> > …1† = < z1 > A1n 7 7 7 . .. 7 > . >: > . 5> > . > > ; : > …1† zn A …1†

A0n

nn

…31†

1318

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

In their formulation, since the analog equation of the governing di€erential equation at s0 is not established, z_ 0 and z0 will not be used. From, Eq. (30) 8 9 2 …1† 38 9 …1† A11


A1n > > = = < z_ 1 > < z1 > 6 . 7 .. .. . ˆ …32† .. 5 . : 4 .. . > > > ; ; : > : …1† …1† z_ n zn An1


Ann Note that Eq. (32) can be obtained from Eq. (28) directly even without modifying the weighting coecient matrix since z0 ˆ 0 anyway. From Eq. (25), 8 9 8 9 8 9 8 9 8 9 8 9 8 9 u_ > z_ > > s1 > z1 > > u 1 > > 1 1> > > > > > = = > = < > < = < 1> < 1> < = < = < = .. .. .. .. .. .. .. ˆ v ; v ˆ u ; Dt Dt 0 0 0 . . > > > > > > : > > > ; > ; > . > > . > > . ; > . ; > . ; : : > : : ; ; : : 1 1 zn un sn z_ n u_ n 8 9 8 9  z > > u > > > = > = < 1> < 1> .. .. ˆ : …33† . . > > > > > ; > ; : > : >  zn un Hence, Eq. (32) becomes 8 9 8 9 2 …1† A11 > = <1= < u_ 1 > .. 6 . .. v ˆ Dt .. 4 0 . :.; > ; : > …1† 1 u_ n A n1

Note that 2 …1† A11 6 . 4 .. …1† An1









3 …1† 8 9 A1n < 1 = .. .. 7 . 5: . ; ˆ 1 A…1† nn









308 9 …1† > A1n = < u1 > B .. .. 7 @ 5 . . > ; : > …1† un Ann

8 …1† 9 > = < A10 > .. > ; : .…1† > An0

Hence, Eq. (34) is equivalent to Eq. (22a). Similarly, from Eq. (31), 8 9 2 …1† 32 …1† …1† A11


A1n > A11


= < z1 > 6 . 6 . .. .. 7 ˆ . .. 4 5 4 . . . > ; : > …1† …1† …1† zn An1


Ann An1




2 and

8 9 <1= . u0 .. : ; 1

…1†

A11 6 . 4 ..

…1† An1









8 91 > = < s1 > . C v0 Dt .. A: > ; : > sn

38 9 8 9 …1† A1n > = <1= < s1 > . .. .. 7 5 . > ˆ : .. ;: . > : ; …1† 1 sn Ann

38 9 …1† > z1 > A1n < = .. .. 7 5 . . > ; : > …1† zn Ann

and making use of Eq. (33), Eq. (36) becomes 8 9 2 …1† 8 91 32 …1† 308 9 8 9 …1† …1† u1 > A11


A1n > > > A11


A1n = = = <1= < < u1 > < s1 > .. B .. 6 . 6 . .. .. C .. 7 .. 7 u ˆ v Dt . . @ 4 5 4 5 0 0 . . > . >A . . . . :.; > > > ; ; ; : > : : …1† …1† 1  un un sn An1


A…1† An1


A…1† nn nn 8 …1† 9 2 …1† 38 9 2 …1† 38 …1† 9 …1† 2 …1† > A11


A1n > A11


A1n > = = = < u1 > < A10 > < A10 > 6 . 7 6 . 7 .. . . . .. ˆ 4 .. ‡ u ‡ v : Dt .. 5 . . . 4 5 0 0 . : . > > . ; > .…1† > > ; ; : : .…1† > …1† …1† …1† …1† un An1


Ann An1


Ann An0 An0

…34†

…35†

…36†

…37†

Hence, Eq. (37) is equivalent to Eq. (22b) as well. In other words, the approach to incorporate initial conditions by modifying the weighting coecient matrices is in fact equivalent to the modi®ed di€erential quadrature rule proposed in [8]. The advantage of the modi®ed di€erential quadrature rule is that the interpolation of v…s† is also speci®ed explicitly.

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1319

Alternatively, the weighting coecient matrices can also be modi®ed di€erently. For example, following the idea presented in [11], since z_ 0 ˆ 0, the ®rst row of the weighting coecient matrix could be modi®ed to zero, i.e. 8 9 2 3 9 0 0


0 8 z_ 0 > > > > > z0 > …1† …1† …1† > = 6 A10 A11


A1n 7> = < z_ 1 > < z1 > 6 7 ˆ6 . …38† 7 .. . . . .. > .. 5> .. > 4 .. . > > > > > > ; > : : ; …1† …1† z_ n zn An0 An1


A…1† nn and hence z0 , z1 ; . . . ; zn are 8 9 2 …1† …1† A00 A01 z0 > > > > > > < z1 = 6 A…1† A…1† 6 11 ˆ 6 10 .. .. > > 4 ... . > > . > ; : > …1† …1† zn An0 An1

related to z0 , z1 ; . . . ; zn by 9 2 …1† …1† 38 A00


A0n > z_ 0 > > > > …1† 7> …1† = < 6


A1n 7 z_ 1 6 A10 ˆ 7 6 . .. 5> . > 4 .. . > . > . > ; : > …1† …1† z _ n An0


Ann

…1†









…1†






A01 …1† A11 .. . An1

…1† 32 0 A0n …1† 76 A…1† A1n 76 10 6 .. .. 7 . 54 .…1† …1† An0 Ann

0 …1† A11 .. . …1†

An1












3 9 0 8 z0 > …1† > > = < > A1n 7 7 z1 .. 7 .. : > . 5> > ; : . > …1† zn Ann …39†

It can be shown that z1 ; . . . ; zn in Eq. (39) are related to z1 ; . . . ; zn by Eq. (36) as well. 4. Accuracy of the approximate solutions To investigate the accuracy of the approximate solutions, a linear single-degree-of-freedom system is considered. The governing di€erential equation is given by d2 u du ‡ 2nx ‡ x2 u ˆ 0 2 dt dt with initial conditions u…t ˆ 0† ˆ u0

and

…40†

du ˆ v0: dt tˆ0

…41†

Eq. (40) can be normalized as  u ‡ 2fxDtu_ ‡ x2 Dt2 u ˆ 0:

…42†

Note that the time derivative of u is de®ned as v and is given by vˆ

du 1 du 1 _ ˆ ˆ u: dt Dt ds Dt

…43†

_ ˆ 0† ˆ v0 Dt. As a result, the initial conditions for Eq. (42) are u…s ˆ 0† ˆ u0 and u…s The accuracy of the approximate solutions can be obtained by comparing the approximate solutions u…s† and v…s† with the exact solutions u…s† and v…s†. If 0 6 n < 1, u…s†, v…s† are given by (e.g. [21])     nx 1 u…s† ˆ u0 exp… nxsDt† cos…xd sDt† ‡ sin…xd sDt† ‡ v0 exp… nxsDt† sin…xd sDt† ; …44a† xd xd q 1 d u and xd ˆ 1 n2 x: …44b† v…s† ˆ Dt ds 4.1. Imposing the initial conditions by the conventional approach It can be veri®ed that in general, if the initial conditions are imposed by the conventional approach presented in Section 2, then for all types of sampling grid points under consideration, u…s†

u…s† ˆ O…Dtn‡2 †

and

v…s†

v…s† ˆ O…Dtn‡1 †:

…45†

1320

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

In other words, the order of accuracy for the approximate solutions within a time interval is n only in general, as limited by the accuracy of the velocity. The order of accuracy for the approximate solutions at the end of a time step (i.e. s ˆ 1 or t ˆ Dt) can be improved if the Legendre or Radau sampling grid points are used. For the Legendre sampling grid points, u…s ˆ 1†

u…s ˆ 1† ˆ O…Dt2n‡2 †

and

v…s ˆ 1†

v…s ˆ 1† ˆ O…Dt2n‡1 †;

…46†

so that the order of accuracy for the approximate solutions at the end of a time step is 2n in general. Similarly, for the Radau sampling grid points, u…s ˆ 1†

u…s ˆ 1† ˆ O…Dt2n‡1 †

and

v…s ˆ 1†

v…s ˆ 1† ˆ O…Dt2n †;

so that the order of accuracy for the approximate solutions at the end of a time step is 2n

…47† 1 in general.

4.2. Imposing initial conditions by the modi®ed di€erential quadrature rule If the initial conditions are imposed by the modi®ed di€erential quadrature rule, it can be veri®ed that the order of accuracy for u…s† and v…s† are both n in general with u…s†

u…s† ˆ O…Dtn‡1 †

and

v…s†

v…s† ˆ O…Dtn‡1 †:

…48†

Furthermore, if the Legendre sampling grid points are used u…s ˆ 1†

u…s ˆ 1† ˆ O…Dt2n‡1 †

and

v…s ˆ 1†

v…s ˆ 1† ˆ O…Dt2n‡1 †

…49†

As a result, the order of accuracy for the approximate solutions at the end of a time step is 2n in general. Similarly, for the Radau sampling grid points, u…s ˆ 1†

u…s ˆ 1† ˆ O…Dt2n †

and

v…s ˆ 1†

v…s ˆ 1† ˆ O…Dt2n †:

…50†

As a result, the order of accuracy for the approximate solutions at the end of a time step is 2n 1 in general. In conclusion, it can be seen that by using n undetermined coef®cients with the conventional approach to impose the initial conditions, the orders of accuracy for displacement and velocity are n ‡ 1 and n, respectively, in general for all types of sampling grid points under consideration. On the other hand, by using n undetermined coef®cients with the modi®ed differential quadrature rules to impose the initial conditions, the orders of accuracy for displacement and velocity are both n, respectively, in general for all types of sampling grid points under consideration. Hence, the orders of accuracy for various types of sampling grid points and approaches to impose the initial conditions are comparable for numerical solutions within a time interval. However, it can be seen that the accuracy at the end of a time interval (i.e. s ˆ 1 or t ˆ Dt) could be improved if the Legendre or Radau sampling grid points are used. Hence, the accuracy of the approximate solution is higher if the Legendre or Radau sampling grid points are used. However, it will be shown in the next section that the algorithms are unconditionally stable only if the modi®ed di€erential quadrature rule is used to impose the initial conditions. Note that although the orders of accuracy discussed above are for linear equations only, it can be shown that the higher orders of accuracy are maintained for non-linear equations as well if the Legendre and Radau sampling grid points are used [22,23].

5. Stability of the di€erential quadrature algorithms To investigate the stability properties of the algorithms, Eq. (40) with n ˆ 0 is considered. The stability property can be determined from the spectral radius of the numerical ampli®cation matrix. The numerical ampli®cation matrix ‰ZŠ relates to the ®nal states at the end of a time interval to the given initial conditions at the beginning of the time interval. In the present situation, the ®nal displacement and the ®nal velocity at

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1321

the end of a time interval can be related to the initial displacement and the initial velocity at the beginning of the time interval by 

u…t ˆ Dt† v…t ˆ Dt†



 ˆ

Z11 Z21

Z12 Z22



u0 v0



 ˆ ‰ZŠ

 u0 : v0

…51†

The eigenvalues k of ‰ZŠ may be real or complex. The largest magnitude of the eigenvalues is de®ned as the spectral radius [24±26]. For unconditionally stable algorithms, the spectral radius should always be less than unity for all time step size Dt. It is well known that the di€erential quadrature method is equivalent to the collocation method [17,27,28]. As the stability characteristics for various collocation methods are very well studied, many existing results could be used to investigate the stability properties of the various di€erential quadrature methods under consideration. Figs. 1±5 show the magnitude of the eigenvalues jkj versus xDt for algorithms with the initial conditions imposed by the conventional approach. Figs. 6±9 show jkj versus xDt for algorithms with the initial conditions imposed by the modi®ed di€erential quadrature rule. 5.1. Imposing the initial conditions by the conventional approach 5.1.1. Legendre sampling grid points Figs. 1(a)±(d) show the magnitudes of the eigenvalues of the numerical ampli®cation matrix jkj versus xDt when the Legendre sampling grid points are used with the initial conditions imposed by the conventional approach. The spectral radius corresponds to the larger magnitude of the two eigenvalues. It can be seen that initially the eigenvalues are complex conjugates with unit magnitude (i.e. jkj ˆ 1). As xDt increases, the eigenvalues may become real with one of the eigenvalues exceeding unity. The algorithms are therefore conditionally stable. For n P 2, the eigenvalues may collapse and become complex conjugates with unit magnitude again. In this case, a bubble is formed with bifurcation and reverse bifurcation in the plot. It can be seen that the number of bubbles increases with n and in general, the bifurcation occurs at xDt  p, 2p; . . . ; np. If the time step Dt is chosen such that xDt falls within the bubbles, the numerical results will be unstable as the errors would grow with the time steps. However, it can be seen that as n increases, the size of the ®rst bubble decreases. For example, when ˆ 3, the bubble at xDt  p is neglipn gible. The stability limit could therefore be assumed to be xDt < 240=7 … 5:8554†. Golley [29] also studied the use of the second-order Gauss quadrature points (i.e. the Legendre points) for collocation in solving dynamic problems. It was also found that the resultant algorithm was fourthorder accurate and conditionally stable only. This algorithm is in fact equivalent to the present algorithm using Legendre points with n ˆ 2 and the algorithm presented in [30] with W13 ˆ 1=6. 5.1.2. Radau sampling grid points Fig. 2(a)±(e) show the magnitudes of the eigenvalues of the numerical ampli®cation matrix jkj versus xDt when the Radau sampling grid points are used with the initial conditions imposed by the conventional approach. It can be veri®ed that the determinant of the numerical ampli®cation matrix decreases from 1 to 0 as xDt increases from zero to in®nity. Hence, the product of the two eigenvalues would be less than unity in magnitude. In other words, jkj 6 1 if the eigenvalues are complex conjugates. The algorithms may still be conditionally stable if the eigenvalues become real with one of them exceeding unity in magnitude. It can be seen that initially, the eigenvalues are complex conjugates with magnitudes less than unity. As xDt increases, the eigenvalues may become real and bifurcate. When n ˆ 1 (in this case s1 ˆ 1), the algorithm is unconditionally stable since the bifurcated eigenvalues do not exceed the unit magnitude. In fact, this method corresponds to the Backward Euler method, which is known to be unconditionally stable but is ®rst-order accurate only. When n increases, some bubbles may be formed along the curves with spectral radii (i.e. the maximum magnitude of the eigenvalues) exceeding unity. For example, the stability limits for n ˆ 2 and 3 are around

1322

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Fig. 1. Legendre sampling grid points with the conventional approach to impose the initial conditions: (a) n ˆ 1 and n ˆ 2; (b) n ˆ 3; (c) n ˆ 4; (d) n ˆ 5.

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1323

Fig. 2. Radau sampling grid points with the conventional approach to impose the initial conditions: (a) n ˆ 1; (b) n ˆ 2; (c) n ˆ 3; (d) n ˆ 4; (e) n ˆ 5.

2.9345 and 3.0910. Similar to the Legendre sampling grid points, the size of the initial bubble decreases as n increases as well. Hence, for n ˆ 4, even though the spectral radius exceeds unity over a very narrow interval around xDt ˆ p, the stability limit of the algorithm could be regarded as 6.0417. 5.1.3. Chebyshev sampling grid points Fig. 3(a)±(d) show the magnitudes of the eigenvalues of the numerical ampli®cation matrix jkj versus xDt when the Chebyshev sampling grid points are used with the initial conditions imposed by the conventional approach. It can be shown that the curves are qualitative similar to those given by the Legendre sampling grid points. In other words, the eigenvalues are complex conjugates with unit magnitude (i.e.

1324

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Fig. 2. (Continued).

jkj ˆ 1) initially. As xt increases, the eigenvalues may become real and bifurcate with one of the eigenvalues exceeding unity. The algorithms are therefore conditionally stable. 5.1.4. Chebyshev±Gauss±Lobatto and equally spaced sampling grid points Figs. 4(a)±(d) and Fig. 5(a)±(c) show the magnitudes of the eigenvalues of the numerical ampli®cation matrix jkj versus xDt when the Chebyshev±Gauss±Lobatto sampling grid points and the equally spaced sampling grid are used, respectively, with the initial conditions imposed by the conventional approach. It can be seen that the curves are qualitatively similar to those given by the Radau sampling grid points. When n ˆ 2, the two algorithms are p the same with s1 ˆ 1=2 and s2 ˆ 1. The algorithm is conditionally stable with a stability limit xDt 6 8. When n P 3, there are bubbles and bifurcation in the curves. However, it can be seen that the spectral radii may exceed unity even when xDt is very small. As a result, the algorithms are conditionally unstable. In other words, to maintain numerical stability, xDt must not be less than a critical value. This is different from the conditionally stable algorithms. Hence, these conditionally unstable algorithms are of limited value, as the time step could not be reduced repeatedly to improve the solution accuracy. In [6], the Chebyshev±Gauss±Lobatto sampling grid points with n ˆ 7 were used to construct time step integration algorithms. The resultant algorithms are in fact conditionally stable only as shown in Fig. 4(e). The stability limit should be xDt < 9:4354. It is interesting to note that when n ˆ 5 for the equally spaced sampling grid points, the denominator of the determinant of the numerical ampli®cation matrix would become zero when xDt  6:31404 and 9.98041. As a result, the spectral radius becomes very large near this region. Hence, the equally spaced sampling grid points are not suitable for the solution of dynamic problems.

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1325

Fig. 3. Chebyshev sampling grid points with the conventional approach to impose the initial conditions: (a) n ˆ 1 and n ˆ 2; (b) n ˆ 3; (c) n ˆ 4; (d) n ˆ 5.

1326

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Fig. 4. Chebyshev±Gauss±Lobatto sampling grid points with the conventional approach to impose the initial conditions: (a) n ˆ 2; (b) n ˆ 3; (c) n ˆ 4; (d) n ˆ 5; (e) n ˆ 7.

In summary, the algorithms obtained from the di€erential quadrature method by using the conventional approach to impose the initial conditions are not satisfactory. They are in general conditionally stable only, as veri®ed by Houwen et al. [22] and Kramarz [31]. However, it is possible to generate unconditionally stable algorithms if the residual at the initial time is also included in the formulation [32]. 5.2. Modi®ed di€erential quadrature rule to impose the initial conditions When the modi®ed di€erential quadrature rule is used to impose the initial conditions, it has been shown in [8] that the eigenvalues of the numerical ampli®cation matrix are complex conjugate pair. As a result, in

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1327

Fig. 4. (Continued).

the curves showing the magnitudes of the eigenvalues of the numerical ampli®cation matrix jkj versus xDt, there is no bifurcation or bubble. However, the spectral radius may still exceed unity if the determinant of the numerical ampli®cation matrix exceeds unity. 5.2.1. Legendre, Radau and Chebyshev sampling grid points Fig. 6 shows the spectral radii for the various algorithms employing the Legendre, Radau and Chebyshev sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions. It has been shown in [8] that the algorithms obtained by using the Legendre sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions are unconditionally stable. In fact, the algorithms are non-dissipative with the spectral radii equal to unity for all xDt. Besides, it is also shown in [8] that the algorithms obtained by using the Radau sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions are unconditionally stable. In fact, the algorithms are asymptotically annihilating with spectral radii approaching zero as xDt approaching in®nity. Although the accuracy of the solutions given by using the Chebyshev sampling grid points are not as high as those given by the Legendre or Radau sampling grid points, the algorithms are unconditionally stable. The method also gives uniformly low truncation errors for interpolated results found within a time step [15]. Furthermore, the algorithms are non-dissipative since the Chebyshev sampling grid points s1 ; . . . ; sn are symmetrically placed with respect to s ˆ 1=2 in the interval. In fact, it can be shown [33,34] that any two, three or four symmetrically placed s1 ; . . . ; sn could produce unconditionally stable algorithms. The conditions of the symmetrically placed s1 ; . . . ; sn for unconditionally stable algorithms with n 6 7 have been discussed in [35]. Besides, it has been veri®ed

1328

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Fig. 5. Equally spacing sampling grid points with the conventional approach to impose the initial conditions: (a) n ˆ 3; (b) n ˆ 4; (c) n ˆ 5.

[36] that if s1 ; . . . ; sn are chosen from the zeros of the shifted Gegenbauer (or ultraspherical) polynomials, Pn…k† …s† ˆ

Cn …s

s2 †k …1=2†

dn  …s dsn

s2 †n‡k

…1=2†

 ;

where Cn is a normalization constant;

…52†

so that s1 ; . . . ; sn will be symmetrically placed in the interval, the resultant algorithms will be non-dissipative and unconditionally stable for 1 6 n 6 100 if 0:5 < k 6 1:5. The Legendre and Chebyshev points correspond to the two special cases with k ˆ 0:5 and 0, respectively.

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1329

Fig. 6. Radau sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions.

Fig. 7. Chebyshev±Gauss±Lobatto sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions.

Fig. 8. Equally spaced sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions.

5.2.2. Chebyshev±Gauss±Lobatto and equally spaced sampling grid points Figs. 7 and 8 show the spectral radii for the various algorithms using the Chebyshev±Gauss±Lobatto and equally spaced sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions. It can be seen that the spectral radii are less than unity when xDt is greater than certain values. However, the algorithms are conditionally unstable with spectral radii exceeding unity when xDt is small. The

1330

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

Fig. 9. Chebyshev±Gauss±Lobatto sampling grid points with the modi®ed di€erential quadrature rule to impose the initial conditions …n ˆ 10†.

maximum spectral radii can be found in Figs. 7 and 8. As a result, the algorithms are not too useful. The algorithms are not unconditionally A-stable even if the Chebyshev±Gauss±Lobatto sampling grid points are used. It should be noted that no differential quadrature analog equation of the differential equations has been established at s0 ˆ 0. In [13], the Chebyshev±Gauss±Lobatto sampling grid points with n ˆ 10 were used. The maximum spectral radius is found to be 1.00156 when xDt ˆ 10:129 as shown in Fig. 9. Strictly speaking, the algorithms are not unconditionally stable. There may be numerical stability problems for very long-term integration. In summary, the most useful algorithms are those generated by using the Legendre, Radau and Chebyshev sampling grid points with the initial conditions imposed by the modi®ed di€erential quadrature rule. The algorithms using the Legendre and Radau sampling grid points have an additional advantage that they are of higher-order accuracy at the end of a time step than the other algorithms. 6. Conclusions In this paper, the in¯uences of the sampling grid points and the ways to impose the initial conditions on the solutions of dynamic problems using the di€erential quadrature method are studied. The following ®ndings can be observed: 1. The algorithms are in general conditionally stable if the initial conditions are imposed by the conventional approach. 2. Unconditionally stable algorithms can be obtained if the initial conditions are imposed by the modi®ed di€erential quadrature rule and using the Legendre, Radau and Chebyshev sampling grid points. 3. The Legendre and Radau sampling grid points could be used to generate higher-order accurate results at the end of a time step. 4. The performance of the commonly used Chebyshev±Gauss±Lobatto sampling grid points is not very satisfactory. The resultant algorithms do not have very high-order accuracy at the end of the time step. Besides, they are in general conditionally stable or conditionally unstable only. 5. The equally spaced sampling grid points are not suitable for the solution of dynamic problems using differential quadrature method. References [1] C.W. Bert, M. Malik, Di€erential quadrature method in computational mechanics: a review, Appl. Mech. Rev. 49 (1996) 1±28. [2] C.W. Bert, S.K. Jang, A.G. Striz, Two new approximate methods for analysing free vibration of structural components, AIAA J. 26 (1988) 612±618.

T.C. Fung / Comput. Methods Appl. Mech. Engrg. 191 (2002) 1311±1331

1331

[3] S.K. Jang, C.W. Bert, A.G. Striz, Application of di€erential quadrature to static analysis of structural components, Int. J. Numer. Meth. Engrg. 28 (1989) 561±577. [4] W.L. Chen, A.G. Striz, C.W. Bert, A new approach to the di€erential quadrature method for fourth-order equations, Int. J. Numer. Meth. Engrg. 40 (1997) 1941±1956. [5] X. Wang, H. Gu, Static analysis of frame structures by the di€erential quadrature element method, Int. J. Numer. Meth. Engrg. 40 (1997) 759±772. [6] T.Y. Wu, G.R. Liu, A di€erential quadrature as a numerical method to solve di€erential equations, Comput. Mech. 24 (1999) 197±205. [7] T.Y. Wu, G.R. Liu, The generalized di€erential quadrature rule for initial-value di€erential equations, J. Sound Vib. 233 (2000) 195±213. [8] T.C. Fung, Solving initial value problems by di€erential quadrature method ± part 2: second- and higher-order algorithms, Int. J. Numer. Meth. Engrg. 50 (2001) 1429±1454. [9] X. Wang, C.W. Bert, A new approach in applying di€erential quadrature to static and free vibrational analysis of beams and plates, J. Sound Vib. 162 (1993) 566±572. [10] X. Wang, C.W. Bert, A.G. Striz, Di€erential quadrature analysis of de¯ection, buckling, and free vibration of beams and rectangular plates, Comput. Struct. 48 (1993) 473±479. [11] M. Malik, C.W. Bert, Implementing multiple boundary conditions in the DQ solution of higher-order PDE's: application to free vibration of plates, Int. J. Numer. Meth. Engrg. 39 (1996) 1237±1258. [12] C. Shu, H. Du, Generalized approach for implementing general boundary conditions in the GDQ free vibration analysis of plates, Int. J. Solids Struct. 34 (1997) 837±846. [13] M. Tanaka, W. Chen, Dual reciprocity BEM applied to transient elastodynamic problems with di€erential quadrature method in time, Comput. Methods Appl. Mech. Engrg. 190 (2001) 2331±2347. [14] T.C. Fung, Solving initial value problems by di€erential quadrature method ± part 1: ®rst-order equations, Int. J. Numer. Meth. Engrg. 50 (2001) 1411±1427. [15] J.C. Butcher, The role of orthogonal polynomials in numerical ordinary di€erential equations, J. Comput. Appl. Math. 43 (1992) 231±242. [16] R. Bellman, B.G. Kashef, J. Casti, Di€erential quadrature: a technique for the rapid solution of nonlinear partial di€erential equations, J. Comput. Phys. 10 (1972) 40±52. [17] J.R. Quan, C.T. Chang, New insights in solving distributed system equations by the quadrature method. I. Analysis, Comput. Chem. Engrg. 13 (1989) 779±788. [18] J.R. Quan, C.T. Chang, New insights in solving distributed system equations by the quadrature method. II. Numerical experiments, Comput. Chem. Engrg. 13 (1989) 1017±1024. [19] C.W. Bert, X. Wang, A.G. Striz, Convergence of the DQ method in the analysis of anisotropic plates, J. Sound Vib. 170 (1994) 140±144. [20] G.R. Liu, T.Y. Wu, Numerical solution for di€erential equations of Dung-type nonlinearity using the generalized di€erential quadrature rule, J. Sound Vib. 237 (2000) 805±817. [21] R.W. Clough, J. Penzien, Dynamics of Structures, McGraw-Hill, New York, 1993. [22] P.J. van der Houwen, B.P. Sommeijer, N.H. Cong, Stability of collocation based Runge±Kutta±Nystr om methods, BIT 31 (1991) 469±481. [23] T.C. Fung, On the equivalence of the time domain di€erential quadrature method and the Runge±Kutta collocation method, Int. J. Numer. Meth. Engrg., to appear. [24] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cli€s, NJ, 1987. [25] W.L. Wood, Practical Time-Stepping Schemes, Clarendon Press, Oxford, 1990. [26] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, fourth ed., McGraw-Hill, New York, 1991. [27] C.W. Bert, X. Wang, A.G. Striz, Di€erential quadrature for static and free vibration analyses of anisotropic plates, Int. J. Solids Struct. 30 (1993) 1737±1744. [28] N. Bellomo, Nonlinear models and problems in applied sciences from di€erential quadrature to generalized collocation methods, Math. Comput. Modelling 26 (1997) 13±34. [29] B.W. Golley, A time-stepping procedure for structural dynamics using Gauss point collocation, Int. J. Numer. Meth. Engrg. 39 (1996) 3985±3998. [30] T.C. Fung, Unconditionally stable higher-order accurate Hermitian time ®nite elements, Int. J. Numer. Meth. Engrg. 39 (1996) 3475±3495. [31] L. Kramarz, Stability of collocation methods for the numerical solution of y 00 ˆ f …x; y†, BIT 20 (1980) 215±222. [32] T.C. Fung, Unconditionally stable collocation algorithms for second-order initial value problems, J. Sound Vib., to appear. [33] M. Crouzeix, P.A. Ruamps, On rational approximations to the exponential, RAIRO Analyse Numerique 11 (1977) 241±243. [34] E. Hairer, G. Wanner, Solving Ordinary Di€erential Equations II, Springer, Berlin, 1991. [35] R. Barrio, Characterization of low degree A-stable symmetric RK collocation methods, J. Comput. Appl. Math, 111 (1999) 1±11. [36] R. Barrio, On the A-stability of Runge±Kutta collocation methods based on orthogonal polynomials, SIAM J. Numer. Anal. 36 (1999) 1291±1303.