python odeint

1#Solving ODE in python using ODEINT 
2#Luv Kumar
3
4import numpy as np
5from scipy.integrate import odeint
6import matplotlib.pyplot as plt
7import math
8#function that return dtheta/dt
9def model(theta,t,b,g,l,m):
10	theta1 = theta[0]
11	theta2 = theta [1]
12	dtheta1_dt = theta2
13	dtheta2_dt = -(b/m)*theta2 -(g/l)*math.sin(theta1)
14	dtheta_dt = [dtheta1_dt,dtheta2_dt]
15	return dtheta_dt
16
17
18
19
20
21b = 0.02
22g = 9.81
23l = 1
24m = 0.1
25
26#initial condition
27theta_0 = [0,5]
28#time points
29t = np.linspace(0,10,150)
30#solve ODE
31theta = odeint(model,theta_0,t,args = (b,g,l,m))
32#plot results
33plt.plot(t,theta[:,0])
34plt.plot(t,theta[:,1])
35plt.ylabel('plot')
36plt.xlabel('time')
37plt.legend()
38plt.show()
39
40ct = 1
41for i in theta[:,0]:
42	x0 = 0
43	y0 = 0
44	x3 = -1.5
45	y3 = 0
46	x4 = 1.5
47	y4 = 0
48	x1 = l*math.sin(i)
49	y1 = -l*math.cos(i)
50	filename = str(ct)+'.png'
51	ct = ct +1
52	plt.figure()
53	plt.plot([x3,x4],[y3,y4])
54	plt.plot([x0,x1],[y0,y1],x1,y1,'o',markersize=20,color='red')
55	plt.xlim([-2,2])
56	plt.ylim([-2,2])
57	plt.title('oscillating pendulum')
58	plt.savefig(filename)
59
60