# Coding a numerical approximation for a damped pendulum

Hi there.

I have a question about the damped pendulum. I am working on an exercise where I have already numerically approximated the solution for a simple pendulum without dampening. Now, the excercise says that I can simply change the code of this simple situation to describe a pendulum with dampening.

For the simple pendulum without dampening I have used Eulers method of approximating the differential equation. Here is the code for the pendulum without damping.

Code:
import numpy as np
from math import *
import matplotlib.pyplot as plt
def simplePendulumSimulation(theta0,omega0,tau,m,g,l,numSteps,plotFlag):
# This function simulates a simple pendulum using the Euler-Cromer method.
# Inputs: theta0 (the initial angle, in rad)
#    omega0 (the initial angular velocity, in rad/s)
#    tau (the time step)
#    m (mass of the pendulum)
#       g (gravitational constant)
#    l (length of the pendulum)
#    numSteps (number of time steps to take, in s)
#       plotFlag (set to 1 if you want plots, 0 otherwise)
# Outputs: t_vec (the time vector)
#     theta_vec (the angle vector)
# initialize vectors
time_vec = [0]*numSteps
theta_vec = [0]*numSteps
omega_vec = [0]*numSteps
KE_vec = [0]*numSteps
PE_vec = [0]*numSteps
# set initial conditions
theta = theta0
omega = omega0
time = 0
# begin time stepping
for i in range(0,numSteps):
omega_old = omega
theta_old = theta
# update the values
omega = omega_old - (g/l)*np.sin(theta_old)*tau
theta = theta_old + omega*tau
# record the values
time_vec[i] = tau*i
theta_vec[i] = theta
omega_vec[i] = omega
KE_vec[i] = (1/2)*m*l**2*omega**2
PE_vec[i] = m*g*l*(1-np.cos(theta))
# make graphs
if plotFlag == 1:
plt.figure(0)
plt.plot(time_vec,theta_vec)
plt.xlabel('Time (s)')
plt.savefig('plot1.png', bbox_inches='tight')
plt.figure(1)
plt.plot(time_vec,KE_vec,label='Kinetic Energy')
plt.plot(time_vec,PE_vec,label='Potential Energy')
plt.plot(time_vec,TE_vec,label='Total Energy')
plt.legend(loc='upper left')
plt.xlabel('Time (s)')
plt.ylabel('Energy (J)')
plt.savefig('plot2.png', bbox_inches='tight')
plt.figure(2)
plt.plot(theta_vec,omega_vec)
plt.savefig('plot3.png', bbox_inches='tight')

plt.show()

# return the vectors

simplePendulumSimulation(0.2,0,0.03,1,1,1,3000,1)
I was wondering if anyone could tell me if there is a simple way to change this code to introduce damping instead of beginning al over again. Thanks in advance.

BvU

Related MATLAB, Maple, Mathematica, LaTeX News on Phys.org
BvU
Homework Helper
2019 Award
For a damped oscillator you have ##\ddot x +\beta \dot x + \omega_0^2 x = 0## instead of ##\ddot x + \omega_0^2 x = 0## so all you have to do is add a term ##\beta \omega## to ##\tt (g/l)*np.sin(theta\_old)## :
Python:
# omega = omega_old - (g/l)*np.sin(theta_old)*tau
omega = omega_old - (0.02 * omega_old + (g / l) * np.sin(theta_old)) * tau
Works instantly: I see your three plots come out as expected

George Jones
Staff Emeritus
Gold Member
I have a question about the damped pendulum. I am working on an exercise where I have already numerically approximated the solution for a simple pendulum without dampening. Now, the excercise says that I can simply change the code of this simple situation to describe a pendulum with dampening.

For the simple pendulum without dampening I have used Eulers method of approximating the differential equation. Here is the code for the pendulum without damping.

I was wondering if anyone could tell me if there is a simple way to change this code to introduce damping instead of beginning al over again. Thanks in advance.
From where does the following line of your code come?
Code:
omega = omega_old - (g/l)*np.sin(theta_old)*tau

For a damped oscillator you have ##\ddot x +\beta \dot x + \omega_0^2 x = 0## instead of ##\ddot x + \omega_0^2 x = 0## so all you have to do is add a term ##\beta \omega## to ##\tt (g/l)*np.sin(theta\_old)## :
Python:
# omega = omega_old - (g/l)*np.sin(theta_old)*tau
omega = omega_old - (0.02 * omega_old + (g / l) * np.sin(theta_old)) * tau
Works instantly: I see your three plots come out as expected
Thank you so much! I had tried almost the same, expect that I didn't use old omega but just omega and it exploded..

BvU
BvU
Homework Helper
2019 Award
Ah, the peculiarities of this Euler-Cromer method.
I didn't know of it, so I'm grateful for your thread ! -- in earlier threads solving Newton I always encouraged the poster to start with simple Euler and postpone changing to a more complicated intgrator, but this is quite elegant and effective !

PS use [code=python] -- it looks a lot better.

PS2

Ah, the peculiarities of this Euler-Cromer method.
I didn't know of it, so I'm grateful for your thread ! -- in earlier threads solving Newton I always encouraged the poster to start with simple Euler and postpone changing to a more complicated intgrator, but this is quite elegant and effective !

PS use [code=python] -- it looks a lot better.

PS2
Thank you very much! And now you are here I am going to bother you just a little bit more..

I now have to introduce a driving term to the motion, namely:
Acos(ωt)

I don't see how I can put this time dependent term into the code within the Euler method. Maybe you have some thoughts to share? Also, and I have already asked this question in the classical physics thread, do you maybe know what a steady state is in the context of the pendulum. I am familiar with the steady state in Quantum Mechanics but can't immediately see it's meaning when applied on a pendulum.

BvU