import numpy as np
import matplotlib.pyplot as plt
# dtheta/dt=w
# dw/dt=bcos(w_0*t)-sin(theta)-q*dtheta/dt
def f(r, t, wqb):
f_theta=r[1]
f_w=wqb[2]*np.cos(wqb[0]*t)-np.sin(r[0])-wqb[1]*r[1]
return np.array([f_theta, f_w])
if __name__=="__main__":
wqb=[[2/3, 1/2, 9/10], [2/3, 1/2, 1.15]]
g=9.8
l=2; theta_0=np.pi/2; m=2; w_0=-1
t=np.linspace(0, 30, 1000)
h=t[1]-t[0]
tt=[]
tw=[]
for _wqb in wqb:
theta=np.zeros_like(t)
theta[0]=theta_0
w=np.zeros_like(t)
w[0]=0
for i in range(len(t)-1):
r=[theta[i], w[i]]
k_1=h*f(r, t[i], _wqb)
k_2=h*f(r+k_1/2, t[i]+h/2, _wqb)
k_3=h*f(r+k_2/2, t[i]+h/2, _wqb)
k_4=h*f(r+k_3, t[i]+h, _wqb)
theta[i+1], w[i+1]=r+(k_1+2*k_2+2*k_3+k_4)/6
tt.append(theta)
tw.append(w)
plt.plot(t, tt[0])
plt.plot(t, tt[1])
plt.plot(t, tw[0])
plt.plot(t, tw[1])
plt.show()