单摆仿真，结果与实际不符，找不到

Question

我在youtube上观看了3Blue1Brown频道（微分方程概述|第1章）。 这段视频讲的是模拟钟摆，我试图编码，但奇怪的结果出现了。 我基本上复制了视频中的代码，只添加了一些东西。你知道吗

===
我试图模拟这个微分方程： 这是数学符号，不是编程：） θ双θ点=-mu*θ点-（g/L）*sin（θ）

含义：
θ：=角，θ点：=角速度，θ双点：=角加速度，μ：=空气摩擦常数，g：=重力常数，L：=摆长

delta_t:=是时间步长

======

初始条件： θ：=60度，θ点：=0速度

问题1:
测试1的值：
g=9.8，L=2，mu=0.09，δt=0.01
结果：加速永动机自由能4:）

测试2的值：
g=9.8，L=2，mu=0.01（空气摩擦力小于上述试验），δt=0.01
结果：这是更现实的

空气阻力越大，系统就越倾向于永动机和自由能的获得？！你知道吗

问题2： 无论我选择了什么值来调用theta(..some value here..)，模拟总是以相同的值结束。你知道吗

问题3：（一些额外的）：我怎样才能使代码更快？ 我认为append的操作非常昂贵，因为每次都要复制数组。你知道吗

import numpy as np
import matplotlib.pyplot as plt

# Physical constants
g = 9.8
L = 2
mu = 0.1

# Initial conditions
THETA_0 = np.pi / 3  # 60 degrees
THETA_DOT_0 = 0  # No initial angual velocity

theta_values = np.array([])
theta_dot_values = np.array([])
theta_double_dot_values = np.array([])

# Definition of ODE
def get_theta_double_dot(theta, theta_dot):
    return -mu * theta_dot - (g / L) * np.sin(theta)


def print_graph(t, delta_t):
    plt.plot(np.arange(0, t, delta_t), theta_values, label='theta')
    plt.plot(np.arange(0, t, delta_t), theta_dot_values, label='theta_dot')
    plt.plot(np.arange(0, t, delta_t), theta_double_dot_values, label='theta_double_dot')
    plt.legend()
    plt.show()

# Solution to the differential equation
def theta(t):
    # Initialize changing values
    global theta_values
    global theta_dot_values
    global theta_double_dot_values

    theta = THETA_0
    theta_dot = THETA_DOT_0
    delta_t = 0.001  # Some time step

    for time in np.arange(0, t, delta_t):
        theta_double_dot = get_theta_double_dot(theta, theta_dot)

        #  saving values into arrays, this is for printing results
        theta_values = np.append(theta, theta_values)
        theta_dot_values = np.append(theta_dot, theta_dot_values)
        theta_double_dot_values = np.append(theta_double_dot, theta_double_dot_values)

        # updating thetas values
        theta += theta_dot * delta_t
        theta_dot += theta_double_dot * delta_t

    print_graph(t, delta_t)
    return theta

theta(15)