numerical prediction of a pendulums motion in python not predicting properly

Question

I'm working on a numerical prediction of a pendulums motion for a project. while running the program I've noticed the time between each peak (when the value reaches the same theta value as the beginning value of theta) is not the same as the following equation predicts(T = 2pisqrt(l/g)). If there's anyone how could take a look I would be so grateful(I know its kind of a mess...lol)

p.s only the blue graph is relevant, ignore the orang on also ignore the path of "graph_pick == "2"".

import math
import matplotlib.pyplot as plt
import tabulate
import pandas as pd

time = list()
amp = list()
velocity = list()
force = list()
acceleration = list()
step = list()
velocity_change = list()
step_change = list()
time_change = list()
amplitude_change = list()
period = list()

def print_plot(x, y):
    t_data = pd.read_excel(r'C:\Users\n\Desktop\עבודת גמר פיזיקה\data for python numeric model.xlsx', usecols="A:A",
                           engine='openpyxl')
    theta_data = pd.read_excel(r'C:\Users\n\Desktop\עבודת גמר פיזיקה\data for python numeric model.xlsx', usecols="B:B",
                               engine='openpyxl')
    plt.figure()
    plt.scatter(x, y)
    plt.scatter(t_data, theta_data)
    plt.xlabel("time")
    plt.ylabel("amplitude")
    plt.show()


def prediction(dt_step, mass, b_const, gravity, length, angle, drag):
    t = 0
    v = 0
    i = 0
    step.append(i)
    if drag == '1':
        f = round(-b_const * length * v - mass * gravity * math.sin(angle), 4)
    if drag == '2':
        f = round(-mass * gravity * length * math.sin(angle), 4)
    a = round(f / mass, 4)

    time.append(t)
    amp.append(angle)
    velocity.append(v)
    force.append(f)
    acceleration.append(a)

    while t <= 30:
        i += 1
        t = round(t + dt_step, 4)
        v = round(v + a * dt_step, 4)
        angle = round(angle + v * dt_step, 4)
        if drag == '1':
            f = round(-b_const * length * v - mass * gravity * length * math.sin(angle), 4)
        if drag == '2':
            f = round(-mass * gravity * length * math.sin(angle), 4)
        a = round(f / mass, 4)

        time.append(t)
        amp.append(angle)
        velocity.append(v)
        force.append(f)
        acceleration.append(a)
        step.append(i)

    return time, amp, velocity, force, acceleration, step


def print_table(t, a, v, f, acc, i):
    print(str(i) + "   " + str(t[i]) + "    " + str(a[i]) + "    " + str(f[i]) + "    " + str(v[i]) + "    " + str(
        acc[i]))

# constants
dt = 0.01
m = 0.056
b = 3.49 * pow(10, -5)
g = 9.81
l = 0.7
theta = 0.172
print("calculate --> 1:with drag  2:no drag")
use_drag = input("enter: ")
print("graph --> 1:amplitude over time 2:period over time")
graph_pick = input("enter: ")

time, amp, velocity, force, acceleration, step = prediction(dt, m, b, g, l, theta, use_drag)

if graph_pick == "1":
    print(tabulate.tabulate(
        {"step": step, "time": time, "amplitude": amp, "velocity": velocity, "force": force, "acceleration": acceleration},
        headers=["step", "time", "amplitude", "velocity", "force", "acceleration"]))
  
    print_plot(time, amp)

if graph_pick == "2":
    velocity_change.append(velocity[0])
    step_change.append(0)
    time_change.append(time[0])
    amplitude_change.append(amp[0])
    period.append(0)
    k = 1
    end = 30/dt
    n = 1
    while k <= end:
        if abs(amp[k]) >= abs(amp[k-1]) and abs(amp[k]) >= abs(amp[k+1]):
            if abs(velocity[k]) < abs(velocity[k-1]) and abs(velocity[k]) < abs(velocity[k+1]):
                velocity_change.append(velocity[k])
                step_change.append(k)
                time_change.append(time[k])
                amplitude_change.append(amp[k])
                n += 1
        k += 1


    print(tabulate.tabulate(
       {"step": step_change, "time": time_change, "amplitude": amplitude_change, "velocity": velocity_change},
       headers=["step", "time", "amplitude", "velocity"]))

""""
 if velocity[k - 1] > 0 and velocity[k + 1] < 0:
    velocity_change.append(velocity[k])
    step_change.append(k)
    time_change.append(time[k])
    amplitude_change.append(amp[k])
    n += 1
if velocity[k - 1] < 0 and velocity[k + 1] > 0:
    velocity_change.append(velocity[k])
    step_change.append(k)
    time_change.append(time[k])
    amplitude_change.append(amp[k])
    n += 1
k += 1

Answer 1

Well, one explanation might be that you have modeled the non-linear pendulum (which is a bit more realistic than the linear one), while the formula you are comparing to is the period for the linear pendulum. The calculation for the period of the non-linear pendulum is much more involved and is not given by such a simple formula, but in terms of an elliptic function. In fact, the period for the non linear pendulum is different for different starting positions. Moreover, you are also modelling the pendulum with velocity dependent damping, so there you cannot talk about periods there.

Also, you may have errors in the way you are modelling the force. Make sure you have the following expressions:

f = - b*l*v - m*g*sin(angle)
a = - (b/m)*l*v - g*sin(angle)

and when b = 0:

f = - m*g*sin(angle)
a = - g*sin(angle)

I made some corrections to your script, look at them and verify the formulas. I do not guarantee there are other issues.

    import math
    import matplotlib.pyplot as plt
    import tabulate
    import pandas as pd
    
    time = list()
    amp = list()
    velocity = list()
    force = list()
    acceleration = list()
    step = list()
    velocity_change = list()
    step_change = list()
    time_change = list()
    amplitude_change = list()
    period = list()
    
    def print_plot(x, y):
        t_data = pd.read_excel(r'C:\Users\n\Desktop\עבודת גמר פיזיקה\data for python numeric model.xlsx', usecols="A:A",
                               engine='openpyxl')
        theta_data = pd.read_excel(r'C:\Users\n\Desktop\עבודת גמר פיזיקה\data for python numeric model.xlsx', usecols="B:B",
                                   engine='openpyxl')
        plt.figure()
        plt.scatter(x, y)
        plt.scatter(t_data, theta_data)
        plt.xlabel("time")
        plt.ylabel("amplitude")
        plt.show()
    
    
    def prediction(dt_step, mass, b_const, gravity, length, angle, drag):
        t = 0
        v = 0
        i = 0
        step.append(i)
        if drag == '1':
            f = round( - b_const * length * v - mass * gravity * math.sin(angle), 4)
        if drag == '2':
            # old version was:
            #f = round(-mass * gravity * length * math.sin(angle), 4)
            # new version is:
            f = round(- mass * gravity * math.sin(angle), 4)
        a = round(f / mass, 4)
    
        time.append(t)
        amp.append(angle)
        velocity.append(v)
        force.append(f)
        acceleration.append(a)
    
        while t <= 30:
            i += 1
            t = round(t + dt_step, 4)
            v = round(v + a * dt_step, 4)
            # old version:
            # angle = round(angle + v * dt_step, 4)
            # new version:
            angle = round(angle + (v/length) * dt_step, 4)
            if drag == '1':
                # old version:          
                # f = round(-b_const * length * v - mass * gravity * length * math.sin(angle), 4)
                # new version:
                f = round( - b_const * length * v - mass * gravity * math.sin(angle), 4)
            if drag == '2':
                # old version was:
                #f = round(-mass * gravity * length * math.sin(angle), 4)
                # new version is:
                f = round(- mass * gravity * math.sin(angle), 4)
            a = round(f / mass, 4)
    
            time.append(t)
            amp.append(angle)
            velocity.append(v)
            force.append(f)
            acceleration.append(a)
            step.append(i)
    
        return time, amp, velocity, force, acceleration, step
    
    
    def print_table(t, a, v, f, acc, i):
        print(str(i) + "   " + str(t[i]) + "    " + str(a[i]) + "    " + str(f[i]) + "    " + str(v[i]) + "    " + str(
            acc[i]))
    
    # constants
    dt = 0.01
    m = 0.056
    b = 3.49 * pow(10, -5)
    g = 9.81
    l = 0.7
    theta = 0.172
    print("calculate --> 1:with drag  2:no drag")
    use_drag = input("enter: ")
    print("graph --> 1:amplitude over time 2:period over time")
    graph_pick = input("enter: ")
    
    time, amp, velocity, force, acceleration, step = prediction(dt, m, b, g, l, theta, use_drag)
    
    if graph_pick == "1":
        print(tabulate.tabulate(
            {"step": step, "time": time, "amplitude": amp, "velocity": velocity, "force": force, "acceleration": acceleration},
            headers=["step", "time", "amplitude", "velocity", "force", "acceleration"]))
      
        print_plot(time, amp)
    
    if graph_pick == "2":
        velocity_change.append(velocity[0])
        step_change.append(0)
        time_change.append(time[0])
        amplitude_change.append(amp[0])
        period.append(0)
        k = 1
        end = 30/dt
        n = 1
        while k <= end:
            if abs(amp[k]) >= abs(amp[k-1]) and abs(amp[k]) >= abs(amp[k+1]):
                if abs(velocity[k]) < abs(velocity[k-1]) and abs(velocity[k]) < abs(velocity[k+1]):
                    velocity_change.append(velocity[k])
                    step_change.append(k)
                    time_change.append(time[k])
                    amplitude_change.append(amp[k])
                    n += 1
            k += 1
    
    
        print(tabulate.tabulate(
           {"step": step_change, "time": time_change, "amplitude": amplitude_change, "velocity": velocity_change},
           headers=["step", "time", "amplitude", "velocity"]))
    
    """"
     if velocity[k - 1] > 0 and velocity[k + 1] < 0:
        velocity_change.append(velocity[k])
        step_change.append(k)
        time_change.append(time[k])
        amplitude_change.append(amp[k])
        n += 1
    if velocity[k - 1] < 0 and velocity[k + 1] > 0:
        velocity_change.append(velocity[k])
        step_change.append(k)
        time_change.append(time[k])
        amplitude_change.append(amp[k])
        n += 1
    k += 1