Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download

Jupyter notebook test-zhidkov.ipynb

Project: zhidkov
Views: 191
Kernel: Python 2 (system-wide)
from numpy import sin, cos import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate from ipywidgets import *
G = 9.8 # acceleration due to gravity, in m/s^2 def derivs(state, t,params): L1, L2, M1, M2 = params dydx = np.zeros_like(state) dydx[0] = state[1] del_ = state[2] - state[0] den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_) dydx[1] = (M2*L1*state[1]*state[1]*sin(del_)*cos(del_) + M2*G*sin(state[2])*cos(del_) + M2*L2*state[3]*state[3]*sin(del_) - (M1 + M2)*G*sin(state[0]))/den1 dydx[2] = state[3] den2 = (L2/L1)*den1 dydx[3] = (-M2*L2*state[3]*state[3]*sin(del_)*cos(del_) + (M1 + M2)*G*sin(state[0])*cos(del_) - (M1 + M2)*L1*state[1]*state[1]*sin(del_) - (M1 + M2)*G*sin(state[2]))/den2 return dydx
def example(th1,w1,th2,w2,l1,l2,m1,m2, Tmax): pars = (l1, l2, m1, m2) ders = lambda state, t: derivs(state, t, pars) dt = 0.001 t = np.arange(0.0, Tmax, dt) state = np.radians([th1, w1, th2, w2]) y = integrate.odeint(ders, state, t) x1 = l1*sin(y[:, 0]) y1 = -l1*cos(y[:, 0]) x2 = l2*sin(y[:, 2]) + x1 y2 = -l2*cos(y[:, 2]) + y1 plt.plot(x2,y2, label = 'bottom pendulum') plt.legend() plt.scatter([0.0],[0.0],c='orange',s=36.0) plt.scatter([x2[0]],[y2[0]],c='green',s=36.0) plt.scatter([x2[-1]],[y2[-1]],c='red',s=36.0)
interact(example, th1 = BoundedFloatText(value=120.0, min = - 180.0, max = 180.0, description='Angle 1 (degrees):'), w1 = FloatText(value=0.0, description='Angular velocity 1:'), th2 = BoundedFloatText(value=-10.0, min = - 180.0, max = 180.0, description='Angle 2 (degrees):'), w2 = FloatText(value=0.0, description='Angular velocity 2:'), l1 = BoundedFloatText(value=1.0, min = 0.01, max = 5.0, description='L1:'), l2 = BoundedFloatText(value=1.0, min = 0.01, max = 5.0, description='L2:'), m1 = BoundedFloatText(value=1.0, min = 0.01, max = 10.0, description='M1:'), m2 = BoundedFloatText(value=1.0, min = 0.01, max = 10.0, description='M2:'), Tmax = BoundedFloatText(value=20.0, min = 0.01, max = 1000.0, description='Max time:'), __manual=True)
<function __main__.example>