REDUCING THE COMPUTATIONAL REQUIREMENTS OF THE
DIFFERENTIAL QUADRATURE METHOD

Wen Chena , Xinwei Wangb and Yongxi Yua
c

Mechanical Engineering Department, Shanghai Jiao Tong University, Shanghai 200030,

P.R.China
b

Department of Aircraft Engineering, Nanjing University of Aeronautics and Astronautics,

Nanjing 210016, P. R. China.

The mail address for Wen Chen for proofs is P.O.Box 93B1, Shanghai Jiao Tong
University, Shanghai 200030, P.R.China.

Abstract
This paper shows that the weighting coefficient matrices of the differential quadrature
method (DQM ) are centrosymmetric or skew-centrosymmetric if the grid spacings are
symmetric irrespective of whether they are equal or unequal. A new skew centrosymmetric
matrix is also discussed. The application of the properties of centrosymmetric and skew
centrosymmetric matrix can reduce the computational effort of the DQM for calculations of
the inverse, determinant, eigenvectors and eigenvalues by 75%. This computational
advantage are also demonstrated via several numerical examples.

1

1. Introduction
The differential quadrature method (DQM ) is a rather efficient numerical method for the
rapid solution of linear and nonlinear partial differential equations involving one dimensions
or multiple dimensions[1-9]. Compared with the standard methods such as the finite
element and finite difference methods, the DQM requires less computer time and storage.
The essence of the DQM is that a partial derivative of a function is approximated by a
weighted linear sum of the function values at given discrete points. Its weighting coefficients
do not relate to any special problem and only depend on the grid spacing. Thus, any partial
differential equation can be easily reduced to a set of algebraic equations using these
coefficients. The DQM coefficient matrices play a prominent role in the application of the
DQM . Quan and Chang [10] pointed out that the skew centrosymmetric and
centrosymmetric structures of the DQM weighting coefficient matrices for the 1st and 2nd
order derivatives, and this properties can be used to decompose the DQM coefficient matrix
for ordinary differential equations of systemic models into two smaller size matrices. But in
that reference they did not use the centrosymmetric and skew centrosymmetric
terminologies and thus can not be fully aware of the possible advantages in application of
these properties. Chen and Yu [10] pointed out clearly the relation between the DQM
weighting coefficient matrices and centrosymmetric or skew centrosymmetric matrices.

2

The present study is a further development of the above work. Firstly, we proves that the
structure of the weighting coefficient matrices for any order derivative in the differential
quadrature method are centrosymmetric or skew centrosymmetric if the grid spacing is
symmetric with respect to the center point whether the grid spacings are equal or nonequal.
It is known that, in the evaluation of the determinant, inverse and eigenvalue and
eigenvectors,

a centrosymmetric matrix can be factorized into two or three products,

depending on whether the order N of the matrix is even or odd. By applying these
properties, the multiplication complexity can be reduced by 75% and the efficiency of the
DQM can be significantly increased, which is also demonstrated by solving the free
vibration of beams and plates in this paper.

Secondly, the skew centrosymmetric matrix is defined and its properties are found to be
similar to those of the centrosymmetric matrices. It is shown that this kind of matrix is
related to the DQM weighting coefficient matrices for all odd order derivatives.

Next, the computational reduction by using the centrosymmetric properties in the DQM is
extended to non-systemic problems.

Finally, based on the fact that DQM is a polynomial approach [8], a new truncation error
estimation and its distribution in the entire variable domain are given in appendix 1.

3

2. Formulas for Direct Computing the Weighting Coefficients of the DQM
The differential quadrature method approximates the derivative of a function with respect
to a variable at a given discrete point as a weighted linear sum of the function values at all
discrete points. The general differential quadrature approximation at the ith discrete point is
given by (Civan and Sliepcevich [4,5])
m
∂ f ( x)
∂ xm

N

= ∑ w ij f ( x j ),
j =1

x =x i

i = 1, 2, . . . , N. ,

(1)

where xj are the discrete points in the variable domain. f(xj) and wij are the function values at
these points and the weighting coefficients, respectively.

In order to determine the weighting coefficients wij, equation (1) must be exact for all
polynomials

of degree less than or equal to (N-1). To avoid the ill-conditioning

Vandermonde matrix (see Bellman[12]) when N is relatively large, the Lagrange
interpolation formula is used as a test function [6, 9, 10, 11], namely,
f

k

( x) =

x − xi
k − xi

N

∏x

i ≠k

.

(2)

Substituting Eq. (2) into Eq. (1) yields the following two formulae to compute directly the
weighting coefficient of the 1st order derivative [6, 10], i.e.,
Aij =

1
x j −xi

N

xi −x k
,
j −x k

∏x

k ≠i , j
i ≠j

i = 1,2,K , N and j = 1,2,K , N .

4

(3a)

and
N

Aii = − ∑
k ≠i

1
,
xi −xk

i = 1,2,K, N .

(3b)

For higher order derivatives, the weighting coefficients can be generated by two recursion
formulae [6]
wij

( m+1)

= m( Aij wij

(m )

−

( m −1)

wij
),
xi −xj

i ≠j

(4a)

and
wii

( m+1)

N

= − ∑ wij
j ≠i

( m+1)

,

(4b)

where the superscript (m+1) denotes the order of the derivative.

3. The centrosymmetricity of weighting coefficient matrices in DQM
First, the definitions about centrosymmetric and skew centrosymmetric matrices are given
Definition 1: A N×N matrix Q=[qij] is centrosymmetric if
qij= = -qN+1-i,N+1-j,

i, j=1, 2, ...., N.

(5)

Then, Q can be characterized by
JQJ=Q,

(6)

where J denotes the contra-identity matrix given by

5

K 0 1
L 1 0
 ,
N M M

L 0 0

0
0
J =
M

1

(7)

which has unit elements along the secondary diagonal and zeros elsewhere, noting that JT=J
and J2=I, the unit matrix, The effect of premultiplying any matrix by J reverses the order of
its rows, and postmultipling reverses the order of its columns.

Definition 2: A new, skew centrosymmetric matrix R=[ rij]N×N with the following property,
is introduced:
rij =−rN+1-i,N+1-j, i,j=1,2,....,N.

(8)

and
R=−JRJ.

(9)

To our knowledge, this is a new form of matrix, which has not been reported previously. In
the following, the structure of weighting coefficient matrices of the DQM is analyzed. By
applying Eq. (3a) and (3b), it can be deduced that
Aij

=

x

j

1
− xi

N

xi − xk
,
j − x k

∏x

k ≠i , j

= 1, 2, K , N

i, j

(10a)

i ≠j

and
Aii

=

N

∑(1 −

k ≠i

1
,
x i) − ( 1 − x k)

i = 1, 2, K , N. .

6

(10b)

The grid spacings often used in DQM are symmetric such as equally spaced grids and the
roots of the shifted Legendre or the shifted Chebyshev polynomials. For symmetric grid
spacing, over a domain 0≤ x ≤1, it is true that:
x

N + 1− i

= 1 − xi

(11)

where N is the number of grid points. Thus, substituting Eq. (11) into (10a) and (10b) yield
=

Aij

x

N + 1− j

1
− x N + 1 −i

N

∏

k ≠ N +1 − i , N + 1 − j

x N + 1− i − x k
,
x N + 1− j − x k

i, j

= 1, 2, K , N (12a)

i ≠j

and
=

Aii

N

∑

k ≠ N + 1− i

1
x

N + 1− i

− xk

,

i = 1, 2, K , N. .

(12b)

Thus,
Aij = − A N +1−i ,N +1 −j .

(13)

According to the definition given by Eq. (8), the weighting coefficient matrix for the 1st
order derivative can be concluded to be skew centrosymmetric if the grid spacing is
symmetric. Next, consider the coefficient matrix for the high order derivative. Using Eqs.
(13) and the recursion formulas (4a) and (4b) and after some manipulations, it can be shown
that
wij

( m+1)

= ( −1) ( m+1) wN + 1 − i , N + 1 − j ( m+1 ) .

(14)

where the wij(m+1) denotes the DQM coefficient for the (m+1)th order derivative. Therefore,
the DQM weighting coefficient matrices are skew centrosymmetric for odd order

7

derivatives and centrosymmetric for even order derivatives when the grid spacing is
symmetric.

4. Some Properties of the Centrosymmetric and Skew Centrosymmetric Matrices
4.1 Centrosymmetric Matrices
For completeness, the interesting structure of the centrosymmetric matrix of even degree
[13,14] is described briefly. The corresponding similar properties exist for odd order
centrosymmetric matrix. Here CN×N denotes the set of N×N centrosymmetric matrices.
Lemma 1. If Q 1, Q2∈CN×N and Q1-1 exists, then
(1) Q=Q1Q2∈ CN×N
(2) Q=Q1+Q2∈CN×N

(15)

(3) Q1-1∈ CN×N

Lemma 2. If Q∈ CN×N and N=2M, Q can be written as
 A JCJ 
Q = 

C JAJ  ,

(16)

where A and C are arbitrary M ×M matrices. The determinant of matrix Q can be evaluated
by
Q = A + JC A − JC

(17)

and the inverse of the matrix is given as

8

P
Q=
R

JRJ 
,
JPJ 

(18)

where 2P=(A+JC)-1+(A-JC)-1 and 2R=(A+JC)-1 - (A-JC)-1.

Applying Eqs. (17) and (18), the calculation effort of the inverse and determinant of a
centrosymmetric matrix can be reduced

by 75%. Note that premultiplying or

postmultiplying any matrix by J only moves the elements of the matrix and requires a very
little computational time.

Lemma 3. If Q∈CN×N and N=2M, then
 A JCJ 
Q=

C JAJ 

and

 A − JC
 0



A + JC 
0

(19)

are orthogonally similar. Thus, the evaluation of the eigenvectors and eigenvalues of Q is
equivalent to the that of two M ×M matrices A-JC and A+JC.

Lemma 4. The eigenvector of centrosymmetric matrix Q in Lemma 3 is the sum of a
symmetric vector and a skew-symmetric vector. The N/2 skew symmetric orthonormal
eigenvectors Ui and the corresponding eigenvalues λ i can be determined by
( A − JC )v i = λ v

i

(20a)

and

9

Ui =

1
2

[

vi

− (Jv i )

T

T

].
T

(20b)

The N/2 symmetric eigenvectors Wi and the corresponding eigenvalues Pi are found by
solving the following equation:

(A + JC ) y i

=P y

(21a)

i

and
Wi =

1
2

[

yi

T

(Jy i )

T

]
T

(21b)

4.2 Skew Centrosymmetric Matrices
Next, the structure of skew centrosymmetric matrix is discussed. The proofs for these
properties are similar to those for centrosymmetric matrices [13, 14], so they are omitted
here for brevity. We herein analysis only the even order skew centrosymmetric matrix, but
it should be pointed out that there are similar properties for odd order skew
centrosymmetric matrix. In the following, NCN×N is the set of N×N skew centrosymmetric
matrices.
Lemma 5. If Q 1, Q2∈NCN×N and Q1-1 exists, then
1) Q=Q1Q2∈ CN×N
(2) Q=Q1+Q2∈NCN×N

(22)

(3) Q1-1∈NCN×N

10

Lemma 6. If Q ∈NCN×N and N=2M, then Q can be stated as
 A −JCJ 
Q=

C −JAJ 

(23)

where A and C are arbitrary M ×M matrices. The determinant of matrix Q can be expressed
as
Q = A + JC A − JC

(24)

and the inversion of the matrix is given as
P
Q=
R

− JRJ 
− JPJ 

(25)

where 2P=(A+JC)-1+(A-JC)-1 and 2R=(A+JC)-1 - (A-JC)-1.

Lemma 7. If Q∈NCN×N and N=2M, then
 A −JCJ 
Q=

C −JAJ 

and

 0
 A + JC


A − JC 
0 

(26)

are orthogonally similar. Thus, the evaluation of the eigenvectors and eigenvalues of Q is
equivalent to the that of an M ×M matrices.

5. Application
In this section, several examples are solved by the DQM and the aforementioned properties
of the DQM coefficient matrix are applied to reduce the computational effort. There are
several approaches available in applying the DQM for high order boundary value problems.

11

In this paper, according to Wang and Bert [8], the weighting coefficient matrices are
modified in terms of the specific boundary conditions. It is obvious that such coefficient
matrices are also centrosymmetric when the boundary conditions are symmetric. In the
applications, equally spaced points and the roots of the shifted Chebyshev polynomials
given below are adopted as the grid points.
x i = cos

( 2i − 1) π
,
2N

i = 1, 2, K , N. .

(27)

For the later, the following formulae can be derived.
Aii

=

xi
,
2
( 1 - xi )

i = 1, 2, K , N. .

(28a)

and
Aij

=

( −1)

( i −j )

( xi − x

j )

2

(1 − xj )
2

( 1 − xi )

.

i ≠ j ,

(28b)

where Aij Ôs represent the weighting coefficients of the 1st order derivative. The weighting
coefficients for higher order derivatives can be easily computed from Eqs ( 4a) and (4b).

Example 1: Flexural Vibration of a Simply Supported Beam
The governing differential equation for this example can be expressed as:
W

iv

= ϖ 2 W( x ) ,

(29a)

where the nondimensionalized frequency is ϖ 2=ρA0L4ω2/EI. A0, L and ρ are the constant
cross-sectional area, the length of the beam, the density and I the constant area moment of

12

inertia about the neutral axis, respectively. The boundary conditions at the simply
supported ends are given by:
w( 0) = 0 ,
w( 1) = 0 ,

w ′( 0) = 0 ,

(29b)

w ′′( 1) = 0 .

(29c)

In terms of differential quadrature, Eq. (29a) is expressed as:
N −1

∑D

j =2

ij

2
W j = ϖ Wi ,

i = 2, K , ( N − 1).

(29d)

Note that the boundary conditions have been applied in the formulation of the weighting
coefficient matrix D ij , which is a centrosymmetric matrix as shown in the appendix 2. Eq.
(29d) can be written in the matrix form as:
 D 22
 w2 
D 23 L D 2, N − 1   w 2 





D 32
D 33 L D 3, N −1
w
w
 3  = ϖ 2  3  .

M
M
M
 M
 M 
 M 
 D



w 
D N −1, 3 L D N −1, N −1  w N − 1
N − 1, 2
N −1

(29e)

According to Lemma 2, we have

{D ij }

( N − 2 ) ×( N − 2 )

P
= 
R

JRJ 

JPJ  .

(29f)

In the present study, eight equally spaced grid points are used. Therefore, the order of the
matrices P and R is 3. In addition, it is known that the eigenvectors corresponding to the
fundamental frequency of a simply supported beam is symmetric. Thus, according to
lemma 4, the fundamental frequency and the corresponding eigenvector can be computed
by:

13

( P + JR) ui = ϖ 2 ui .

(29g)

Similarly, the 2nd order frequency with skew-symmetric eigenvector can be calculated by:

(P − JR ) v i = ϖ 2 v i

(29h)

The fundamental frequency obtained by the present method is 9.8683, comparing well with
the exact solution of 9.8696. The relative error is -0.006%. The 2nd order frequency is
39.2411, and the exact solution 39.4784, hence the relative error is 0.6%. If a nonuniform
grid formed from the Chebyshev polynomials are used, the errors are further reduced to
0.000% for the fundamental frequency and to -0.17%for the 2nd frequency. In contrast
with Wang and Bert [8], the present DQM reduces the computational effort by 87.5%
because of the use of the centrosymmetric property of weighting coefficients matrix.

Example 2: Transverse vibration of a thin rectangular plate simply supported on all four
edges.
Applying the DQM , the governing equation [7] for this case is given by:
N −1

∑D

k =2

ik

2
Wkj + ( 2 α )

N −1

N −1

∑ B ∑B

m=2

jm

k =2

ik

4
Wkm + ( α )

N −1

∑D

k =2

jk

2
Wik = ϖ Wij ,

(30)

where ϖ 2=ρha4ω2/D, and α=a/b=1.5, and D and B are the modified weighting coefficient
matrices for the fourth and second derivatives, respectively.

14

The DQM with a rather small number of grid points can produce very good results. In this
case, 8×8 equally spaced grids are used. The DQM results for the fundamental of
ϖ=32.0721 and for the 2nd frequency of ϖ=61.4449 agree well with the exact solutions of
32.0762 and 61.6850. The relative error is -0.01% and -0.39%. If the roots of the shifted
Chebyshev polynomials are used as the grid points, the DQM solutions for 1st and 2nd
frequency are calculated as 32.0761 and 61.6159, and the errors are further reduced to 0.0003% and 0.11%. For the same problem, Wang and Bert [8] need compute a (N-2)2×(N2)2 order matrix eigenvalue problem. In contrast, the present approach only requires solving
a (N-2)2 /2 ×(N-2)2 /2 matrix due to the use of the DQM centrosymmetric properties.
Therefore, the computational effort are only 12.5% of that of Wang and Bert [8].

Example 3. The transverse vibration of a thin, isotropic, skew plate vibration
This problem have been handled by the DQ method in reference [15]. But since its
governing and boundary condition equations contain cross derivative term, and thus its DQ
formulation exists the centrosymmetric and skew centrosymmetric matrices simultaneously,
we analyze this problem again to demonstrate the computational reduction by the
centrosymmetric properties for this sort of problem. The governing equation for this case is
given by
w , ξξξξ − ( 4 cos θ ) w, ξξξη + 2( 1 + 2 cos
+ w, ηηηη = ( ρh ω /
2

D) sin

4

2

θ ) w , ξξηη − ( 4 cos θ ) w , ξηηη
(31a)

θ

15

where θ=skew angle. ρ density, ω circular natural frequency and D flexural rigidity. The
simply supported boundary conditions are
w = 0,

w, ξξ − 2 w , ξn cos θ = 0,

at

ξ = 0, a

w = 0,

w , ηη − 2 w , ξn cos θ = 0,

at

η = 0, b

(31b)

The clamped boundary conditions are
w = 0,

w, ξ = 0,

at

ξ = 0, a

w = 0,

w , η = 0,

at

η = 0, b

(31c)

let x=(2ζ-a)/a, y=(2η-b)/b. The matrix formulations for derivatives of a two-variable
function are presented herein, i.e.

∂4 w
= D x w,
∂x 4
∂4w
3

∂x ∂y

∂4 w
2

∂x ∂y

∂4 w

T

2

= B x w By ,

∂4 w

T

= C x w Ay ,

∂x ∂y

3

= A x w Cy

∂y

4

= w Dy

T

,

(31d)

T

Note that W is (n-2)×(n-2) matrix. The DQM formulation for the governing equation (31a)
can be stated as
Dx W − (4 β cos θ ) C x w A y

(4 β

cos θ ) A x w C y

T

T

+ 2 β

4

+ β w Dy

T

=

(1
(ω

2

2

+ 2 cos

2

16)w

θ ) B x w By

T

−

(31f)

Applying the kronecker product of matrix to the above equation, we have

{[I n − 2] ⊗
(4 β

cos

[D ] − (4 β cos θ )[A ] ⊗ [C ] + 2 β (1 + 2 cos θ )[B ] ⊗ [B ] −
θ )[C ] ⊗ [A ] + β [D ] ⊗ [I
vec {W} = (ω 16)vec {W}
]}
2

x

y

4

y

x

2

x

y

y

x

2

n −2

(31g)

where ⊗ denotes the kronecker product in matrix computation (the DQ formulation for this
case in Eq. (10) of the reference [15] was incorrectly expressed). [I

16

n−2

]

represents a unite

matrix, whose n is the number of grid points. vec(w) is the vector-function of a rectangular
matrix

formed

by

4
2
ϖ 2=(ρh a ω D) sin

stacking
4

θ .

the

columns

of

matrix

into

one

long

vector.

β and D are aspect ratio and flexural rigidity, respectively.

Weighting coefficient matrices A ,

C , D and B are modified by the boundary conditions.

For clamped boundary conditions, the D and B are the same as the rectangular plate. But
the problem is more complex for simply supported boundary conditions because it contain
cross

derivative

w , xy = A x w A y

T

term

w , xy .

Using

the

DQ

method,

we

have

= ( A y ⊗ Ax ) vec ( w ) . It is easily proved that Ax ⊗ A y

is a

centrosymmetric matrix, since Ax and Ay are the skew centrosymmetric matrices in using
symmetric distributed grid points as in reference [15]. Therefore, D
centrosymmetric matrices and A ,

and B

will be

C the skew centrosymmetric matrix. Obviously, the

resulting matrix for the governing equation (31g) is a centrosymmetric matrix. Thus, we can
apply the computational reduction of centrosymmetric matrix for this case similar to
example 2. Namely, the computational effort is reduced by 87.5%.

Example 4. The convection-diffusion problem
Civan and Sliepcevich [16] computed the convection-diffusion problems by the differential
quadrature method. We analyze this example again to show that the computational
reduction in centrosymmetric matrix can be used for those problems with non-symmetric

17

boundary conditions. We consider the steady-state case. In terms of the DQ method, the
governing equation [16] can be expressed as
−

φ ij

N

x

+ α ∑ b ik
4a
k =1

x

φ kj + β

N

y

∑b

k =1

y
jk

φ jk = 0 ,

(32a)

where α and β are constant quantities. bikÔs are the DQ weighting coefficient for the second
derivative and have been not modified by boundary conditions, which is different from
example 1-3. The boundary conditions are approximated by the DQ method as
φij=prescribed value, j=1,2, ..., Ny

(32b)

φiNy = 0, j=1,2, ..., Nx

(32c)

∂φ ij
∂y

Ny

=

φNx j =0,

∑a

k =1

y
jk

φ ik = 0 , at j=1.

(32d)

j=1,2, ..., N x

(32c)

Since the boundary conditions in this problem are not systemic, we need handle with Eq.
(32a) and Eq. (32b-c), respectively. Obviously, the resulting coefficient matrix in Eq. (32a)
is a centrosymmetric matrix irrespective of whether the boundary conditions are symmetric
or non-symmetric. Solving the Eq. (32a) by using the centrosymmetric properties, the
function values at inner grid points are expressed by the side lateral function values using a
matrix formula. Substituting this formula into the boundary condition equations (32b-c), we
can obtain the algebraic equations only containing the side lateral function values. The order
of this equation is much smaller than that of Eq. (32a). The computational effort for solving
it will be little. The main computational effort in this procedure is solving Eq. (32a), while

18

the computational reduction by using the centrosymmetric matrix can be achieved in this
step. So the present computational effort is nearly only 25% as that in reference [16].

For the deflection and vibrational analysis of plates and beams with non-systemic boundary
conditions, the present reduction approach is also applicable, but the approach presented
by Wang and Bert [8] in applying the DQM to structural problems can be not used for the
cases. Another two different approaches proposed separately by Bert et. al. [7, 17] and by
Chen et. al [18] should be used, in which the governing and boundary condition equations
can be handled with separately.

6. Conclusion
This paper proves generally the centrosymmetricity of the weighting coefficient matrices in
the DQM when the grid spacing is symmetric. Also a new type of matrix, called as the
skew centrosymmetric matrix, is defined and its main properties are discussed. The
structure of such matrices allows for factorization of the determinant and the characteristic
equations into two, smaller and nearly equal size matrices. Therefore, the computational
complexity can be reduced by 75%. Furthermore, it is shown that these matrices possess
either symmetric or skew-symmetric eigenvectors. In many problems, the need of
symmetric or skew-symmetric eigenvectors can be applied as a constraint. Therefore, the
computational effort can be reduced further, as demonstrated by the foregoing examples.

19

The present work emphasizes the importance of the analysis of the weighting coefficient
matrices in the DQM and makes the DQM even more attractive for the practical engineering
analysis.

Reference
1. R. E. Bellman, B. G. Kashef and J. Casti, ÒDifferential quadrature: A technique for the
rapid solution of nonlinear partial differential equations,Ó J. Comput. Phys. 10, 4052(1972).
2. L. C. Hu and C. Hu, ÒIdentification of rate constants by differential quadrature in partly
measurable compartmental models,Ó Math. Biosci. 21, 71-76(1974).
3. F. Civan and C. M . Sliepcevich, ÒApplication of differential quadrature to transport
processes,Ò J. Math. Anal. Appl. 93, 206-221(1983).
4. F. Civan and C. M . Sliepcevich, ÒSolving integro-differential equations by the quadrature
method, Ò Proc. of the 1st Int. Conference on Integral M ethods in Sci. & Engrg., UT
Arlington, Mar. 18-21, 1985, Hemisphere Publ. Co., New York, 106-111, 1986.
5. F. Civan and C. M . Sliepcevich, ÒOn the solution of the Thomas-Fermi equation by
differential quadrature,Ó J. Comput. Phys., 56, 343-348(1984).
6. C. Shu and B. E. Richards, ÒParallel simulation of incompressible viscous flows by
generalized differential quadrature,Ó Computing Systems in Engineering, 3, 271-281(1992).

20

7. Yusheng. Feng and C. W. Bert, Ò Application of the quadrature to flexural vibration
analysis of a geometrically nonlinear beam,Ó Nonlinear Dynamics 3, 13-18(1992).
8. X. Wang and C. W. Bert, Ò A new approach in applying differential quadrature to static
and free vibrational analysis of beams and plates, Ò J. Sound & Vibration, 162, 566-572
(1993).
9. C. W. Bert, Xinwei Wang and A. G. Striz, ÒDifferential quadrature for static and free
vibrational analysis of anisotropic plates, Ò Int. J. Solids Structures 30, 1737-1744 (1993).
10. J. R. Quan and C. T. Chang, ÒNew insights in solving distributed system equations by
the quadrature methods - 1,Ó Computers Chem. Engrg 13, 779-788(1989).
11. Wen Chen and Yongxi Yu, ÒCalculation and analysis of weighting coefficient matrices in
differential quadrature,Ó Proc. of the 1st Pan-Pacific Conference on Computational
Engineering, Seoul, Korea, 1-5 Nov., 1993, Elsevier Science Publ. B. V. , Netherlands, 1993,
157-162.
12. R. E. Bellman, M ethods of nonlinear analysis, Vol. 2, chap. 16, Academic Press. New
York, 1973.
13. I. J. Good, ÒThe inverse of a centrosymmetric matrix,Ó Technometrics 12, 925-928
(1970).
14. L. Datta and S. D. M orgera, ÒOn the reducibility of centrosymmetric matrices
applications in engineering problems,Ó Circuits Systems Signal Process 8, 71-96 (1989).

21

15. X. Wang, A. G. Striz and C. W. Bert, ÒBucking and vibration analysis of skew plates by
the DQ method,Ó AIAA J. 37, 886-889 (1994).
16. F.Civan and C. M. Sliepcevich, ÒDifferential quadrature for multidimensional
problems,Ò J. Math. Anal. Appl., 101, 423-443(1984).
17. C. W. Bert, S. K. Jang and A. G. Striz, ÒTwo new methods for analyzing free vibration
of structure components,Ó AIAA Journal, 26, 612-618(1988).
18. Wen Chen and Yongxi Yu, ÒDifferential quadrature method for high order boundary
value problems,Ó Proc. of the 1st Pan-Pacific Conference on Computational Engineering,
Seoul, Korea, 1-5 Nov., 1993, Elsevier Science Publ. B. V. , Netherlands, 1993, 163-168.

Appendix 1
The DQM Truncation Error
If a function can be approximated by a Lagrangian polynomial [1, 9], then
f ( x ) = ψ j ( x ) f ( x j ) + R( x ) ,

j = 1, 2, K , N.

(A1)

where summation convection is applied on the repeated index j, ψj(x) are the Lagrangian
polynomials, and R(x) is the truncation error. Hence,
f (′ x i ) = ψ ′ j ( x i ) f ( x j ) + R (′ x i )
= Aij f ( x j ) + R (′ x i ) ,

I, j=1,2,...,N.

(A2)

where Aij are the DQM weighting coefficient of the 1st order derivative. xi and xj are the
discrete interpolation points, RÕ(xi) is the truncation error of the DQM to approximate 1st

22

order derivative. The truncation error for higher order derivatives can be obtained in a similar
way. Note that
( N)

f

R( x ) =

where W( x ) =
f

R( x i ) =

( ζ ) W( x )
,
N!
N

∏( x

i =1

( N)

− x i ).

(A3)

Thus,

( ζ ) W′( x i )
,
N!

i=1,2,...,N.

(A4)

It should be pointed out that the truncation error formula (A4) is different from that of
Bellman [1]. Considering equally spaced grid points and assuming f
R ′( x i )

≤

K
N ( N − 1)

N −1

and

R ′( x N)

≤

( N)

K
N ( N − 1)

N −1

( ζ)

≤ K. Then,

(A5)

for both ends and
R ′( x N / 2)

≤

( N / 2 − 1) !( N /
N ! ( N − 1)

2) ! K

N −1

(A6)

for the center points if N is an even number.

Comparing the truncation error formulae (A5) for the ends with that for the center point
(A6), it may be concluded that largest truncation error could possibly occur in the vicinity
of the ends if equally spaced grid points are used in the DQM . The two examples as
accuracy tests (Eqs. (4) and (5) of ref. 4) can be as the numerical computation
demonstration for the present error analysis (see fig. 1 and fig. 2 of ref. 4). Also a similar

23

error distribution is expected for unequally spaced grids such as defined by the roots of the
shifted Chebyshev polynomials.

Appendix 2

Application of the boundary conditions

in the formulation of the weighting

coefficients
Reference [8[ provided a new approach in applying the boundary conditions in the DQM ,
which applied the boundary conditions in formulation of the weighting coefficients. For a
simply supported beam, reference [8] gave a detailed description. When the symmetric grid
spacing is used, obviously, the modified 1st order derivative weighting coefficient matrix

[A ] (Eq. (3) of ref. 8) is a skew centrosymmetric matrix like [A] (Eq. (2) of ref. 8). The
modified weighting coefficient matrix for the 2nd order derivative is given by (Eq. (4) of ref.
8)

[B ] = [A][A ] .

(B1)

According to Lemma 5, [B ] is a centrosymmetric matrix. The modified weighting coefficient
matrix for the 4th order derivative [D] is decided by (Eq. (6) of ref. 8)

[D] = [B ][B ] .

(B2)

According to Lemma 1, [D] is also a centrosymmetric matrix.

24

