Page 336 - MATLAB an introduction with applications

Direct Numerical Integration Methods ——— 321

6.2.2 Central Difference Method
Let the duration over which the numerical solution of Eq. (6.1) is required be divided into n equal parts of

( =
interval h =∆t each. The initial conditions are assumed to be given by Xt 0) = X and Xt = 0) = X .
(
0
0
The accuracy of the solution always depends on the size of the time step. The critical time step is given by
∆ cri = τt n /π , where τ is the natural period of the system. If ∆ t is selected to be larger than ∆ t , the
n
cri
method becomes unstable, that is the truncation of higher order terms or rounding-off in the computer
causes errors to grow and makes the dynamic response calculations meaningless. Numerical methods, which
require the use of time, step ∆t smaller than the critical time step ∆ t are said to be conditionally stable.
cri

A safe rule to use is to choose ≤τ n /10 . Substituting Eqs. (6.5) and (6.4) for X and X respectively in
h
i
i
Eq. (6.1) at mesh or grid point i give

 X − 2X + X 1   X − X 1 

M  i+ 1 i 2 i−  + C  i+ 1 i−  + KX = F i ...(6.6)
i
  ( ) t ∆    2 t ∆ 
() and F =
where X = X t i i F () t . Solving Eq. (6.6) for X i+ 1 gives
i
i
 
 1   2M   C M  


X
X i+ 1 =    − K  { } +  −  X i− 1 + F i  ...(6.7)
i
t ∆
t ∆
 M + C    () 2   2 t ∆ () 2   
 () 2 2 t ∆ 
t ∆


which is known as the recurrence formula.
Thus, the displacement of the mass, X i +1 , can be calculated using Eq. (6.7) if we know the previous
displacements, X , X i +1 and the present external force F . Repeated application of Eq. (6.7) gives us the
i
i
response time history of the behaviour of the system. Since the solution of X i+1 given in Eq.(6.7) is based
on the previous displacement X given in Eq. (6.6), the integration procedure is known as an explicit
i
integration method. Note that in order to compute X , both X and X are required but the initial conditions
–1
0
1

provide only the values of X and X . Therefore, the central difference method is not self-starting. However,
0
0

the value of X can be obtained from Eqs. (6.4) and (6.5) as follows. First calculate the value of X by
–1
0

substituting the given values of X 0 and X in Eq. (6.1) as follows,
0

X = 1   F (t = 0) – CX 0 – KX  ...(6.8)
0 
0
M
The value of X is then obtained by the application of Eqs. (6.4) and (6.5) at i = 0.
–1

X –1 = X − t X∆ 0 + ( t∆ ) 2 X 0 ...(6.9)
0
2
6.2.3 Runge-Kutta Method
In the Runge-Kutta method, an approximation to X t+∆ t is obtained from X in such a way that the power
t
series expansion of the approximation coincides, up to terms of a certain order () t∆ N in the time interval t ∆ ,
with the actual Taylor series expansion of (t +∆ in powers of t∆ . The method is based on the assumption
) t
that the higher derivatives exist at points required.