Simulation of Chaotic Operation Of A Damped Driven Pendulum Using

Python

JOAN JANI

Department of Engineering Physics

Polytechnic University of Tirana

Rr. S. Delvina, 1001

ALBANIA

Abstract: In this paper, we are presenting a new pedagogical method for the introduction of the study of non-

linear systems. Our approach is based on the use of open-source software which is publicly available. In response

to this motivation, we are using the Python programming language which offers a holistic approach to scientific

research. We will present the analysis of the pendulum motion under the influence of an external force. The

differential equation governing the system will be presented and solved using numerical methods. Moreover, the

phase diagram of the system will be presented for various system parameters. We will describe the transition to

a chaotic operation and the key factors of this procedure. The chaotic behaviour is verified by calculating the

maximum Lyapunov exponent. This pedagogical approach emerging is based on the physical properties of the

system and not on the numerical methods used so that the student can understand the dynamics of the system more

comprehensively.

Key-Words: Chaos, simple harmonic oscillator, damped driven pendulum, phase space, Lyapunov exponent

Received: October 19, 2021. Revised: October 23, 2022. Accepted: November 25, 2022. Published: January 25, 2023.

Introduction

Computer modeling is an essential part of engineer-

ing education, as it provides students with the op-

portunity to develop and apply their knowledge in a

practical setting. By utilizing computer models, en-

gineers can gain a better understanding of how dif-

ferent components interact within complex systems,

[1]. For this purpose, the knowledge of computer

programming should be an integral part of the edu-

cation of engineers and scientists. In the past, lan-

guages such as Fortran, C, C++, etc. have been used

for expressing the mathematical models which are de-

scribing the system, [2]. Using these languages re-

quires experience and a deep understanding of how

computers work which should not be a prerequisite

for a new scientist to get an introduction to the field

of dynamical system study. On the other hand, the

Python programming language offers a user-friendly

working environment, offers fairly good documenta-

tion, and can be used even by the new scientist, [3].

Python is an ideal language for physics simula-

tion and modeling due to its ease of use, flexibil-

ity, and availability of a wide range of libraries for

physics. Python allows for fast prototyping of simu-

lations and models and can be used to quickly develop

user interfaces for visualization and control. Addi-

tionally, Python can be used to integrate data from

different sources, such as sensor data, to develop com-

plex models and simulations. Python is also well-

suited for developing physics-based models and sim-

ulations due to its powerful graphics capabilities. Fi-

nally, Python is a great choice for physics research

and engineering, as it can be used to automate data

analysis and the development of models and simula-

tions. The Python programming language has been

generally accepted as a scientific tool in the field of

the study of dynamical systems, [4].

There are several publicly available packages for

the Python programming language that offer ready-

to-use numerical methods, [5]. The use of these

methods focuses students’ interest in understanding

the model and the main parameters and constants

that govern it. As an example of the application of

these methods, we will present the use of the function

odeint() for solving differential the system of dif-

ferential equations which describe the behavior of the

system could be found in the library scipy which is

a scientific computation library that apart from ODE

(Ordinary Differential Equations) solvers offers mod-

ules for optimization, integration, interpolation, linear

algebra, FFT (Fast Fourier transform ), signal, and im-

age processing, [6]. The use of Python for teaching

purposes in the context of courses such as Computa-

tional Physics and Dynamical Systems is constantly

gaining ground over other languages, [7].

Chaos is the appearance of a lack of order in a

system that nevertheless follows laws or rules. [8]

Chaos is present in systems with non-linear behavior,

WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION

DOI: 10.37394/232010.2023.20.1

Joan Jani

E-ISSN: 2224-3410

Volume 20, 2023

the non-linearity in a system means that the measured

values in the output of the system are not propor-

tional to the values of the input. The presence of non-

linearity in a system does not mean that the system

will behave chaotically but a form of non-linearity if

required to achieve chaotic behavior, [9]. A system

showing chaotic behavior undergoes transitions be-

tween non-chaotic and chaotic states in general, [10].

A similar approach could be found in [11] where the

transition to chaos is made by changing the angular

frequency of external force, in our approach the tran-

sition to chaos is made by changing the amplitude of

the external force applied to the system.

1 Simple oscillator

We begin our pedagogical method by using the inte-

gration method on a well-known problem such as that

of the oscillator.

Fig. 1: Simple Harmonic Oscillator with a spring

mass system

In figure 1 a simple harmonic oscillator with a

spring mass system is presented. Since the constant

of elasticity of the spring is denoted with kand the

mass of the object with m, we can express the differ-

ential equation of the motion of the system as follows:

md2x

dt2=−kx (1)

the system is governed by a second-order linear dif-

ferential equation, [12]. To apply a numerical method

of integration for equation (1) we rewrite it as a pair of

first-order differential equations, which are presented

at equation (2).

dt =−k

mxand dx

dt =v(2)

We suppose the system with parameters k=

9N/mand m= 1 kg and with initial conditions

x= 5 mand v= 0 m/s. The simulation will cal-

culate the position and velocity for the time interval

of the oscillator from t= 0 sto t= 10 s. Will be cal-

culated N= 1000 points of the trajectory. Increas-

ing the sampling (N>1000) of the trajectory would

require more computing power to achieve the solu-

tion, while reducing it would result in an increase

in the computational error of the numerical method.

For our purpose, we will use the numpy which is a

Python library used for working with arrays, and the

matplotlib which is a plotting library for the Python

programming language.

# Oscilator

from s c i p y . i n t e g r a t e import odeint

import numpy a s np

import m a t p l o t l i b . p y p l o t a s p l t

# S y s t e m P a r a m e t e r s

[ k , m] = [ 9 , 1 ]

# I n i t i a l C o n d i t i o n s

x = [ 5 , 0 ]

# S i m u l a t i o n p a r a m e t e r s

N = 1000

t = np . l i n s p a c e ( 0 , 1 0 ,N)

d e f f ( x , t ) :

d v d t = −( k /m)*x [ 0 ]

d x d t = x [ 1 ]

r e t u r n [ d xd t , d v d t ]

s o l u t i o n = o d e i n t ( f , x , t )

The above code implements the solution. The

code can be easily understood by someone who has

no programming knowledge. In the beginning, the

necessary packages are loaded into the program. The

commands for the system parameters and the ini-

tial conditions are given following with the function

f(x,t) with the differential equations of the system.

The function odeint solves the system and the results

are saved to variable solution

p l t . p l o t ( t , s o l u t i o n [ : , 0 ] )

p l t . x l a b e l ( ’ Time t [ s ] ’ )

p l t . y l a b e l ( ’ x [m] ’ )

p l t . show ( )

To plot the solution which we are getting from

odeint, the above code is used. The presented code

uses the library matplotlib which has a very similar

syntax to MATLAB. The result is presented in figure (2)

from were the sinusoidal correlation of the position of

the oscillator in relation to time can be seen.

The phase space of the system is created with some

minor changes to the above code. (fig. 3). The phase

diagram gives us an overview of the evolution of the

system. Τhe position of the oscillator is shown on the

horizontal axis, while its speed is on the vertical axis.

Furthermore, from the Fig. 3 we observe that the pe-

riod of the system is T= 2.09 sas we expect from

the relation T= 2π/√k/m

2 Driven damped pendulum

Considering a simple pendulum in a constant gravita-

tional field as it is presented in fig. 4. The equation of

motion is given by the differential equation 3 where

its current angle is denoted by θand angular velocity

WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION

DOI: 10.37394/232010.2023.20.1

Joan Jani

E-ISSN: 2224-3410

Volume 20, 2023

0 2 4 6 8 10

Time t[s]

x[m]

Fig. 2: Position of the oscillator vs time

4 2 0 2 4

x[m]

v[m/s]

Phase Space

Fig. 3: Phase space of the harmonic oscillator

by ˙

θ.

−ml¨

θ=mg sin(θ)(3)

The reason that the small-angle approximation is

made is that otherwise the differential equation which

describes the system could not be solved using analyt-

ical methods. The solution of the equation could eas-

ily get using the code below. Regardless of the initial

conditions, the system will go to a steady state. [11]

Equation (4) shows the differential equation of a

damped-driven pendulum. The motion of the pendu-

lum is damped by a force proportional to its velocity

which is expressed in the equation with the term −β˙

θ.

Moreover, equation (4) represents a pendulum with a

time-varying external force with amplitude Aand fre-

quency ωexpressed by the term Asin(ωt).

θ=−g

Lsin θ−β˙

θ+Asin(ωt)(4)

For solving the differential equation (4) we con-

sider a system with parameters, g= 10 m/s2,L=

10 m,β= 0.5Ns/m,A= 1 Nand ω= 2,3rad/s.

The simulation starts from initial position θ= 0.5rad

Fig. 4: A simple pendulum

and ˙

θ= 0 rad/s. The simulation will be performed

for 0< t < 200 susing N= 20000 time samples.

# D r i v e n dumped p endul um

from s c i p y . i n t e g r a t e import odeint

import numpy a s np

import m a t p l o t l i b . p y p l o t a s p l t

# S y s t e m P a r a m e t e r s

[ g , L , b ]= [ 1 0 , 10 , 0 . 5 ]

[w , A] = [ 2 / 3 , 1 ]

# I n i t i a l C o n d i t i o n s

θ= [ 0 . 5 , 0 ]

# S i m u l a t i o n p a r a m e t e r s

N = 20000

t = np . l i n s p a c e ( 0 , 2 0 0 ,N)

d e f f ( θ, t ) :

d v d t = −( g / L)*np . s i n ( θ[0])...

. . . − b*θ[1]+A*np . s i n (w*t )

dθd t = θ[1]

r e t u r n [ dθd t , d v d t ]

s o l u t i o n = o d e i n t ( f , θ, t )

We run the code and the results are saved in the

variable “solution”. The time needed for the solu-

tion is relatively short although the number of cal-

culation points is significantly larger. Then we pro-

ceed with the creation of the corresponding graphical

representations of the angular displacement with time

and phase space.

The graph is shown in Fig. 5 where we see that

after the time of 100 seconds the system performs

harmonic oscillation of constant amplitude and fre-

quency. At this point we let the students experiment

with the various initial conditions of the system. By

choosing different initial conditions we see that the

system exhibits the same behavior and finally per-

forms harmonic oscillations.

To construct the phase diagram we will use the

states of the system after harmonic oscillation is

reached. We should also ensure that the value of the

WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION

DOI: 10.37394/232010.2023.20.1

Joan Jani

E-ISSN: 2224-3410

Volume 20, 2023

0 25 50 75 100 125 150 175 200

Time t[s]

[rad]

=1.0, = 0.5, A=1, =0.667

Fig. 5: Position of the oscillator vs time

angle should not exceed the range −π≤θ≤π.

210123

Position [rad]

Speed [rad/s]

=1.0, = 0.5, A=1, =0.667

Fig. 6: Phase space of the system with A= 1 N

In order to bring the system described by the equa-

tion (4) in chaotic operation we use the method of

period doubling by changing the amplitude of the

external force. We have changed the angular fre-

quency of the external force to ω= 2/3 rad/sand

we change the amplitude of the external force to A=

[1.35,1.45,1.47,and 1.5]

The transition to chaos is shown in figure 7 where

the phase space of the system is presented when the

amplitude of the driving force is A= 1.35. The tra-

jectory of the system in the phase space diagram is a

closed loop as expected. Increasing the amplitude to

A= 1.45 the trajectory in the phase space is chang-

ing, and a second path is appearing denoting the pe-

riod doubling of the system as it is shown in fig. 8.

By increasing the amplitude to A= 1.47 the phase

space diagram of the system changes again as it is pre-

sented in figure 9. The system doubles its period, but

its behavior is still predictable.

By choosing the amplitude of the external force

A= 1.5the operation of the system becomes chaotic

3 2 1 0 1 2 3

Position [rad]

Speed [rad/s]

=1.0, = 0.5, A=1.35, =0.667

Fig. 7: Phase space of the system with A= 1.35 N

3 2 1 0 1 2 3

Position [rad]

Speed [rad/s]

=1.0, = 0.5, A=1.45, =0.667

Fig. 8: Phase space of the system with A= 1.45 N

as it shown in Fig. 10. It has to be mentioned that

the motion of the system is not random. There is

a periodicity but the period is large enough. This

makes a small observation of the operation look to-

tally stochastic. In this state, the system is called a

strange attractor. In this situation, a small change in

the initial conditions of the system will cause a big

difference in its motion, regardless that the shape of

the phase space will not change.

The chaotic operation could be characterized by

the rate of separation of infinitesimally close trajecto-

ries. this is done by the calculation of the Lyapunov

exponent described below. The presence of chaos is

verified also in electronic circuits [13].

3 Lyapunov exponent calculation

The chaotic behavior of a systems is verified by the

calculation of the Lyapunov exponent, [14]. The cal-

culation of this exponent is done by equation (5):

δt∼

=δ0eλt λ=1

∑

j=1

log δt,j

δ0,j

(5)

WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION

DOI: 10.37394/232010.2023.20.1

Joan Jani

E-ISSN: 2224-3410

Volume 20, 2023

3 2 1 0 1 2 3

Position [rad]

Speed [rad/s]

=1.0, = 0.5, A=1.47, =0.667

Fig. 9: Phase space of the system with A= 1.47 N

3 2 1 0 1 2 3

Position [rad]

Speed [rad/s]

=1.0, = 0.5, A=1.5, =0.667

Fig. 10: Phase space of the system with A= 1.5

were λis the Lyapunov exponent δtis the separation

of trajectories after ttime, δ0is the initial separation.

For discrete time systems, like in our case, N is the

time samples needed for the calculation. In order to

make the calculation, we simulate two trajectories of

the system whose initial conditions differ very little.

The initial conditions θ1oand θ2owhere θ2o=θ1o+ϵ,

where ϵis value very close to the machine epsilon, for

our case ϵ= 2.220e−16. so the initial separation of

the trajectories in the phase space is:

δ0=|θ2o−θ1o|

After the calculation of each trajectory for the cor-

responding initial condition we calculate their differ-

ence through the relation.

δ0= ∆θ=|θ2−θ1|

the difference is growing exponentially with time as

it is shown in fig. 11. The graph is assumed to be

linear when we create a log-linear plot. From our

calculations, the Lyapunov exponent of this system

is λ= 0.06132158503950399. The positive value

verifies the chaotic operation as was expected.

0 25 50 75 100 125 150 175 200

t [s]

100

102

104

106

108

[radians]

Fig. 11: Exponential increase of the difference be-

tween two very close orbits with respect to time. The

red line represents the linear regression.

4 Discussion

This approach is applied in the context of the Compu-

tational Physics course taught to second-year students

in the Applied Physics department of the Polytechnic

University of Tirana in Albania. After the explanation

of the basic principles governing the phenomenon, the

students, with the help of the lecturer, construct the

mathematical model that describes it. This model is

then run on the computer and the graphs describing

its evolution are created.

Our didactic approach presents clearly and com-

prehensively the transition of a nonlinear system from

order to chaos. The demonstration could be presented

in the lecture or a laboratory environment giving stu-

dents a hands-on experience in the field of dynamic

systems and chaos.

In the future, we will try to make a create the code

that implements the Poincare section of the system.

5 Conclusion

In this paper, we demonstrate the use of Python pro-

gramming language in simulating a simple linear sys-

tem. We showed that this methodology can be used in

simulations of a non-linear system. This demonstra-

tion can be used as a didactic tool for introduction to

dynamical system analysis of its elegance and ease of

use.

References:

[1] E. A. Gago, C. L. Brstilo, and N. de Brito,

“Computational simulation: Multidisciplinary

teaching of dynamic models from the linear

algebra perspective,” WSEAS Transactions On

Advances In Engineering Education, 2022.

[2] P. W. Gaffney, “A performance evaluation of

some fortran subroutines for the solution of

WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION

DOI: 10.37394/232010.2023.20.1

Joan Jani

E-ISSN: 2224-3410

Volume 20, 2023

stiff oscillatory ordinary differential equations,”

ACM Transactions on Mathematical Software

(TOMS), vol. 10, no. 1, pp. 58–72, 1984.

[3] P. Borcherds, “Python: a language for computa-

tional physics,” Computer Physics Communica-

tions, vol. 177, no. 1, pp. 199–201, 2007. Pro-

ceedings of the Conference on Computational

Physics 2006.

[4] B. W. l. Margolis, “Simupy: A python frame-

work for modeling and simulating dynamical

systems,” Journal of Open Source Software,

vol. 2, no. 17, p. 396, 2017.

[5] D. R. Gwynllyw, K. L. Henderson, E. G. Guil-

lot, et al., “Using python in the teaching of nu-

merical analysis.,” MSOR Connections, vol. 18,

no. 2, 2020.

[6] J. Nunez-Iglesias, S. Van Der Walt, and H. Dash-

now, Elegant SciPy: The art of scientific python.

” O’Reilly Media, Inc.”, 2017.

[7] J. Kusalaas, “Numerical methods in engineering

with python 3,” 2013.

[8] A. BenSaïda, “A practical test for noisy chaotic

dynamics,” SoftwareX, vol. 3-4, pp. 1–5, 2015.

[9] S. Strogatz, Nonlinear Dynamics and Chaos:

With Applications to Physics, Biology, Chem-

istry, and Engineering. CRC Press, 2018.

[10] H. Nagashima, Y. Baba, and M. Nakahara, In-

troduction to chaos: physics and mathematics

of chaotic phenomena. CRC Press, 2019.

[11] E. Ayars, “Computational physics with python,”

California State University, 2013.

[12] K. N. Anagnostopoulos, Computational

Physics, Vol I: A Practical Introduction to Com-

putational Physics and Scientific Computing.

Konstantinos Anagnostopoulos, 2014.

[13] M. Hanias, I. Giannis, and G. Tombras, “Chaotic

operation by a single transistor circuit in the

reverse active region,” Chaos: An Interdisci-

plinary Journal of Nonlinear Science, vol. 20,

no. 1, p. 013105, 2010.

[14] J. Jani and P. Malkaj, “Numerical calculation of

lyapunov exponents in various nonlinear chaotic

systems,” International Journal of Scientific &

Technology Research, vol. 3, no. 7, pp. 87–90,

2014.

WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION

DOI: 10.37394/232010.2023.20.1

Joan Jani

E-ISSN: 2224-3410

Volume 20, 2023

Contribution of Individual Authors to the

Creation of a Scientific Article (Ghostwriting

Policy)

The author contributed in the present research, at all

stages from the formulation of the problem to the

final findings and solution.

Sources of Funding for Research Presented in a

Scientific Article or Scientific Article Itself

No funding was received for conducting this study.

Conflict of Interest

The author has no conflict of interest to declare that

is relevant to the content of this article.

Creative Commons Attribution License 4.0

(Attribution 4.0 International, CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en

_US