# Coding a numerical approximation for a damped pendulum

• Marchionni
In summary, the conversation is about changing a code for a simple pendulum without damping to one with damping. The solution is to add a term for damping and use the Euler-Cromer method. The code provided for the damped pendulum is shown to be effective and produces expected results. The use of the Euler-Cromer method is praised for its elegance and effectiveness.
Marchionni
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
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

Marchionni said:
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

BvU said:
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
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

BvU said:
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.

Marchionni said:
I don't see how I can put this time dependent term into the code within the Euler method
If you are brave you try something similar as before ...

You want to change your damped pendulum to a driven damped pendulum. The differential equation is second order and the homogeneous equation ( 0 instead of ##A \cos (\omega t)## ) has a solution that dampens out. As you have found. That's called the transient solution.

Any single so called particular solution to the http://www.math.psu.edu/tseng/class/Math251/Notes-2nd%20order%20ODE%20pt2.pdf ##\ddot x +\beta \dot x + \omega_0^2 x = A\cos(\omega t)## is then enough to provide the complete solution: transient + particular. The latter is often called the steady-state solution when in fact it's periodic.

## 1. How does a damped pendulum work?

A damped pendulum is a system where the motion of a pendulum is gradually slowed down due to friction or resistance. This causes the pendulum to eventually come to rest at its equilibrium position.

## 2. What is the formula for calculating the period of a damped pendulum?

The formula for calculating the period of a damped pendulum is T = 2π√(L/g), where T is the period, L is the length of the pendulum, and g is the acceleration due to gravity.

## 3. How is the damping coefficient incorporated into the numerical approximation for a damped pendulum?

The damping coefficient is incorporated into the numerical approximation for a damped pendulum by adding a term to the equation for the angular velocity. This term is proportional to the angular velocity and represents the force due to damping.

## 4. What are some common methods for numerically approximating a damped pendulum?

Some common methods for numerically approximating a damped pendulum include Euler's method, Runge-Kutta methods, and the Verlet algorithm.

## 5. How accurate are numerical approximations for a damped pendulum?

The accuracy of a numerical approximation for a damped pendulum depends on the method used and the step size chosen. In general, the smaller the step size, the more accurate the approximation will be. However, the numerical approximation will never be exact due to the simplifications and assumptions made in the model.

• Other Physics Topics
Replies
1
Views
2K
Replies
7
Views
1K
• Programming and Computer Science
Replies
1
Views
1K