EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 17, No. 2, 2024, 1254-1264
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Two- dimensional coupled asymmetric van der Pol
oscillator

V. S. Erturk1, B. Rath2, Taqwa M. Al-Khader3, Noorhan Alshaikh4,
P. Mallick2, Jihad Asad4,∗

1 Department of Mathematics, Faculty of Arts and Sciences, Ondokuz Mayis University,
55139 Samsun, Turkey
2 Department of Physics, Maharaja Sriram Chandra Bhanja Deo University, Baripada-
757003, Odisha, India
3 Department of Applied Mathematics, Faculty of Applied Sciences, Palestine Technical
University- Kadoorie, Tulkarm P 305 , Palestine
4 Department of Physics, Faculty of Applied Sciences, Palestine Technical University- Kadoorie,
Tulkarm P 305 , Palestine

Abstract. A reliable algorithm based on an adaptation of the standard differential transform
method (DTM) is presented, which is the multi-step differential transform method (Ms-DTM)
since it may be difficult to directly apply differential transform method (DTM) to obtain the
series solutions for the present two- dimensional coupled asymmetric van der Pol oscillator. The
solutions of a two- dimensional coupled asymmetric van der Pol oscillator were obtained by Ms-
DTM. Figurative comparisons between the Ms-DTM and the classical fourth order Runge-Kutta
method (RK4) are given. The obtained results reveal that the proposed technique is a promising
tool to solve the considered van der Pol oscillator and yield same information on the phase portrait
confirming the stability of the system, effectively. It can be said that the considered approach
can be easily extended to other nonlinear van der Pol oscillator systems and therefore is widely
applicable in engineering and other sciences.

2020 Mathematics Subject Classifications: 34A45, 34C05, 65L05, 34C15, 34C25

Key Words and Phrases: Coupled van der Pol, two- dimension, phase portrait, semi analytical

1. Introduction

The study of nonlinear oscillations is a crucial area of research in various fields of
mechanical structures, physical science, and other mathematical sciences. Since most
real-world systems are described by nonlinear differential equations, mechanical systems,

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v17i2.5140

Email addresses: vserturk@omu.edu.tr (V. S. Erturk ), biswanathrath10@gmail.com (B. Rath),
t.alkhader@ptuk.edu.ps (T. M. Al-Khader), pravanjanphy@gmail.com (P. Mallick),
j.asad@ptuk.edu.ps (J. Asad)

https://www.ejpam.com 1254 © 2024 EJPAM All rights reserved.



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1255

mathematical physics, and engineering face several difficulties. Understanding the charac-
teristics and methods for resolving nonlinear equations in mechanical systems has drawn
more attention in recent years, with the aim of developing more accurate and efficient
methods for analyzing and predicting their behavior. Nonlinear oscillator problems have
recently been solved using a variety of analytical and numerical techniques, including
homotopy perturbation method [7−−10], energy balance method [2, 16, 27], frequency-
amplitude formulation [5, 12, 13], parameter expanding method [11, 25, 26], and vari-
ational iteration method [14, 15], and recently accurate technique is used to find an
approximate solution to the fractional-order Duffing-Van der Pol oscillators [18].

In chaos theory, the van der Pol equation is a typical reference. Early on in the model’s
creation, researchers in the biological and physical sciences made use of it. This formula
is considered to be an example of an oscillator that damps nonlinearly. Several authors
and scholars set out to find an analytical answer to the van der pol dilemma. The first
person to simulate electrical circuits in a vacuum was Van der Pol [6]. The equation is often
researched as a template for testing numerical techniques and for its extensive collection of
chaotic solutions [3]. Deriving the closed-form solutions of a nonlinear differential equation
lacks a general framework. As a result, there has been a growing emphasis on numerical
techniques to get solutions (some of which are provided in [19]), perturbative approaches,
approximated analytical approaches based on linearization [20], or adomian decomposition
[17].

In order to examine the oscillatory behavior, the parameter is crucial. These pa-
rameters can be computed in a variety of ways [27]. Among of these techniques, the
approximation Lie theorem was given by Baikov [5, 12, 13]. For the perturbed ODEs,
Naeem and Mahomed [11, 26] proposed the partial Lagrangian technique [14, 25]. Many
research works have been developed for the systems of approximation Hamiltonians.

The mathematical model for the system is a well- known second order ordinary differen-
tial equation with cubic nonlinearity – the Van der Pol equation. Since then, thousands of
papers have been published achieving better approximations to the solutions occurring in
such non- linear systems. The Van der Pol oscillator is a classic example of self-oscillatory
system and is now considered as very useful mathematical model that can be used in
much more complicated and modified systems. The Van der pol equation is so impor-
tant to mathematicians, physicists and engineers and is still being extensively studied
[4, 21 − −24]. Since its introduction in the 1920’s, the Van der Pol equation has been a
prototype for systems with self-excited limit cycle oscillations. The classical experimental
setup of the system is the oscillator with vacuum triode. The investigations of the forced
Van der Pol oscillator behaviour have carried out by many researchers. The equation
has been studied over wide parameter regimes, from perturbations of harmonic motion to
relaxation oscillations [1].

Let us discuss the exact solution of a Van der Pol like equation in two dimensions [8]
as:

x′′ + ω2x+ ε
(
αx2 + βy2 − 1

)
x′ = 0 (1)

y′′ + ω2y − ε
(
αx2 + βy2 − 1

)
y′ = 0 (2)



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1256

For α = β = 1, the above equations reduce to:

x′′ + ω2x+ ε
(
x2 + y2 − 1

)
x′ = 0. (3)

y′′ + ω2y − ε
(
x2 + y2 − 1

)
y′ = 0. (4)

Using x = A sinωt and y = A cosωt, we have

x′′ + ω2x = 0. (5)

y′′ + ω2y = 0. (6)

The solution of ( 5&6 ) remains the same as our assumptions.
However, for the case α ̸= β ̸= 1 one has to use a specific analysis as given section 2

below.

2. Strong Coupling Limit for ω = 1, α ̸= β, and ε = 1

In this section, we present some numerical analysis for the system given in ( 1 & 2) for
the following case

Consider α = 1.546, β = 1.551, ω = 1 with ε = 1. In this case system ( 1&2) become:

x′′ + x+
(
1.546x2 + 1.551y2 − 1

)
x′ = 0. (7)

y′′ + y −
(
1.546x2 + 1.551y2 − 1

)
y′ = 0. (8)

The above system is a coupled nonlinear one. Numerical technique has to be used to
solve it. Below we use fourth- order Runge- Kutta technique to simulate x(t), y(t), x′(t)
and y′(t) vs time as shown in Figure 1 (a- c).

Figure 1: Simulation of system (7 & 8) using forth- order Runge- Kutta method (a) x(t) vs time, (b) y(t) vs
time (c) x′(t) vs time, and (d) y′(t) vs time



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1257

The time-series of the dynamical variables (x(t) and y(t)) is shown in figure 1 (a and b
). As one can see, these solutions are periodic. It may be noted that the time evolution of
the dynamical variables (x(t) and y(t)) with the same initial conditions nearly show similar
oscillatory behaviour, there are also minute changes in amplitudes. The same conclusion
can be drawn for figure 1 (c and d) for the dynamical variables (x′(t) and y′(t)).

Furthermore, In Figure 2( a& b) the phase portrait of x′(t)vsx(t), and y′(t) vs y(t) has
been plotted. These phase portraits are closed and show nearly rectangular shape. This
implies that the suggested system a stable one.

Figure 2: Phase portrait of system (7 & 8): (a) x′(t) vs x(t), and (b) y′(t) vs y(t).

3. The Ms- DTM method

The multi-step differential transform technique, or Ms-DTM for short, is the method
used to solve ODEs numerically, and it is the sole topic of this section. Let [0, T ] be the
interval for the nonlinear initial value issue, where a finite series may be used to represent
f
(
t, x, x′, . . . , x(r)

)
= 0

x(t) =
N∑

n=0

ant
n, t ∈ [0, T ] (9)

with an initial condition x(r)(0) = ck for k = 0, 1, . . . , r − 1.

Using the nodes tm = mh, we suppose that interval [0, T ] is split into M subintervals
[tm−1, tm] with m = 1, 2, . . . ,M and step size h = T/M . Ms-DTM is a method for
computational performance for values of h. The main idea of Ms-DTM is we first apply
the differential transform method(in short, DTM) to ODE with initial value problem



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1258

f
(
t, x, x′, . . . , x(r)

)
= 0 over the interval [0, t1]. For initial condition x

(k)
1 (0) = ck, following

approximate solution denoted by

x1(t) =
K∑

n=0

a1nt
n, t ∈ [0, t1] . (10)

For m ≥ 2 and subinterval [tm−1, tm] with initial conditions x
(k)
m (tm−1) = x

(k)
m−1 (tm−1),

we apply the DTM to f
(
t, x, x′, . . . , x(r)

)
= 0, where t0 in Y (k) = 1

k!
dky(t)
dtk

]
t=t0

, which is

called the differential transformation(DT, in short) of a function y(t) is replaced by tm−1.
The process is repeated and generates a sequence of approximate solutions xm(t),m =
1, 2, . . . ,M with N = K.M for the solution x(t) and denoted by

xm(t) =
K∑

n=0

amn (t− tm−1)
n , t ∈ [tm, tm+1] . (11)

Finally, the Ms-DTM assumes the following solution denoted by x(t)

x(t) =



x1(t), t ∈ [0, t1]
x2(t), t ∈ [t1, t2]

·
·
·

xm(t), t ∈ [tM−1, tM ] .

(12)

4. Solution of the system via Ms-DTM

In this section, we apply the Ms-DTM to Eqs. (7)-(8) to demonstrate the effectiveness
of Ms-DTM as an approximate tool for solving the ODE. Let x = u1, x

′ = u2, y = u3 and
y′ = u4. Substituting back into the ODE (keeping in mind that x′′ = u′2 and y′′ = u′4 ) we
get:

u′2 + u2
(
1.546u21 + 1.551u23 − 1

)
+ u1 = 0. (13)

and
u′3 − u4

(
1.546u21 + 1.551u23 − 1

)
+ u3 = 0. (14)

Thus Eqs. (7)-(8) can be expressed in the form of four simultaneous first-order differential
equations in terms of u1, u2, u3 and u4, i.e.

u′1 = u2,

u′2 = −u2
(
1.546u21 + 1.551u23 − 1

)
− u1

u′3 = u4

u′4 = u4
(
1.546u21 + 1.551u23 − 1

)
− u3. (15)



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1259

Taking the DT of Eq. (15) by using the well-known differential transform formulas,
we obtain

U1(k + 1) =
1

k + 1
U2(k),

U2(k + 1) =
1

k + 1−1.546
k∑

k2=0

k2∑
k1=0

U2 (k1)U1 (k2 − k1)U1 (k − k2) − 1.551
k∑

k2=0

k2∑
k1=0

U2 (k1)U3 (k2 − k1)U3 (k − k2) + U2(k) − U1(k)


U3(k + 1) =

1

k + 1
U4(k)

U4(k + 1) =
1

k + 1
,1.546

k∑
k2=0

k2∑
k1=0

U4 (k1)U1 (k2 − k1)U1 (k − k2) + 1.551
k∑

k2=0

k2∑
k1=0

U4 (k1)U3 (k2 − k1)U3 (k − k2) − U4(k) − U3(k)



(16)

where U1(k), U2(k), U3(k) and U4(k) are the differential transforms of u1, u2, u3 and u4,
respectively. The DT of the initial conditions are given by U1(0) = 1, U2(0) = 0, U3(0) = 1,
and U4(0) = 0. In view of y(t) =

∑N
k=0 Y (k) (t− t0)

n, which is called the differential
inverse transform, DTM series solution for Eqs. (7)-(8) can be obtained as,

u1(t) =
N∑

n=0

U1(n)t
n,

u2(t) =
N∑

n=0

U2(n)t
n,

u3(t) =
N∑

n=0

U3(n)t
n,

u4(t) =
N∑

n=0

U4(n)t
n.

(17)

Now, according to the Ms-DTM, taking N = K.M , the series solution for Eqs. (7)-(8)
is given by,

u1(t) =



∑K
n=0 U11(n)t

n, t ∈ [0, t1]∑K
n=0 U21(n)t

n, t ∈ [t1, t2]
·
·∑K

n=0 UM1(n)t
n, t ∈ [tM−1, tM ]

(18)



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1260

u2(t) =



∑K
n=0 U12(n)t

n, t ∈ [0, t1]∑K
n=0 U22(n)t

n, t ∈ [t1, t2]
·
·∑K

n=0 UM2(n)t
n, t ∈ [tM−1, tM ]

(19)

u3(t) =



∑K
n=0 U13(n)t

n, t ∈ [0, t1]∑K
n=0 U23(n)t

n, t ∈ [t1, t2]
...
·∑K

n=0 UM3(n)t
n, t ∈ [tM−1, tM ]

(20)

and

u4(t) =



∑k
n=0 U14(n)t

n, t ∈ [0, t1]∑K
n=0 U24(n)t

n, t ∈ [t1, t2]
...
·∑K

n=0 UM4(n)t
n, t ∈ [tM−1, tM ]

(21)

where Un1, Un2, Un3 and Un4, for n = 1, 2, . . . ,M , satisfy the following recurrence
relations,

Un1(k + 1) =
1

k + 1
Un2(k),

Un2(k + 1) =
1

k + 1−1.546
k∑

k2=0

k2∑
k1=0

Un2 (k1)Un1 (k2 − k1)Un1 (k − k2) − 1.551
k∑

k2=0

k2∑
k1=0

Un2 (k1)Un3 (k2 − k1)Un3 (k − k2) + Un2(k) − Un1(k)


Un3(k + 1) =

1

k + 1
Un4(k),

U4(k + 1) =
1

k + 11.546
k∑

k2=0

k2∑
k1=0

Un4 (k1)Un1 (k2 − k1)Un1 (k − k2) + 1.551
k∑

k2=0

k2∑
k1=0

Un4 (k1)Un3 (k2 − k1)Un3 (k − k2) − Un4(k) − Un3(k)


(22)

such that Un1(0) = U(n−1)1(0), Un2(0) = U(n−1)2(0), Un3(0) = U(n−1)3(0) and Un4(0) =
U(n−1)4(0).



J. Asad et al. / Eur. J. Pure Appl. Math, 17 (2) (2024), 1254-1264 1261

Finally, if we start with U01(0) = 1, U02(0) = 0, U03(0) = 1 and U04(0) = 0, using the
recurrence relations given in Eq. (11), then we can obtain the Ms-DTM solution given
in Eqs. (7), (8), (9) and (10). Figure 3 shows the phase-portrait as obtained from the
approximate solutions for Eqs. (7)-(8) obtained using Ms-DTM. From Figure 2 (a-b) and
Figure 3 (a-b), we can observe that the phase-portrait trajectories obtained using Ms-
DTM are in high agreement with the phase-portrait trajectories obtained using fourth-
order Runge- Kutta method.

Figure 3: Phase portrait of system (7 & 8): (a) x(t) vs x′(t), and (b) y(t) vs y′(t) using Ms- DTM approach.

Figure 4 (a-d) below show the approximate solutions for Eqs. (7)-(8) obtained using
Ms-DTM. One can see that there is a high agreement with the time-series of the dynamical
variables obtained using fourth- order Runge- Kutta method shown in figure 1 (a-d).

Figure 4: Simulation of system ( 7 & 8) using Ms-DTM (a) x(t) vs time, (b) y(t) vs time (c) x′(t) vs time,
and (d) y′(t) vs time



REFERENCES 1262

5. Conclusion

Van der Pol oscillator is one of the paradigms of nonlinear dynamics and a common
model for nonlinear phenomena in science and engineering. In this study, an algorithm
for solving a two-dimensional coupled van der Pol like oscillator was introduced via Ms-
DTM. A reliable accuracy solution was obtained via this algorithm. The considered two-
dimensional coupled asymmetric van der Pol oscillator are discussed and focused it’s quasi
analytical solution of a coupled two-dimensional oscillator in strong coupling regime ω =
ε = 1, α ̸= β. For α = β = 1, solution is a sinusoidal nature having elliptical nature
phase portrait. However, here α ̸= β ̸= 1 we get nearly rectangular type in both the
cases, rectangular nature is asymmetric. Similarly, we also plot displacement vs time
graph (x ∼ t, y ∼ t) and velocity vs time graph (x′ ∼ t, y′ ∼ t , )We believe present two-
dimensional van der Pol model will generate interest among people to look into coupling
regime. The comparison between Ms-DTM solution and the fourth order Runge-Kutta
method are discussed and given graphically. The solution via Ms-DTM is continuous on
this domain and analytical at each subdomain. In this work all calculations are made by
Mathematica Software.

Acknowledgements

Taqwa M. Al-Khader, Noorhan Alshaikh and Jihad Asad would like to thank Palestine
Technical University- Kadoorie for funding this research financially.

References

[1] Ali Al-Qahtani, Houari B. Khenous, and Shaban Aly. Synchronization of impulsive
real and complex van der pol oscillators. Applied Mathematics, 6:6, 2015.

[2] H Askari, M K Yazdi, and Z Saadatnia. Frequency analysis of nonlinear oscillators
with rational restoring force via he’s energy balance method and he’s variational
approach. Nonlinear Science Letters A, (1):425, 2010.

[3] F M Atay. Van der pol’s oscillator under delayed feedback. Journal of Sound and
Vibration, 218(2):333, 1998.

[4] M. L. Cartwright. Balthazar van der pol. J.London Mathematical Society, 35:367–376,
1960.

[5] A G Davodi, D D Gangi, R Azami, and H Babazadeh. Application of improved
amplitude-frequency formulation to nonlinear differential equation of motion equa-
tions. Modern Physics Letters B, 23:3427–3436, 2009.

[6] Serge D’Alessio. Solutions of the van der pol equation. The College Mathematics
Journal, 24(2):90, 2023.



REFERENCES 1263

[7] D D Ganji. The application of he’s homotopy perturbation method to nonlinear
equations arising in heat transfer. Physics Letters A, 355:337, 2006.

[8] J H He. A coupling method of a homotopy technique and a perturbation technique
for non-linear problems. International Journal of Non-Linear Mechanics, 35(1):37,
2000.

[9] J H He. Application of homotopy perturbation method to nonlinear wave equations.
Chaos Solitons and Fractals, 26:695, 2005.

[10] J H He. New interpretation of homotopy perturbation method. International Journal
of Modern Physics B, 20:2561, 2006.

[11] J H He. Some asymptotic methods for strongly nonlinear equations. International
Journal of Modern Physics B, 20:1141, 2006.

[12] J H He. Comment on ’he’s frequency formulation for nonlinear oscillators. European
Journal of Physics, 29:19, 2008.

[13] J H He. Max-min approach to nonlinear oscillators. International Journal of Non-
linear Science and Numerical Simulation, 9:211, 2008.

[14] J H He, G C Wu, and F Austin. The variational iteration method which should be
followed. Nonlinear Science Letters A, 1:1, 2010.

[15] N Herisanu and V Marinca. A modified variational iteration method for strongly
nonlinear problems. Nonlinear Science Letters A, 1:183, 2010.

[16] I Mehdipour, D D Ganji, and M Mozaffari. Application of the energy balance method
to nonlinear vibrating equations. Current Applied Physics, 10:104, 2010.

[17] S S. Motsa and P Sibanda. A note on the solutions of the van der pol and duff-
ing equations using a linearisation method. Mathematical Problems in Engineering,
2012:693453, 2012.

[18] Attia N., Seba D., Akgul A., and NourA. Solving duffing-van der pol oscillator
equations of fractional order by an accurate technique. Journal of Applied and Com-
putational Mechanics, 7(3):1480–1487, 2021.

[19] P V Ramana and P B K Raghu. Modified adomian decomposition method for van
der pol equations. International Journal of Non Linear Mechanics, 65:121, 2014.

[20] A S Soomro, G A Tularam, and M M S Shaikh. A comparison of numerical methods
for solving the unforced van der pol’s equation. Mathematical Theory and Modeling,
3(2):66, 2013.

[21] B. van der Pol. A theory of the amplitude of free and forced vibrations. Radio Review,
1:701–710, 754–762, 1920.



REFERENCES 1264

[22] B. van der Pol. On “relaxation-oscillations”. Phil .Mag., 2:978–992, 1926.

[23] B. van der Pol. The nonlinear theory of electric oscillations. Proc.IRE, 22:1051–1086,
1934.

[24] B. van der Pol and J. van der Mark. Frequency demultiplication. Nature, 120:363–364,
1927.

[25] L Xu. Determination of limit cycle by he’s parameter-expanding method for strongly
nonlinear oscillators. Journal of Sound and Vibration, 302:178, 2007.

[26] L Xu. Variational iteration method-reality, potential, and challenges. Jownal of
Computational and Applied Mathematics, 207:148, 2007.

[27] D Younesian, H Askari, Z Saadatnia, and M Kalami Yazdi. Frequency analysis of
strongly nonlinear generalized duffing oscillators using he’s frequency–amplitude for-
mulation and he’s energy balance method. Computers and Mathematics with Appli-
cations, 59:3222, 2010.