# Characteristics of the Differential Quadrature Method and Its Improvement

Authors: Wang Fangzong; Liao Xiaobing; Xie Xiong
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The differential quadrature method has been widely used in scientific and engineering computation. However, for the basic characteristics of time domain differential quadrature method, such as numerical stability and calculation accuracy or order, it is still lack of systematic analysis conclusions. In this paper, according to the principle of differential quadrature method, it has been derived and proved that the weighting coefficients matrix of differential quadrature method meets the importantV-transformation feature. Through the equivalence of the differential quadrature method and the implicit Runge-Kutta method, it has been proved that the differential quadrature method is A-stable and s-stage s-order method. On this basis, in order to further improve the accuracy of the time domain differential quadrature method, a class of improved differential quadrature method of s-stage 2s-order have been proposed by using undetermined coefficients method and Padé approximations. The numerical results show that the improved differential quadrature method is more precise than the traditional differential quadrature method.


---

## Body

## 1. Introduction

The differential quadrature method (DQM) was first proposed by Bellman and his associates in the early 1970s [1, 2], which is usually used for solving ordinary and partial differential equations. As an analogous extension of the quadrature for integrals, it can be essentially expressed as the values of the derivatives at each grid point as weighted linear sums approximately of the function values at all grid points within the domain under consideration.The differential quadrature method is conceptually simple and the implementation is straightforward. It has been recognized that the differential quadrature method has the capability of producing highly accurate solutions with minimal computational effort [3, 4] when the method is applied to problems with globally smooth solutions. So far, the differential quadrature method has been widely applied to boundary-value problems in many areas of engineering and science, such as structural mechanics [5–8], transport process [9], dynamic systems [10–12], and calculation of transmission line transient response [13, 14]. A comprehensive review of the chronological development of the differential quadrature method can be found in [4]. Although the differential quadrature method has been successfully applied in so many fields, for the basic characteristics of the method, such as numerical stability and calculation accuracy or order, not much work about them has been done in this area for the differential quadrature method. According to Fung [15], using Lagrange interpolation functions as test functions, the differential quadrature in time domain was shown to be equivalent to the recast implicit Runge-Kutta method [16–18]; besides, some low-order algorithms were discussed in detail. However, the method used by Fung is not the traditional sense of differential quadrature method but involved postprocessing (i.e., numerical solution at the end of grid points adopts polynomial extrapolation).In this paper, using general polynomial as test functions [19], the weighting coefficients matrix of differential quadrature method is proved to satisfy V-transformation [17, 20]. The equivalent implicit Runge-Kutta method is constructed through the differential quadrature method. Hence, making use of Butcher fundamental order theorem and the method of linear stability analysis [17, 18], the basic characteristics of the differential quadrature method can be systematically analysed. Unfortunately, the differential quadrature method is only a method of s-stage s-order and A-stable. Consequently, the differential quadrature method cannot yield higher accurate solutions to the boundary-value problems with fewer computational efforts. Based on above deduction, the method of undetermined coefficients is used to make the stability function of the equivalent Runge-Kutta method become the diagonal Padé approximations to the exponential function [17, 18]. Therefore, a class of improved differential quadrature method of s-stage 2s-order is derived.The paper is arranged as follows. In Section2, the weighting coefficients matrix of traditional differential quadrature method using general polynomial as test functions is briefly discussed. In Section 3, the equivalent relationship between the differential quadrature method and the Runge-Kutta methods is deduced. In Section 4, the stability and accuracy characteristics of the differential quadrature method are studied. A class of improved differential quadrature method of s-stage 2s-order and A-stable is proposed in Section 5. In Section 6, the transient response of a double-degree-of-freedom system is computed, which is given to verify the computational accuracy with the defined three grid points. Conclusions are then given in Section 7.


## 2. Traditional Differential Quadrature Method

Suppose functionf(x) is sufficiently smooth in the whole interval; there are (s+1) grid points with coordinates as ci,i∈(0,s). The first order derivative f1(ci) at each grid point ci,i∈(1,s), is approximated by a linear sum of all the function values in the whole domain; that is,(1)f1ci=∑j=0sgijfcj,i∈1,s,where f(ci) represent function values at a grid point ci, and gij is the weighting coefficients.In order to compute the weighting coefficientsgij in (1), the test functions can be chosen as(2)rkx=xk,k=0,1,…,s.Substituting (2) into (1) gives(3)k=0,0=∑j=0sgij,i∈1,s,(4)∑j=0sgijcjk=k·cik-1,i,k∈1,s.Equation (3) can be expanded into matrix form as(5)g11g12⋯g1sg21g22⋯g2s⋮⋮⋱⋮gs1gs2⋯gss11⋮1+g10g20⋮gs0=00⋮0.Let(6)G0=g10g20⋮gs0,G=g11g12⋯g1sg21g22⋯g2s⋮⋮⋱⋮gs1gs2⋯gss,e=11⋮1.Using (6), (5) can be simplified to(7)G0≡-Ge.From k=1,2,…,s, (4) can be expanded as(8)gi0c0+gi1c1+gi2c2+⋯+giscs=1gi0c02+gi1c12+gi2c22+⋯+giscs2=2ci⋮gi0c0s+gi1c1s+gi2c2s+⋯+giscss=scis-1.Since initial grid point c0 is usually defined as 0, (8) reduces to(9)gi1c1+gi2c2+⋯+giscs=1gi1c12+gi2c22+⋯+giscs2=2ci⋮gi1c1s+gi2c2s+⋯+giscss=scis-1.From i=1,2,…,s, (9) can be also expanded into matrix form as(10)Gc1c12⋯c1sc2c22⋯c2s⋮⋮⋱⋮cscs2css=12c1⋯sc1s-112c2⋯sc2s-1⋮⋮⋱⋮12cs⋯scss-1=1c1⋯c1s-11c2⋯c2s-1⋮⋮⋱⋮1cs⋯css-1.12⋱s.Vandermonde matrix V is defined as follows:(11)V=1c1⋯c1s-11c2⋯c2s-1⋮⋮⋱⋮1cs⋯css-1.Making use of (11), (10) can be expressed as(12)G-1V=c1c12⋯c1sc2c22⋯c2s⋮⋮⋱⋮cscs2⋯css112⋱1s=1c1c12⋯c1s1c2c22⋯c2s⋮⋮⋮⋱⋮1cscs2⋯css00⋯0112⋱1s.Finally, it can be inferred that(13)V-1G-1V=As,where As is(14)As=000⋯α1100⋯α2012⋮⋮⋮⋮⋮⋮⋱⋮00⋯1s-1αs,with(15)αs=α1,α2,…,αsT=1sV-1cs,cs=c1s,c2s,…,cssT.Equation (13), that is, G=VAs-1V-1, is called the implicit expression of the weighting coefficients matrix of the differential quadrature method and is also called V-transformation.When the grid points have been selected, the weighting coefficients matricesG and G0 are easy to calculate with the above formula. Obviously, the weighting coefficients of the differential quadrature method depend on the test functions and distribution of grid points but are independent of some specific problems. There are four typical grid points’ distributions: Legendre grid points, Chebyshev grid points, Chebyshev-Gauss-Lobatto grid points, and Uniform grid points (also called equally spaced grid points) [18]. This paper will focus on the latter three kinds of commonly used grid points, which are defined as follows:(1)
Chebyshev grid points:(16)ck=121-cos2k-12s-2π,k∈1,s-1,c0=0,cs=1;(2)
Chebyshev-Gauss-Lobatto grid points:(17)ck=121-cosksπ,k∈0,s;(3)
Uniform grid points:(18)ck=ks,k∈0,s.


## 3. The Equivalence of Differential Quadrature Method and Runge-Kutta Method

In order to analyse the numerical stability and order of the differential quadrature method, the differential quadrature method in time domain can be transformed into equivalent implicit Runge-Kutta  method. Consider the following ordinary differential equation(19)dxdt=ft,x,0<t≤Txt=0=x0.In the following, tn, tn+1 represent, respectively, the beginning and the end points at each step. h=tn+1-tn will be used to denote the step size. The time interval [tn,tn+1] will be normalized. That is, c=(t-tn)/h,t∈tn,tn+1,c∈[0,1]. At the same time, (19) can be made in the standard normalized form(20)ddcx=hftn+ch,x~,x~=xtn+ch;then, using s-stage differential quadrature method to solve (20) yields(21)Gx~1⋮x~s+G0xn=hftn+c1h,x~1⋮ftn+csh,x~s,where x~i=x(tn+cih),i∈(1,s). Since G0≡-Ge, (21) reduces to(22)x~1⋮x~s=xn1⋮1+hG-1ftn+c1h,x~1⋮ftn+csh,x~s.Let(23)G-1=A=aij,i,j∈1,s.Clearly, making use of (13) and (23) leads to(24)A=VAsV-1.Therefore, the weighting coefficients matrix A also satisfies V-transformation. It can be inferred from (22) that(25)x~i=xn+h∑j=1saijftn+cjh,x~j,i∈1,s.Since cs=1, tn+csh=tn+h=tn+1; therefore, x~i(i=s) is the approximate solution at the end of the step. Then, x~s can be rewritten as the following form:(26)x~s=xn+1=xn+h∑j=1sasjftn+cjh,x~j=xn+h∑j=1sbjftn+cjh,x~j,where bj=asj,j∈(1,s). It can be seen that (25) and (26) are the standard forms for an s-stage Runge-Kutta method. Since x~s=xn+1, the equivalent Runge-Kutta method is a reducible method [20]. In fact, the traditional differential quadrature method generally does not involve postprocessing, so the Runge-Kutta method converted from traditional differential quadrature method will naturally become a reducible method. The Runge-Kutta method can be conveniently summarized in the Butcher tableau [18] as(27)cAbTz21=c1a11⋯a1s⋮⋮⋮csas1⋯assb1⋯bs,where bT=(b1,b2,…,bs)=(as1,as2,…,ass).


## 4. Analysis of the Basic Characteristics of the Differential Quadrature Method

The stability and accuracy characteristics of the equivalent Runge-Kutta method will be investigated next. From (13) and (23), it can be inferred that(28)a11a12⋯a1sa21a22⋯a2s⋮⋮⋱⋮as1as2⋯ass1c1⋯c1s-11c2⋯c2s-1⋮⋮⋱⋮1cs⋯css-1=c1c12⋯c1sc2c22⋯c2s⋮⋮⋱⋮cscs2css112⋱1s.Equation (28) reduces to(29)∑j=1saijcjk-1=cikk,i∈1,s,k∈1,s.On the other hand, since bj=asj,j∈(1,s) and cs=1, from (28), it can be inferred that(30)bTV=1,12,…,1s.Similarly, (30) reduces to(31)∑i=1sbicik-1=1k,k∈1,s.Obviously, it has been shown that the equivalent Runge-Kutta method at least satisfies simplifying assumptions C(s) and B(s) from (29) and (31). Furthermore, it can be verified that the equivalent Runge-Kutta method only satisfies simplifying assumptions D(0). From Theorem 5.1 on page 71 in [17], it can be concluded that the implicit Runge-Kutta method or the corresponding differential quadrature method is s-stage s-order.The stability functionR(z) of the equivalent Runge-Kutta method or the corresponding differential quadrature method is given by the formula [16–18](32)Rz=det⁡I+zebT-Adet⁡I-zA,where, as usual, I is the identity matrix of dimension s. Due to grid points’ asymmetric distribution, the equivalent implicit Runge-Kutta method is not a symmetric method. As a result, there is a unique adjoint method (also called reflected method) [16, 18], which is defined as c∗A∗b∗T, satisfying(33)c∗=e-Pc,PA∗P=ebT-A,b∗=Pb,where P is(34)P=11⋰1∈Rs×s.Furthermore, from Theorem 343B on page 221 in [18], if the original method satisfies the simplifying assumptions C(s) and B(s), then the adjoint method also satisfies the same simplifying assumptions. Hence, the adjoint method enjoys V-transformation:(35)A∗=V∗As∗V∗-1,V∗=1c1∗⋯c1∗s-11c2∗⋯c2∗s-1⋮⋮⋱⋮1cs∗⋯cs∗s-1,where As∗ is(36)As∗=000⋯β1100⋯β2012⋮⋮⋮⋮⋮⋮⋱⋮00⋯1s-1βs,with(37)βs=β1,β2,…,βsT=1sV∗-1c∗s,c∗s=c1∗s,c2∗s,…,cs∗sT.Equation (32) can be reduced to(38)Rz=detI+zPA∗P-1detI-zVAsV-1=detI+zA∗detI-zAs=det⁡I+zV∗AsV∗-1det⁡I-zAs=det⁡I+zAs∗det⁡I-zAs.Because As and As∗ are a class of special matrices, (38) can be evaluated as(39)Rz=detI+zAs∗detI-zAs=1-∑k=s1βkk-1!/s-1!-zs-k+11-∑k=s1αkk-1!/s-1!zs-k+1,since cs=1, resulting in c1∗=0, it implies β1=0. It can be verified that the stability function of equivalent Runge-Kutta method is A-acceptability of p-order (p≥s-1) rational approximation to exponential function [21]. Therefore, the corresponding differential quadrature method is A-stable.In the following, the three-stage differential quadrature method using Uniform grid points will be given as an example. Whens=3, c1, c2, and c3 are given by 1/3, 2/3, 1. It can be worked out that matrices G0 and G are given by(40)G0=-112-1,G=-323-12-332192-9112,and the Butcher tableau of equivalent Runge-Kutta method is(41)cAbTz22=132336-495362379-29191340141z223401z224.V-transformation of matrix Α is(42)A=11319123491110022710-1127012231131912349111-1,the stability function of the equivalent Runge-Kutta method is(43)Rz=1+1/3z+1/27z21-2/3z+11/54z2-1/27z3,and the Butcher tableau of the adjoint method is(44)c∗A∗b∗Tz222=000013-136295362319491414z221403z224.V-transformation of matrix Α∗ is(45)A∗=100113191234900010-227012131001131912349-1.


## 5. Improved Differential Quadrature Method

Based on the above deduction, the traditional differential quadrature method is a method ofs-stage s-order. Compared with the multistage high-order Runge-Kutta method, for example, Gauss method (s-stage 2s-order method), it has the disadvantage of lower precision. As it is well-known that if the stability function of a numerical method is diagonal Padé approximations to the exponential function, then this method is the method of A-stable and s-stage 2s-order [22]. Inspired by this idea, the stability function of new Runge-Kutta method or new differential quadrature method have been converted into diagonal Padé approximation to the exponential function by using undetermined coefficients method. From (38), it can be seen that the stability function of the equivalent Runge-Kutta method will be determined by As and As∗. Suppose As=As∗=A~s, without changing b and V, a new Runge-Kutta method cA~bT is redefined as(46)A~=G~-1=V000⋯γ1100⋯γ2012⋮⋮⋮⋮⋮⋮⋱⋮00⋯1s-1γsV-1=VA~sV-1.Then, the stability function of this new Runge-Kutta method becomes(47)Rz=detI+zA~sdetI-zA~s=1-∑k=s1γkk-1!/s-1!-zs-k+11-∑k=s1γkk-1!/s-1!zs-k+1.From (46) and (47), it can be inferred that the last column elements in A~s determine the stability function. To improve the order of new Runge-Kutta method, undetermined coefficients γ can be selected so that the stability function of new Runge-Kutta method becomes the diagonal Padé approximations to the exponential function (defined by ezss):(48)Rz=1-∑k=s1γkk-1!/s-1!-zs-k+11-∑k=s1γkk-1!/s-1!zs-k+1=ezss.By comparing the coefficients on both sides of (48), undetermined coefficients γ can be conveniently obtained as(49)s=2,γ2=-112,12,s=3,γ3=160,-15,12,s=4,γ4=-1280,114,-928,12,⋮After getting γs=γ1,…,γsT, coefficients matrix A~ or G~ can also be easily computed through using (46). Therefore, a class of new Runge-Kutta method of s-stage 2s-order has been successfully constructed. In other words, a class of improved differential quadrature method of s-stage 2s-order has been derived. Besides, the adjoint method of new Runge-Kutta method is also s-stage 2s-order.Take the same as above, the improved differential quadrature method using Uniform grid points will be given as an example. It can be worked out that the new matricesG0 and G are given by(50)G~0=4343-12,G~=-17210-176-112416752-42332,and the Butcher tableau of new Runge-Kutta method is(51)cA~bTz22=1373120-2360131202397120-17601712012740320740740z223401z224.V-transformation of matrix Α~ is(52)A~=11319123491110016010-15012121131912349111-1and the Butcher tableau of the adjoint method is(53)cA~∗bTz22=0340-32034013131201760-7120231712023601712014z221403z224.V-transformation of matrix A~∗ is(54)A~∗=10011319123490016010-15012121001131912349-1.


## 6. Numerical Examples

Consider a two-degree-of-freedom system governed by(55)2001x¨1x¨2+6-2-24x1x2=010,with initial condition(56)x1=0,x1˙=0x2=0,x2˙=0.The exact solution of this problem is(57)x1=1-53cos2t+23cos5t,x2=3-53cos2t-43cos5t.The differential quadrature method can be used to find the numerical solutions to transient response directly. The detailed calculation steps of solving second-order differential equations in [11] can be seen. Figures 1, 2, and 3 show the displacement error trajectories comparison of traditional differential quadrature method and improved differential quadrature method with the same step size h=0.5s. In these Figures, the exact solution at each step is used for comparison. From Figures 1–3, it is evident that improved differential quadrature method is two orders of magnitude higher than traditional differential quadrature method. The error of improved differential quadrature method ranges between 10−5 and 10−4, even with a large step size h=0.5s.Figure 1
Error trajectories comparison of different DQM using Chebyshev grid (s=3, h=0.5s).
(a)
Error trajectories ofx1(t)
(b)
Error trajectories ofx2(t)Figure 2
Error trajectories comparison of different DQM using Chebyshev-Gauss-Lobatto grid (s=3, h=0.5s).
(a)
Error trajectories ofx1(t)
(b)
Error trajectories ofx2(t)Figure 3
Error trajectories comparison of different DQM using Uniform grid (s=3, h=0.5s).
(a)
Error trajectories ofx1(t)
(b)
Error trajectories ofx2(t)In order to compare the calculation precision of the traditional differential quadrature method and the improved differential quadrature method better, dynamic equations (55) have been computed for different time steps. Table 1 shows the computational results at t=60 s. Comparing with the exact solution, the numerical results using traditional differential quadrature method are the same until the first decimal place when s=5, h=1s, while the numerical results using improved differential quadrature method are the same until the third decimal place. When s=10, h=2s, the numerical results using traditional differential quadrature method are the same until the fifth decimal place, while the numerical results using improved differential quadrature method are the same until the seventh decimal place (especially using Chebyshev grid points and Chebyshev-Gauss-Lobatto grid points, the numerical results are almost the exact solution).Table 1
Computational results of the displacement (t = 60 s).
t = 60 s
x
1
(
t
)
x
2
(
t
)
Exact solution
2.264386552
5.469005157
s = 5 h = 1s
Traditional DQM-Chebyshev
2.253975383
5.491424542
Improved DQM-Chebyshev
2.264406270
5.468965760
Traditional DQM-Chebyshev-Gauss-Lobatto
2.276613635
5.439926064
Improved DQM-Chebyshev-Gauss-Lobatto
2.264406270
5.468965760
Traditional DQM-Uniform
2.342513616
5.278336836
Improved DQM-Uniform
2.264406270
5.468965760
s = 10 h = 2s
Traditional DQM-Chebyshev
2.264383332
5.469011768
Improved DQM-Chebyshev
2.264386552
5.469005157
Traditional DQM-Chebyshev-Gauss-Lobatto
2.264399479
5.468978613
Improved DQM-Chebyshev-Gauss-Lobatto
2.264386552
5.469005157
Traditional DQM-Uniform
2.265087814
5.467565205
Improved DQM-Uniform
2.264386561
5.469005173


## 7. Conclusion

In this paper, the linear stability and the order of differential quadrature method in time domain are systematically studied in detail and a class of new differential quadrature method ofs-stage 2s-order is proposed. From the above analysis and derivation, the following conclusions can be made.(1)
Based on general polynomial as test functions, the weighting coefficients matrix of the differential quadrature method satisfies theV-transformation. It plays an extremely important role in the analysis of basic characteristics of differential quadrature method and its improvement.(2)
The traditional differential quadrature method can be converted into equivalent Runge-Kutta method of A-stable ands-stage s-order. Therefore, compared with the commonly used single-stage low-order numerical integral methods, even with a small number of grid points, the traditional differential quadrature method can also give more accurate solutions. This is the main mechanism that differential quadrature method has been successfully applied in many fields.(3)
Finally, by making the stability function of equivalent Runge-Kutta method become the diagonal Padé approximations to the exponential function, a class of improved differential quadrature method ofs-stage 2s-order and A-stable is proposed. Therefore, the improved differential quadrature method can be extended to multi-degree-of-freedom time domain dynamic systems, which can produce higher accurate solutions at lower computational cost.



---

*Source: 657494-2015-06-01.xml*