interactive_predator_prey.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "%matplotlib inline\n",
    "from ipywidgets import interactive\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sympy import symbols,Matrix,solve,diff,init_printing\n",
    "from scipy.integrate import odeint\n",
    "init_printing()\n",
    "\n",
    "#Always use a large enough font size!\n",
    "plt.rcParams['axes.labelsize'] = 16\n",
    "plt.rcParams['xtick.labelsize'] = 14\n",
    "plt.rcParams['ytick.labelsize'] = 14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predator prey model:\n",
    "\n",
    "x: Prey (Rabbits) <br>\n",
    "y: Predator (Foxes)\n",
    "\n",
    "$F(x,y) = (\\underbrace{a x}_{\\text{Rabbits reproduce}} - \\underbrace{b x y}_{\\text{Rabbits get eaten by fox}}, \\underbrace{-c y}_{\\text{Foxes die out without rabbits}}+\\underbrace{d x y}_{\\text{Foxes thrive by eating rabbits}})$\n",
    "\n",
    "$a,b,c,d > 0$ are parameters of the model. In simple terms, $a$ describes the reproduction rate of rabbits in the absence of foxes, $b$ describes the reduction of rabbit population if a rabbit meets a fox, $c$ describes the decay rate of the fox population if there are not rabbits (starving foxes), and finally $d$ describes the increase of fox population if a fox meets a rabbit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#symbols is a function from the sympy package, defines symbolic variables\n",
    "x,y,a,b,c,d = symbols('x y a b c d')\n",
    "\n",
    "def F(x,y,a,b,c,d):\n",
    "    return [a*x - b*x*y, -c *y+d*x*y]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get fixpoints\n",
    "Solve $F(x,y)=(0,0)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAAAmCAYAAADeK5lgAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAIa0lEQVR4Ae2c25HVRhCGz7o2AC4RGDLgEgGQAV5HAGRgiid4oyADQwRcMgBHsEAGxhEYyAD/n3ZapSPNpUcaHS+Fukoraaanu6fnn56LZs/R48eP7+52uze6hvTgyZMnL4YJ2/M8D8iPf6rkO93fzpNwfkqpDhdkzV+6bun5W8myUPf7I76LRwF0v4vht1Hm9rrQA8Hplzy+Fc8zqfs3qKRBPyjt00ITisWl44qY0P1ZF/rpHHeV/lz3CSn9mhIJUtf1XATeUID4v+r911+GidtzOw/IwfTwG7pnO7Pyr+j6W7ynuj/n0vN1XQBhVZKu21LwUddDPXOhm8j8TM9EtQkpnY6AbS8nmc6EDXROR9WwqWGIBjTMLUe5d+J5qzLD4fdEaaSvRgFURCzARpQz4vmz0pJRTHlMvS7oztSsmo6rS2wFPB4gCrzINRxClE9UYXh7yruR0i/a84p3OgXAGc/diX7DDpAygQj+j8q/15UEaKzwFuliXlmQpgag9xPp9oCUEElEq260hKza5E73sJBsZ0ilExSjbADaa/HSwapoA12Vu1zMRBDmZtneHxqYRi42sEtrBVNGN0DcKf+9Uxx1ZdEBUN20gc7tqjKjnE+EowGKUU68Bkq77ylQPsPc2jScy6Hrjq5uxSz9RTCJh/IA9IEuN1WDToqqUO22pAHjObDtkapRM1wyn5qsblUPIsiXBi6JipB8gA5Y+rZUmk0LPoRCd5Q2BmVMHouR8V5cjK9Pq1pIyIg/VLJb3fQSztcD2w/JPaYDmErDuXu9bGUTnu0JQMa2CcSQ+1Rp0QjYcbT5A9hfSg9telnXK11EOuwhjUWOhwBvt5JVOc8CZOfeHJZAHHpT94djS5SG0yA2F6/qwnBPL6GMmzx6xIPDvuk+XpW59cxhlD6GVva8ruq5ed3n2HSoMqovHYYIn+1wyvdvDouZ3vdI9xjgcDQbm7a5CA+fffrQrffFJHkuPeJjg5MIgs2HpBsok96fCnDBwUS7rv7hPXvzzumIZJNwKwczlhNa+7CqZ4YF3if8SptFM/Sg26LvLJ0zCvWT8Bllf/QiRDoivYu8oDtRw8eGK+YFse+Dp0q/rTKtok2VnmArNrfS73Emq82fMcrhm67e8rcLeEXQSRBzuZQzcXRslWX8rZb9c/Rgw4muQxEAt3ofSud50WP1vuIxqAg6CWHYYMzeI2cUubRXaMbLAj3YjO2HJFuBHlLnedBloHO1twd0TBBjzjQFzOFS1GJ4m6sHm109L2W8N10dw/TEor5XzA/Lp/obBlzt7QEdguY6k/2fQ1BMDzYbGNa2weXstY04B/Jj7TAx63iSMk0g0hiSh7k5IFp0Yt9uKc3VQ8g/FBisvjE/TeqvyPB9kphIEO9RImt28or6Xf72gC5aOUKqLvJiiizNxvqoDE/iAj2pzuJRW8tjHcPqnS2vOjUHUlbhKPP/1u8ZXnFoyplM1mNDmPX8yQJkVH/v6xw92Gxg8OrZ+JZ5IDb3n0j0gI5oFQMWwvjYG9uJvq70T0QpmCA9p4B7xpD/69YzEAPwi5F2oV2mzuppnc3SV7k3snkN28wPWdke0H2ShJsxKao8G8ZfdGcvr6PgEPbH7oWkXUj7qjufsqpJ5Vx6RoI74I/S9l6X2jUQZhF1SccaiEs/NrQ5raQyRzZZUDI/ZCV4QMfpg2sZKTQux2DsdAInSfkXNcDakZ7pAd3EXs9zG6aoJ6izGxvK2QOSjezaBTnUkcMOq1IrmxsbaW2KD4pUXEiokgyTO13XeB5LVBqKsqcLKCM+Tl/0EXEsp/Tu1RN0dT1PZYpzyqV2DezOTUMGbPlH2UMHp+NSh9d6n/i2lc2SA1hs6sKINbfTWKSzs3gSmyZPpKM0H88nlU+LTeZwNMrVG5ISfBmcdKn54N/CLhwem9/6LA5c8g+dnKgOIHKRerHNtIUuvtrQYYodVDwpYvqV/Q+yYUEX6GQYcyoOSBqihzJczyqLE1vs22X1BRuxNXZAYVK2oV3MVzlxQz0XkWQwNYCiQGho85mWs+lTDuDGl7oTnSejYIrZBbpQmJMeS44r3ZezOOu2NmFjTVRuZRf/GQUZYM7e5v3too/8lRoVWtm8kw6zNwpwp/nIcIPWDbrgAA5HcjK3mlRudcAF27CxuFViFWhll+QAEBquxSEDGjEJglY2Bx+UAB7Y4jfZUg3a47ioeKoU0JirgyeuvZzauDHKCqccTMqr5r+ymSkLZfAt0w/AxnD1VNchCNAwj7Rgwvzsld77g7kFIwAt5d0d3R3pCoq37DMPMMQyr7Pen/VL4GMuOPwtEVsAJSNdVmhFpvQz/wTgAJ9fJCCg3NP1Rs+keYjT41UdZAOdx61OHjUUQywNZ8BJlgwNTmSM/ZYIkSM1n0vKnJFhnaP/6a+gF93F7S3xdtFRd29U7EzcQDejpXJF1ABs17B6JoLkCGASFcerbICwepQLhtnQOAb4BeVfDjy52yNlEhmraANdlbvczACPDd4c8alwD1wCII3NsOZeCeYUOPImAJcNFv1Oc+XFR5Rjb64qyiFzA13OszPz1BBd9NI9+p/vSgdcXGNwAcSd8vfASNpKBMDH4GJrDEraIPsoR5QzXvjdtIHO7apqRn6bjvkaDZSi8YqvG+5gVrnib4mkhFam9zZIJx0B4LPt9C0mJ/DQWe7puS8b402lHYcMKjg8zYrS8VwjJWNLj3hA/rNPTKwE+4k6rCGPSNIDUmlM3JkHWoThEEX10KXyNYR8Pt3Z1wQWNnzvzbU9PBzuKNomHjbqJ9H+6Pv3IdbEslFTD8jxAIsFgzVsJ1/vRBXmfQxvTNpf6SK6sMAgjZ9p2CujtKY0sIGIhT0flZYEXOC/obt1jFn2/Ad1KQnG863vtwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[ \\left( 0, \\  0\\right), \\  \\left( \\frac{c}{d}, \\  \\frac{a}{b}\\right)\\right]$"
      ],
      "text/plain": [
       "⎡        ⎛c  a⎞⎤\n",
       "⎢(0, 0), ⎜─, ─⎟⎥\n",
       "⎣        ⎝d  b⎠⎦"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#solve is a function from the sympy package\n",
    "fps=solve(F(x,y,a,b,c,d),[x,y]) #solves a*x - b*x*y == 0 && -c *y+d*x*y == 0, solving for x and y\n",
    "fps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fix(i,a1,b1,c1,d1):\n",
    "    vl=list(zip([a,b,c,d],[a1,b1,c1,d1]))\n",
    "    return [z.subs(vl) for z in fps[i]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobian Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "i,j = symbols('i j')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a vector with all state variables, easier to generalize to more dimensions\n",
    "v=[x,y]\n",
    "def jij(i,j):\n",
    "    #diff is a function from the sympy package\n",
    "    return diff(F(x,y,a,b,c,d)[i],v[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKIAAAAzCAYAAAAUyNP9AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHJ0lEQVR4Ae2d7XUUNxSGlxwXQEwHTgchVJClAwwVEDqAk1/2Pw50AK6AmA4IFWDoIHQQZztI3mdWUrTyfGh3Rtr50D1H6HNGV++9utJodc29i4uL+6saury83NQUl6KCQC8EpFe1+vaD3vpV4Z8gvOnVW3m4INCMwLWqQn27vieL+JcqXkhT/2x+ttQUBNIgIL37TW8+xyIWKggcHYGTo3OwBwNmf4Fp/0XhVvmf9nh8UU2nhtWkLKLA3Sg8lkZ9VyhbiZapNTWsJqWIHu4/K/3Jy5dkMwKTwGpyiqiZvjaYF4vYrHxVzZSwmpwiCuFqaWbp6ZBDqZ4QVlEfKxL6maTK2SJ7s78VPio8UflbxbkJi/hNfb80HT9S/EF5eFop5sCUIwE+ZL4q/15xRabuSvG5LZt53IoVYzeYgBcEls8VkPczBeiL2lTYbrNp/u20iGKCwXDo/UppAsr3TuGN0rWn5KpLQqY/9jwA9d7wAnDXSlMG/W7K2UOGB/NPVfaERnOnSKyAATm+NZh9Uf5KYa38K8XIN8RQRcNTqyKawXBcggJiDS2R/q6y3MsjkwL61fZtYvjAQqOkgAmxhN9Wqf//oWwpe8tWrIBEeGEJfUWrcFSZXUVOg3pl09BJx2th8r4YtozZ5gyy01zrOWbUZ95hH4yIz/Xct4Z2KBLLMoD5xPsfKDA57LNYv9d+I6XhOywLmowjOwB2XVgx0Bv14xuYhypz+Kou2xamSxER5o4FMQCxDHYen6gtCsPghqI7E0B9UAaxl6kU1JShnG4CqQxrSdnOeJQfJZmx9MGuFSsGrT7spLUY1E1eW5c0blyaxSRCI4QKB7MM4hgCZQLYpRc2IDtrfX4oczO7arW1hhyIh+Cb6tlFsVhVAxcuKC7ydiudylgNKUtOJxE9+Kab5pXJJyEm2fgj8LAN1anI9WVAYmJwaWPjdYgQXDtTvqT9oYXCYRBiZfLs//lYYRK/UGCiumeU58OPj5bk1GYRESwMItSKxBSKxxJ3sy1ZPQ4YN8XJImYrvzNbAsg/xINbgk2FD+ZK9cx2QmjdTfNZRl1YWUxuhQ9W79ZHwWAWrj5+k0HTXRaRJY5zt5eK+Rj4oIBlYRZRxjFOTnquzuCHM0LA49gmVEL4YRbTjo8tzj3t5Qh/+VbxrKkLK7AAOxRyJaxYVT4pIFOO61BQt0zTJiUt4j6iAEUh14r7bP5TymGx75ZM5nkfEaVTYEZXpDSWk8FO4thmy/Xy/m3cI04YCvaxbCEs8UsBv8JkW2ZsxyWOR6Brjxj/pvG0ZH94KsVjD8ve8J3SS9objkcSe3AyO0Uslm8P6Y+o6RyX5hHBW1iJRaAoYixSpV1SBIoiJoW3vDwWgaKIsUiVdkkRKIqYFN7y8lgEiiLGIlXaJUXgaMc3Ombh8gSHzWcKXFzg9kehkSOQSm7HVETuBT7UwP5VvKRbMSNXtXb2JK8kcjvq0qxBrc2wyy8f7fI/uBaMFayX3sHv8R9MIbejKqIGx5WyYzhh+bjOPc2lD8KQNLjcjq2IWMRiDYdUkTzvGlxu2faIMudnwoh7gdye5rIqCsgHi7uepTb2ylZxjhcwwsNe7FVutVHA6y6rz43665QbzKmdlR3ZRwpczOXZZwpQq6N+FosoJplB3BH0nfQBGfItYnGOFyDC60yBP6CK8KzzO5d6LWZKpif1HSs3mOnlqJ9cETUYZgq+JXVO+s7TTu2wjtZHgj3Ijg+F8pT5SqvsbIlThI/CxL9DiZNYttOFWLkhAbXt7aifY2lmFsc46c/COd4I8OA/KqDn8RlhSXNbFqUR9o/ETWSew4KFdEqB6uvOaTEE1h03fC5WbjzX21E/hyIyk3csmQaPlQRsN8NVxh5opRgwqXdOUSrDWlK28x7lR0dmHH18Yyq8zHuix6f2dYq2Ujk31lnq9/2DWVFyg0G9O9y38uzORKJdGyVdmsUgykNwCmeYgVEGUKdYi3WOb8HLwJYnauGjTW4Vc3rWGhK3reB95p2NA8hhEel8x89YefZ71SwSg8xYlgjbBktp00pWtIj9oTDAwZ0Bb6pRB/+obq1QN3mDloNl6+SwIzf1xF6eb4BejvqpLSKAAhzKVZGARPFYam+2JXec9HcGD/hqRwitqnl8dhFbkjv7NuHAng2hJyf1tY/crHzwg2b12+HRyM9+hDbyftJYM1wFoO7jpL9o53gJDkd3rAuKxxEOhIBfqwwFyUWxcsPQ9HbUH72DvREIS1KfD4BcwhtdP8Lt0I+VLGMRf+NzsEfpFIpz/LAqgBXNaUkP4j7H0rwPY8xe/+/pXClfnOP3QTBoq4md8+Mm6D0+OzZFLM7x8bKbVctRKaJmrzt7mhXKZTCdCCQ9vunsvTQoCBgErEXkJyCOCBwpP/oNrmO2JCaDQKhnYrzSO2sR+UDw/zNnzrAKFQRSIMCvMHd07T8hFPjf8NlIIgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a - b y & - b x\\\\d y & - c + d x\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡a - b⋅y    -b⋅x  ⎤\n",
       "⎢                 ⎥\n",
       "⎣  d⋅y    -c + d⋅x⎦"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Marix function from sympy:\n",
    "J = Matrix(len(v),len(v),jij)\n",
    "#J now is the Jacobian matrix in symbolic form\n",
    "J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #Explicit way\n",
    "# J = Matrix(np.empty([len(v),len(v)]))\n",
    "# for i in range(0,len(v)):\n",
    "#     for j in range(0,len(v)):\n",
    "#         J[i,j]=diff(F(x,y,a,b,c,d)[i],v[j])\n",
    "\n",
    "# #J now is the Jacobian matrix in symbolic form\n",
    "# J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stability of first fixpoint\n",
    "Calculate Eigenvalues of J at first fixpoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADwAAAAVCAYAAAAaX42MAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACD0lEQVRYCe2Y7VHCQBCGD8cCUDvQDmTsADoAS9AO4Cf8c7QDsQJHOgArEO1AOvCjBJ83uWNyySXA8BUz7Myyd3t7u/veXjbRWr/f7xpjbmFRZzAYfMTDavyCx8NXA/A90MYsTKoBMYwCfBHOo/BydbUHwNWtbYzsUOGqV/i4CCCd7Zx1dbcZ/AWP4Db6B2TpyHZi5Sn6hafovNdsLmAMm2x4gRuMBdggx4gmcgjLYSmIXFQY5dZjrKIY5COiA7c0dxR8hjGuYyCwchCBtRs0nqErDVibl8COyCsCa3XXSOk9yquwrnEdB0PPmuoyTzpNLa8+JYYO9xWWXJbmX4TsVyVV4bvkZvQnybkb5wHW6XhfXjYxd3Xc/rUlfnVbGms4inK1fha6yVxpC0ynnb4OcmxY9w5iYYQtGhTkmhs1AzhhmXx2pdbDH3U8AqlTq9oR2cBuujNJXN0OkZPxzP6yrkfQowxg60RVTAJqM7+Ep3Z3CzvXuXUbfpi/27VdC/UZdWOPyEd96NtTMgn+tYSxQDzBb/AZ/AzrFOVEugk28/cb4090Ir3CgqcdL2/nl5jKS+TyUP7eq9PajINNyyadOTWchHQG+wtYt2AvROzesoEzV3rZjQG7KwLvvLqBPApVGwEMUF0h90lXGHDfixsBDIgbQJfy+zp9wK5puedv/gWTNvyvcwrRJXf9z+4U7vwBabCuS2Ht+s8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[ a, \\  - c\\right]$"
      ],
      "text/plain": [
       "[a, -c]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f1l=J.subs(list(zip([x,y],fps[0]))).eigenvals()\n",
    "#.subs replaces x and y with the values in fps[0], hence x=0,y=0. \n",
    "#.eigenvals returns diccionary with eigenvalues and their multiplicity.\n",
    "#following line extracts only Eigenvalues\n",
    "[evs for evs,multiplicity in f1l.items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One Eigenvalue is positive, hence fixpoint is unstable. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stability of second fixpoint\n",
    "Calcualte Eigenvalues of J at second fixpoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJUAAAAZCAYAAAA40GzsAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD40lEQVRoBe2a7VXUQBSGFw4FAHaAHYhWIHQA2gF0wB5/wT+PdiBW4JEOxAoEOpAOVErweeIkm4Rs9Gwy+VhzzxkymZncj3fe3JmdsHF+fv55Nps9o6RyeHFxcZ/e9HHFvv7c9mG7ziZ+bdT1x+wLmLzGxllMOyvqzvNlvoWSPRx+uqKyWI+9Q/E+ft3FMjBCvWIyB5P5UH3HtxN8290cmoM4lmTNiVCLmRkbJoMjFVC+ofhWTrJAYFSYDIpUvJF74OhyfL3A8/+ujRGTQZEK+rhfGOyeoSd6jw6TwZCKN3KbSXs+ZakFdceKyWBIBZTuo94uIJ1qY8XEI4XeJbyRB1xPYziDXvdqktbzlB+UK8oR7e+5PhLaHes45YFyQ1unxxvYM3NHw8TAsBEFl0GQivii/boBuAP0e8DruVdySMf1C/dO2CVF0iRCXZDt8zxI4s24fuByTDn0vkOJhokxEFc0XDZjg4TzZgQDqBPHXNYNWKUPnb7tEkqS5E99rd/TlhEq6JdQV7QnhAptr7ja3pqgvzdMDAL7UXFZmqmC4a/4oAP/Ksc8ly0T1P2k4BvnJO5XKQljzAaPpAUfXMa20VMmrCTPE2fGGH0wUxX2dbTv0NaaoK8RJjqCDuekydxExaWOVL7FlUQwsL8JgXsy7mQ+oZx5T8kIl3v+lPbKz0S0N/IBG2aZwpkXOp2QdJnLufFnbLCZb2+tju7GmOhM8HHluUFFVFyiLX8EfheCT998M1ZB6PdbUWWWKgxc4QbdksdSXroEdEZ/RraasQ5tTbDTKyYGUhNra7hEI1U6EwRhtjFjuY8wQ+TFg73y0pTvb6Oe30upzw13kjHxJ/Ep+Gifvj4S+rM9IXWJ2kiCvT4x0f9WcckDEp1UwZhruJKdlgPsEfduiisnMhnd4E/QazbKiBxsugTdBNX5f/Nxkv2VVxCe0fefNlKXUL+4tvFvOZ1jEmIQ71ZxUW9eNvh/qltAarI+5/UtrWPDZcg3fof6Q5iYl9aXPtSwA92S4CPlG8W93SeK9pxQ264Zk+3zqKcT/Z0+xefLxw5pn0cUjXzn+c4xMSjsto5L0Ot2ZtYlqSSUIJqtTL1miSiHneiOKvhtlpWQTUm1NpgIOHgkpOpq+dOgKdes4IbdkmYFqqOTF8TTiFBGvGaYZJPYGamCRX8Jmno9eCxvFDOnhlzBb/1PP+G04eroMSmD0CmpmBAPHM1W2Ya97NAI7k+Io/Kb4Sq+rwkmhdA721MVrE43a4kAL0i3e6q1RHEKqhKBTpe/Sg+mxrVDYMuISFvp2Yu3+QNB7yeZEFiKANw5ozM9GtqlPv8NJvrozd9hSuwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[ - \\sqrt{- a c}, \\  \\sqrt{- a c}\\right]$"
      ],
      "text/plain": [
       "⎡   ______    ______⎤\n",
       "⎣-╲╱ -a⋅c , ╲╱ -a⋅c ⎦"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f2l=J.subs(list(zip([x,y],fps[1]))).eigenvals()\n",
    "[evs for evs,multiplicity in f2l.items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since a,b,c,d > 0, both Eigenvalues are purely imaginary. We cannot say anything about the stability of the fixpoint from this (more advanced techniques would have to be applied), but the imaginary Eigenvalues suggest oscillating behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Streamplot and numerical integration\n",
    "\n",
    "Indeed we find an oscillating behavior by looking at the streamplot as well as the numerical solution of the differential equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e3d759e913f249bdbdee4178911e7102",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='a', max=2.0), FloatSlider(value=1.0, description='b'…"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r=5\n",
    "x1 = np.linspace(0,r,101)\n",
    "X,Y = np.meshgrid(x1,x1)\n",
    "\n",
    "#initial condition for integration\n",
    "x0 = [1.5,0.8]\n",
    "\n",
    "#define function in vector like input form. \n",
    "def F2(vx,t,a,b,c,d):\n",
    "    x,y=vx\n",
    "    return F(x,y,a,b,c,d)\n",
    "\n",
    "t = np.linspace(0, 40, 201)\n",
    "    \n",
    "def plotF(a,b,c,d):\n",
    "    #Evaluate streamplot data\n",
    "    U, V = F(X,Y,a,b,c,d)\n",
    "    #Numberical integration\n",
    "    sol = odeint(F2, x0, t, args=(a,b,c,d))\n",
    "    \n",
    "    f, ax = plt.subplots(1,2,figsize=(12,6))\n",
    "    \n",
    "    ax[0].streamplot(X, Y, U, V,density=1.5,linewidth=0.5,color='black')\n",
    "    #hightlight stream lines originating from initial conditions:\n",
    "    ax[0].streamplot(X, Y, U, V,density=1.5,linewidth=1.5,color='g',start_points=[x0])\n",
    "    ax[0].plot(x0[0], x0[1], 'bo',color='g')\n",
    "    #plot fixpoint in black:\n",
    "    myfp=fix(1,a,b,c,d)\n",
    "    ax[0].plot(myfp[0], myfp[1], 'bo',color='black')\n",
    "    ax[0].axis((0, r, -0, r))\n",
    "    ax[0].set_xlabel(\"x\")\n",
    "    ax[0].set_ylabel(\"y\")\n",
    "    \n",
    "    \n",
    "    ax[1].plot(t, sol[:, 0],  label='x(t)')\n",
    "    ax[1].plot(t, sol[:, 1],  label='y(t)')\n",
    "    ax[1].legend(loc=1,fontsize=14)\n",
    "    ax[1].set_xlabel('t')\n",
    "    ax[1].set_ylabel('Population')\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "interactive_plot = interactive(plotF,  a=(0,2,0.1), b=(0,2, 0.1),c=(0,2,0.1),d=(0,2,0.1))    \n",
    "\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}