python中钟摆运动的数值预测不能正确预测

Question

我正在对一个项目的钟摆运动进行数值预测。在运行程序时，我注意到每个峰值之间的时间（当值达到与 theta 的起始值相同的 theta 值时）与以下等式预测（T = 2pi平方（升/克））。如果有人能看一下我会很感激（我知道它有点乱......哈哈）

p.s 只有蓝色图是相关的，忽略上的猩猩也忽略"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

嗯，一种解释可能是您已经对非线性钟摆建模（这比线性钟摆更现实一点），而您要比较的公式是线性钟摆的周期。非线性摆的周期计算要复杂得多，不是用这么简单的公式给出的，而是用椭圆函数给出的。事实上，非线性摆的周期对于不同的起始位置是不同的。此外，您还在使用与速度相关的阻尼对钟摆进行建模，因此您不能在那里谈论周期。

此外，您在建模力的方式上可能存在错误。确保您具有以下表达式：

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

当b = 0:

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

我对您的脚本进行了一些更正，查看它们并验证公式。我不保证还有其他问题。

    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