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)