# 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

BvU
Homework Helper
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
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.