在单摆的微分方程中遇到一个奇怪的问题

在3b1b 的微分方程概论第一章中介绍了一个单摆系统的案例，一个单摆在理想情况下做简谐振动。但在实际情况中，需要考虑偏角大小，还有空气阻力。通过受力分析，可以得出角度相关的微分方程。

image.png image.png
在理想模型下的求解方程是比较简单的：

到曲线是一条光滑的余弦曲线。

但在实际模型中通过微分方程

来获得表达式和曲线是十分困难的，所以视频中提供了一种数值逼近的方法来求解任意时刻t的值。

我在视频中提供的代码的基础上加上了一些修改，以便绘制从0～t 的时间内，、、的变化曲线,并与理想模型下对比。

# 单摆系统的微分方程
import numpy as np
import matplotlib.pyplot as plt

# 重力加速度
g = 0.98
# 单摆长度
L = 2
# 空气阻力系数
mu = 0.1
# 初始角度
THEAT_0 = np.pi/3
# 初始角速度
THEAT_DOT_0 = 0


def get_theta_double_dot(theat, theat_dot):
    """[求角速度的加速度]

    Args:
        theat ([type]): [角度]
        theat_dot ([type]): [角速度]

    Returns:
        [type]: [加速度]
    """
    return -mu*theat_dot-(g/L)*np.sin(theat)


def theat(t):
    """[求角度]

    Args:
        t ([type]): [时间t]

    Returns:
        [type]: [时间序列、理想模型下角度序列，实际角度序列、实际角速度序列、实际加速度序列]
    """
    theat = THEAT_0  # 角度 初始化
    theat_dot = THEAT_DOT_0  # 角速度 初始化
    delta_t = 0.01  # 微分时间
    x_data = []  # 时间序列
    y_data_0 = []
    y_data_1 = []
    y_data_2 = []
    y_data_3 = []

    for i in np.arange(0, t, delta_t):
        theat_double_dot = get_theta_double_dot(theat, theat_dot)
        theat += theat_dot*delta_t
        theat_dot += theat_double_dot*delta_t
        print(i, theat)
        x_data.append(i)
        y_data_0.append(THEAT_0*np.cos(np.sqrt(g/L)*i))
        y_data_1.append(theat)
        y_data_2.append(theat_dot)
        y_data_3.append(theat_double_dot)
    return x_data, y_data_0, y_data_1, y_data_2, y_data_3


# 测试。。。。。。。。。。。。
t = 100  # 时间
x, y0, y1, y2, y3 = theat(t)
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.subplot(2, 2, 1)
plt.xlabel('时间')
plt.title('理想角度')
plt.plot(x, y0, color='red', label='theta-temp')


plt.subplot(2, 2, 2)
plt.xlabel('时间')
plt.title('实际角度 theta')
plt.plot(x, y1, color='red')
plt.subplot(2, 2, 3)
plt.xlabel('时间')
plt.title('实际角速度 theta_dot')
plt.plot(x, y2, color='green')

plt.subplot(2, 2, 4)
plt.xlabel('时间')
plt.title('实际加速度 theta_double_dot')

plt.plot(x, y3, color='blue')
plt.tight_layout()
plt.show()

运行绘制得到如下曲线:

Figure_1.png

到这里还是正常的，可以看到在实际情况下，随着时间的流逝，角度会在震荡中不断衰减，直至无限到达零点的稳定位置。而角速度、和加速度也符合预期。所以我开始测试一些其他初始条件下的曲线，比如初始角度为，初始角速度不为0，重力加速度增加或减少，空气阻力系数增大或减少，都是正常的。但当我想测试真空下，也就是空气阻力为0时，函数曲线开始出现一些奇怪的事情。

image.png

表面上看上去很正常，单摆的动能和重力势能相互转化，但仔细看却发现的峰值在不断上升！！角速度也是。也就是说，在没有空气阻力的情况下，初始速度为零，也没有其他外界对其做功，这个单摆反而越摆越高！这完全不合常理。我开始以为是图像渲染的错误问题，所以我把时间t拉长到t=500，并打印出每个时刻的参数数值，却发现依旧是如此，图像渲染并没问题。在其他参数不变，而时，单摆的角度、角速度波峰居然在递增？

image.png

而当我把时间继续拉大，t=1000时，更加奇怪的事情发生了。居然在某个临界值突然放大，然后一路狂奔。也就是说，这个单摆越转越高，越转越快，最后彻底变成了一个轮子停不下来了，我造出了一个永动机！

image.png

我尝试把微分时间减少，但结果依旧是如此。如果你看到这篇文章，请看看我的代码哪里有问题，或者这个微分方程的数值求解哪里有问题，谢谢。

原视频链接 https://www.bilibili.com/video/BV1tb411G72z