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)