Простое моделирование маятника, результат не соответствует действительности, не может найти ошибку

Question

Я смотрел на YouTube канал 3Blue1Brown (Обзор дифференциальных уравнений | Глава 1). В видео рассказывается об имитации маятника, я пытался его кодировать, но появляются странные результаты. Я в основном скопировал код из видео и добавил еще немного.

====
Я пытаюсь смоделировать это дифференциальное уравнение: Это математическая запись, а не программирование :) theta_double_dot = -mu * theta_dot - (г / л) * грех (тета)

Значения:
theta: = угол, theta_dot: = угловая скорость, theta_double dot: = угловое ускорение, mu: = постоянная трения в воздухе, g: = гравитационная постоянная, L: = длина маятника

delta_t: = шаг по времени

====

Начальные условия: тета: = 60 градусов, тета_dot: = 0 скорость

Задача 1:
test1 со значениями:
g = 9,8, L = 2, mu = 0,09, delta_t = 0,01
результат: ускорение вечного движения свободной энергии 4 всем:)

test2 со значениями:
g = 9,8, L = 2, mu = 0,01 (меньшее трение воздуха, чем при испытании выше), delta_t = 0,01
результат: это более реалистично

чем больше сопротивление воздуха, тем больше система стремится к постоянному движению и получению свободной энергии?!

Проблема 2: все значения, которые я выбрал для вызова theta(..some value here..), симуляция ВСЕГДА заканчивается одинаковыми значениями.

Проблема 3: (немного больше): Как я могу сделать код быстрее? Я думаю, что операция добавления действительно дорогостоящая, потому что она должна копировать каждый раз массив.

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)