NORTH.HOILAND
Analysis of the Response of a Strongly Nonlinear D a m p e d System using a Differential Transformation Technique Ming-Jyi Jang
Department of Mechanical Engineering Far East College Tainan, Taiwan, Republic of China and Chieh-Li Chen
Institute of Aeronautics and Astronautics National Cheng-Kung University Tainan, Taiwan, Republic of China
ABSTRACT
The differential transformation technique is applied to solve a nonlinear differential equation that describes the underdamped and the overdamped motion of a system subjected to external excitations. Numerical results are compared to those obtained by the Runge-Kutta method to illustrate the preciseness and effectiveness of the proposed method. © Elsevier Science Inc., 1997
1.
INTRODUCTION
The classical perturbation methods such as the Lindstedt-Poincare method, the multiple time scale method, the Krylov-Bogoliubov-Mitropolski method and averaging method [1] are widely used to determining the approximate solutions of weakly nonlinear systems described by
~ + ~2 x = ~ f( x, 3c),
(1)
where s is a small parameter. Many research works [2-5] have been devoted to the extension of the classical method to solve strongly nonlinear systems
APPLIED MATHEMATICSAND COMPUTATION88:137-151 (1997) © Elsevier Science Inc., 1997 655 Avenue of the Americas, New York, NY 10010
0096-3003/97/$17.00 PII S0096-3003(96)00308-6
138
M.-J. JANG AND C.-L. CHEN
of the form ~ + ax + bx 3 = 6 f ( x, ~).
(2)
Among these approaches, the extended Krylov-Bogoliubov-Mitropolski method has been applied to solve a nonlinear damped system subjected to different function excitations. Recently, an approximate technique to solve this type of problem was proposed by E1-Asir et al. [6] based on a method introduced by Mansour and Hussein [7]. In this paper, a differential transformation technique [8], which can be considered as an extended Taylor series method of order k, is applied to a strongly nonlinear problem. The method can be used to obtain the limit cycles of the free nonlinear oscillation systems as well as the response of the forced nonlinear oscillation systems. Numerical results are given to illustrate the effectiveness of the proposed method. 2.
PROBLEM FORMULATION
The equation of motion of a nonlinear, damped system of a single degree of freedom subjected to external excitation can be described as +
+ a2x +
= E( t).
(3)
dx
Where ~ = - ~ , 2y is the linear damping parameter, 6 is the nonlinearity parameter, 1] is the frequency of underdamped motion, f( x, 5) is a nonlinear function of x and its derivative and E(t) is the external excitation. In this paper, (3) with f( x, 5) = 71 x2 x + x3 is considered, i.e., + (2y+ sylx2)}+122x+
3.
ex 3= E(t).
(4)
DIFFERENTIAL TRANSFORMATION TECHNIQUE
Differential transformation is one of the useful techniques to solve the ordinary differential equations with fast convergence rate and small calculation error. One advantage of this method over the integral transformation methods is its ability to handle nonlinear differential equations, either initial value problems or boundary value problems. The basic definitions and operations of differential transformation [8, 10] are introduced as follows.
Differential Transformation Technique
139
DEFINITION 1. If x(t) is analytic in the time domain T then it will be differentiated continuously with respect to time t.
~(t) 8t k
~( t, k) Vt e T.
(5)
For t = t~ then ~(t, k) = ~(t~, k), where k belongs to the set of nonnegarive integers, denoted as the K domain. Therefore, (5) can be rewritten as
X ( k ) = ~( ti, k) = [
Ot k ]t=&Vk ~- K
(6)
where X(k) is called the spectrum of x(t) at t = t~ in the K domain.
DEFINITION 2. represented as
If x(t) can be expressed by Taylor series then x(t) can be
( t - t,) k k! x(k)
x(t)= E
k=O
(7)
Eq. (7) is called the inverse transformation of X(k). Using the symbol " D " denoting the differential transformation process and combining (6) and (7), it is obtained that
x(t)=
E
k=O
( t - t,) k k! x(k)
--- D-l{ X ( k ) } .
(8)
If the X(k) is defined as
• (k) = M(k)
akq( t) x( t) ] where k = 0, 1 , 2 , . . . , ~ c~tk
(9)
] t=to
then the function x(t) can be described as
1 ~ (t-k/o)k x(k) x( t) = q( t) E M( k) k=O
(10)
140
M.-J. JANG AND C.-L. CHEN
where M(k) ~ O, q(t) ~ O. M(k) is called the weighting factor and q(t) is regarded as a kernel corresponding to x(t). If M(k) = 1 and q(t) = 1 then (6) and (9) are equivalent. In this way, (7) can be treated as a special case of (10). In this paper, the transformation with M(k) = Hk/k! and q(t) = 1 is applied, where H is the time horizon of interest and the resulting transformation is called the Taylor transformation. The fundamental properties of the Taylor transformation are given in the Appendix. Using the differential transformation, a differential equation in the time domain can be transformed to be an algebraic equation in the K domain and x(t) can be obtained by finite-term Taylor series plus a remainder, as
1 ~ l ( t - _ k ~ O ) k x(k) x(t) = q(t) k=o M(k----) + RN( t)
(11)
where
1
~ (t--/°)kk
RN( t) = q(t) •
-
k=N
X(k)
M(k)
and RN(t) --* O, as N--* ~ within the time interval of interest, say t [0, H]. If RN(t) is small enough, then x(t) can be represented by finite terms. Using this method, differential equations, either linear or nonlinear, can be solved. For practical problems, the desired H is not necessarily small, therefore, it is required to divide the entire domain H into many sub-domains for good approximation. In this section, the domain split process is described as follows. Suppose the time domain of interest H is split into n sub-domains, as shown in Figure 1. In every sub-domain, there exists a solution function x~(ti), i = 1, 2 , . . . , n, to specify the solution of x(t) for t [y:}~l hi, F.}= 0 hi] with h 0 = 0. Therefore, the Taylor series of x(t) described using the local time scale t, in the first time interval can be shown as x ( t ) = xl( t 1) =
~
)(1(k)
for
t 1 E [0, hi] and t ~ [0, hi]
k=0
(12) where the initial condition is Xl(0) = XI(0 ). The value of xl(t 1) at t 1 = h 1 can be obtained from (12) as N-1
z1(4) = E z,(k) k=O
(13)
Differential Transformation Technique
141
x(t)
v
I
x l(t 1)
x(t )
j/
Time t P local time scale t n local time scale t 2 local time scale t 1
FlG. 1. Domain split of time interval H.
which is the initial value of x2(t 2) in the second sub-domain, i.e., x2(0) = Xl( h 1). The initial values of the i-th sub-domain can be obtained as N-1 xi(O)
~- Xi-l(O) "[- E
Xi-l(k).
(14)
k=l
In this way, xn(t) n can be obtained sequentially and x(t) will be composed of these functions. However, if there exists a calculation error in the previous domain interval then this error will influence the current interval. To prevent the calculation error being accumulated during the calculation process, the accuracy check is performed as follows. In the i-th interval,
xi( t,) = xi(O ) + ~_,
k=l ~ hi]
X~( k).
(15)
142
M.-J. JANG AND C.-L. CHEN
Therefore, • ,_ 1(o) -- ~,( - a,_ 1)
N-l[ --hz_i t k = xi(O ) + ~., Xi(k ). k=i~ h, ]
(16)
The above equation indicates that the initial condition of the previous interval can be expressed by the spectrum of the current interval. After the spectrum in the i-th sub-domain is obtained, the right hand side of (16) can be determined. If the result is not close enough to the left hand side of (16), the transformation is redone by increasing the number power series terms or readjusting the interval length h i. 4. ANALYSIS RESULTS AND DISCUSSIONS Using the transformation technique introduced in section 3, the Taylor transform of (4) with initial conditions x(0)= 0 and k(0)= 0 can be obtained as A 2 x ( k ) "~- 2")/ A X ( k ) + a 2 x ( k ) + 8'~lX(2)(k) $ ^ x(k)
+ s x(3)(k) = E ( k )
(17)
with X(0) = x(0) = 0
(18)
X(1) = hx(0) = 0.
(19)
and
From the operation properties given in the appendix, the corresponding spectrums can be determined as X( K + 2) =
h2
( k + 1 ) ( k + 2)
- 8"/1 {l__~o X ( 2 ) ( l ) ( k - l +h
( k + 1) E ( k ) - 2T----~---- X ( k ) - a 2 X ( k )
1) X(k - l + 1 ) } - 6X(3)(k)}. (20)
143
Differential Transformation Technique
Using (20), the solution x(t) corresponding to different cases is obtained by the spectrum X(k).
x
(tj =
+
E
l/t,t
In the following numerical studies, h = 0.3333 and N = 4. Case 1: Step function external excitation For E(t) = U/-/(t), where U is a constant and H ( t ) is the heaviside step function, then E ( k ) = U~(k) in (20). The time responses of the underdamped and overdamped systems are considered as follows. U n d e r d a r a p e d c a s e s (27 = O. I)
For 71 = 0.1, 6 = 1, ~ = 1 and U = 2, the response of x(t) is shown in Figure 2. For 71 = 0 , 8 = 2, ~ = 1 and U = 2, the response of x(t) is shown in Figure 3. Compare with the results obtained using the fourth order Runge-Kutta method, the proposed results axe very satisfactory. O v e r d a m p e d cases (27 = 2.5)
For "~1 ~- 0, E = 1 , ~-~ = 1 and U = 2, the response of x(t) is shown in Figure 4. In this case the result proposed in [6] is also given. The results show that the proposed method and the Runge-Kutta method achieve the
1
o/ 0
-
-
Runge-Kutta o Taylor Transformation
5
10
15
Time (sec) F I G . 2.
Time response for ~ + 0.1(1 + x2)k + x + x 3 = 2 H(t).
144
M.-J. JANG AND C.-L. CHEN
,.,,.,,
0.5
--
Runge-Kutta
o Taylor Transformation i
5
10 Time (sec)
FIG. 3. Time response for
~+
0.15+ x + 2x 3 = 2H(t).
s a m e result. F o r 7 1 = 2.5, s = 1, ~ = 1 a n d U = 2, t h e response of x ( t ) is shown in F i g u r e 5. C o m p a r e d w i t h t h e n u m e r i c a l i n t e g r a t i o n results obt a i n e d using t h e fourth order R u n g e - K u t t a m e t h o d , t h e p r o p o s e d results are also v e r y s a t i s f a c t o r y w i t h good c o m p u t a t i o n a l efficiency, i.e., w i t h large h a n d small N.
1.2 1 - - Runge-Kutta A
0.8
o Taylor Transformation
v X
- - EI-Asir et al. 0.6
0.4 0.2
i
0
5
10
T i m e (sec)
FIG. 4. Time response for ~ + 2.5~+ x + x3 = 2H(t).
15
Differential Transformation Technique
145
--
A
Rumba
o TaylorT r a ~ I ~ i ~
0.8
0.6
X
0.5
0.4
0
i
0
i
5
10
15
Time (sec)
FIG. 5. Time responsefor ~ + 2.5(1 + x2) J: + x +
X3 = 2 H ( t ) .
Case 2: Periodical function external excitation For
E(t) =
U
7~/,.
cos{~-+
cos(t),
where
\
(i-1)h)and/denotes
\ - -
U
is
a
constant,
hk E(k) = -k~
the i-th split domain in (20). The
/
x - x' phase plane plot of the underdamped and overda.mped systems are considered.
Underdamped cases (2y = O.1) For y, = 0.1, 8 = 1, ~ = 1 and U = 1,4,7 and 10, the corresponding phase plane plot axe shown in Figure 6. Overdamped cases (2~ = 2.5) For " / l = 0 , e = 1, ~ = 1 and U = 1,4,7 and 10, the corresponding phase plane plot is shown in Figure 7, where the limit cycle motion of the system has also been determined. 5.
CONCLUSIONS
The differential transformation technique is applied to solve a strongly nonlinear differential equation that describes the underdamped and the overdamped motion of a system subjected to external excitations. Numerical
146
M.-J. JANG AND C.-L. CHEN
U= 1 x'(t)
x'(t)0
0
-1
-2 -2
-1
0
1
2
-3
-2
x(t)
-1
0
1
2
x(t)
6
101 =
U=IO
4
2 x'(t)
0 -2
-3
-2
-1
o
,
0
1
x(t)
101
2
3
-4
,
-2
0
2
4
x(t)
FIG. 6. Phase plot for ~ + 0.1(1 + x2)~ + x + x3 = Ucos(t).
results are c o m p a r e d to those o b t a i n e d b y t h e R u n g e - K u t t a m e t h o d to illustrate t h e effectiveness of t h e p r o p o s e d m e t h o d . It i s shown t h a t t h e differential t r a n s f o r m a t i o n m e t h o d is v e r y p r o m i s i n g in nonlinear s y s t e m analysis. APPENDIX T a b l e 1 lists some differential t r a n s f o r m . T h e s y m b o l " ^ " denotes t h e differential o p e r a t o r a n d " * " denotes t h e convolution o p e r a t o r in K domain.
Differential Transformation Technique
147 1.5
0.6
U=I 1
0.4
0.5 0.2 x'(t)
x'(t)
o
-0. -0.2
-1
-0.4 -0.5
0 x(t)
0.5
-1.5 -1.5
-0.5
.5
1
1
2
1.5
x(t) 3
3
U
U=7 2
2
1
1
x'(t)
-1
x'(t)
F
0
0
-1 -2
-2 -3 -2
-3
|
-1
o
x(t)
1
-3
2
-2
-1
0
x(t)
FIG. 7. Phase plane plot for ~+ 2.5~+ x + x3 = Ucos(t).
F r o m Table 1, the simplest integral calculation in differential transformation is M(k) = Hk/k! and q(t) = 1, because the binomial coefficient C~ in integral operator is vanished and its differential calculation can be obtained easier c o m p a r e d with the other forms.
,Qk=
k~
t (k- t)!
Basic function operators in K domain are listed in Table 2.
148
M.-J. J A N G
-~
I-.
8 ~
II
~
II
~
z I
© g~
II
II
a:
8~,,,]~
iw
H
tl ~
Z Z
Z
l
© o
~1 II
II
©
©
II
.o
II
.~e I N
<
AND
C.-L. C H E N
Differential Transformation Technique
149
TABLE 2 TIME DOMAINOPERATORAND K DOMAINOPERATOR Time domain Coperator
±
X
d d-~
( )m
K domain operator
±
*
A
()(m)
The Operation Properties of Differential-Transform with q(t) = 1 M(k) = H k / k! f ( t ) and g(t) are two uncorrelated functions with time t and F(k), G(k) are the transformed functions corresponding to ](t), g(t) and the basic properties are shown as follows. 1. Linearity If F ( k ) = D [ f( t)], G(k) = D[ g( t)], and c 1 and c~ arc independent of t and k, then
D [ e l f ( t ) + c2g(t)] = elF(k ) + c2G(k ). 2. Convolution If z(t) = f ( t ) g ( t ) , f ( t ) = D - l [ F ( k ) ] and g(t) = D - l [ G ( k ) ] then D[ z ( t ) ] = D[ f ( t ) g ( t ) ] = F(k) * G(k) k
Y'. F ( l ) G( k - 1).
(A.1)
l=O Therefore, the transform of fro(t), m is a positive integer, can be obtained as
D[ fro( t)] = = F(m- 1)(k) * F ( k ) k
F(m-1)( l) F( k - l). /=0
(A.2)
150
M.-J. JANG AND C.-L. CHEN
3. Integration of ](t) If z(t) = ]~ f(t) dt and 2(k) = /9[ z( t)] then D[ z(t)] = A - 1 F ( k ) H
= -~F(k-
1) + Cog(k )
(A.3)
where co is obtained by boundary condition. 4. Derivative If f ( t ) and its derivatives f'(t), f ' ( t ) , . . . , f(")(t) be continuous functions for time interval [0, H], then D[
d~f(t) I (k+ 1)(k+2)'"(k+ n) H~
F( k + n).
(A.4)
5. Shifting If F~(k) = D{f(t + 7)} and z > 0, then
f,(k)
= m=O
n=O
whereS(k-
n) =
1
0
k=n
kq: n.
(A.5)
Therefore, if n > m, the binomial coefficient C~ is zero and F~(k) can be rewritten as
Ft( k) = ~,
C~'F( m).
(A.6)
m=k
From (A.1)-(A.6), it is found that the linear properties are preserved in the differential-transformation operators. Integral or differential calculations in time domain can be done by algebraic operations in K domain. In this way, computation time can be reduced significantly and accuracy and convergence speed for solving the differential equations can be improved significantly. See Table 3.
Differential Transformation Technique
151
TABLE 3 TRANSFORMED FUNCTIONS
Time function, ~t)
Transformed function F(k)
u(t) t t TM
~(k) H 6 ( k - 1) H ' B ( k - m) (AH) k
e At
(1 + t) ~
k~ Hk __k-ira(m- 1).- ( m - k + 1)
sin(wt + a)
k!
cos(wt + ~)
k------if-,cos T
sin
+ a + a
REFERENCES 1 A.H. Nayfeh Perturbation Method~ John Wiley, New York, 1963. 2 P . G . D . Barham and A. C. Soudack, International Journal of Control 10:377, (1969). 3 A. C. Soudack and P. G. D. Barham, International Journal of Control 13:767, (1971). 4 S.B. Yuste and J. D. Bejarano, Journal of Sound and Vibration 139:151, (1990). 5 Z. Xu and Y. K. Cheung, Journal of Sound and Vibration 174:563, (1994). 6 B. E1-Asir, A. Mansour, A. Abu-Mowis, and N. Kassab, Journal of Sound and Vibration 174:115, (1994). 7 A.R. Mansour and A. M. Hussein, International Communications on Heat and Mass Transfer 17:823, (1990). 8 J. K. Zhou, Differential Transformation and Its Applications for Electrical Circuit~ Wuhahn Huarjung University Press, China (in Chinese), 1986.