In [1]:

    

# Basic libraries import
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import animation
from pathlib import Path
from datetime import datetime
import cv2
from tqdm import tqdm

# Plotting
%matplotlib inline

sns.set_context("paper")
sns.set_style("darkgrid")


    

In [2]:

    

def logistic(r: float, x):
    """
    Logistic function
    :param r: logistic coefficient
    :param x: input
    """
    return r * x * (1 - x)


    

In [3]:

    

# plost sample graph of the logistic function
r = 2
x = np.linspace(0, 1)
ax = sns.lineplot(x, logistic(r, x))
ax.set(xlabel='x', ylabel='logistic(x)')
plt.show()


    

In [4]:

    

def plot_logistic_map(r, x0, n):
    """
    Plot iteration over logistic map
    :param r: logistic coefficient
    :param x: initial input value
    :param n: number of iterations
    """
    # plot logistic function over fixed linespace
    x = np.linspace(0, 1)
    ax = sns.lineplot(x, logistic(r, x))
    
    # iteratively apply logistic from initial value
    # and plot directions
    # (x, x) -> (x, y)
    # (x, y) -> (y, y)
    x = x0
    for i in range(n):
        y = logistic(r, x)
        # Plot the two lines.
        ax.plot([x, x], [x, y], 'k', lw=1)
        ax.plot([x, y], [y, y], 'k', lw=1)
        # Plot the positions with increasing
        # opacity.
        ax.plot([x], [y], 'ok', ms=10, alpha=(i + 1) / n)
        x = y
        
    ax.set_title(f"r={r:.1f}, x_0={x0:.1f}")


    

In [5]:

    

plot_logistic_map(3.5, .1, 20)


    

In [6]:

    

def plot_bifurcation_diagram(x0, min_r, max_r, nb_r_vals, nb_iter, nb_last_iter):
    """
    Plot bifurcation diagram by simulating logistic map runs for different coefficient values.
    For each plot results for the last nb_last_iter results
    :param x0: initial input value
    :param min_r: min value for logistic coefficient
    :param max_r: max value for logistic coefficient
    :param nb_r_vals: number of values on which to run the simulation
    :param nb_iter: number of iterations for each logistic run
    :param nb_last_iter: number of last iterations to plot
    """
    # setup plot
    fig, ax = plt.subplots(1, 1, figsize=(4, 4))
    ax.set_xlim(min_r, max_r)
    ax.set_title("Bifurcation diagram")
    
    # range of logistic coefficient values over which we run the simulation
    r = np.linspace(min_r, max_r, nb_r_vals)
    # initial condition (for all simulations)
    x = x0 * np.ones(nb_r_vals)
    
    # run simulation
    for i in range(nb_iter):
        x = logistic(r, x)
        # plot values if last iterations
        if i >= (nb_iter - nb_last_iter):
            ax.plot(r, x, ',k', alpha=.25)


    

In [7]:

    

plot_bifurcation_diagram(1e-5, 1.5, 4, 10000, 1000, 100)


    

In [8]:

    

import scipy.integrate as spi


    

In [9]:

    

# model params
m = 1.   # particle's mass
k = 1.   # drag coefficient
g = 9.81 # gravity accelleration


    

In [10]:

    

p0 = (0, 0)  # initial position
v0 = (4, 10) # initial speed vector


    

In [11]:

    

# encode everything in single vector to use scipy solver
v0 = np.zeros(4)
v0[2] = 4.
v0[3] = 10.


    

In [12]:

    

def derive_velocity(v, t0, k):
    u, udot = v[:2], v[2:]
    # we compute the second derivative of p
    udotdot = -k / m * udot
    udotdot[1] -= g
    
    return np.r_[udot, udotdot]


    

In [13]:

    

def plot_system_simulation():
    fig, ax = plt.subplots(1, 1, figsize=(8, 4))
    
    # simulate system on 30 linearly spaced times between t=0 and t=3.
    t = np.linspace(0., 3., 30)
    
    # We simulate the system for different values of k.
    for k in np.linspace(0., 1., 5):
        # We simulate the system and evaluate $v$ on the
        # given times.
        v = spi.odeint(derive_velocity, v0, t, args=(k,))
        # We plot the particle's trajectory.
        ax.plot(v[:, 0], v[:, 1], 'o-', mew=1, ms=8, mec='w', label=f'k={k:.1f}')
    ax.legend()
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 6)


    

In [14]:

    

plot_system_simulation()