0

I am having a trouble in providing a graphical representation of my system, which happens to be a harmonically driven pendulum. The problem is shown below for reference. Problem

The source code I used is shown below using the Verlet Scheme.

#Import needed modules
import numpy as np
import matplotlib.pyplot as plt

#Initialize variables (Initial conditions)
g = 9.8 #Gravitational Acceleration
L = 2.0 #Length of the Pendulum
A0 = 3.0 #Initial amplitude of the driving acceleration
v0 = 0.0 #Initial velocity
theta0 = 90*np.pi/180 #Initial Angle
drivingPeriod = 20.0 #Driving Period


#Setting time array for graph visualization
tau = 0.1 #Time Step
tStop = 10.0 #Maximum time for graph visualization derived from Kinematics
t = np.arange(0., tStop+tau, tau) #Array of time
theta = np.zeros(len(t))
v = np.zeros(len(t))


#Verlet Method
theta[0] = theta0
v[0] = v0
for i in range(len(t)-1):
    accel = -((g + (A0*np.sin((2*np.pi*t) / drivingPeriod)))/L) * np.sin(theta[i])
    theta[i+1] = theta[i] + tau*v[i] + 0.5*tau**2*accel[i]
    v[i+1] = v[i] + 0.5*tau*(accel[i] + accel[i+1])


#Plotting and saving the resulting graph
fig, ax1 = plt.subplots(figsize=(7.5,4.5))
ax1.plot(t,theta*(180/np.pi))
ax1.set_xlabel("Time (t)")
ax1.set_ylabel("Theta")

plt.show()

A sample output is shown. Output

The pendulum should just go back to its initial angle. How can I solve this issue? Notice that as time evolves, my angle measure (degrees) also increases. I want it to only have a domain of 0 degrees to 360 degrees.

Rafael Gomez
  • 11
  • 1

1 Answers1

0

Please change the array computation

    accel = -((g + (A0*np.sin((2*np.pi*t) / drivingPeriod)))/L) * np.sin(theta[i])
    theta[i+1] = theta[i] + tau*v[i] + 0.5*tau**2*accel[i]

into the proper element computation at the correct place

    theta[i+1] = theta[i] + tau*v[i] + 0.5*tau**2*accel[i]
    accel[i+1] = -((g + (A0*np.sin((2*np.pi*t[i+1]) / drivingPeriod)))/L) * np.sin(theta[i+1])

Note that you need to compute accel[0] separately outside the loop.

It makes the code more readable if you separate out the specifics of the physical model and declare at the start

def acceleration(t,theta):
    return -((g + (A0*np.sin((2*np.pi*t) / drivingPeriod)))/L) * np.sin(theta)

so that later you just call

accel[i+1]=acceleration(t[i+1],theta[i+1])

And even then, with a forced oscillation your system is open, it is possible that the forcing action pumps energy into the pendulum, causing it to start to rotate. This is what your graph shows.

The Verlet method like any symplectic method only promises to somewhat have a constant energy if the system is closed and conservative, that is, in the most common cases, no outside influence and all forces are gradient forces.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • I think there is no need to put `accel[i+1]` inside the for loop as it will violate the Verlet Scheme that I am using. – Rafael Gomez Mar 21 '21 at 09:13
  • Edit: I now know what the graph looks like! I made used of the scipy.integrate module of python and offers way better readability. – Rafael Gomez Mar 21 '21 at 10:38
  • Your first comment is contradicted by the formula for `v[i+1]` in the loop. My remark is just that you need to compute that value properly. // And yes, it is always a good idea to test ones own attempt against established solutions. What you end up with depends on your goals, if it is just to get a good solution or if it is to learn about integration methods. – Lutz Lehmann Mar 21 '21 at 13:02