{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "514dc19a",
   "metadata": {},
   "source": [
    "# Discrete dynamical systems: practice and discovery of different dynamical regimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65ef3c4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from math import *\n",
    "import numpy as np\n",
    "from numpy.linalg import eig\n",
    "from scipy.integrate import solve_ivp,odeint\n",
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3414f023",
   "metadata": {},
   "source": [
    "## Try out exercise 1.1 from Chapter 1 of the book by Kaplan & Glass\n",
    "\n",
    "We will simulate the discrete dynamical system defined by iterating the function $f(x)=Rx-\\frac{R}{2000}x^2$, for different values of $R$. We will observe different dynamical regimes as the number of iterates increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e74f409",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(i, x, R):\n",
    "    i_new = i+1\n",
    "    x_new = R*x-(R/2000)*x**2\n",
    "    return i_new, x_new\n",
    "\n",
    "def iterate_map(f, x0, R, N):\n",
    "    points = [(0, x0)]\n",
    "    i, x = 0, x0\n",
    "    for _ in range(N):\n",
    "        i, x = f(i, x, R)\n",
    "        points.append((i, x))\n",
    "    return points\n",
    "\n",
    "# \n",
    "x0= 200\n",
    "N = 30\n",
    "R = 1.5\n",
    "points = iterate_map(f, x0, R, N)\n",
    "# \n",
    "#x0= 200\n",
    "#N = 30\n",
    "#R = 3.2\n",
    "#points = iterate_map(f, x0, R, N)\n",
    "# \n",
    "#x0= 200\n",
    "#N = 50\n",
    "#R = 3.7\n",
    "#points = iterate_map(f, x0, R, N)\n",
    "\n",
    "#\n",
    "fig, ax = plt.subplots()\n",
    "line, = ax.plot([], [], 'bo-')\n",
    "ax.set_title(\"Iterates of function f\")\n",
    "ax.set_xlabel(\"i\")\n",
    "ax.set_ylabel(\"x\",rotation=0)\n",
    "ax.grid(True)\n",
    "\n",
    "xs_all, ys_all = zip(*points)\n",
    "ax.set_xlim(min(xs_all) - 5, max(xs_all) + 5)\n",
    "ax.set_ylim(min(ys_all) - 100, max(ys_all) + 100)\n",
    "\n",
    "ax.plot(xs_all,ys_all,'bo')    \n",
    "ax.plot(xs_all,ys_all,'b')    \n",
    "# \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea00c6c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return R*x-(R/2000)*x**2\n",
    "\n",
    "# \n",
    "R=1.5\n",
    "\n",
    "def cobweb_plot(f, x0, iterations, x_range=(0, 2000), resolution=500):\n",
    "    # \n",
    "    x = np.linspace(x_range[0], x_range[1], resolution)\n",
    "    y = f(x)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    ax.plot(x, y, 'k', lw=2, label='f(x)')\n",
    "    ax.plot(x, x, 'r--', label='y = x')\n",
    "\n",
    "    # \n",
    "    xn = x0\n",
    "    for _ in range(iterations):\n",
    "        x_next = f(xn)\n",
    "        # \n",
    "        ax.plot([xn, xn], [xn, x_next], 'b', lw=1)\n",
    "        # \n",
    "        ax.plot([xn, x_next], [x_next, x_next], 'b', lw=1)\n",
    "        xn = x_next\n",
    "\n",
    "    ax.set_title(f'Cobweb Plot for f(x), starting at x₀ = {x0}')\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$f(x)$')\n",
    "    ax.legend()\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "\n",
    "# \n",
    "cobweb_plot(f, x0=200, iterations=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd5d27a0",
   "metadata": {},
   "source": [
    "## Second example with a cubic function\n",
    "Try to modify the initial condition $x_0$ and the parameter value $R$ in order to unveil different dynamical regimes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a302d08",
   "metadata": {},
   "outputs": [],
   "source": [
    "def g(i, x, R):\n",
    "    i_new = i+1\n",
    "    x_new = -x**3/3+x+R\n",
    "    return i_new, x_new\n",
    "\n",
    "def iterate_map(g, x0, R, N):\n",
    "    points = [(0, x0)]\n",
    "    i, x = 0, x0\n",
    "    for _ in range(N):\n",
    "        i, x = g(i, x, R)\n",
    "        points.append((i, x))\n",
    "    return points\n",
    "\n",
    "# \n",
    "x0= 1\n",
    "R=-1.2\n",
    "N = 20\n",
    "points = iterate_map(g, x0, R, N)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots()\n",
    "line, = ax.plot([], [], 'bo-')\n",
    "ax.set_title(\"Iterates of function f\")\n",
    "ax.set_xlabel(\"i\")\n",
    "ax.set_ylabel(\"x\",rotation=0)\n",
    "ax.grid(True)\n",
    "\n",
    "xs_all, ys_all = zip(*points)\n",
    "ax.set_xlim(min(xs_all) - 5, max(xs_all) + 5)\n",
    "ax.set_ylim(min(ys_all) - 5, max(ys_all) + 5)\n",
    "\n",
    "ax.plot(xs_all,ys_all,'bo')    \n",
    "ax.plot(xs_all,ys_all,'b')    \n",
    "# Display\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "426fbf34",
   "metadata": {},
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    return -x**3/3+x+R\n",
    "\n",
    "# \n",
    "R=-1.2\n",
    "\n",
    "def cobweb_plot(g, x0, iterations, x_range=(-3, 3), resolution=500):\n",
    "    # \n",
    "    x = np.linspace(x_range[0], x_range[1], resolution)\n",
    "    y = g(x)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    ax.plot(x, y, 'k', lw=2, label='g(x)')\n",
    "    ax.plot(x, x, 'r--', label='y = x')\n",
    "\n",
    "    # \n",
    "    xn = x0\n",
    "    for _ in range(iterations):\n",
    "        x_next = g(xn)\n",
    "        # \n",
    "        ax.plot([xn, xn], [xn, x_next], 'b', lw=1)\n",
    "        # \n",
    "        ax.plot([xn, x_next], [x_next, x_next], 'b', lw=1)\n",
    "        xn = x_next\n",
    "\n",
    "    ax.set_title(f'Cobweb Plot for g(x), starting at x₀ = {x0}')\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$g(x)$')\n",
    "    ax.legend()\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "\n",
    "# \n",
    "cobweb_plot(g, x0=1, iterations=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a20d3f4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}