In [1]:
%run Steppers.ipynb #Import stepper functions (Euler, RK4, etc)
%run ODEs.ipynb #Import ODEs (Pendulum etc)
%run Tools.ipynb
import matplotlib
matplotlib.use('nbagg')
import matplotlib.pyplot as plt
In [2]:
ncycles = 5
t0 = 0
tf = 5*(2*np.pi)
h = 0.05
t = np.linspace(t0,tf,int((tf-t0)/h)) #Generate timesteps
x0=(0,1)
omega=1
b=0
vectorPoints=np.array([(0,1),(0,-1),(1,0),(2,0),(-2,0),(-1,0),(1,1),(-1,-1)])
vectors=vectorPoints.copy()
for i, x in enumerate(vectorPoints):
vectors[i]=(SimplePendulum(x,t,b,omega))
sol = RK4(SimplePendulum, x0, t, args=(b,omega,))
x, v = np.hsplit(sol, len(sol[0]))
plt.plot(x,v, label=r"$x_0,y_0={}$".format(x0))
plt.xlabel("x")
plt.ylabel("v")
x,y=np.hsplit(vectorPoints,2)
u,v=np.hsplit(vectors,2)
plt.quiver(x,y,u,v,)
plt.axis("scaled")
plt.title("Phase Portrait")
plt.legend()
plt.show()
In [3]:
ncycles = 5
t0 = 0
tf = 5*(2*np.pi)
h = 0.05
t = np.linspace(t0,tf,int((tf-t0)/h)) #Generate timesteps
x0=(0,2)
omega=1
b=0.5
vectorPoints=np.array([(0,1),(0,-1),(1,0),(2,0),(-2,0),(-1,0),(1,1),(-1,-1)])
vectors=vectorPoints.copy()
for i, x in enumerate(vectorPoints):
vectors[i]=SimplePendulum(x,t,b,omega)
sol = RK4(SimplePendulum, x0, t, args=(b,omega,))
x, v = np.hsplit(sol, len(sol[0]))
plt.plot(x,v, label=r"$x_0,y_0={}$".format(x0))
plt.xlabel("x")
plt.ylabel("v")
x,y=np.hsplit(vectorPoints,2)
u,v=np.hsplit(vectors,2)
plt.quiver(x,y,u,v)
plt.axis("scaled")
plt.title(r"Phase Portrait with damping $b={}$".format(b))
plt.legend()
plt.show()
In [4]:
t = np.arange(0,30,0.01)
f0=(-1,-2)
sol=RK4(Nonlinear4, f0, t)
x,y=np.hsplit(sol,len(f0))
lim=10
y[x>lim]=np.inf
x[x>lim]=np.inf
y[x<-lim]=-np.inf
x[x<-lim]=-np.inf
plt.plot(t,x)
plt.plot(t,y)
plt.show()
for f0 in ((1.99,1.99),(1,1),(-3,-5),(-4,-4),(-4,5)):
sol=RK4(Nonlinear4, f0, t)
x,y=np.hsplit(sol,len(f0))
plt.plot(x,y, color="C0")
plt.axis("scaled")
plt.xlim((-10,10))
plt.ylim((-10,10))
plt.show()
In [5]:
t = np.arange(0,30,0.01)
xLims=(-6,5)
yLims=(-6,5)
f0=(-1,-2)
sol=RK4(Nonlinear5, f0, t)
x,y=np.hsplit(sol,len(f0))
lim=10
y[x>lim]=np.inf
x[x>lim]=np.inf
y[x<-lim]=-np.inf
x[x<-lim]=-np.inf
plt.plot(t,x)
plt.plot(t,y)
plt.show()
for f0 in ((1.99,1.99),(1,1),(-3,-5),(-4,-4),(-4,5)):
sol=RK4(Nonlinear5, f0, t)
x,y=np.hsplit(sol,len(f0))
plt.plot(x,y, color="C0")
plt.axis("scaled")
PlotQuiver(Nonlinear5,xLims=(xLims[0],xLims[1]),yLims=(yLims[0],yLims[1]),spacing=0.3)
plt.xlim(xLims)
plt.ylim(yLims)
plt.show()
In [ ]:
In [ ]: