Period off for simple pendulum using Runge Kutta

Question

I'm trying to model a simple pendulum using Runge-Kutta methods in python, and the period should be 2 seconds. Not only that, my values for angular velocity and angular displacement are also incorrect. I'm not sure what my mistake is in my code.

This is the code I tried - v1 is the velocity step, and t1 is the angular displacement step.

g = 9.81
pi = math.pi
l = g/pi

t0 = 0
t = t0
t_list = [t0]

v0 = 0
v = v0
v_list = [v0]

theta0 = 0.1
theta = theta0
theta_list = [theta0]

dt = 0.02
tfinal = 20

def theta2(t, v, theta):
  return -(g/l)*math.sin(theta)

def theta1(t, v, theta):
  return v

while t <= tfinal:
  if t == 0:
    v1 = ((dt)/2)*theta2(t, v, theta)
  else:
    v1 = theta1(t-(dt)/2, v, theta) + (dt)*theta2(t, v, theta)
  t1 = theta + (dt)*v1
  
  v = v1
  theta = t1
  t = t + dt

  t_list.append(round(t, 2))
  v_list.append(v)
  theta_list.append(theta)

print(t_list)
print(v_list)
print(theta_list)
ax = plt.plot(t_list, theta_list)
ax = plt.xlabel('time')
ax = plt.ylabel('theta')
#ax = plt.xlim(0,3.5)
ax = plt.show()

This is the resulting graph. Period should be 2 seconds.

Answer 1

You have a few issues. One of the main reasons is you used an incorrect length for the pendulum.

To get the pendulum length you need to solve for:

T = 2 * pi * sqrt(l/g) ==> l=g/pi^2

Its also been awhile since I took a class in numerical dynamical systems but this does not look like Runge-Kutta so i coded the 4th order Runge-Kutta method. This seemed to provide the correct results.

import math
import matplotlib.pyplot as plt

# Constants
g = 9.81
pi = math.pi
l = g / (pi ** 2)  # Corrected length for a period of 2 seconds

# Initial conditions
t0 = 0
v0 = 0
theta0 = 0.1  # radians

dt = 0.02
tfinal = 20

# Initialize variables
t = t0
v = v0
theta = theta0

# Initialize lists to store data
t_list = [t0]
v_list = [v0]
theta_list = [theta0]

# Define derivative functions
def theta_dot(t, theta, v):
    return v

def v_dot(t, theta, v):
    return - (g / l) * math.sin(theta)

# Runge-Kutta 4th order method
while t <= tfinal:
    # Compute k1
    k1_theta = dt * theta_dot(t, theta, v)
    k1_v = dt * v_dot(t, theta, v)
    
    # Compute k2
    k2_theta = dt * theta_dot(t + dt/2, theta + k1_theta/2, v + k1_v/2)
    k2_v = dt * v_dot(t + dt/2, theta + k1_theta/2, v + k1_v/2)
    
    # Compute k3
    k3_theta = dt * theta_dot(t + dt/2, theta + k2_theta/2, v + k2_v/2)
    k3_v = dt * v_dot(t + dt/2, theta + k2_theta/2, v + k2_v/2)
    
    # Compute k4
    k4_theta = dt * theta_dot(t + dt, theta + k3_theta, v + k3_v)
    k4_v = dt * v_dot(t + dt, theta + k3_theta, v + k3_v)
    
    # Update theta and v
    theta += (k1_theta + 2*k2_theta + 2*k3_theta + k4_theta) / 6
    v += (k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
    
    # Update time
    t += dt
    
    # Append data to lists
    t_list.append(t)
    theta_list.append(theta)
    v_list.append(v)

# Plotting after simulation
plt.plot(t_list, theta_list)
plt.xlabel('Time (s)')
plt.ylabel('Theta (rad)')
plt.title('Simple Pendulum Motion')
plt.show()