# Application of Hybrid Functions for Solving Duffing-Harmonic Oscillator

Authors: Mohammad Heydari; Ghasem Brid Loghmani; Seyed Mohammad. Hosseini; Seyed Mehdi Karbassi
Journal: Journal of Difference Equations
 (2014)

---

## Abstract

A numerical method for finding the solution of Duffing-harmonic oscillator is proposed. The approach is based on hybrid functions approximation. The properties of hybrid functions that consist of block-pulse and Chebyshev cardinal functions are discussed. The associated operational matrices of integration and product are then utilized to reduce the solution of a strongly nonlinear oscillator to the solution of a system of algebraic equations. The method is easy to implement and computationally very attractive. The results are compared with the exact solution and results from several recently published methods, and the comparisons showed proper accuracy of this method.


---

## Body

## 1. Introduction

Most phenomena in our world are essentially nonlinear and are described by nonlinear ordinary differential equations. Nonlinear oscillation in mechanics, physics, and applied mathematics has been a topic of intensive research for many years. Difficulty of solving the nonlinear problems or getting an analytic solution leads one to use numerical methods. Several methods have been used to find approximate solutions to these nonlinear problems. Some of these well-known methods are harmonic balance method [1], multiple scales method [2], Krylov-Bogoliubov-Mitropolsky method [3, 4], modified Lindstedt-Poincare method [5], linearized perturbation method [6], energy balance method [7], iteration perturbation method [8], bookkeeping parameter perturbation method [9], amplitude frequency formulation [10], maximum approach [11], Mickens iteration procedure [12], rational harmonic balance method [13], Adomian decomposition method [14], variational iteration method [15], modified variational iteration method [16, 17], homotopy perturbation method [18], modified differential transform method [19], and modified homotopy perturbation method [20].Recently, hybrid functions have been applied extensively for solving differential equations or systems and proved to be a useful mathematical tool. The pioneering work in the solution of linear systems with inequality constraints via hybrid of block-pulse functions and Legendre polynomials was led in [21] that first derived an operational matrix for the integrals of the hybrid function vector. Razzaghi and Marzban in [22] the variational problems are solved using hybrid of block-pulse and Chebyshev functions. Razzaghi and Marzban [23] applied the hybrid of block-pulse and Chebyshev functions to find approximate solution of systems with delays in state and control. Solution of time-varying delay systems is approximated using hybrid of block-pulse functions and Legendre polynomials in [24]. Maleknejad and Tavassoli Kajani in [25] introduced a Galerkin method based on hybrid Legendre and block-pulse functions on interval [
0,1
) to solve the linear integrodifferential equation system. Razzaghi and Marzban in [26], a direct method for solving multidelay systems using hybrid of block-pulse functions and Taylor series is presented. Marzban et al. [27] implemented hybrid of block-pulse functions and Lagrange-interpolating polynomials to find approximate solution of Volterra’s population model. The Lane-Emden type equations are solved in [28] using hybrid functions of block-pulse and Lagrange-interpolating polynomials. The hybrid of block-pulse functions and Taylor series is employed in [29] to solve the linear quadratic optimal control with delay systems. Application of hybrid of block-pulse functions and Lagrange polynomials for solving the nonlinear mixed Volterra-Fredholm-Hammerstein integral equations is investigated in [30].In this study, we consider the following nonlinear Duffing-harmonic oscillation [31–35]:
(1)
u
′′
+
u
3
1
+
u
2
=
0
,
u
(
0
)
=
A
,
u
′
(
0
)
=
0
,
where is an example of conservative nonlinear oscillatory systems having a rational form for the restoring force. Note that, for small values of u, (1) is that of a Duffing-type nonlinear oscillator; that is,
(2)
u
′′
+
u
3
≃
0
,
while for large values of u the equation approximates that of a linear harmonic oscillator; that is,
(3)
u
′′
+
u
≃
0
.
Hence, (1) is called the Duffing-harmonic oscillator [31]. The system will oscillate between symmetric bounds [
-
A
,
A
], and the frequency and corresponding periodic solution of the nonlinear oscillator are dependent on the amplitude A [33].In this paper, we introduce an alternative numerical method to solve Duffing-harmonic oscillator. The method consists of reducing this equation to a set of algebraic equations by first expanding the candidate function as a hybrid function with unknown coefficients. These hybrid functions, which consist of block-pulse functions plus Chebyshev cardinal functions, are first introduced. The operational matrices of integration and product are given. These matrices are then used to evaluate the coefficients of the hybrid function for the solution of strongly nonlinear oscillators.The outline of this paper is as follows. In Section2, the basic properties of hybrid block-pulse functions and Chebyshev cardinal functions required for subsequent development are described. In Section 3, we apply the proposed numerical method to the Duffing-harmonic oscillation. Results and comparisons with existing methods in the literature are presented in Section 4 and finally conclusions are drawn in Section 5.


## 2. Properties of Hybrid Functions

Marzban et al. in [27, 29, 30] used the hybrid of block-pulse functions and Lagrange-interpolating polynomials based on zeros of the Legendre polynomials. But no explicit formulas are known for the zeros of the Legendre polynomials. In this study we used Chebyshev cardinal functions which are special cases of Lagrange-interpolating polynomials based on zeros of the Chebyshev polynomials of the first kind to overcome this problem. In this paper, we present the properties of hybrid functions which consist of block-pulse functions plus Chebyshev cardinal functions similar to [27, 29, 30]. The hybrid functions are first introduced, and the operational matrices of integration and product are then derived.


### 2.1. Hybrid Functions of Block-Pulse and Chebyshev Cardinal Functions

Hybrid functionsθ
n
,
m
(
t
), n
=
1,2
,
…
N, m
=
0,1
,
…
,
M
-
1, are defined on the interval [
0
,
t
f
) as
(4)
θ
n
,
m
(
t
)
=
{
C
m
(
2
N
t
f
t
-
2
n
+
1
)
,
t
∈
[
(
n
-
1
N
)
t
f
,
n
N
t
f
)
,
0
,
otherwise
,
where n and m are the order of block-pulse functions and Chebyshev cardinal functions, respectively. Here, C
m
(
t
) are defined as [36–38]
(5)
C
m
(
t
)
=
T
M
(
t
)
T
M
,
t
(
t
m
)
(
t
-
t
m
)
,
m
=
0,1
,
…
,
M
-
1
,
where T
M
(
t
) is the first kind Chebyshev polynomial of order M in [
-
1,1
] defined by
(6)
T
M
(
t
)
=
cos
⁡
(
M
arccos
(
t
)
)
,
subscript t denotes t-differentiation, and t
m, m
=
0,1
,
…
,
M
-
1, are the zeros of T
M
(
t
) defined by cos
⁡
(
(
2
m
+
1
)
π
)
/
2
M
), m
=
0,1
,
…
,
M
-
1, with the Kronecker property
(7)
C
m
(
t
i
)
=
δ
m
i
=
{
1
,
if
i
=
m
,
0
,
if
i
≠
m
,
where δ
m
i is the Kronecker delta function.


### 2.2. Function Approximation

A functionf
(
t
), defined over the interval [
0
,
t
f
), may be expanded as
(8)
f
(
t
)
=
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
c
n
m
θ
n
m
(
t
)
.
If the infinite series in (8) is truncated, then (8) can be expressed as
(9)
f
(
t
)
≃
∑
n
=
1
N
∑
m
=
0
M
-
1
‍
c
n
m
θ
n
m
(
t
)
=
C
T
Θ
(
t
)
,
where
(10)
C
=
[
c
10
,
…
,
c
1
M
-
1
,
c
20
,
…
,
c
2
M
-
1
,
…
,
c
N
0
,
…
,
c
N
M
-
1
]
T
,
(11)
Θ
(
t
)
=
[
θ
10
(
t
)
,
…
,
θ
1
M
-
1
(
t
)
,
θ
20
(
t
)
,
…
,
θ
2
M
-
1
(
t
)
,
…
,
θ
N
0
(
t
)
,
…
,
θ
N
M
-
1
(
t
)
]
T
.
In (10) and (11), c
n
m, n
=
1,2
,
…
,
N, m
=
0,1
,
…
,
M
-
1, are the expansion coefficients of the function f
(
t
) in the nth subinterval [
(
(
n
-
1
)
/
N
)
t
f
,
(
n
/
N
)
t
f
) and θ
n
m
(
t
), n
=
1,2
,
…
,
N, m
=
0,1
,
…
,
M
-
1, are defined as in (4). With the aid of (7), the coefficients c
n
m can be obtained as
(12)
c
n
m
=
f
(
t
f
2
N
(
t
m
+
2
n
-
1
)
)
.


### 2.3. The Operational Matrix of Integration

In this section, the operational matrix of integration is derived. The integration of the vectorΘ
(
t
) defined in (11) can be approximated as
(13)
∫
0
t
‍
Θ
(
s
)
d
s
≃
P
Θ
(
t
)
,
where P is the M
N
×
M
N operational matrix of integration for Chebyshev cardinal functions. The matrix P can be obtained by the following process. Let
(14)
∫
0
t
‍
Θ
(
s
)
d
s
=
[
∫
0
t
‍
θ
10
(
s
)
d
s
,
…
,
∫
0
t
‍
θ
1
M
-
1
(
s
)
d
s
,
…
,
m
∫
0
t
‍
θ
N
0
(
s
)
d
s
,
…
,
∫
0
t
‍
θ
N
M
-
1
(
s
)
d
s
]
T
.Using (9), any function ∫
0
t
‍
θ
l
k
(
s
)
d
s, l
=
1,2
,
…
,
N, k
=
0,1
,
…
,
M
-
1, can be approximated as
(15)
∫
0
t
‍
θ
l
k
(
s
)
d
s
≃
∑
n
=
1
N
‍
∑
m
=
0
M
-
1
‍
c
n
m
θ
n
m
(
t
)
.
From (12) we can get
(16)
c
n
m
=
∫
0
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
‍
θ
l
k
(
s
)
d
s
.
With the aid of (4), we consider the following cases.Case 1.
Ifn
<
l, then (
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
<
(
(
l
-
1
)
/
N
)
t
f and we obtain c
n
m
=
0.Case 2.
Ifn
=
l, then (
(
l
-
1
)
/
N
)
t
f
<
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
<
(
l
/
N
)
t
f and we obtain
(17)
c
n
m
=
∫
(
(
l
-
1
)
/
N
)
t
f
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
‍
C
k
(
2
N
t
f
s
-
2
l
+
1
)
d
s
=
t
f
2
N
∫
-
1
t
m
‍
C
k
(
v
)
d
v
.Case 3.
Ifn
>
l, then (
l
/
N
)
t
f
<
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
) and we obtain
(18)
c
n
m
=
∫
(
(
l
-
1
)
/
N
)
t
f
(
l
/
N
)
t
f
‍
C
k
(
2
N
t
f
s
-
2
l
+
1
)
d
s
=
t
f
2
N
∫
-
1
1
‍
C
k
(
v
)
d
v
.
Comparing (13), (14), and (15), we obtain
(19)
P
=
(
E
H
H
⋯
H
0
E
H
⋯
H
0
0
E
⋯
H
⋮
⋮
⋮
⋮
0
0
0
⋯
E
)
,
where E and H are M
×
M matrices that can be obtained as follows.
Let(20)
E
=
(
e
i
j
)
,
H
=
(
h
i
j
)
;
then, for i
,
j
=
0,1
,
…
,
M
-
1, we have
(21)
e
i
j
=
t
f
2
N
∫
-
1
t
j
‍
C
i
(
v
)
d
v
,
h
i
j
=
t
f
2
N
∫
-
1
1
‍
C
i
(
v
)
d
v
,
where t
j, j
=
0,1
,
…
,
M
-
1, are the zeros of the first kind Chebyshev polynomial of order M. It is noted that E is the operational matrix of integration for Chebyshev cardinal functions over interval [
(
(
n
-
1
)
/
N
)
t
f
,
(
n
/
N
)
t
f
).Remark 1.
To calculate the entriese
i
j and h
i
j, i
,
j
=
0,1
,
…
,
M
-
1, we have
(22)
C
m
(
v
)
=
T
M
(
v
)
T
M
,
v
(
t
m
)
(
v
-
t
m
)
=
β
T
M
,
v
(
t
m
)
×
∏
k
=
0
,
k
≠
m
M
-
1
‍
(
v
-
t
k
)
,
where β
=
2
M
-
1 is the coefficient of t
M in the Chebyshev polynomial function T
M
(
v
). Using (22) we get
(23)
e
i
j
=
t
f
β
2
N
T
M
,
t
(
t
i
)
∫
-
1
t
j
‍
∏
k
=
0
,
k
≠
i
M
-
1
‍
(
v
-
t
k
)
d
v
,
h
i
j
=
t
f
β
2
N
T
M
,
t
(
t
i
)
∫
-
1
1
‍
∏
k
=
0
,
k
≠
i
M
-
1
‍
(
v
-
t
k
)
d
v
,
for i
,
j
=
0,1
,
…
,
M
-
1.


### 2.4. The Operational Matrix of Product

The following property of the product of two hybrid function vectors will also be used. Let(24)
Θ
(
t
)
Θ
T
(
t
)
F
≃
F
~
Θ
(
t
)
,
where
(25)
F
=
[
f
10
,
…
,
f
1
M
-
1
,
f
20
,
…
,
f
2
M
-
1
,
…
,
f
N
0
,
…
,
f
N
M
-
1
]
T
,
and F
~ is an M
N
×
M
N product operational matrix. To find F
~, we apply the following procedure. First, by using (11) and (25) we obtain
(26)
Θ
(
t
)
Θ
T
(
t
)
F
=
(
θ
10
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
θ
11
(
t
)
∑
n
=
1
∞
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
⋮
θ
1
M
-
1
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
⋮
θ
N
0
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
θ
N
1
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
⋮
θ
N
M
-
1
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
)
.
Using (9), any function θ
i
j
(
t
)
θ
l
k
(
t
),  i
,
l
=
1,2
,
…
,
N, j
,
k
=
0,1
,
…
,
M
-
1, can be approximated as
(27)
θ
i
j
(
t
)
θ
l
k
(
t
)
≃
∑
n
=
1
N
‍
∑
m
=
0
M
-
1
‍
c
n
m
θ
n
m
(
t
)
,
where
(28)
c
n
m
=
θ
i
j
(
t
f
2
N
(
t
m
+
2
n
-
1
)
)
θ
l
k
(
t
f
2
N
(
t
m
+
2
n
-
1
)
)
=
δ
i
n
δ
j
m
δ
l
n
δ
k
m
.
So, from (26) and (28), we have
(29)
Θ
(
t
)
Θ
T
(
t
)
F
≃
[
f
10
θ
10
(
t
)
,
…
,
f
1
M
-
1
θ
1
M
-
1
(
t
)
,
…
,
m
f
N
0
θ
N
0
(
t
)
,
…
,
f
N
M
-
1
θ
N
M
-
1
(
t
)
]
T
;
therefore, we find the M
N
×
M
N matrix F
~ as
(30)
F
~
=
diag
⁡
[
f
10
,
…
,
f
1
M
-
1
,
f
20
,
…
,
hhh
f
2
M
-
1
,
…
,
f
N
0
,
…
,
f
N
M
-
1
]
.Lemma 2.
The functionsC
j
(
t
), j
=
0,1
,
…
,
M
-
1, are orthogonal with respect to w
(
t
)
=
1
/
1
-
t
2 on [
-
1,1
] and satisfy the orthogonality condition
(31)
〈
C
i
(
t
)
,
C
j
(
t
)
〉
w
=
∫
-
1
1
‍
C
i
(
t
)
C
j
(
t
)
1
-
t
2
d
t
=
{
π
M
,
i
f
j
=
i
,
0
,
i
f
j
≠
i
.
The proof of this lemma is presented in [39].Remark 3.
Sinceθ
n
,
m
(
t
) consists of block-pulse functions and Chebyshev cardinal functions, which are both complete and orthogonal, the set of hybrid of block-pulse functions and Chebyshev cardinal functions is a complete orthogonal set in the Hilbert space L
2
[
0
,
t
f
).Lemma 4.
Let MN vectorsC and C
p be hybrid functions coefficients of u
(
t
) and u
p
(
t
), respectively. If
(32)
C
=
[
c
10
,
…
,
c
1
M
-
1
,
c
20
,
…
,
c
2
M
-
1
,
…
,
c
N
0
,
…
,
c
N
M
-
1
]
T
,
then
(33)
C
p
≃
[
c
10
p
,
…
,
c
1
M
-
1
p
,
c
20
p
,
…
,
c
2
M
-
1
p
,
…
,
c
N
0
p
,
…
,
c
N
M
-
1
p
]
T
,
where p
≥
1 is a positive integer.Proof.
Whenp
=
1, (33) follows at once from u
p
(
t
)
=
u
(
t
). Suppose that (33) holds for p; we will deduce it for p
+
1. Since u
p
+
1
(
t
)
=
u
(
t
)
u
p
(
t
), from (24) and (30) we have:
(34)
u
p
+
1
(
t
)
≃
(
C
T
Θ
(
t
)
)
(
C
p
T
Θ
(
t
)
)
mi
=
C
T
Θ
(
t
)
Θ
T
(
t
)
C
p
≃
C
T
C
~
p
Θ
(
t
)
,
where C
~
p can be calculated in a similar way to matrix F
~ in (24). Now, using (33) we obtain
(35)
C
p
+
1
=
C
T
C
~
p
=
[
c
10
p
+
1
,
…
,
c
1
M
-
1
p
+
1
,
c
20
p
+
1
,
…
,
mmmii
c
2
M
-
1
p
+
1
,
…
,
c
N
0
p
+
1
,
…
,
c
N
M
-
1
p
+
1
]
T
.
Therefore, (33) holds for p
+
1, and the lemma is established.


## 2.1. Hybrid Functions of Block-Pulse and Chebyshev Cardinal Functions

Hybrid functionsθ
n
,
m
(
t
), n
=
1,2
,
…
N, m
=
0,1
,
…
,
M
-
1, are defined on the interval [
0
,
t
f
) as
(4)
θ
n
,
m
(
t
)
=
{
C
m
(
2
N
t
f
t
-
2
n
+
1
)
,
t
∈
[
(
n
-
1
N
)
t
f
,
n
N
t
f
)
,
0
,
otherwise
,
where n and m are the order of block-pulse functions and Chebyshev cardinal functions, respectively. Here, C
m
(
t
) are defined as [36–38]
(5)
C
m
(
t
)
=
T
M
(
t
)
T
M
,
t
(
t
m
)
(
t
-
t
m
)
,
m
=
0,1
,
…
,
M
-
1
,
where T
M
(
t
) is the first kind Chebyshev polynomial of order M in [
-
1,1
] defined by
(6)
T
M
(
t
)
=
cos
⁡
(
M
arccos
(
t
)
)
,
subscript t denotes t-differentiation, and t
m, m
=
0,1
,
…
,
M
-
1, are the zeros of T
M
(
t
) defined by cos
⁡
(
(
2
m
+
1
)
π
)
/
2
M
), m
=
0,1
,
…
,
M
-
1, with the Kronecker property
(7)
C
m
(
t
i
)
=
δ
m
i
=
{
1
,
if
i
=
m
,
0
,
if
i
≠
m
,
where δ
m
i is the Kronecker delta function.


## 2.2. Function Approximation

A functionf
(
t
), defined over the interval [
0
,
t
f
), may be expanded as
(8)
f
(
t
)
=
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
c
n
m
θ
n
m
(
t
)
.
If the infinite series in (8) is truncated, then (8) can be expressed as
(9)
f
(
t
)
≃
∑
n
=
1
N
∑
m
=
0
M
-
1
‍
c
n
m
θ
n
m
(
t
)
=
C
T
Θ
(
t
)
,
where
(10)
C
=
[
c
10
,
…
,
c
1
M
-
1
,
c
20
,
…
,
c
2
M
-
1
,
…
,
c
N
0
,
…
,
c
N
M
-
1
]
T
,
(11)
Θ
(
t
)
=
[
θ
10
(
t
)
,
…
,
θ
1
M
-
1
(
t
)
,
θ
20
(
t
)
,
…
,
θ
2
M
-
1
(
t
)
,
…
,
θ
N
0
(
t
)
,
…
,
θ
N
M
-
1
(
t
)
]
T
.
In (10) and (11), c
n
m, n
=
1,2
,
…
,
N, m
=
0,1
,
…
,
M
-
1, are the expansion coefficients of the function f
(
t
) in the nth subinterval [
(
(
n
-
1
)
/
N
)
t
f
,
(
n
/
N
)
t
f
) and θ
n
m
(
t
), n
=
1,2
,
…
,
N, m
=
0,1
,
…
,
M
-
1, are defined as in (4). With the aid of (7), the coefficients c
n
m can be obtained as
(12)
c
n
m
=
f
(
t
f
2
N
(
t
m
+
2
n
-
1
)
)
.


## 2.3. The Operational Matrix of Integration

In this section, the operational matrix of integration is derived. The integration of the vectorΘ
(
t
) defined in (11) can be approximated as
(13)
∫
0
t
‍
Θ
(
s
)
d
s
≃
P
Θ
(
t
)
,
where P is the M
N
×
M
N operational matrix of integration for Chebyshev cardinal functions. The matrix P can be obtained by the following process. Let
(14)
∫
0
t
‍
Θ
(
s
)
d
s
=
[
∫
0
t
‍
θ
10
(
s
)
d
s
,
…
,
∫
0
t
‍
θ
1
M
-
1
(
s
)
d
s
,
…
,
m
∫
0
t
‍
θ
N
0
(
s
)
d
s
,
…
,
∫
0
t
‍
θ
N
M
-
1
(
s
)
d
s
]
T
.Using (9), any function ∫
0
t
‍
θ
l
k
(
s
)
d
s, l
=
1,2
,
…
,
N, k
=
0,1
,
…
,
M
-
1, can be approximated as
(15)
∫
0
t
‍
θ
l
k
(
s
)
d
s
≃
∑
n
=
1
N
‍
∑
m
=
0
M
-
1
‍
c
n
m
θ
n
m
(
t
)
.
From (12) we can get
(16)
c
n
m
=
∫
0
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
‍
θ
l
k
(
s
)
d
s
.
With the aid of (4), we consider the following cases.Case 1.
Ifn
<
l, then (
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
<
(
(
l
-
1
)
/
N
)
t
f and we obtain c
n
m
=
0.Case 2.
Ifn
=
l, then (
(
l
-
1
)
/
N
)
t
f
<
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
<
(
l
/
N
)
t
f and we obtain
(17)
c
n
m
=
∫
(
(
l
-
1
)
/
N
)
t
f
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
)
‍
C
k
(
2
N
t
f
s
-
2
l
+
1
)
d
s
=
t
f
2
N
∫
-
1
t
m
‍
C
k
(
v
)
d
v
.Case 3.
Ifn
>
l, then (
l
/
N
)
t
f
<
(
t
f
/
2
N
)
(
t
m
+
2
n
-
1
) and we obtain
(18)
c
n
m
=
∫
(
(
l
-
1
)
/
N
)
t
f
(
l
/
N
)
t
f
‍
C
k
(
2
N
t
f
s
-
2
l
+
1
)
d
s
=
t
f
2
N
∫
-
1
1
‍
C
k
(
v
)
d
v
.
Comparing (13), (14), and (15), we obtain
(19)
P
=
(
E
H
H
⋯
H
0
E
H
⋯
H
0
0
E
⋯
H
⋮
⋮
⋮
⋮
0
0
0
⋯
E
)
,
where E and H are M
×
M matrices that can be obtained as follows.
Let(20)
E
=
(
e
i
j
)
,
H
=
(
h
i
j
)
;
then, for i
,
j
=
0,1
,
…
,
M
-
1, we have
(21)
e
i
j
=
t
f
2
N
∫
-
1
t
j
‍
C
i
(
v
)
d
v
,
h
i
j
=
t
f
2
N
∫
-
1
1
‍
C
i
(
v
)
d
v
,
where t
j, j
=
0,1
,
…
,
M
-
1, are the zeros of the first kind Chebyshev polynomial of order M. It is noted that E is the operational matrix of integration for Chebyshev cardinal functions over interval [
(
(
n
-
1
)
/
N
)
t
f
,
(
n
/
N
)
t
f
).Remark 1.
To calculate the entriese
i
j and h
i
j, i
,
j
=
0,1
,
…
,
M
-
1, we have
(22)
C
m
(
v
)
=
T
M
(
v
)
T
M
,
v
(
t
m
)
(
v
-
t
m
)
=
β
T
M
,
v
(
t
m
)
×
∏
k
=
0
,
k
≠
m
M
-
1
‍
(
v
-
t
k
)
,
where β
=
2
M
-
1 is the coefficient of t
M in the Chebyshev polynomial function T
M
(
v
). Using (22) we get
(23)
e
i
j
=
t
f
β
2
N
T
M
,
t
(
t
i
)
∫
-
1
t
j
‍
∏
k
=
0
,
k
≠
i
M
-
1
‍
(
v
-
t
k
)
d
v
,
h
i
j
=
t
f
β
2
N
T
M
,
t
(
t
i
)
∫
-
1
1
‍
∏
k
=
0
,
k
≠
i
M
-
1
‍
(
v
-
t
k
)
d
v
,
for i
,
j
=
0,1
,
…
,
M
-
1.


## 2.4. The Operational Matrix of Product

The following property of the product of two hybrid function vectors will also be used. Let(24)
Θ
(
t
)
Θ
T
(
t
)
F
≃
F
~
Θ
(
t
)
,
where
(25)
F
=
[
f
10
,
…
,
f
1
M
-
1
,
f
20
,
…
,
f
2
M
-
1
,
…
,
f
N
0
,
…
,
f
N
M
-
1
]
T
,
and F
~ is an M
N
×
M
N product operational matrix. To find F
~, we apply the following procedure. First, by using (11) and (25) we obtain
(26)
Θ
(
t
)
Θ
T
(
t
)
F
=
(
θ
10
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
θ
11
(
t
)
∑
n
=
1
∞
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
⋮
θ
1
M
-
1
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
⋮
θ
N
0
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
θ
N
1
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
⋮
θ
N
M
-
1
(
t
)
∑
n
=
1
∞
‍
∑
m
=
0
∞
‍
f
n
m
θ
n
m
(
t
)
)
.
Using (9), any function θ
i
j
(
t
)
θ
l
k
(
t
),  i
,
l
=
1,2
,
…
,
N, j
,
k
=
0,1
,
…
,
M
-
1, can be approximated as
(27)
θ
i
j
(
t
)
θ
l
k
(
t
)
≃
∑
n
=
1
N
‍
∑
m
=
0
M
-
1
‍
c
n
m
θ
n
m
(
t
)
,
where
(28)
c
n
m
=
θ
i
j
(
t
f
2
N
(
t
m
+
2
n
-
1
)
)
θ
l
k
(
t
f
2
N
(
t
m
+
2
n
-
1
)
)
=
δ
i
n
δ
j
m
δ
l
n
δ
k
m
.
So, from (26) and (28), we have
(29)
Θ
(
t
)
Θ
T
(
t
)
F
≃
[
f
10
θ
10
(
t
)
,
…
,
f
1
M
-
1
θ
1
M
-
1
(
t
)
,
…
,
m
f
N
0
θ
N
0
(
t
)
,
…
,
f
N
M
-
1
θ
N
M
-
1
(
t
)
]
T
;
therefore, we find the M
N
×
M
N matrix F
~ as
(30)
F
~
=
diag
⁡
[
f
10
,
…
,
f
1
M
-
1
,
f
20
,
…
,
hhh
f
2
M
-
1
,
…
,
f
N
0
,
…
,
f
N
M
-
1
]
.Lemma 2.
The functionsC
j
(
t
), j
=
0,1
,
…
,
M
-
1, are orthogonal with respect to w
(
t
)
=
1
/
1
-
t
2 on [
-
1,1
] and satisfy the orthogonality condition
(31)
〈
C
i
(
t
)
,
C
j
(
t
)
〉
w
=
∫
-
1
1
‍
C
i
(
t
)
C
j
(
t
)
1
-
t
2
d
t
=
{
π
M
,
i
f
j
=
i
,
0
,
i
f
j
≠
i
.
The proof of this lemma is presented in [39].Remark 3.
Sinceθ
n
,
m
(
t
) consists of block-pulse functions and Chebyshev cardinal functions, which are both complete and orthogonal, the set of hybrid of block-pulse functions and Chebyshev cardinal functions is a complete orthogonal set in the Hilbert space L
2
[
0
,
t
f
).Lemma 4.
Let MN vectorsC and C
p be hybrid functions coefficients of u
(
t
) and u
p
(
t
), respectively. If
(32)
C
=
[
c
10
,
…
,
c
1
M
-
1
,
c
20
,
…
,
c
2
M
-
1
,
…
,
c
N
0
,
…
,
c
N
M
-
1
]
T
,
then
(33)
C
p
≃
[
c
10
p
,
…
,
c
1
M
-
1
p
,
c
20
p
,
…
,
c
2
M
-
1
p
,
…
,
c
N
0
p
,
…
,
c
N
M
-
1
p
]
T
,
where p
≥
1 is a positive integer.Proof.
Whenp
=
1, (33) follows at once from u
p
(
t
)
=
u
(
t
). Suppose that (33) holds for p; we will deduce it for p
+
1. Since u
p
+
1
(
t
)
=
u
(
t
)
u
p
(
t
), from (24) and (30) we have:
(34)
u
p
+
1
(
t
)
≃
(
C
T
Θ
(
t
)
)
(
C
p
T
Θ
(
t
)
)
mi
=
C
T
Θ
(
t
)
Θ
T
(
t
)
C
p
≃
C
T
C
~
p
Θ
(
t
)
,
where C
~
p can be calculated in a similar way to matrix F
~ in (24). Now, using (33) we obtain
(35)
C
p
+
1
=
C
T
C
~
p
=
[
c
10
p
+
1
,
…
,
c
1
M
-
1
p
+
1
,
c
20
p
+
1
,
…
,
mmmii
c
2
M
-
1
p
+
1
,
…
,
c
N
0
p
+
1
,
…
,
c
N
M
-
1
p
+
1
]
T
.
Therefore, (33) holds for p
+
1, and the lemma is established.


## 3. Hybrid Functions Method to Solve Duffing-Harmonic Oscillator

In this section, by using the results obtained in the previous section about hybrid functions, an effective and accurate method for solving Duffing-harmonic oscillator (1) is presented.Consider the following nonlinear Duffing-harmonic oscillator:(36)
u
′′
+
u
3
1
+
u
2
=
0
,
t
∈
[
0
,
t
f
)
,
with the initial conditions
(37)
u
(
0
)
=
A
,
u
′
(
0
)
=
0
.At first, we write (36) in the following form:
(38)
(
1
+
u
2
)
u
′′
+
u
3
=
0
.Let(39)
u
′′
(
t
)
≃
U
T
Θ
(
t
)
,
where Θ
(
t
) is defined in (11) and U is a vector with M
N unknowns as follows:
(40)
U
=
[
u
10
,
…
,
u
1
M
-
1
,
u
20
,
…
,
u
2
M
-
1
,
…
,
u
N
0
,
…
,
u
N
M
-
1
]
T
.
By expanding u
(
0
)
=
A and f
(
t
)
=
1 in terms of hybrid functions we get
(41)
u
(
0
)
=
A
=
A
T
Θ
(
t
)
,
(42)
f
(
t
)
=
1
=
E
T
Θ
(
t
)
,
where A
=
[
A
,
A
,
…
,
A
︷
M
N
]
T and E
=
[
1
,
1
,
…
,
1
︷
M
N
]
T. Integrating (39) from 0 to t and using (41), we obtain
(43)
u
′
(
t
)
=
∫
0
t
‍
U
T
Θ
(
s
)
d
s
+
u
′
(
0
)
≃
U
T
P
Θ
(
t
)
,
u
(
t
)
=
∫
0
t
‍
u
′
(
s
)
d
s
+
u
(
0
)
≃
∫
0
t
‍
U
T
P
Θ
(
s
)
d
s
+
A
T
Θ
(
t
)
≃
U
T
P
2
Θ
(
t
)
+
A
T
Θ
(
t
)
=
U
T
Θ
(
t
)
,
where U
=
(
U
T
P
2
+
A
T
)
T and P is the operational matrix of integration given in (13). Using Lemma 4 the functions u
2
(
t
) and u
3
(
t
) can be expanded as
(44)
u
2
(
t
)
≃
U
2
T
Θ
(
t
)
,
(45)
u
3
(
t
)
≃
U
3
T
Θ
(
t
)
.
Therefore, by using (39) and (42)–(45), the right side of (38) can be approximated as
(46)
(
1
+
u
2
(
t
)
)
u
′′
(
t
)
+
u
3
(
t
)
≃
(
E
T
+
U
2
T
)
Θ
(
t
)
Θ
T
(
t
)
U
︸
U
~
Θ
(
t
)
+
U
3
T
Θ
(
t
)
≃
[
(
E
T
+
U
2
T
)
U
~
+
U
3
T
]
Θ
(
t
)
=
0
,
where U
~ can be calculated in a similar way to matrix F
~ in (24). Since the above equation is satisfied for every t
∈
[
0
,
t
f
), we can get
(47)
(
E
T
+
U
2
T
)
U
~
+
U
3
T
=
0
.
This is a system of algebraic equations with M
N equations and M
N unknowns, which can be solved by Newton’s iteration method to obtain the unknown vector U.Remark 5.
The approximate period and frequency of the hybrid functions method (HFM) can be obtained as follows:(48)
T
HFM
=
α
,
ω
HFM
=
2
π
α
,
where α is the first positive root of equation u
(
t
)
-
A
=
0. Here, we use the famous Newton’s iteration method for finding a proper approximation α of nonlinear equation u
(
t
)
-
A
=
0, in the following form:
(49)
α
n
+
1
=
α
n
-
F
(
α
n
)
F
′
(
α
n
)
,
n
=
0,1
,
…
,
where F
(
t
)
=
u
(
t
)
-
A
=
(
U
T
-
A
T
)
Θ
(
t
) and α
0 is initial approximation.


## 4. Results and Discussions

In this section, we illustrate the accuracy of the hybrid functions method (HFM) by comparing the approximate solutions previously obtained with the exact angular frequencyω
ex. All the results obtained here are computed using the Intel Pentium 5, 2.2 GHz processor and using Maple 17 with 64-digit precision.The exact angular frequency,ω
ex, of the Duffing-harmonic oscillator was found by Lim and Wu in [32] as
(50)
ω
ex
=
π
2
(
∫
0
π
/
2
‍
A
cos
⁡
t
d
t
A
2
cos
⁡
2
t
+
ln
⁡
[
1
-
A
2
cos
⁡
2
t
/
(
1
+
A
2
)
]
)
-
1
.
By using alternative form (38) and applying the harmonic balance method (HBM) [13], Mickens [31] obtained the first approximate angular frequency
(51)
ω
HBM
=
(
3
A
2
4
+
3
A
2
)
1
/
2
.
Ozis and Yildirim [35] obtained the angular frequency using the energy balance method (EBM) in the following form:
(52)
ω
EBM
=
(
1
-
2
A
2
ln
⁡
(
1
+
A
2
1
+
(
A
2
/
2
)
)
)
1
/
2
.
Ganji et al. in [34] obtained the same approximation as that in (52). Using a single-term approximate solution u
(
t
)
=
A
cos
⁡
(
ω
t
) to (36) and the Ritz procedure [40], Tiwari et al. [41] obtained an approximate angular frequency as follows:
(53)
ω
Tiw
=
(
1
+
2
A
2
(
1
1
+
A
2
-
1
)
)
1
/
2
.The computed results for the HFM frequencyω
HFM with exact frequency ω
ex [32], HBM frequency ω
HBM [31], EBM frequency ω
EBM [35], and Tiwari’s frequency ω
Tiw [41] are listed in Tables 1 and 2. Table 2 shows that the maximum percentage error between ω
HFM and exact frequency ω
ex is 0.118%. Comparison of the exact frequency ω
ex obtained by (50) with ω
HBM, ω
EBM, ω
Tiw, and ω
HFM is shown in Figure 1 for 0
≤
A
≤
10, 3.0
≤
A
≤
3.25, 5.0
≤
A
≤
5.25, and 8.0
≤
A
≤
8.25.Table 1
Comparison of various approximate angular frequencies with exact angular frequency.
A
ω
HBM
ω
EBM
ω
Tiw
ω
ex
HFM
N
M
t
f
ω
HFM
0.01
0.00866
0.00866
0.00866
0.00847
10
5
750
0.00846
0.05
0.04326
0.04326
0.04327
0.04232
10
5
150
0.04236
0.1
0.08628
0.08627
0.08624
0.08439
10
5
75
0.08446
0.5
0.39736
0.39638
0.39423
0.38737
10
5
16.5
0.38757
1.0
0.65465
0.65164
0.64359
0.63678
15
5
10
0.63673
5.0
0.97435
0.97343
0.96731
0.96698
14
10
6.5
0.96699
10.0
0.99340
0.99314
0.99095
0.99092
14
10
6.5
0.99093Table 2
Error percentage comparison between previous results and HFM with variousA.
A
Error percentage
HBM
EBM
Tiw
HFM
0.01
2.242
2.242
2.242
0.118
0.05
2.223
2.219
2.212
0.095
0.1
2.239
2.225
2.197
0.083
0.5
2.579
2.326
1.771
0.052
1.0
2.807
2.334
1.070
0.008
5.0
0.763
0.667
0.034
0.001
10.0
0.250
0.224
0.003
0.001Figure 1
Comparison of the approximate frequencies with corresponding exact frequency.From Figure1 and Tables 1 and 2, it can be observed that the hybrid functions method (HFM) yields excellent approximate frequencies for both small and large amplitudes.Figures2 and 3 give a comparison between the present HFM results and the numerical results obtained by using the 4th order Runge-Kutta method (RK4) with time step Δ
t
=
0.001, where
(54)
Δ
(
τ
i
)
=
|
u
HFM
(
τ
i
)
-
u
RK
4
(
τ
i
)
|
,
i
=
0,1
,
…
,
t
f
β
Δ
t
,
τ
i
=
i
β
Δ
t, i
=
0,1
,
…
,
t
f
/
(
β
Δ
t
), and β
∈
N. It can be seen from these figures that the solutions obtained by the proposed procedure are in good agreement with the RK4 based solutions.Figure 2
Plots of displacementu versus time t. Solid line: HFM; solid circle: RK4.Figure 3
Plots of error valuesΔ
(
r
i
) with β
=
50.


## 5. Conclusion

In this paper, we presented a numerical scheme based on hybrid block-pulse functions and Chebyshev cardinal functions for solving Duffing-harmonic oscillator. This algorithm reduces the solution of Duffing-harmonic oscillator differential equation to the solution of a system of algebraic equations in matrix form. The merit of this method is that the system of equations obtained for the solution does not need to consider collocation points; this means that the system of equations is obtained directly. A comparative study between HBM [31], EBM [35], Tiwari’s method [41], and the proposed method was discussed in Section 4. The obtained results showed that the HFM is accurate, capable, and effective technique for the solution of the Duffing-harmonic oscillator. Further research can concentrate on other strongly nonlinear oscillators and more complicated cases.



---

*Source: 210754-2014-08-14.xml*