在3b1b 的微分方程概論第一章中介紹了一個(gè)單擺系統(tǒng)的案例乾胶，一個(gè)單擺在理想情況下做簡諧振動(dòng)佛点。但在實(shí)際情況中，需要考慮偏角大小挂洛，還有空氣阻力礼预。通過受力分析，可以得出角度相關(guān)的微分方程虏劲。
在理想模型下的求解
到
但在實(shí)際模型中通過微分方程
來獲得表達(dá)式和曲線是十分困難的，所以視頻中提供了一種數(shù)值逼近的方法來求解任意時(shí)刻t的
值柒巫。
我在視頻中提供的代碼的基礎(chǔ)上加上了一些修改励堡，以便繪制從0～t 的時(shí)間內(nèi)，堡掏、
应结、
的變化曲線,并與理想模型下對比。
# 單擺系統(tǒng)的微分方程
import numpy as np
import matplotlib.pyplot as plt
# 重力加速度
g = 0.98
# 單擺長度
L = 2
# 空氣阻力系數(shù)
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]): [時(shí)間t]
Returns:
[type]: [時(shí)間序列泉唁、理想模型下角度序列鹅龄，實(shí)際角度序列、實(shí)際角速度序列亭畜、實(shí)際加速度序列]
"""
theat = THEAT_0 # 角度 初始化
theat_dot = THEAT_DOT_0 # 角速度 初始化
delta_t = 0.01 # 微分時(shí)間
x_data = [] # 時(shí)間序列
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 # 時(shí)間
x, y0, y1, y2, y3 = theat(t)
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.subplot(2, 2, 1)
plt.xlabel('時(shí)間')
plt.title('理想角度')
plt.plot(x, y0, color='red', label='theta-temp')
plt.subplot(2, 2, 2)
plt.xlabel('時(shí)間')
plt.title('實(shí)際角度 theta')
plt.plot(x, y1, color='red')
plt.subplot(2, 2, 3)
plt.xlabel('時(shí)間')
plt.title('實(shí)際角速度 theta_dot')
plt.plot(x, y2, color='green')
plt.subplot(2, 2, 4)
plt.xlabel('時(shí)間')
plt.title('實(shí)際加速度 theta_double_dot')
plt.plot(x, y3, color='blue')
plt.tight_layout()
plt.show()
運(yùn)行繪制得到如下曲線:
到這里還是正常的溶耘，可以看到在實(shí)際情況下二拐，隨著時(shí)間的流逝，角度會(huì)在震蕩中不斷衰減凳兵，直至無限到達(dá)零點(diǎn)的穩(wěn)定位置百新。而角速度、和加速度也符合預(yù)期庐扫。所以我開始測試一些其他初始條件下的曲線饭望，比如初始角度為
铅辞，初始角速度不為0，重力加速度增加或減少萨醒，空氣阻力系數(shù)增大或減少斟珊，都是正常的。但當(dāng)我想測試真空下富纸，也就是空氣阻力為0時(shí)囤踩，
函數(shù)曲線開始出現(xiàn)一些奇怪的事情。
表面上看上去很正常胜嗓，單擺的動(dòng)能和重力勢能相互轉(zhuǎn)化高职，但仔細(xì)看卻發(fā)現(xiàn)的峰值在不斷上升９痴А辞州！角速度也是。也就是說寥粹，在沒有空氣阻力的情況下变过，初始速度為零，也沒有其他外界對其做功涝涤，這個(gè)單擺反而越擺越高媚狰！這完全不合常理。我開始以為是圖像渲染的錯(cuò)誤問題阔拳，所以我把時(shí)間t拉長到t=500崭孤，并打印出每個(gè)時(shí)刻的參數(shù)數(shù)值，卻發(fā)現(xiàn)依舊是如此，圖像渲染并沒問題辨宠。在其他參數(shù)不變遗锣，而
時(shí)，單擺的角度嗤形、角速度波峰居然在遞增精偿？
而當(dāng)我把時(shí)間繼續(xù)拉大，t=1000時(shí)赋兵，更加奇怪的事情發(fā)生了笔咽。居然在某個(gè)臨界值突然放大，然后一路狂奔霹期。也就是說叶组，這個(gè)單擺越轉(zhuǎn)越高，越轉(zhuǎn)越快历造，最后徹底變成了一個(gè)輪子停不下來了扶叉，我造出了一個(gè)永動(dòng)機(jī)！
我嘗試把微分時(shí)間減少帕膜，但結(jié)果依舊是如此枣氧。如果你看到這篇文章，請看看我的代碼哪里有問題垮刹，或者這個(gè)微分方程的數(shù)值求解哪里有問題达吞，謝謝。
