Vietnam Journal of Mechanics, VAST, Vol. 43, No. 2 (2021), pp.

171 – 182
DOI: https://doi.org/10.15625/0866-7136/15758

ON TWO IMPROVED NUMERICAL ALGORITHMS FOR

VIBRATION ANALYSIS OF SYSTEMS INVOLVING
FRACTIONAL DERIVATIVES

Nguyen Van Khang1,∗ , Duong Van Lac1 , Pham Thanh Chung1

1
Hanoi University of Science and Technology, Vietnam
∗
E-mail: [email protected]

Received: 15 December 2020 / Published online: 21 June 2021

Abstract. Zhang and Shimizu (1998) proposed a numerical algorithm based on Newmark
method to calculate the dynamic response of mechanical systems involving fractional
derivatives. On the basis of Runge–Kutta–Nyström method and Newmark method, the
present study proposes two new numerical algorithms, namely, the improved Newmark
algorithm using the second order derivative and the improved Runge–Kutta–Nyström al-
gorithm using the second order derivative to solve the fractional differential equations
of vibration systems. The accuracy of new algorithms is investigated in detail by nu-
merical simulation. The simulation result demonstrated that the Runge–Kutta–Nyström
algorithm using the second order derivative for the vibration analysis of systems involv-
ing fractional derivatives is more effective than the Newmark algorithm of Zhang and
Shimizu.

Keywords: vibration, fractional differential equation, numerical algorithm, dynamical sys-

tems.

1. INTRODUCTION
A differential equation is called the fractional differential equation if it includes at
least one fractional-order derivative in the expression. Ordinary differential equations
involving fractional differential operators of Riemann–Liouville’s type or of Caputor’s
type are known to have many potential applications in mathematical modeling, in areas
like mechanics, and the life sciences [1–9].
Among the approximate methods for finding a solution of nonlinear fractional dif-
ferential equations, the decomposition method and the numerical method are often used.
The decomposition method is a nonnumerical method for solving nonlinear differential
equations [10–15]. The method was developed by George Adomian in 1984. Essentially,
it approximates the solution of a non-linear differential equation with a series of func-
tions. This method is getting into use for the solution of fractional differential equations.
By using the decomposition method, one needs to express nonlinear terms in the form
© 2021 Vietnam Academy of Science and Technology
172 Nguyen Van Khang, Duong Van Lac, Pham Thanh Chung

of power series. That is a difficult problem for many nonlinear fractional differential
equations.
The use of numerical algorithms for the solution of differential equations involving
fractional derivatives has been discussed in several works [16–28]. Yuan and Agrawal
[22] have rewritten the definition of a fractional derivative and turned a fractional dif-
ferential equation into a system of linear differential equations. However, Schmidt and
Gaul [23] have shown that in some cases, this method loses the advantages of fractional
calculus over integer-order calculus.
Zhang and Shimizu [28] presented the numerical method for dynamic problems in-
volving fractional operators. Using the idea of Zhang and Shimizu, a new algorithm is
developed by incorporating one-step schemes of well-known Newmark types [29] into
its formula. Further, based on Riemann–Liouville’s definition of fractional derivatives
and the well-known Runge–Kutta–Nyström numerical method for calculating the solu-
tion of differential equations [30], we present a new algorithm for calculating nonlinear
fractional differential equations. It is shown that the proposed algorithm is very efficient
in many cases.
This study is organized into four sections. In Section 2, we present three numerical al-
gorithms for solving fractional differential equations, including a well-known algorithm
and two new algorithms. In Section 3, the effectiveness of the numerical algorithms is
studied in detail. Section 4 includes some concluding remarks of the study.

2. SOME NUMERICAL ALGORITHMS FOR CALCULATING RESPONSES OF

MECHANICAL SYSTEMS INVOLVING FRACTIONAL DERIVATIVES
2.1. Preliminaries
Fractional integrals and derivatives are deduced from the generalization of the integer-
order operations. It is usual to define the integral operator D −q as

Zt
−q 1
a Dt x ( t ) = (t − τ )q−1 x (τ )dτ, (1)
Γ(q)
a

where q > 0 and Γ( x ) is the Gamma function

Z∞
Γ( x ) = e−z z x−1 dz. (2)
0

For a continuous x (t),

D − p D −q x (t) = D −( p+q) x (t), (3)
as given in [3] (if both p and q are non-negative).
With the fractional integral operator, fractional derivatives are easily introduced. For
a real α > 0, D α is defined by the Riemann–Liouville definition [3], using the above
 On two improved numerical algorithms for vibration analysis of systems involving fractional derivatives 173

fractional integral operator

Zt
!
dn d−(n−α) x (t) 1 dn
α
a Dt x ( t ) = n = (t − τ )n−α−1 x (τ )dτ. (4)
dt dt−(n−α) Γ(n − α) dtn
a

Another choice is the Caputo definition

Zx
dn
 
C α 1 n − α −1
a Dt x ( t ) = (t − τ ) x (τ ) dτ. (5)
Γ(n − q) dτ n
a

In both cases (n − 1) < α < n.

Actually, the two definitions only differ in the consideration of conditions at the start
of the interval
n −1
1 Γ(n − α)
α C α
a Dt x ( t ) = a Dt x ( t ) + ∑
Γ ( n − α ) k =0 Γ ( k − α + 1 )
( t − a ) k − α x ( k ) ( a ). (6)

The differential equation to be solved is the vibration equation with fractional damp-
ing, with one degree of freedom
m ẍ (t) + b ẋ (t) + µc( x ) D α x (t) + g( x ) = f (t), 0 < α < 1, (7)
with the initial conditions
x (0) = x0 , ẋ (0) = ẋ0 . (8)
The existence and uniqueness of the solutions of Eq. (7) are presented in [4]. Note
that this study focuses on developing numerical algorithms for solving this equation.
In the applications, D practically always means 0 Dt , and most authors use the
Riemann–Liouville, or the mathematically equivalent Gruenwald–Letnikov definition
(see [3] for precise conditions of equivalence). Also, since the Riemann–Liouville def-
inition has a singularity for non-zero initial conditions, the initial conditions are often
considered zero. For a physical interpretation of this singularity, see [5–7].
Using the step-size
h = ∆t = ti − ti−1 , (9)
we have
tn = t0 + nh = t0 + n∆t, n = 1, 2, 3, . . . (10)
Using the notations x (ti ) = xi , from Eq. (7) we have the following iterative compu-
tational scheme at the time tn as follows
m ẍn + b ẋn + µc( xn ) D q xn + kxn = f (tn ), n = 1, 2, 3, . . . (11)
2.2. The Newmark-based algorithm proposed by Zhang and Shimizu: A review
Based on the single-step integration method by Newmark [29], Zhang and Shimizu
(1998) have obtained the following approximation formulas [28]
 
1 1 1
ẍn = ( x n − x n −1 ) − ẋn−1 − − 1 ẍn−1 = ψ2 ( β, ẋn−1 , ẍn−1 , xn−1 , xn ) , (12)
β∆t2 β∆t 2β
174 Nguyen Van Khang, Duong Van Lac, Pham Thanh Chung

   
α α
ẋn = ẋn−1 + (1 − α)∆t ẍn−1 + α∆t ẍn = 1− ẋn−1 + 1 − ∆t ẍn−1
β 2β
(13)
α
+ ( xn − xn−1 ) = ψ1 (α, β, ẋn−1 , ẍn−1 , xn−1 , xn ) .
β∆t
Constants α, β are parameters associated with the quadrature scheme [28, 29]. The
numerical algorithm to calculate the fractional derivative at t = tn of Eq. (4) is
t 
Zn−1 Ztn
x ( t 0 ) −q 1 ẋ ( τ ) dτ ẋ ( τ ) dτ
D q x (tn ) = tn + +
Γ (1 − q ) Γ (1 − q ) (tn − τ )q (tn − τ )q
 
0 t n −1 (14)
1
= ( I0 + In−1 + ∆In ) ,
Γ (1 − q )
where we denote
tZn−1 " #
n −2
ẋ (τ )dτ h ẋ0 ẋn−1 ẋ (ih)
In−1 = q ≈ q + +2 ∑ q , (15)
(tn − τ ) 2 tn hq i =1 ( tn − ih )
0
Ztn
ẋ (τ )dτ ∆t1−q ∆t2−q ∆t2−q
∆In = q = ẋn−1 + (1 − α) ẍn−1 + α ẍn .
(tn − τ ) 1−q (1 − q)(2 − q) (1 − q)(2 − q)
t n −1
(16)
By substituting ẍn in Eq. (12) into Eq. (16) we obtain
Ztn
ẋ (τ )dτ
∆In =
(tn − τ )q
t n −1 (17)
∆t1−q
     
α α α
= ( x n − x n −1 ) + 2 − q − ẋn−1 + 1 − ∆t ẍn−1 .
(1 − q)(2 − q) β∆t β 2β
Then, from Eq. (14) we have
1
D q x (tn ) = ( I0 + In−1 + ∆In ) = ψq (α, β, ẋ0 , ẋ1 , ẋ2 , . . . , ẋn−1 , ẍn−1 , xn−1 , xn ) .
Γ (1 − p )
(18)
From Eqs. (12), (13) and (18) we can rewrite Eq. (11) in the following form
mψ2 ( β, ẋn−1 , ẍn−1 , xn−1 , xn ) + bψ1 (α, β, ẋn−1 , ẍn−1 , xn−1 , xn )
(19)
+ µc( xn )ψq (α, β, ẋ0 , ẋ1 , ẋ2 , . . . , ẋn−1 , ẍn−1 , xn−1 , xn ) + kxn = f (tn ).
Eq. (19) is a nonlinear algebraic equation to find unknown xn . We can then calculate
ẋn , ẍn according to the following formulas
ẋn = ẋn−1 + (1 − α)∆t ẍn−1 + α∆t ẍn ,
(20)
 
1 1 1
ẍn = 2
( x n − x n −1 ) − ẋn−1 − − 1 ẍn−1 .
β∆t β∆t 2β
 On two improved numerical algorithms for vibration analysis of systems involving fractional derivatives 175

2.3. The improved Runge–Kutta–Nyström (RKN) algorithm

From the definition of the Liouville–Riemann’s fractional derivative, Eq. (4), we can
apply the composition rule to D q x (t) [1–3], that is
Ztn Ztn
q 1 d x (τ ) x ( t0 ) − q 1 ẋ (τ )
D x (tn ) = q dτ = tn + dτ
Γ (1 − q) dt (tn − τ ) Γ (1 − q ) Γ (1 − q ) (tn − τ )q
t0 0
 
t n
x ( t0 ) − q 1 tn
Z
= tn + − ẋ (τ ) (tn − τ )1−q + ẍ (τ ) (tn − τ )1−q dτ 
Γ (1 − q ) Γ (2 − q ) 0
0
 
t n
x ( t0 ) − q 1
Z
= tn +  ẋ (t0 )t1n−q + ẍ (τ )(tn − τ )1−q dτ 
Γ (1 − q ) Γ (2 − q )
0
1 1
= I0 + ( J0 + J (tn )) ,
Γ (1 − q ) Γ (2 − q )
(21)
where we denote
Ztn Ztn
1− q 1− q
J0 = ẋ (t0 )tn , J (tn ) = ẍ (τ )(tn − τ ) dτ = ytn (τ )dτ. (22)
0 0

We approximate the integrals according to Eq. (22) for every instance tn by trapezoid
numerical integration with an accuracy of O(t3 ) as follows
J (t0 ) = 0,
n −1  n −2 h 
h  h
J (tn ) ≈ ∑ 2
ytn (τj ) + ytn (τj+1 ) = ∑
2
ytn (τj ) + ytn (τj+1 ) + ytn (tn−1 ),
2
( n ≥ 1),
j =0 j =0
  n −1
h h  h
J tn +
2
≈ ∑ 2
ytn +h/2 (τj ) + ytn +h/2 (τj+1 ) + ytn +h/2 (tn ),
4
( n ≥ 0),
j =0
(23)
Thus, the formulas for determining the level fractional derivatives at tn , tn+h2 and tn +h
have the following forms
D q x (tn ) = ψ1 ( x0 , ẋ0 , ẍ0 , ẍ1 , ẍ2 , . . . , ẍn−1 ) ,
D q x (tn + h/2) = ψ2 ( x0 , ẋ0 , ẍ0 , ẍ1 , ẍ2 , . . . , ẍn−1 ) , (24)
q q
D x (tn + h) = D x (tn+1 ) = ψ3 ( x0 , ẋ0 , ẍ0 , ẍ1 , ẍ2 , . . . , ẍn ) .
From the approximation formula determining D q x (tn ) at step t = tn , Eq. (7) can be
rewritten in the following form
m ẍn + b ẋn + µc( xn )ψ1 ( x0 , ẋ0 , ẍ0 , ẍ1 , ẍ2 , . . . , ẍn−1 ) + kxn = f (tn ),
1 (25)
ẍn = ( f (tn ) − b ẋn − µc( xn )ψ1 ( x0 , ẋ0 , ẍ0 , ẍ1 , ẍ2 , . . . , ẍn−1 ) − kxn ) = g(tn , xn , ẋn ),
m
176 Nguyen Van Khang, Duong Van Lac, Pham Thanh Chung

It should be noted that when implementing expansion (21) we have used the assumption
Zt
ẋ (τ )
on convergence of integral I = dτ if 0 < q < 1. Indeed, since the function ẋ (t)
(t − τ )q
0
is continuous on each finite time interval, so we have
Zt Zt
ẋ (τ ) dτ
| ẋ (τ )| ≤ M ⇒ I = dτ ≤ M . (26)
(t − τ )q (t − τ )q
0 0

Zt
ẋ (τ )
From the above inequality we deduce that the integral I = dτ converges if
(t − τ )q
0
0 < q < 1. Then we can develop
 
Ztn Ztn
ẋ (τ ) 1  tn
I= dτ = − ẋ (τ ) (tn − τ )1−q + ẍ (τ ) (tn − τ )1−q dτ 
(tn − τ )q 1−q 0
0 0
  (27)
t n
1 
Z
1− q 1− q
= ẋ (t0 )tn + ẍ (τ )(tn − τ ) dτ  .
1−q
0

Applying the Runge–Kutta–Nyström algorithm with an accuracy of O(t4 ) to the differ-

ential equation (25), we have a straightforward schema [30] as below.
ẍn = g(tn , xn , ẋn ), x ( t0 ) = x0 , ẋ (t0 ) = ẋ0 , (28)

h
xn+1 = xn + hẋn + (k1 + k2 + k3 ) ,
3
1 (29)
ẋn+1 = ẋn + (k1 + 2k2 + 2k3 + k4 ) ,
3
ẍn+1 = g(tn + h, xn+1 , ẋn+1 ).
where
h
k1 =
g (tn , xn , ẋn ) ,
2  
h h h h
k2 = g tn + , xn + ẋn + k1 , ẋn + k1 ,
2 2 2 4
  (30)
h h h h
k3 = g tn + , xn + ẋn + k1 , ẋn + k2 ,
2 2 2 4
h
k4 = g (tn + h, xn + hẋn + hk3 , ẋn + 2k3 ) .
2
2.4. The improved Newmark algorithm
Using the second-order derivative by numerical integral we propose an algorithm
based on the well-known Newmark algorithm [28, 29] to find the solution of Eq. (7).
 On two improved numerical algorithms for vibration analysis of systems involving fractional derivatives 177

Firstly, we have rewritten the Eqs. (12) and (13) as follows

 
1 1 1
ẍn = ( x n − x n −1 ) − ẋn−1 − − 1 ẍn−1 = ψ2 ( β, ẋn−1 , ẍn−1 , xn−1 , xn ) , (31)
β∆t2 β∆t 2β
ẋn = ẋn−1 + (1 − α)∆t ẍn−1 + α∆t ẍn
   
α α α
= 1− ẋn−1 + 1 − ∆t ẍn−1 + ( x n − x n −1 ) (32)
β 2β β∆t
= ψ1 (α, β, ẋn−1 , ẍn−1 , xn−1 , xn ) .
The numerical algorithm to calculate the fractional derivative at t = tn of Eq. (7) is
 
Ztn
x ( t0 ) − q 1  ẋ (t0 )t1n−q + ẍ (τ )(tn − τ )1−q dτ 
D q x (tn ) = tn +
Γ (1 − q ) (1 − q ) Γ (1 − q )
0 (33)
1 1
= I0 + ( J0 + J (tn )) ,
Γ (1 − q ) (1 − q ) Γ (1 − q )
where
Ztn Ztn
1− q 1− q
J0 = ẋ (t0 )tn , J (tn ) = ẍ (τ )(tn − τ ) dτ = ytn (τ )dτ. (34)
0 0
Similarly, we approximate the integrals in Eq. (34) for every instance tn by trapezoid
numerical integration as follows
J (t0 ) = 0,
n −1  n −2 h 
h  h
J (tn ) ≈ ∑ 2
ytn (τj ) + ytn (τj+1 ) = ∑
2
ytn (τj ) + ytn (τj+1 ) + ytn (tn−1 ),
2
( n ≥ 1)
j =0 j =0

⇒ D q x (tn ) = ψq ( x0 , ẋ0 , ẍ0 , ẍ1 , ẍ2 , . . . , ẍn−1 ) .

(35)
From Eqs. (31), (32) and (35), we can rewrite Eq. (7) in the following form
mψ2 ( β, ẋn−1 , ẍn−1 , xn−1 , xn ) + bψ1 (α, β, ẋn−1 , ẍn−1 , xn−1 , xn )
(36)
+µc( xn )ψq ( ẋ0 , ẋ1 , ẋ2 , . . . , ẋn−1 , xn−1 ) + kxn = f (tn ).
Eq. (36), a nonlinear algebraic equation of an unknown xn , can be solved by the Newton–
Raphson method of iteration. The variables ẋn , ẍn can then be determined by
ẋn = ẋn−1 + (1 − α)∆t ẍn−1 + α∆t ẍn ,
(37)
 
1 1 1
ẍn = 2
( x n − x n −1 ) − ẋn−1 − − 1 ẍn−1 .
β∆t β∆t 2β

3. NUMERICAL RESULTS
To compare the accuracy of the numerical algorithms, these motion equations of a
one-degree-of-freedom oscillator would be considered to evaluate.
 Raphson method of iteration. The variables x!n , !!xn can then be determined by
x!n = x!n -1 + (1 - a )Dtx!!n -1 + aDtx!!n
1 1 æ 1 ö (37)
2 ( n
xn =
!! x - xn -1 ) - x!n -1 - ç - 1÷ !!
xn -1
bDt bDt è 2b ø
3. NUMERICAL RESULTS
178 To compare the accuracy
NguyenofVan
theKhang,
numerical algorithms,
Duong Van these Chung
Lac, Pham Thanh motion equations of a one-
degree-of-freedom oscillator would be considered to evaluate.
Example 1: Consider the following system
Example 1: Consider 0.5the following system
x ( t ) + 0.8 D
!! 3
x (t ) + x = f (t ) (38)
ẍ (t) + 0.8D0.5 x (t) + x3 = f (t) , (38)
where
where æ 9 öæ
f 9( t )= t - 7 ö+ 4t æ t - 7 ö + 4t æ
2 ç t - 7 ÷ç
9 ö 2
÷ ç 7 ÷ ç t - ÷9+ 2t

f (t) = 2 t− tè − 10 øè +104tø t −è 10 ø+ 4tè t 10 −ø + 2t2
(39)
10 10 10 10 3
8 æ 128 7 128 5 42 3 ö  é 2 æ  9öæ 7 öù  (39)
+ √ t - √t + 42 t√÷ + êt ç t - ÷ç t - 9 ÷ ú
) çè 35 t7 − 128
3
8 10G ( 0.5128

25 t5 +25 øt3 ë +è t210 øè 7
+ t − 10 ø û t − ,
10Γ 0.5) 35are
and the initial( conditions
25 25 10 10
x ( 0 ) = 0, are
and the initial conditions x! ( 0 ) = 0 (40)
x (0) = 0, ẋ (0) = 0. (40)
The exact solution of this equation is known as [11, 13, 25] (see Fig. 1).
The exact solution of this equation is known as [11, 13, 25] (see Fig. 1).
æ 9 öæ 7ö
xexact = t 2 ç t - ÷ ç t - ÷  9

7
 (41)
è 10
xexact 10
ø è = øt 2
t− t− . (41)
Fig. 1 shows the exact solution x of Eq. (38). 10 10

Fig. 1. Exact solution of Eq.(38)

Fig. 1. Exact solution x of Eq. (38)

Three numerical algorithms considered 7in Section 2 are then applied to find the ap-
proximate solution of Eq. (38) to evaluate their accuracy. Table 1 and Table 2 show numer-
ical values of the solution to Eq. (38) by our algorithms, exact solution, and the solution
obtained by Ray et al. [13], and Atanackovic et al. [25]. Table 1 shows the five numerical
solutions in comparison with the exact solution and their relative errors in percentage. It
can be seen that the results of the improved Newmark and Newmark method in [28] are
quite similar. However, the improved Runge–Kutta–Nyström algorithm shows the finest
results with the highest accuracy over other algorithms.
 On two improved numerical algorithms for vibration analysis of systems involving fractional derivatives 179

Table 1. The exact solution, numerical solutions, and error in percentage of Eq. (38)
(time step size ∆t = 0.001)

Time xexact xAtanackovic [25] xRay [13] xZhang–Shimizu [28] ximproved Newmark ximproved RKN
0.018253191 0.0182813 0.018507464 0.018508641 0.018281311
0.25 0.01828125
(0.153486%) (0.000274%) (1.237407%) (1.243849%) (0.000332%)
0.019851524 0.0200026 0.020643625 0.020647787 0.020000016
0.5 0.02
(0.742382%) (0.013000%) (3.218126%) (3.238937%) (0.000080%)
−0.004492587 −0.00419593 −0.003348193 −0.003343554 −0.004218918
0.75 −0.00421875
(6.490951%) (0.540919%) (20.635433%) (20.745378%) (0.003975%)
0.029380011 0.0300995 0.030526185 0.030530651 0.029999869
1 0.03
(2.066632%) (0.331667%) (1.753949%) (1.768835%) (0.000437%)
| xmethod − xexact |
Noted: The parentheses (.) indicate the relative errors (i.e., × 100%) of the numerical results, in which
xexact
the lowest error (highest accuracy) is underlined.

Example 2: To further evaluate these methods, let’s consider the below system [31]
ẍ (t) + 0.8D0.5 x (t) + x (t)2 = f (t) , (42)
where
      
3 8 3 8
f (t) = 2 t − t− + 4t t − + 4t t − + 2t2
10 10 10 10
(43)
128 √ 7 88 √ 5 16 √ 3
     2
8 2 3 8
+ √ t − t + t + t t− t− ,
10 π 35 25 25 10 10
and the initial conditions are
x (0) = 0, ẋ (0) = 0. (44)

Fig. 2. Exact solution of Eq. (42)

Fig. 2. Exact solution x of Eq. (42)
Table 2. The exact solution, numerical solutions, and error in percentage of Eq. (42) (time
step size ∆t = 0.001).

Time xexact xZhang -Shimizu ximproved Newmark ximproved RKN

0.25 0.00171875 0.001858277 0.001858493 0.001718750

180 Nguyen Van Khang, Duong Van Lac, Pham Thanh Chung

The exact solution of this equation as indicated in [31] is (see Fig. 2)

  
2 3 8
xexact = t t − t− . (45)
10 10

Table 2. The exact solution, numerical solutions, and error in percentage of Eq. (42)
(time step size ∆t = 0.001)

Time xexact xZhang–Shimizu ximproved Newmark ximproved RKN

0.001858277 0.001858493 0.001718750
0.25 0.00171875
(8.1179263%) (8.1305000%) (0.0000077%)
−0.014694001 −0.014694632 −0.015000092
0.5 −0.015
(2.0399913%) (2.0357850%) (0.0006150%)
−0.012534713 −0.012537936 −0.012656366
0.75 −0.01265625
(0.9602941%) (0.9348281%) (0.0009155%)
0.139197963 0.139202096 0.140000451
1 0.14
(0.5728838%) (0.5699317%) (0.0009155%)

Similar to the results of Example 1, Table 2 shows the five numerical solutions in
comparison with the exact solution and the relative error between them. The improved
RKN also indicates robustness with the highest calculating accuracy. Therefore, the above
comparison of the calculation accuracy leads to the following remark: the accuracy of the
Runge–Kutta–Nyström algorithm using the second-order derivative is very good and
better than the other aforementioned methods.

4. CONCLUSIONS
Using the idea of Zhang and Shimizu [28], two new numerical algorithms for find-
ing the solution of nonlinear fractional differential equations are introduced in this pa-
per. Based on Liouville–Riemann definition of fractional derivatives, and using the well-
known Newmark numerical integration method [29], the well-known Runge–Kutta
–Nyström integration method [30] for differential equations, we proposed two new nu-
merical algorithms for solving the second-order systems containing fractional deriva-
tive components, namely, the improved Newmark algorithm and the improved Runge–
Kutta–Nyström algorithm.
Compared to the aforementioned algorithms, two new algorithms have advantages
in simplicity in the approximation of the fractional derivative components (see Eqs. (23)
and (35)) and calculating scheme (see Eqs. (28)–(30)). The accuracy of these methods
was verified by two examples in Section 3. It reveals that the improved Runge–Kutta–
Nyström is very accurate for the vibration analysis of systems involving fractional deriva-
tives. Noted that the improved RKN method could be easily extended for other systems
with higher degrees of freedom as indicated in [32], and Duffing and Vander Pol systems
as implemented in [33]. However, the problem of convergence and error in the calculat-
ing procedure is required for further investigation.
 On two improved numerical algorithms for vibration analysis of systems involving fractional derivatives 181

ACKNOWLEDGMENTS
This paper was completed with the financial support of the Vietnam National
Foundation for Science and Technology Development (NAFOSTED) under grant number
107.04-2020.28.
REFERENCES
[1] K. B. Oldham and J. Spanier. The fractional calculus. Academic Press, Boston, New York,
(1974).
[2] K. S. Miller and B. Ross. An introduction to the fractional calculus and fractional differential equa-
tions. John Wiley & Sons, New York, (1993).
[3] I. Podlubnv. Fractional differential equations. Academic Press, Boston, New York, (1999).
[4] D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo. Fractional calculus, models and numerical
methods. World Scientific Publishing, Singapore, (2011).
[5] A. Kochubei and Y. Luchko. Handbook of fractional calculus with applications, Volume 2: Frac-
tional differential equations. De Gruyter, Berlin/Boston, (2019).
[6] V. E. Tarasov. Handbook of fractional calculus with applications, Volume 4: Applications in physics,
Part A. De Gruyter, Berlin/Boston, (2019).
[7] V. E. Tarasov. Handbook of fractional calculus with applications, Volume 5: Applications in physics,
Part B. De Gruyter, Berlin/Boston, (2019).
[8] D. Baleanu and A. M. Lopes. Handbook of fractional calculus with applications, Volume 7: Appli-
cations in engineering, life and social sciences, Part A. De Gruyter, Berlin/Boston, (2019).
[9] D. Baleanu and A. M. Lopes. Handbook of fractional calculus with applications, Volume 8: Appli-
cations in engineering, life and social sciences, Part B. De Gruyter, Berlin/Boston, (2019).
[10] G. Adomian. A new approach to nonlinear partial differential equations. Journal of Math-
ematical Analysis and Applications, 102, (1984), pp. 420–434. https://doi.org/10.1016/0022-
247x(84)90182-3.
[11] G. Adomian. A review of the decomposition method and some recent results for
nonlinear equations. Mathematical and Computer Modelling, 13, (7), (1990), pp. 17–43.
https://doi.org/10.1016/0895-7177(90)90125-7.
[12] S. S. Ray, B. P. Poddar, and R. K. Bera. Analytical solution of a dynamic system containing
fractional derivative of order one-half by adomian decomposition method. Journal of Applied
Mechanics, 72, (2005), pp. 290–295. https://doi.org/10.1115/1.1839184.
[13] S. S. Ray, K. S. Chaudhuri, and R. K. Bera. Analytical approximate solution
of nonlinear dynamic system containing fractional derivative by modified decom-
position method. Applied Mathematics and Computation, 182, (2006), pp. 544–552.
https://doi.org/10.1016/j.amc.2006.04.016.
[14] H. Jafari and V. Daftardar-Gejji. Revised Adomian decomposition method for solving sys-
tems of ordinary and fractional differential equations. Applied Mathematics and Computation,
181, (2006), pp. 598–608. https://doi.org/10.1016/j.amc.2005.12.049.
[15] Q. Wang. Numerical solutions for fractional KdV–Burgers equation by Adomian de-
composition method. Applied Mathematics and Computation, 182, (2006), pp. 1048–1055.
https://doi.org/10.1016/j.amc.2006.05.004.
[16] K. Diethelm. An algorithm for the numerical solution of differential equations of fractional
order. Electronic Transactions on Numerical Analysis, 5, (1), (1997), pp. 1–6.
[17] K. Diethelm and N. J. Ford. Analysis of fractional differential equations.
Journal of Mathematical Analysis and Applications, 265, (2002), pp. 229–248.
https://doi.org/10.1006/jmaa.2000.7194.
182 Nguyen Van Khang, Duong Van Lac, Pham Thanh Chung

[18] K. Diethelm and J. Ford. Numerical solution of the Bagley-Torvik equation. BIT Numerical
Mathematics, 42, (3), (2002), pp. 490–507.
[19] K. Diethelm, N. J. Ford, and A. D. Freed. A predictor-corrector approach for thenumerical
solution of fractional differential equations. Nonlinear Dynamics, 29, (1/4), (2002), pp. 3–22.
https://doi.org/10.1023/a:1016592219341.
[20] K. Diethelm and N. J. Ford. Numerical solution of linear and non-linear fractional differential
equations involving fractional derivatives of several orders. Numerical Analysis Report No. 379,
Manchester Center for Computational Mathematics, Manchester, England, (2003).
[21] J. T. Edwards, N. J. Ford, and A. C. Simpson. The numerical solution of linear multi-term
fractional differential equations: Systems of equations. Journal of Computational and Applied
Mathematics, 148, (2002), pp. 401–418. https://doi.org/10.1016/s0377-0427(02)00558-7.
[22] L. Yuan and O. P. Agrawal. A numerical scheme for dynamic systems contain-
ing fractional derivatives. Journal of Vibration and Acoustics, 124, (2002), pp. 321–324.
https://doi.org/10.1115/1.1448322.
[23] A. Schmidt and L. Gaul. On a critique of a numerical scheme for the calculation of fraction-
ally damped dynamical systems. Mechanics Research Communications, 33, (2006), pp. 99–107.
https://doi.org/10.1016/j.mechrescom.2005.02.018.
[24] J.-Z. Wang and et al. Coiflets bases method in solution of nonlinear of dynamic systems
containing fractional derivative. In Proceedings of Fourth international Conference on Nonlinear
Mechanics, (2002).
[25] T. M. Atanackovic and B. Stankovic. On a numerical scheme for solving differential
equations of fractional order. Mechanics Research Communications, 35, (2008), pp. 429–438.
https://doi.org/10.1016/j.mechrescom.2008.05.003.
[26] A. Pálfalvi. Efficient solution of a vibration equation involving fractional deriva-
tives. International Journal of Non-Linear Mechanics, 45, (2010), pp. 169–175.
https://doi.org/10.1016/j.ijnonlinmec.2009.10.006.
[27] J. T. Machado. Numerical calculation of the left and right fractional derivatives. Journal of
Computational Physics, 293, (2015), pp. 96–103. https://doi.org/10.1016/j.jcp.2014.05.029.
[28] W. Zhang and N. Shimizu. Numerical algorithm for dynamic problems involving
fractional operators. JSME International Journal Series C, 41, (3), (1998), pp. 364–370.
https://doi.org/10.1299/jsmec.41.364.
[29] N. M. Newmark. A method of computation for structural dynamics. Journal of the Engineering
Mechanics Division, 85, (1959), pp. 67–94. https://doi.org/10.1061/jmcea3.0000098.
[30] L. Collatz. Numerische behandlung von differentialgleichungen. Springer-Verlag, Berlin, (1951).
https://doi.org/10.1007/978-3-662-22248-5.
[31] S. Banerjee, S. Shaw, and B. Mukhopadhyay. A modified series solution method for fractional
integro-differential equations. International Research Journal of Engineering and Technology (IR-
JET), 3, (8), (2016), pp. 1966–1973.
[32] D. V. Lac. Calculating vibration of systems involving fractional derivatives. Engineering Gradua-
tion Project, Hanoi University of Science and Technology, (2014). (in Vietnamese).
[33] D. V. Lac. Development of Runge-Kutta-Nyström method for calculating vibration of systems involv-
ing fractional derivatives. Master Science Thesis, Hanoi University of Science and Technology,
(2016). (in Vietnamese).