In [2]:
%matplotlib inline
import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
plt.style.use('bmh')
In [3]:
# 4th-order Runge-Kutta
def rk4(x, t, h, f):
# x is coordinates (as a vector)
# h is timestep
# f(x) is a function that returns the derivative
# "Slopes"
k1 = f(x, t)
k2 = f(x + k1*h/2, t + h/2)
k3 = f(x + k2*h/2, t + h/2)
k4 = f(x + k3*h, t + h)
# Update time and position
x_ = x + h*(k1 + 2*k2 + 2*k3 + k4)/6
return x_
def streamline(x, T, f, h = 0.05):
# Return a trajectory of length T (in time units)
# starting at position given by the vector x
# Calculated by Runge-Kutta, with timestep h
# Number of timesteps
Nt = int(T / h)
traj = np.zeros((2, Nt))
traj[:,0] = x
for i in range(Nt-1):
t = i*h
traj[:,i+1] = rk4(traj[:,i], h, t, f)
return traj
def get_streamlines(f, T, h, N, xlims, ylims):
Nt = int(T/h)
trajectories = np.zeros((2, Nt, N))
trajectories[0,0,:] = np.random.uniform(low = xlims[0], high = xlims[1], size = N)
trajectories[1,0,:] = np.random.uniform(low = ylims[0], high = ylims[1], size = N)
for i in range(1, Nt):
t = i*h
trajectories[:,i,:] = rk4(trajectories[:,i-1,:], t, h, f)
return trajectories
def plot_streamlines(f, T, h, N, xlims, ylims, figsize=(9,9)):
# Calculate N streamlines from random initial positions
lines = get_streamlines(f, T, h, N, xlims, ylims)
# Plot all the lines
fig = plt.figure(figsize = figsize)
plt.plot(lines[0,:,:], lines[1,:,:], lw = 1, c = '#A60628', alpha = 0.25)
# Set plot limits
plt.xlim(*xlims)
plt.ylim(*ylims)
In [4]:
def f(x, t):
dx = x[0,:] + np.exp(-x[1])
dy = -x[1]
return np.array([dx, dy])
T = 0.5
h = 0.005
N = 1000
xlims = (-3, 3)
ylims = (-2, 2)
plot_streamlines(f, T, h, N, xlims, ylims)
In [5]:
def f(x, t):
a = 0
dx = -x[1] + a*x[0]*(x[0]**2 + x[1]**2)
dy = x[0] + a*x[1]*(x[0]**2 + x[1]**2)
return np.array([dx, dy])
T = 0.5
h = 0.005
N = 1000
xlims = (-3, 3)
ylims = (-3, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [6]:
def f(X, t):
x = X[0]
y = X[1]
dx = x*(3 - x - 2*y)
dy = y*(2 - x - y)
return np.array([dx, dy])
T = 2
h = 0.002
N = 2000
xlims = (0, 4)
ylims = (0, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [7]:
def f(X, t):
x = X[0]
y = X[1]
dx = y
dy = x - x**3
return np.array([dx, dy])
T = 5
h = 0.005
N = 200
xlims = (-3, 3)
ylims = (-3, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [8]:
def f(X, t):
x = X[0]
y = X[1]
dx = -2*np.cos(x) - np.cos(y)
dy = -2*np.cos(y) - np.cos(x)
return np.array([dx, dy])
T = 5
h = 0.005
N = 2000
xlims = (-10, 10)
ylims = (-10, 10)
plot_streamlines(f, T, h, N, xlims, ylims)
In [9]:
def f(X, t):
x = X[0]
y = X[1]
dx = y
dy = -np.sin(x)
return np.array([dx, dy])
T = 1
h = 0.001
N = 1000
xlims = (-7, 7)
ylims = (-3.5, 3.5)
plot_streamlines(f, T, h, N, xlims, ylims)
In [10]:
def f(X, t):
x = X[0]
y = X[1]
dx = x**2*y
dy = x**2 - y**2
return np.array([dx, dy])
T = 1.0
h = 0.001
N = 1000
xlims = (-2, 2)
ylims = (-2, 2)
plot_streamlines(f, T, h, N, xlims, ylims)
In [11]:
def f(X, t):
x = X[0]
y = X[1]
r = np.sqrt(x**2 + y**2)
t = np.arctan2(y, x)
dr = np.cos(2*2*np.pi*r)
dt = 1/2
dx = dr*np.cos(t) - r*dt*np.sin(t)
dy = dr*np.sin(t) + r*dt*np.cos(t)
return np.array([dx, dy])
T = 0.5
h = 0.005
N = 5000
xlims = (-2, 2)
ylims = (-2, 2)
plot_streamlines(f, T, h, N, xlims, ylims)
In [12]:
def f(X, t):
x = X[0]
y = X[1]
r = np.sqrt(x**2 + y**2)
t = np.arctan2(y, x)
dr = 1 + np.cos(2*np.pi*r)
dt = 1/2
dx = dr*np.cos(t) - r*dt*np.sin(t)
dy = dr*np.sin(t) + r*dt*np.cos(t)
return np.array([dx, dy])
T = 0.5
h = 0.005
N = 2000
xlims = (-2, 2)
ylims = (-2, 2)
plot_streamlines(f, T, h, N, xlims, ylims)
In [13]:
def f(X, t):
x = X[0]
y = X[1]
r = np.sqrt(x**2 + y**2)
t = np.arctan2(y, x)
dr = 1.1 + np.cos(2*np.pi*r)
dt = 1/2
dx = dr*np.cos(t) - r*dt*np.sin(t)
dy = dr*np.sin(t) + r*dt*np.cos(t)
return np.array([dx, dy])
T = 0.5
h = 0.005
N = 2000
xlims = (-2, 2)
ylims = (-2, 2)
plot_streamlines(f, T, h, N, xlims, ylims)
In [14]:
def f(X, t):
x = X[0]
y = X[1]
mu = 1
dx = y
dy = -mu*(x**2 - 1)*y - x
return np.array([dx, dy])
T = 5
h = 0.005
N = 200
lims = 6
xlims = (-lims, lims)
ylims = (-lims, lims)
plot_streamlines(f, T, h, N, xlims, ylims)
In [15]:
def f(X, t):
x = X[0]
y = X[1]
mu = 1
dx = y
dy = -mu*(x**2 - 1)*y - x
return np.array([dx, dy])
T = 5
h = 0.005
N = 5
lims = 6
xlims = (-lims, lims)
ylims = (-lims, lims)
plot_streamlines(f, T, h, N, xlims, ylims)
In [16]:
def f(X, t):
x = X[0]
y = X[1]
mu = 3
dx = y
dy = -mu*(x**2 - 1)*y - x
return np.array([dx, dy])
T = 5
h = 0.001
N = 500
lims = 6
xlims = (-lims, lims)
ylims = (-lims, lims)
plot_streamlines(f, T, h, N, xlims, ylims)
In [17]:
def f(X, t):
x = X[0]
y = X[1]
a = 0.08
b = 0.6
dx = -x + a*y + x**2*y
dy = b - a*y - x**2*y
return np.array([dx, dy])
T = 20
h = 0.01
N = 100
xlims = (0, 3)
ylims = (0, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [18]:
def f(X, t):
x = X[0]
y = X[1]
a = 0.13
b = 0.6
dx = -x + a*y + x**2*y
dy = b - a*y - x**2*y
return np.array([dx, dy])
T = 10
h = 0.01
N = 500
xlims = (0, 3)
ylims = (0, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [19]:
def f(X, t):
x = X[0]
y = X[1]
a = 0.4
b = 0.6
dx = -x + a*y + x**2*y
dy = b - a*y - x**2*y
return np.array([dx, dy])
T = 20
h = 0.02
N = 100
xlims = (0, 3)
ylims = (0, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [20]:
def f(X, t):
x = X[0]
y = X[1]
mu = 5
dx = mu*(y - (x**3/3 - x))
dy = -x/mu
return np.array([dx, dy])
T = 10
h = 0.01
N = 500
xlims = (-3, 3)
ylims = (-3, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [21]:
def f(X, t):
x = X[0]
y = X[1]
dx = x*y
dy = - x**2
return np.array([dx, dy])
T = 20
h = 0.02
N = 200
xlims = (-3, 3)
ylims = (-3, 3)
plot_streamlines(f, T, h, N, xlims, ylims)
In [22]:
def f(X, t):
x = X[0]
y = X[1]
mu = 1.0
dx = mu - x**2
dy = -y
return np.array([dx, dy])
T = 1
h = 0.001
N = 1000
w = 1.2
xlims = (-w, w)
ylims = (-w, w)
plot_streamlines(f, T, h, N, xlims, ylims)
In [23]:
def f(X, t):
x = X[0]
y = X[1]
mu = -1.9
dx = mu*x + y + np.sin(x)
dy = x - y
return np.array([dx, dy])
T = 1
h = 0.001
N = 1000
w = 1
xlims = (-w, w)
ylims = (-w, w)
plot_streamlines(f, T, h, N, xlims, ylims)
In [24]:
def f(X, t):
x = X[0]
y = X[1]
omega = 0.25
b = 0.5
mu = -0.2
r = np.sqrt(x**2 + y**2)
t = np.arctan2(y, x)
dr = mu*r + r**3 - r**5
dt = omega + b*r**2
dx = dr*np.cos(t) - r*dt*np.sin(t)
dy = dr*np.sin(t) + r*dt*np.cos(t)
return np.array([dx, dy])
T = 20
h = 0.02
N = 50
w = 1.5
xlims = (-w, w)
ylims = (-w, w)
plot_streamlines(f, T, h, N, xlims, ylims)
In [25]:
def f(r, mu):
return mu*r + r**3 - r**5
x = np.linspace(0, 1.1, 100)
plt.plot(x, f(x, 0.0), label = 0.0)
plt.plot(x, f(x, -0.25), label = -0.25)
plt.plot(x, f(x, -0.5), label = -0.5)
plt.legend(loc = 'best')
Out[25]:
In [ ]:
##########################################################
#### ####
#### Running this cell will create a series of images ####
#### which illustrate the transition between a ####
#### homoclinic and hetroclinic orbit ####
#### ####
##########################################################
def f(X, t):
x = X[0]
y = X[1]
dx = y
dy = mu*y + x - x**2 + x*y
return np.array([dx, dy])
T = 10
h = 0.01
N = 200
w = 2.0
xlims = (-w, w)
ylims = (-w, w)
#plot_streamlines(f, T, h, N, xlims, ylims)
fig = plt.figure(figsize = (12, 12))
x = np.random.random(size = N) * (xlims[1] - xlims[0]) + xlims[0]
y = np.random.random(size = N) * (ylims[1] - ylims[0]) + ylims[0]
for j, mu in enumerate(np.linspace(-0.95, -0.75, 40)):
plot_streamlines(f, T, h, N, xlims, ylims)
for i in range(N):
# Get random initial conditions by scaling a random number between 0 and 1
line = streamline([x[i], y[i]], T, f)
plt.plot(line[0,:], line[1,:], lw = 1, c = '#A60628', alpha = 0.25)
# r = np.linspace(0, 1.25, 100)
# plt.plot(r, mu*r + r**3 - r**5)
# plt.plot(r, np.zeros(len(r)), c = 'k', alpha = 0.8, lw = 1)
plt.xlim(*xlims)
plt.ylim(*ylims)
plt.savefig('homoclinic_mu=%.04f.png' % (mu + 1))
plt.clf()
In [27]:
def f(X, t):
x = X[0]
y = X[1]
mu = -0.3
dx = y
dy = -mu*y - np.sin(x)
return np.array([dx, dy])
T = 10
h = 0.01
N = 200
w = 0.5
xlims = (-w, w)
ylims = (-w, w)
plot_streamlines(f, T, h, N, xlims, ylims)
In [28]:
def f(X, t):
x = X[0]
y = X[1]
mu = -0.2
dx = mu*x - y + x*y**2
dy = x + mu*y + y**3
return np.array([dx, dy])
T = 500
h = 0.05
N = 100
w = 1.0
xlims = (-w, w)
ylims = (-w, w)
fig = plt.figure(figsize = (12, 12))
plot_streamlines(f, T, h, N, xlims, ylims)
In [29]:
def f(X, t):
x = X[0]
y = X[1]
mu = 1.2
r = np.sqrt(x**2 + y**2)
t = np.arctan2(y, x)
dr = r*(1-r**2)
dt = mu - np.sin(t)
dx = dr*np.cos(t) - r*dt*np.sin(t)
dy = dr*np.sin(t) + r*dt*np.cos(t)
return np.array([dx, dy])
T = 2
h =0.02
N = 1000
w = 1.75
xlims = (-w, w)
ylims = (-w, w)
plot_streamlines(f, T, h, N, xlims, ylims)
In [30]:
def f(X, t):
x = X[0]
y = X[1]
mu = -0.92
dx = y
dy = mu*y + x - x**2 + x*y
return np.array([dx, dy])
T = 10
h = 0.005
N = 500
w = 2.0
xlims = (-w, w)
ylims = (-w, w)
plot_streamlines(f, T, h, N, xlims, ylims)
In [ ]: