Видео с пояснениями теории и кода, файлы с исходным кодом и теорией по маятнику можно скачать из тг-канала: https://t.me/oldfiz/9?single.
Для копирования кода щёлкните на него и нажмите на появившийся справа значок со стрелкой, направленной вниз.
import numpy as np
import scipy
import matplotlib.pyplot as plt
def sys_x(z,t):
phi,omega = z
dphi_dt=omega
domega_dt=-g*np.sin(phi)/L-rd*omega
return [dphi_dt, domega_dt]
g=9.81
L=1
m=1
V0=0
w0=V0/L
gamma_=90
gamma=gamma_*np.pi/180
alpha=0.00
rd=alpha/m
z0x=[gamma,w0] #начальные условия
td=5 #длительность колебаний
t_eval=np.linspace(0,td,1000)
sol_phi=scipy.integrate.odeint(sys_x,z0x,t_eval)
phi=sol_phi[:,0]
omega=sol_phi[:,1]
phi_acc=gamma*np.cos(np.power(g/L,0.5)*t_eval)
diff_phi=phi_acc-phi
plt.title('Зависимость угла от времени')
plt.xlabel('Время, с')
plt.ylabel('Угол, радианы')
#plt.plot(phi,omega,label='Фазовый портрет',linewidth=3)
plt.plot(t_eval,phi,label='Точное решение',linewidth=3)
plt.plot(t_eval,phi_acc,label='Приближеное решение',linewidth=3,linestyle='dashed')
plt.plot(t_eval,diff_phi,label='Разница решений',linewidth=3,linestyle='dotted')
#plt.plot(t_eval,omega,label='Точное решение угловой скорости',linewidth=3)
plt.grid(True)
plt.legend(loc=8)
plt.show()