Kernel: Python 3 (system-wide)
In [3]:
#driver: Brandon Soto #navigator: Jasmine Hadad
In [9]:
#exercise: 8.4
In [15]:
from numpy import array,arange,sin from pylab import plot,xlabel,show g = 9.81 l = 0.1 def f(r,t): theta = r[0] omega = r[1] ftheta = omega fomega = -(g/l)*sin(theta) return array([ftheta,fomega],float) omega_0 = 0 #setting the initial conditions theta_0 = 3.124 a = 0.0 #a and b are the time intervals b = 50.0 N = 1000 #steps h = (b-a)/N tpoints = arange(a,b,h) thetapoints = [] omegapoints = [] r = array([theta_0,omega_0],float) #setting the initial conditions for t in tpoints: thetapoints.append(r[0]) omegapoints.append(r[1]) k1 = h*f(r,t) k2 = h*f(r+0.5*k1,t+0.5*h) k3 = h*f(r+0.5*k2,t+0.5*h) k4 = h*f(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 plot(tpoints,thetapoints) # to plot xlabel("t") show()
Out[15]:
In [12]:
#exercise 8.6 part a
In [13]:
from numpy import array,arange,sin from pylab import plot,xlabel,show def f(r,t): x = r[0] y = r[1] fx = y fy = (-omega**2)*x return array([fx,fy],float) x_0 = 1 #intial conditions y_0 = 0 omega = 1 a = 0.0 #a and b are the time intervals b = 50.0 N = 1000 #steps h = (b-a)/N tpoints = arange(a,b,h) xpoints = [] ypoints = [] r = array([x_0,y_0],float) #setting the initial conditions for t in tpoints: xpoints.append(r[0]) ypoints.append(r[1]) k1 = h*f(r,t) k2 = h*f(r+0.5*k1,t+0.5*h) k3 = h*f(r+0.5*k2,t+0.5*h) k4 = h*f(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 plot(tpoints,xpoints) # to plot xlabel("t") show()
Out[13]:
In [14]:
#exercise 8.6 part b
In [16]:
from numpy import array,arange,sin from pylab import plot,xlabel,show def f(r,t): x = r[0] y = r[1] fx = y fy = (-omega**2)*x return array([fx,fy],float) x_0 = 2 #intial conditions y_0 = 0 omega = 1 a = 0.0 #a and b are the time intervals b = 50.0 N = 1000 #steps h = (b-a)/N tpoints = arange(a,b,h) xpoints = [] ypoints = [] r = array([x_0,y_0],float) #setting the initial conditions for t in tpoints: xpoints.append(r[0]) ypoints.append(r[1]) k1 = h*f(r,t) k2 = h*f(r+0.5*k1,t+0.5*h) k3 = h*f(r+0.5*k2,t+0.5*h) k4 = h*f(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 plot(tpoints,xpoints) # to plot xlabel("t") show()
Out[16]:
In [17]:
#exercise 8.6 part c
In [21]:
from numpy import array,arange,sin from pylab import plot,xlabel,show def f(r,t): x = r[0] y = r[1] fx = y fy = (-omega**2)*x**3 return array([fx,fy],float) x_0 = 3 #intial conditions y_0 = 0 omega = 1 a = 0.0 #a and b are the time intervals b = 50.0 N = 1000 #steps h = (b-a)/N tpoints = arange(a,b,h) xpoints = [] ypoints = [] r = array([x_0,y_0],float) #setting the initial conditions for t in tpoints: xpoints.append(r[0]) ypoints.append(r[1]) k1 = h*f(r,t) k2 = h*f(r+0.5*k1,t+0.5*h) k3 = h*f(r+0.5*k2,t+0.5*h) k4 = h*f(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 plot(tpoints,xpoints) # to plot xlabel("t") show()
Out[21]:
In [19]:
#exercise 8.6 part d
In [31]:
from numpy import array,arange,sin from pylab import plot,xlabel,show def f(r,t): x = r[0] y = r[1] fx = y fy = (-omega**2)*x return array([fx,fy],float) x_0 = 1 #intial conditions y_0 = 0 omega = 1 a = 0.0 #a and b are the time intervals b = 50.0 N = 1000 #steps h = (b-a)/N tpoints = arange(a,b,h) xpoints = [] ypoints = [] r = array([x_0,y_0],float) #setting the initial conditions for t in tpoints: xpoints.append(r[0]) ypoints.append(r[1]) k1 = h*f(r,t) k2 = h*f(r+0.5*k1,t+0.5*h) k3 = h*f(r+0.5*k2,t+0.5*h) k4 = h*f(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 plot(xpoints,ypoints) # to plot show() #when x is cubed in fy then we get an oval rather than a circle
Out[31]:
In [23]:
#exercise 8.6 part e
In [29]:
from numpy import array,arange,sin from pylab import plot,xlabel,show def f(r,t): x = r[0] y = r[1] fx = y fy = mu*(1-x**2)*y-omega**2*x return array([fx,fy],float) x_0 = 1 #intial conditions y_0 = 0 omega = 1 mu = 1 a = 0.0 #a and b are the time intervals b = 20.0 N = 1000 #steps h = (b-a)/N tpoints = arange(a,b,h) xpoints = [] ypoints = [] r = array([x_0,y_0],float) #setting the initial conditions for t in tpoints: xpoints.append(r[0]) ypoints.append(r[1]) k1 = h*f(r,t) k2 = h*f(r+0.5*k1,t+0.5*h) k3 = h*f(r+0.5*k2,t+0.5*h) k4 = h*f(r+k3,t+h) r += (k1+2*k2+2*k3+k4)/6 plot(xpoints,ypoints) # to plot show() #as mu increases than the space between the osscillations become closer together
Out[29]:
In [0]: