Coding a numerical approximation for a damped pendulum

[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.

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))
 TE_vec = np.add(KE_vec,PE_vec)
 # make graphs
 if plotFlag == 1:
  plt.figure(0)
  plt.plot(time_vec,theta_vec)
  plt.xlabel('Time (s)')
  plt.ylabel('Displacement (rad)')
  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.xlabel('Displacement (rad)')
  plt.ylabel('Velocity (rad/s)')
  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)$ :

# 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]

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?

omega = omega_old - (g/l)*np.sin(theta_old)*tau

 

[Marchionni]

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)$ :

# 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

 

[Marchionni]

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.

Thank you again in advance, your help is greatly appreciated!

 

[BvU]

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.