Plot_phase_portraits.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()\n",
    "import numpy as np\n",
    "%matplotlib widget\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8b5cbe907a644090a14188bff3cb74c9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.streamplot.StreamplotSet at 0x7f8f6dd053d0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y, X = np.mgrid[-7:7:500j, -7:7:500j]\n",
    "\n",
    "U = np.sin(X + Y)\n",
    "V = X - Y\n",
    "plt.streamplot(X, Y, U, V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e748d54c64fa42cfb370b1cb043baee1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "\n",
    "def lorenz(t, X):\n",
    "    # Unpack `X` into 3 coordinates\n",
    "    x, y, z = X\n",
    "    # Formulas for `\\dot x`, `\\dot y`, `\\dot z`\n",
    "    dx = 10 * (y - x)\n",
    "    dy = x * (28 - z) - y\n",
    "    dz = x * y - (8 / 3) * z\n",
    "    return (dx, dy, dz)\n",
    "\n",
    "for pt in [(0.01, 0, 0), (0, 0, 1)]:\n",
    "    r = solve_ivp(lorenz, (0, 200), pt, dense_output=True)\n",
    "    T = np.linspace(50, 80, 10000)\n",
    "    Y = r.sol(T)\n",
    "    ax.plot(Y[0], Y[1], Y[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a75bd517527548d68936bfdea0649086",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "\n",
    "def pendulum2(t, X):\n",
    "    # Unpack `X` into 4 coordinates\n",
    "    x, y, z, u = X\n",
    "    # Formulas for `\\dot x`, `\\dot y`, `\\dot z`, `\\dot u`\n",
    "    dx = y\n",
    "    dy = -x\n",
    "    dz = u\n",
    "    du = -z\n",
    "    return (dx, dy, dz, du)\n",
    "\n",
    "for pt in [(0.5, 0, 0, 0.5), (0, 0.5, 1, 0)]:\n",
    "    r = solve_ivp(pendulum2, (0, 200), pt, dense_output=True)\n",
    "    T = np.linspace(50, 80, 10000)\n",
    "    Y = r.sol(T)\n",
    "    ax.plot(Y[0], Y[1], Y[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2bafbc82af7144efa49b9c2d04560efc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "for pt in [(0.5, 0, 0, 0.5), (0, 0.5, 1, 0)]:\n",
    "    r = solve_ivp(pendulum2, (0, 200), pt, dense_output=True)\n",
    "    T = np.linspace(50, 80, 10000)\n",
    "    Y = r.sol(T)\n",
    "    ax.plot(Y[0], Y[2])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myenv",
   "language": "python",
   "name": "myenv"
  },
  "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}