04d_Newton_arclength.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)\n",
    "\n",
    "Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.\n",
    "Held at Technische Universität Bergakademie Freiberg.\n",
    "\n",
    "Comments to:\n",
    "\n",
    "*Prof. Dr. Thomas Nagel  \n",
    "Chair of Soil Mechanics and Foundation Engineering  \n",
    "Geotechnical Institute  \n",
    "Technische Universität Bergakademie Freiberg.*\n",
    "\n",
    "https://tu-freiberg.de/en/soilmechanics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Softening behavior and its consequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example problem\n",
    "\n",
    "We consider a minimal system to describe buckling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![minimal_buckling_model.png](https://raw.githubusercontent.com/nagelt/Numerical_Methods_Introduction/master/Images/minimal_buckling_model.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use our intuition for two limit cases. For a tiny mass and a stiff spring we expect the bar to stand upright, whereas for a huge mass and a soft spring, we expect it hanging down (either to the left or right). \n",
    "From these limits we can conclude that a stability loss of the upright position occurs. This stability loss, in terms of mechanics, corresponds to passing a state of zero stiffness. Thinking scalarly of F=ku we anticipate that solving for u=F/k fails when $k=0$. Indeed applying Newton's method the iterations do not converge (compare to exercise 04a starting Newton iterations on a point with a horizonal tangent).\n",
    "However non-converging iterations are neither a reliable nor an efficient indicator. So we go for an advanced method. In addition we obtain a glance into the exciting world of nonlinear systems.\n",
    "It goes without saying that you won't find nonlinear effects in a _linear world_, so make sure to enable nonlinear geometry in your FEM-code, if you want to analyse a structure for buckling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np #numerical methods\n",
    "from scipy.optimize import root # Newton's method outsourced\n",
    "import matplotlib.pyplot as plt # plotting\n",
    "\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "#Some plot settings\n",
    "import plot_functions.plot_settings\n",
    "%run plot_functions/Newton_arclength_plots.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principle of Minimum Potential Energy\n",
    "\n",
    "Since nature is lazy, all stable equilibria minimize potential energy (statics, conservative system). Note that maxima describe instable equilibria.\n",
    "For our system we have the potential energy of the rotational spring and of a mass in a gravitational field\n",
    "\\begin{align}\n",
    "\\Pi_\\mathrm{spring} &= \\frac{1}{2}c(\\varphi-\\varphi_0)^2, \\\\\n",
    "\\Pi_\\mathrm{grav} &= mgl\\sin \\varphi.\n",
    "\\end{align}\n",
    "The necessary condition for an equilibrium\n",
    "\\begin{equation}\n",
    "\\frac{\\partial \\Pi}{\\partial \\varphi}=0,\n",
    "\\end{equation}\n",
    "corresponds to an equilibrium of forces.\n",
    "From the second derivative we may identify the equilibrium as minimum $\\left(\\frac{\\partial^2 \\Pi}{\\partial \\varphi^2}>0\\right)$ or maximum of potential energy $\\left(\\frac{\\partial^2 \\Pi}{\\partial \\varphi^2}<0\\right)$, or we have to check higher derivatives $\\left(\\frac{\\partial^2 \\Pi}{\\partial \\varphi^2}=0\\right)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the spring relaxed at standing upright ($\\varphi_0=\\pi/2$) we define functions for potential energy and its derivatives. Note that we make potential energy nondimensional, thus leaving only one parameter\n",
    "\\begin{equation}\n",
    "\\frac{\\Pi}{\\Pi_\\mathrm{ref}} = \\frac{1}{2}(\\varphi-\\varphi_0)^2 + \\mu \\sin \\varphi,\n",
    "\\end{equation}\n",
    "with $\\mu=\\frac{mgl}{c}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi0 = 0.5*np.pi   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Epot = lambda phi, mu: 0.5*(phi-phi0)**2 + mu*np.sin(phi) # energy\n",
    "dEpot_dphi = lambda phi, mu: phi - phi0 + mu*np.cos(phi)  # force\n",
    "d2Epot_dphi2 = lambda phi, mu: 1 - mu*np.sin(phi)         # stability\n",
    "d2Epot_dphi_dmu = lambda phi, mu: np.cos(phi)             # for tangent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "25c652be23b24ce4ab4987ff66292ef1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='$\\\\mu$', max=2.0), Output()), _dom_classes=('widget-…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_energy_landscape()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe how the energy landscape changes from single-well ($\\mu<1$) to a double well ($\\mu>1$). Basically this figure already explains everything, for a small load there is one equilibrium (upright) and for a huge load the upright position gets instable (maximum) and there occur two stable equilibria (minima), one to the left and the other to the right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution with the pseudo arc-length method\n",
    "\n",
    "The equation (equilibrium) to solve is transcendental\n",
    "\\begin{equation}\n",
    "f(\\varphi) = \\varphi-\\varphi_0 + \\mu \\cos \\varphi = 0,\n",
    "\\end{equation}\n",
    "so we resort to numerical methods.\n",
    "Newtons method, however, runs into trouble for $\\mu \\to 1$, because then the tangent \n",
    "\\begin{equation}\n",
    "\\frac{\\mathrm{d}f}{\\mathrm{d}\\varphi} = 1 - \\mu \\sin \\varphi\n",
    "\\end{equation}\n",
    "at the equilibrium ($\\varphi=\\varphi_0=\\frac{\\pi}{2}$) approaches zero.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it is time for the clever move, instead of solving for $\\varphi$ with the fixed parameter $\\mu$, we think of a solution curve in $\\varphi$-$\\mu$ space. Starting from a known point of this curve, e.g. the trivial unloaded case $\\mu=0$, we progress pointwise by imposing a fixed distance (called pseudo arc-length because coordinate $\\varphi$ and parameter $\\mu$ do not measure a spatial distance) from one point to another (index _p_ denotes previous point)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_squared = lambda phi,mu, phi_p,mu_p: 0.1*(phi-phi_p)**2 + (mu-mu_p)**2   # pseudo arc-length squared\n",
    "d_distance_squared_dphi = lambda phi,mu, phi_p,mu_p: 2*0.1*(phi-phi_p)  # for tangent\n",
    "d_distance_squared_dmu = lambda phi,mu, phi_p,mu_p: 2*(mu-mu_p)   # for tangent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the weigthing for the arc-length required some divine inspiration, we are not teaching here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The extended system (two equations for two unknowns), we have to solve, reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sys(z, zp, ds):\n",
    "    phi = z[0]\n",
    "    mu = z[1]\n",
    "    phi_p = zp[0]\n",
    "    mu_p = zp[1]\n",
    "    \n",
    "    r0 = dEpot_dphi(phi, mu)   # equilibrium\n",
    "    r1 = distance_squared(phi, mu, phi_p, mu_p) - ds**2 # distance constraint\n",
    "    \n",
    "    res = np.array([r0, r1])\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And its tangent reads "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsys(z, zp, ds):\n",
    "    phi = z[0]\n",
    "    mu = z[1]\n",
    "    phi_p = zp[0]\n",
    "    mu_p = zp[1]\n",
    "    \n",
    "    J00 = d2Epot_dphi2(phi, mu)\n",
    "    J01 = d2Epot_dphi_dmu(phi, mu)\n",
    "    J10 = d_distance_squared_dphi(phi, mu, phi_p, mu_p)\n",
    "    J11 = d_distance_squared_dmu(phi, mu, phi_p, mu_p)\n",
    "    \n",
    "    J = np.matrix([[J00, J01], [J10, J11]])\n",
    "    return J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from the unloaded initial state we increase the load (direction guided by start value of the iterations). In addition we observe the stability of the calculated states (from nonlinear systems theory it is known that a change of stability typically indicates branching of solutions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvYAAAH3CAYAAADQVRsYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2MElEQVR4nO3deZheZX3/8fc3k4QwArIFRCEJIEsgRJAHRYgtwQ2pdaOKMVrR4Li0adVWUUeL2/hDULGtG1PDz+U3jistaitCIYhBXEYIIRDAQkhQlEWQbQhJJt/fH+cZmUwny+SZec7Mmffruua689znPud8n1y5Tj65c59zIjORJEmSNL5NKrsASZIkSY0z2EuSJEkVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAGTyy5gLDvllFPykksuKbsMSZIkVV80egBn7LfivvvuK7sESZIkabsY7CVJkqQKMNhLkiRJFWCwlyRJkirAYC9JkiRVgMFekiRJqgCDvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJkqQKMNhLkiRJFWCwlyRJkirAYC9JkiRVgMFekiRJqgCDvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJI6K7G+bMgZaWou3uLrsiSZpYJpddgCRp/OvuhvZ2WLIE5s2DZctg0aJi24IF5dYmSROFM/aSpIZ1dBShfv58mDKlaJcsKfolSc0RmVl2DWNWrVbLnp6essuQpDGvpQXWrStCfb8NG2DaNOjrK68uSRpHotEDOGMvSWrY7NnF8puBli0r+iVJzWGwlyQ1rL29WFO/dGkxU790afG5vb3syiRp4vDmWUlSw/pvkF28GFatKmbqOzq8cVaSmsk19lvhGntJkiQ1iWvsJUmSJBnsJUmSpEow2EuSJEkVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAEGe0mSJKkCDPaSJElSBRjsJUmSpAow2EuSJEkVYLCXJEmSKsBgL0mSJFVAqcE+IvaLiEsiIsusQ5IkSRrvJpd14oh4JfBpYMMw9zsJuAhYO8TmZwBnZuaS+tgPAW8C7h807qrM/LvhVSxJkiSNXWXO2J8FvAC4egf2/V5mHj3wB3gbsI4i9A/0T4PHGuolaeR1d8OcOdDSUrTd3WVXJEkTS2kz9sCJmbkxIoa73yrggiH6zwS+k5kPNFyZJGlYuruhvR2WLIF582DZMli0qNi2YEG5tUnSRFHajH1mbtzB/e7OzGsG9kXErsDpQOdI1CZJGp6OjiLUz58PU6YU7ZIlRb8kqTmq8lScBcDazPzJENtOiYjLI+KGiLg2Ij4SEa3NLlCSqmzVqmKmfqB584p+SVJzVCXYn8nQs/W9wKPAaZl5FLAIeB1weURMaWJ9klRps2cXy28GWras6JckNce4D/YRcRRwFPDVwdsy89zMXJSZf6x/vo7ipt3jgVdv4XhtEdETET333nvv6BUuSRXS3l6sqV+6FDZsKNpFi4p+SVJzlHnz7Eh5M/DdzBz8SMst+Xm9PR7oGrwxMzupz/7XajWfry9J26H/BtnFi4vlN7NnF+vrvXFWkppnXAf7iNgJWAi8fAvbp2fm4Gn3vnrbMoqlSdKEs2CBQV6SyjTml+JExF4RMXULm08D7t7CTbMAayJicIA/tt5eOyIFSpIkSWPAmA72EXEg8Fvg4i0MORP4t60cYmfgw/3hPiJmAucAtwBfH8FSJUmSpFKVthQnIs6jePPsjPrn5fVNz8rM9fVfPwbcD9w1xP4HU6yT/6utnGYh8FpgeT3ctwKXAB/MzN4R+BqSJEnSmBCZ3h+6JbVaLXt6esouQ5IkSdUXjR5gTC/FkSRJkrR9DPaSJElSBRjsJUmSpAow2EuSJEkVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAEGe0mSJKkCDPaSJElSBRjsJUmSpAow2EuSRkR3N8yZAy0tRdvdXXZFkjSxTC67AEnS+NfdDe3tsGQJzJsHy5bBokXFtgULyq1NkiYKZ+wlSQ3r6ChC/fz5MGVK0S5ZUvRLkpojMrPsGsasWq2WPT09ZZchSWNeSwusW1eE+n4bNsC0adDXV15dkjSORKMHcMZektSw2bOL5TcDLVtW9EuSmsNgL0lqWHt7saZ+6dJipn7p0uJze3vZlUnSxOHNs5KkhvXfILt4MaxaVczUd3R446wkNZNr7LfCNfaSJElqEtfYS5IkSTLYS5IkSZVgsJckSZIqwGAvSZIkVYDBXpIkSaoAg70kSZJUAQZ7SZIkqQIM9pIkSVIFGOwlSZKkCjDYS5IkSRVgsJckSZIqwGAvSZIkVYDBXpIkSaoAg70kSZJUAQZ7SdKI6O6GOXOgpaVou7vLrkiSJpbJZRcgSRr/uruhvR2WLIF582DZMli0qNi2YEG5tUnSROGMvSSpYR0dRaifPx+mTCnaJUuKfklSc0Rmll3DmFWr1bKnp6fsMiRpzGtpgXXrilDfb8MGmDYN+vrKq0uSxpFo9ADO2EuSGjZ7drH8ZqBly4p+SVJzlBrsI2K/iLgkIvxvA0kax9rbizX1S5cWM/VLlxaf29vLrkySJo7Sbp6NiFcCnwY27MC+VwL7AOsHbfp0Zn510NgXAh8DdgamAF8BPpGZm3agbEnSEPpvkF28GFatKmbqOzq8cVaSmqnMp+KcBbwAaAeevgP7n5qZd2xtQETMA34AvDoz/yMiDgB+DuwGvG8HzilJ2oIFCwzyklSmMpfinJiZvx7lc3wC+Flm/gdAZt4JnA/8Q0Q8dZTPLUmSJDVNacE+MzeO5vEjYj/gBGDpoE1XUCzJeelonl+SJElqpvH8VJx3RcTVEXFzRFwVEW8ctP2oert6UH//57mjW54kSZLUPOM12P8R+B/gJOBI4F+AL0TEJweM2bvePjxo34fq7V6jWJ8kSZLUVGXePLvDMvPlg7q+ExHzgXdGxL9k5todPXZEtAFtADNmzNjxIiVJkqQmGq8z9kP5OcX3Oa7++b56u+ugcbvV2z8MdZDM7MzMWmbWpk+fPvJVSpIkSaNg3AX7iJgaEU8eYlP/S8tb6u0N9XbWoHEH1tsVI1yaJEmSVJoxH+wjYq+ImDqg6wTgW0MMPbbeXgeQmb8DrqFYhz/QfIqXYn1/ZCuVJEmSyjOmg31EHAj8Frh40KbnRcRfDBh3EvAW4GuDno3/HuCEiHhpfdz+wDuBT2Xmb0exdEmSJKmpSrt5NiLOo3jz7Iz65+X1Tc/KzPX1Xz8G3A/cNWDXaykC+/sj4uPAk4D1wMeA8waeIzOXRcRfAh+LiI8BU4HPAeeMxneSJEmSyhKZWXYNY1atVsuenp6yy5AkSVL1RaMHGNNLcSRJkiRtH4O9JEmSVAEGe0mSJKkCDPaSJElSBRjsJUkjorsb5syBlpai7e4uuyJJmlhKe9ylJKk6uruhvR2WLIF582DZMli0qNi2YEG5tUnSROGMvSSpYR0dRaifPx+mTCnaJUuKfklSc/gc+63wOfaStH1aWmDduiLU99uwAaZNg76+8uqSpHHE59hLkso3e3ax/GagZcuKfklScxjsJUkNa28v1tQvXVrM1C9dWnxuby+7MkmaOLx5VpLUsP4bZBcvhlWripn6jg5vnJWkZnKN/Va4xl6SJElN4hp7SZIkSQZ7SZIkqRIM9pIkSVIFGOwlSZKkCjDYS5IkSRVgsJckSZIqwGAvSZIkVYDBXpIkSaoAg70kSZJUAQZ7SZIkqQIM9pIkSVIFGOwlSZKkCjDYS5IkSRVgsJckSZIqwGAvSZIkVYDBXpI0Irq7Yc4caGkp2u7usiuSpIllctkFSJLGv+5uaG+HJUtg3jxYtgwWLSq2LVhQbm2SNFE4Yy9JalhHRxHq58+HKVOKdsmSol+S1ByRmWXXMGbVarXs6ekpuwxJGvNaWmDduiLU99uwAaZNg76+8uqSpHEkGj2AM/aSpIbNnl0svxlo2bKiX5LUHAZ7SVLD2tuLNfVLlxYz9UuXFp/b28uuTJImDm+elSQ1rP8G2cWLYdWqYqa+o8MbZyWpmVxjvxWusZckSVKTuMZekiRJksFekiRJqgSDvSRJklQBpQb7iNgvIi6JCBf6S5IkSQ0o7ak4EfFK4NPAhmHutwewCHgNsDPFP07uBD6SmcsGjf0yMA94ZNBhvpGZ5+xY5ZIkSdLYU+bjLs8CXgC0A08fxn5vAf4ROCkzV0ZEC/Ap4McR8eLMvHTQ+DMz88qRKFiSJEkaq8pcinNiZv56B/f9YmauBMjMPuC9QB+weKSKkyRJksaT0mbsM3PjDu563hDHWhcRDwB7NFaVJEmSND6Nu6fiZGZffZb+T+rr7vcGrhxil9dHxFURcVNE/Cwi3hURvnFXkiRJlVKVgHsmcA9w/qD+h+s/i4HHgPnAt4ETgL9qZoGSJEnSaBr3wT4ijqRYY39aZv5h4LbMHLzm/oqIOAc4NyJOzMyrhzheG9AGMGPGjFGqWpIkSRpZ424pzkARMQP4PvDmYTz55uf19vihNmZmZ2bWMrM2ffr0EahSkiRJGn3jNthHxAHAZcB7MvOiIba3RMSeQ+zavz6/ZTTrkyRJkpppzAf7iNgrIqYO6usP9e2Z+Z0B/d8bMOwA4JohDnlsvb12pGuVJEmSyjKm19hHxIHAKmAp8OJ63/71z8uBaRHxugG7zB10iEMjoi0zO+v7zgXeB1wFXD661UuSJEnNU9qMfUScFxHLgZfWPy+v/wycnX8MuB+4a0Dfu4CDgdOArw36GeguirfUnh4RKyLiFuAi4KvAqZmZI/+tJGni6u6GOXOgpaVou7vLrkiSJpYw325ZrVbLnp6essuQpDGvuxva22HJEpg3D5Ytg0WLoKMDFiwouzpJGhei0QOM+TX2kqSxr6OjCPXz58OUKUW7ZEnRL0lqDmfst8IZe0naPi0tsG5dEer7bdgA06ZBX9+W95Mk/Ykz9pKk8s2eXSy/GWjZsqJfktQcBntJUsPa24s19UuXFjP1S5cWn9vby65MkiaOMf24S0nS+NB/g+zixbBqVTFT742zktRcrrHfCtfYS5IkqUlcYy9JkiTJYC9JkiRVgsFekiRJqgCDvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJkqQKMNhLkiRJFWCwlyRJkirAYC9JkiRVgMFekiRJqgCDvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJkqQKMNhLkkZEdzfMmQMtLUXb3V12RZI0sUwuuwBJ0vjX3Q3t7bBkCcybB8uWwaJFxbYFC8qtTZImCmfsJUkN6+goQv38+TBlStEuWVL0S5KaIzKz7BrGrFqtlj09PWWXIUljXksLrFtXhPp+GzbAtGnQ11deXZI0jkSjB3DGXpLUsNmzi+U3Ay1bVvRLkprDYC9Jalh7e7GmfunSYqZ+6dLic3t72ZVJ0sThzbOSpIb13yC7eDGsWlXM1Hd0eOOsJDWTa+y3wjX2kiRJahLX2EuSJEky2EuSJEmVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAEGe0mSJKkCDPaSJElSBRjsJUmSpAooNdhHxH4RcUlE+PpbSZIkqQGlBfuIeCVwDXDwDu7/joi4KSJWRMS1EfHyLYxbGBHX18etjIg3N1C2JEmSNCZNLvHcZwEvANqBpw9nx4h4L/CPwLMz87aIeAHwXxHx0sz84YBxrwEuBJ6bmb+IiLnANRExKTMvGLFvIkmSJJWszKU4J2bmr4e7U0TsDnwQ+Hxm3gaQmZcBlwKfHDBuEnAe8M3M/EV93Argy8DHI2KnRr+AJEmSNFaUFuwzc+MO7noK0AosHdR/BXBERBxe//wsYP8tjNsTmL+D55ckSZLGnPH4VJy59Xb1oP7Vg7Zv7zhJkiRp3BuPwX7vevvwoP6H6u1ewxwnSRoBXV0waxZMmlS0XV1lVyRJE0uZN8+OSRHRBrQBzJgxo+RqJGl86OqCtjbo7S0+r1lTfAZYuLC8uiRpIhmPM/b31dtdB/XvVm//MMxxm8nMzsysZWZt+vTpDRUqSRNFe/sTob5fb2/RL0lqjvEY7FfU21mD+g8ctH17x0mSGrR27fD6JUkjb8wH+4jYKyKmDui6BOgFTho0dD5wU2beXP/8C+A3Wxh3P3DlSNcqSRPVllYuuqJRkppnTAf7iDgQ+C1wcX9fZv4R+CjwNxFxUH3c84EXUby0qn/cJuDdwGsiolYfdxRwBvD+zFzXnG8hSdXX0QGtrZv3tbYW/ZKk5ijt5tmIOI/izbMz6p+X1zc9KzPX13/9GMXs+l0D983McyJiHfCDiNgI9AGvGvjW2fq4b0REC3BhREDxD5l3ZWbn6HwrSZqY+m+QbW8vlt/MmFGEem+claTmicwsu4Yxq1arZU9PT9llSJIkqfqi0QOM6aU4kiRJkraPwV6SJEmqAIO9JEmSVAEGe0mSJKkCtivYR8TrI2JORPyv8RGx08iXJUmSJGk4tvdxl18BEng8Im4Erqv/LAdqEXF6Zs4bnRIlSZIkbcv2Bvs9gWPqP88ETgTexBOP5Xlo5EuTJEmStL22K9jX3/a6tP4DQERMB/4OeCewYDSKkyRJkrR9dvjNs5l5L/DBiNiTIthfMmJVSZIkSRqWkXgqzqXAy0fgOJIkSZJ20PY+FeeDEXFqROw3xOZ9gEdHtixJkiRJw7G9S3HeDewCZETcA1xL8VScPuBM4KzRKU+SJEnS9tjeYP9k4FCKJ+I8EzgWeDuwe337pyPidcCvgJ7MvGiE65QkSZK0FZGZO75zxEEUIf+ZA372zMyWkSmvXLVaLXt6esouQ5IkSdUX2x6ydTv8VByAzLwduB349p8qipjRaFGSJEmShmcknoqzmcxcO9LHlCSNfV1dMGsWTJpUtF1dZVckSRNLQzP2kiRBEeLb2qC3t/i8Zk3xGWDhwvLqkqSJZMRn7CVJE097+xOhvl9vb9EvSWoOg70kqWFrt7AIc0v9kqSRZ7CXJDVsxhYem7ClfknSyDPYS5Ia1tEBra2b97W2Fv2SpOYw2EuSGrZwIXR2wsyZEFG0nZ3eOCtJzdTQC6qqzhdUSZIkqUkafkGVM/aSJElSBRjsJUmSpAow2EuSJEkVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAEGe0mSJKkCDPaSJElSBRjsJUmSpAow2EuSJEkVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAEGe0nSiOjqglmzYNKkou3qKrsiSZpYJpd14ojYBzgfqNW7bgDekZm/2cZ+ZwDnAL8ftGkKcATw/My8vD72y8A84JFBY7+Rmec0Ur8k6QldXdDWBr29xec1a4rPAAsXlleXJE0kpczYR8RU4DJgKnAkRSB/FFgaEbtsxyG+mJlHD/wBPgbcBVw5aOyZg8ca6iVpZLW3PxHq+/X2Fv2SpOYoaynOG4C5wFmZuTEz+4CzgIOAt21j36uAi4boPxO4sH4sSVITrV07vH5J0sgrK9ifBqzNzNv7OzLz98BN9W1blJm3Z+aKgX0RcRDw58CXRqFWSdI2zJgxvH5J0sgrK9jPBVYP0b8aOGoHjncmcGlmrhli2+sj4qqIuCkifhYR74qI0u4tkKQq6uiA1tbN+1pbi35JUnOUFez3Bh4eov8hoDUidt7eA0VEC8XSns4hNj8M3A2cQrGW//1AO/CN4RYsSdqyhQuhsxNmzoSIou3s9MZZSWqmKsxc/0W9/cHgDZm5eFDXFRFxDnBuRJyYmVcP3ici2oA2gBn+H7IkbbeFCw3yklSmsmbs7wN2HaJ/N6A3Mx8bxrHOBP5vZm7czvE/r7fHD7UxMzszs5aZtenTpw+jDEmSJKk8ZQX7FcCsIfoPpHie/XaJiP2AFwH/NsS2lojYc4jd+p+a07K955EkSZLGurKC/UXAzIiY1d8REfsCs4HvDhwYEftGxJbqfCNwxRZumj0AuGaI/mPr7bXDLVqSJEkaq8oK9l+mmJn/RERMrgf3cyieivOF/kERcSLFS6c+N/gAERHAmxj6ptl+h9bXzPfvMxd4H8Wz8C9v/GtIkiRJY0MpwT4z1wMvoFgWcxOwimJ9/cmZ+ciAoY8ADwK/G+Iw84FW4PtbOM1dwFuA0yNiRUTcQvE/BV8FTs3MHInvIkmSJI0FYb7dslqtlj09PWWXIUmSpOqLRg9Q1lIcSZIkSSPIYC9JkiRVgMFekiRJqgCDvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJkqQKMNhLkiRJFWCwlyRJkirAYC9JGhFdXTBrFkyaVLRdXWVXJEkTy+SyC5AkjX9dXdDWBr29xec1a4rPAAsXlleXJE0kzthLkhrW3v5EqO/X21v0S5Kaw2AvSWrY2rXD65ckjTyDvSSpYTNmDK9fkjTyDPaSpIZ1dEBr6+Z9ra1FvySpOQz2kqSGLVwInZ0wcyZEFG1npzfOSlIzRWaWXcOYVavVsqenp+wyJEmSVH3R6AGcsZckSZIqwGAvSZIkVYDBXpIkSaoAg70kSZJUAQZ7SZIkqQIM9pIkSVIFGOwlSZKkCjDYS5IkSRVgsJckSZIqwGAvSZIkVYDBXpIkSaoAg70kSZJUAQZ7SZIkqQIM9pIkSVIFGOwlSSOiqwtmzYJJk4q2q6vsiiRpYplcdgGSpPGvqwva2qC3t/i8Zk3xGWDhwvLqkqSJxBl7SVLD2tufCPX9enuLfklScxjsJUkNW7t2eP2SpJFnsJckNWzGjOH1S5JGnsFektSwjg5obd28r7W16JckNYfBXpLUsIULobMTZs6EiKLt7PTGWUlqpsjMsmsYs2q1Wvb09JRdhiRJkqovGj1AaTP2EbFPRHRFxC31n+9ExP7bue8dEbF8iJ/nDzF2YURcHxErImJlRLx55L+NJEmSVK5SnmMfEVOBy4BbgSOBBC4ElkbEMZn5yLaOkZlHb8d5XlM/7nMz8xcRMRe4JiImZeYFjXwHSZIkaSwpa8b+DcBc4KzM3JiZfcBZwEHA20biBBExCTgP+GZm/gIgM1cAXwY+HhE7jcR5JEmSpLGgrGB/GrA2M2/v78jM3wM31beNhGcB+wNLB/VfAewJzB+h80iSJEmlKyvYzwVWD9G/Gjhqew4QEedGxE8j4taIuDQiXjrEOfqPOfgcA7dLkiRJ415ZwX5v4OEh+h8CWiNi523sfw9wLTCPYo3+xcDFEfG3g87BEOd5qN7uNayKJUmSpDFsXD7HPjOflZnfyMxNmbkhMz8H/BfF2vlpjRw7Itoioicieu69996RKViSJEkaZWUF+/uAXYfo3w3ozczHduCYP68f88gB52CI8+xWb/8w1EEyszMza5lZmz59+g6UIUmSJDVfWcF+BTBriP4DgRu2tmNE7BwRuwyxqa/etgw4B0Oc58BB2yVJkqRxr6xgfxEwMyJm9XdExL7AbOC7AwdGxL71R1f2Ox341BDHPBZ4nOLJOgC/AH4DnDRo3HzgfuDKHa5ekiRJGmPKCvZfppiZ/0RETK4H93Monljzhf5BEXEicBfwuUH7L4iI4waMOx14OXBu/8utMnMT8G7gNRFRq487CjgDeH9mrhuVbyZJkiSVoJQ3z2bm+oh4AXA+xQx7AiuBkwe9dfYR4EHgdwP6fkjx4qnPR8QUYHfgAeCtmdk56DzfiIgW4MKIgOIfMu8aPE6SJEka7yIzy65hzKrVatnT01N2GZIkSaq+aPQA4/Jxl5IkSZI2Z7CXJEmSKsBgL0kaEV1dMGsWTJpUtF1dZVckSRNLKTfPSpKqpasL2tqgt7f4vGZN8Rlg4cLy6pKkicQZe0lSw9rbnwj1/Xp7i35JUnMY7CVJDVu7dnj9kqSRZ7CXJDVsxozh9UuSRp7BXpLUsI4OaG3dvK+1teiXJDWHwV6S1LCFC6GzE2bOhIii7ez0xllJaibfPLsVvnlWkiRJTeKbZyVJkiQZ7CVJkqRKMNhLkiRJFWCwlyRJkirAYC9JkiRVgMFekiRJqgCDvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJkqQKMNhLkiRJFWCwlyRJkirAYC9JkiRVgMFekiRJqgCDvSRJklQBBntJ0ojo6oJZs2DSpKLt6iq7IkmaWCaXXYAkafzr6oK2NujtLT6vWVN8Bli4sLy6JGkiccZektSw9vYnQn2/3t6iX5LUHAZ7SVLD1q4dXr8kaeQZ7CVJDZsxY3j9kqSRZ7CXJDWsowNaWzfva20t+iVJzWGwlyQ1bOFC6OyEmTMhomg7O71xVpKaKTKz7BrGrFqtlj09PWWXIUmSpOqLRg/gjL0kSZJUAQZ7SZIkqQIM9pIkSVIFGOwlSZKkCjDYS5IkSRVgsJckSZIqYHJZJ46IfYDzgVq96wbgHZn5m23stx/wVuAVFP8wmQzcBJydmTcMGnslsA+wftBhPp2ZX230O0iSJEljRSnBPiKmApcBtwJHAglcCCyNiGMy85Gt7H42cDLwvMy8MyKmAV8Dfh4Rzx4c7oFTM/OOEf8SkiRJ0hhS1lKcNwBzgbMyc2Nm9gFnAQcBb9uO/c/NzDsBMnMd8F5gZ6BtlOqVJEmSxrSygv1pwNrMvL2/IzN/T7Gk5rRt7Pu3FLP7A91Vb/cYsQolSZKkcaSsYD8XWD1E/2rgqK3tWJ/h3zSo+9B6e+UQu7wrIq6OiJsj4qqIeOOwq5UkSZLGuLKC/d7Aw0P0PwS0RsTOwzxeG3AjxVr7gf4I/A9wEsVa/n8BvhARnxzm8SVJkqQxrbSn4oyUiHgecDrwZ5n5+MBtmfnyQcO/ExHzgXdGxL9k5tohjtdGfa3+jBkzRqdoSZIkaYSVNWN/H7DrEP27Ab2Z+dj2HCQingF8BXhpZt60nef+OcX3Pm6ojZnZmZm1zKxNnz59Ow8pSZIklausYL8CmDVE/4EUz7PfpoiYC/wH8JrM/OkQ26dGxJOH2LWv3rZsV6WSJEnSOFBWsL8ImBkRs/o7ImJfYDbw3YEDI2LfiJg0qG8ucDHw+sxcVu/bLyIuGDDsBOBbQ5z72Hp7XaNfQpL0hK4umDULJk0q2q6usiuSpImlrGD/ZYqZ+U9ExOR6cD+H4qk4X+gfFBEnUjzK8nMD+o4CLgd+BMyKiNdFxOso1tkfNug8z4uIvxiw70nAW4CvZeavR/5rSdLE1NUFbW2wZg1kFm1bm+FekpopMrOcExcz9OcDNYo3z64E3tH/4qn6mGcAS4HPZOZH6n0XAa/YwmF/nJkn1cftBpxJ8Vz8XYAnAespnpxzXmZu3FaNtVote3p6duj7SdJEMmtWEeYHmzkT7rij2dVI0rgUDR+grGA/HhjsJWn7TJpUzNQPFgGbBr95RJI0lIaDfVlLcSRJFbKlpwP71GBJah6DvSSpYR0d0Nq6eV9ra9EvSWoOg70kqWELF0JnZ7GmPqJoOzuLfklSc7jGfitcYy9JkqQmcY29JEmSJIO9JEmSVAkGe0mSJKkCDPaSJElSBRjsJUmSpAow2EuSJEkVYLCXJEmSKsBgL0mSJFWAwV6SJEmqAIO9JEmSVAEGe0mSJKkCDPaSJElSBRjsJUmSpAow2EuSJEkVYLCXJEmSKsBgL0kaEV1dMGsWTJpUtF1dZVckSRPL5LILkCSNf11d0NYGvb3F5zVris8ACxeWV5ckTSQGe0nSdrln7b2svPoWVv3sVh64+4+bbfvk98+kt3e3zfp6e+Fv2h5i9fe/BEBE8LSn78fs4w/hiBMOY9c9dmlW6ZI0IRjsJUn/y6ZNm7hj5Z2sXHYzK69excplN3PvnX/Y4vgH8x1D9/fuwo+/dc2Q22bNOYA5Jx7OnHmzOeq5h7PPjOkjUbokTViRmWXXMGbVarXs6ekpuwxJGnWZyZ233MXyK1ay/MqVXL/0Rh76w8Objdl1jydxxAmHMfv4Q3nqwU8h4oltr/nbY7n7vmn/67j77r2Ob3z2VwBs3NDH6hvWctM1t3DLL29jw+MbNhv71IP35RknzeHo+UdyzPOOYo99dx/x7ylJY1hse8g2DmCw3zKDvaQqW79uPb+8ZDm/vGQ5PT9azt1r7t1s+/T992Lunx/BnHmzmTPvcGbMfhqTJg39zIXBa+wBWluhs3PoNfbrH9/Ar391e/E/AstWseKqm+h96LHNxjz9mAOpvehonn3qMRx54uFENPx3niSNZQb70WSwl1RF6x/fwI8uvIKvf/wi7vvt/X/qf/Leu3LM847i6PlzOPrkOfVZ+e3/e6arC9rbYe1amDEDOjq2/8bZvr4+/ue6O1h+xUquu+IGbrjqJtave2JG/9Dawbzhw6dz3ClHG/AlVZXBfjQZ7CVVSWZy1bev4Uvv6+L3q+8BinXuJ736RGqnHM0hzzxwizPyzfb4Y4+z4qpV/OpHy7n868v44z0PAjD3z4+g7by/5rDawSVXKEkjzmA/mgz2kqrituvv4LOLl7By2c0AzJj9NN7w4dOZ98pnj5kwvyWPPbqOH3zhUrrP+Xcevv8RAE5543wWnbOQ3ac/ueTqJGnEGOxHk8Fe0ni3ft16vvbhb/OtT36PTX2b2H2fJ3PGR07nlDedTMvklrLLG5ZH/vgo3R+/iIv++T/ZuKGPXfd4Em/7zBt5/uv+zOU5kqrAYD+aDPaSxrPVN6zh46/9Z+648U4igpe+/UWc8dHXsMvuTyq7tIb85ta7+NfFS7j2shUAnPiKZ/Guf3sru+25a8mVSVJDDPajyWAvaby67Ks/5jNvvYD16zZwwGFP5R8vfDtHPOewsssaMZnJpV+5ks+/4//S+9BjTD9gLz500bs59FjX3ksatxoO9mN7YaUkaVgyky+99/9x7hmfZf26DZzyxvl8rucTlQr1ULzF9kVnzOeL153H4c8+hHvv/APvfO4H+clFPy+7NEkqjcFekipi06ZNnN92Ad8892JaJrfw919o4x+WvJ2dn/S/XxxVFfsduC+fuvLDnPKmk1m/bgMfe/WnuPQrV5ZdliSVwmAvSRWQmXzu7y7kh0suZ6edp/LR753FS97ygrLLaoqpO03hXf/2Vv76Q69m06bkk2/6PD/+1k/LLkuSms5gL0kV8O///F987/M/YspOU/jI997LcaccU3ZJTRURvP6fXsUZH30Nmckn3vBZVv3812WXJUlNZbCXpHHulp7b6HzP1wB4z5f/hmc+76iSKyrPa9//Sv7yrS9kw+Mb+Njpn+bRh3rLLkmSmsZgL0nj2KZNm/jnt15A38Y+XvF3p3LS6SeWXVKpIoK3feYMDq0dzD1r7+NrH/pW2SVJUtMY7CVpHFvafTW/vnY10/ffizd2LCi7nDFhytQpvPOCtxARXPy5S/jd6rvLLkmSmsJgL0njVGbyzXP/A4DXn/2qSj/9ZriefsyBnPzaeWzc0MdF5/9n2eVIUlOUGuwjYp+I6IqIW+o/34mI/bdz3ykR8dGIuDkiVkbETyNi3hbGviMiboqIFRFxbUS8fES/iCSVYNXPbmX1DWvZ8ym78/zX/1nZ5Yw5r373ywC47Gs/5vHHHi+5GkkafaUF+4iYClwGTAWOBI4AHgWWRsQu23GIfwVOB56bmXOAC4FLI+LoQed5L/AB4C8zcy5wFvDtiHjxSH0XSSrDVd++BoCTX/tcpkydUnI1Y89Bc2dyyDMP5NEHe7n2v28ouxxJGnVlzti/AZgLnJWZGzOzjyJ0HwS8bWs7RsRhQBtwTmbeC5CZXwJWAx0Dxu0OfBD4fGbeVh93GXAp8MmR/kKS1EzXXbESgONfcmzJlYxdx7+kBsB1lxvsJVVfmcH+NGBtZt7e35GZvwduqm/bmlcAASwd1H8F8MIBM/6nAK1bGHdERBy+g7VLUqk2rN/AHTfeSURw2LOeXnY5AHR1waxZMGlS0XZ1lV0RzD7+EAD+Z/nqkiuRpNFXZrCfSzHDPthqYFsPYZ4LbALWDrHvZIplPf3j+vsHjxu4XZLGlft+cz+b+jax99P2ZFrrTmWXQ1cXtLXBmjWQWbRtbeWH+6cdsh8Ad99xb7mFSFITlBns9wYeHqL/IaA1Inbexr699eU7g/cF2GvAOIY4z+BxkjSu9L94aZc9nlRyJYX2dugd9C6o3t6iv0z9vz+PPuiLqiRVn4+7HCQi2iKiJyJ67r3XGR5JY1NLS3H53tS3qeRKCmsH///pNvqbpf/3Z1KLf91Jqr4yr3T3AbsO0b8bxWz8Y9vYtzUiWobYF+APA8YxxHkGj/uTzOzMzFpm1qZPn76VEiSpPLvtXVzG7v/9H8stpG7GjOH1N8sD9d+f3afvtvWBklQBZQb7FcCsIfoPBLb1+IIVFLUfMMS+GyluwO0fxxDnOXDQdkkaV/Z8yu7svMs0Hr7/Ee7//QNll0NHB7S2bt7X2lr0l2n1yjuBJ9baS1KVlRnsLwJmRsSs/o6I2BeYDXx34MCI2DciBtb670ACJw065nzg0sx8pP75EqB3C+NuysybG/sKklSOiODwZxdPfFnx45u2MXr0LVwInZ0wcyZEFG1nZ9Ffpv7fm8OfdUi5hUhSE5QZ7L9MMTP/iYiYXA/u51A8seYL/YMi4kTgLuBz/X2ZeQvQCbwvIvauj3sjcDDQPmDcH4GPAn8TEQfVxz0feBHwj6P43SRp1B13yjEAXPXdn5VcSWHhQrjjDti0qWjLDvV9G/v46cW/AOC4Fx9dbjGS1ASlBfvMXA+8AOijWDqzimLt+8kDZtwBHgEeBH436BCLgW8DV0fESooXVr0wM5cPOs85FC+t+kFErADOA16VmT8c8S8lSU100uknMGlScM3Fv+S+u+4vu5wx5+r/+AUP3P0gBxz2VA555kFllyNJo25ymSfPzLuB125jzPXAnkP0bwA+UP/Z1nk+A3xmh4qUpDFq+v57ceIrnsVPvvtzvnXuxbz9M28su6QxY9OmTXz94xcB8NK/OYWIKLkiSRp9Pv9LksaxhR/4KyKC73/hR6xZ9ZuyyxkzLv3yldy2/A6m778XL150ctnlSFJTGOwlaRw7+BmzePGik9m4oY/zzvgsG9ZvKLuk0v1u9d188R++AsCZ5yxkp53LfzOvJDWDwV6SxrkzP/E69p05nVt+eRvnv+UCMrPskkrz6IOPcvbLz+XRB3s54WXHMX/BvLJLkqSmMdhL0ji36x678E/f+Qemte7EZV/5MV9811cmZLh/7JHHaH/J/2H1DWs54PCn8Y8Xvt219ZImFIO9JFXAoccezAe+9S4mT2nhon/+T85vu4C+jX1ll9U0f7z3Qd79vA9z49W3MP2Avfj4f72fXffYpeyyJKmpDPaSVBHPPvWZnP3ddzN12hR+uORy3n9qBw/94eGyyxp1t11/B4uf/T5u+eVtPOXAffjkFR/iKbP2KbssSWo6g70kVcjxLzmWc//7bHbf58lc+9838NZj3s2Kq8p/M+1oyEy+/4Uf8XfPeT+/v+NeDjvuYD6z7GM89eCnlF2aJJUiJuI6zO1Vq9Wyp6en7DIkadjuufM+Pnb6p1n1s18TEbzy70/ljI8tYFprNZ4Qc8/aezn/LRfQ86PrAXjxoufxt//6JqZOm1pyZZK0wxq+KchgvxUGe0nj2cYNG/l/H/kO3ef8O5v6NrHvzOm87fwzOOFlx43bm0o3rN/Ad8//T7o++h3W9T7Ornvuwt9//s38+atPKLs0SWqUwX40GewlVcGtv7qNT535BW6/fg0AR88/krbz/ppDnnlQyZVtv8zkJ9/9GUve18Vdt90NwHP/6nj+9l/exJ5P2aPk6iRpRBjsR5PBXlJV9G3s4/tfvJSvnv1NHn7gUQCee9qzef3Zr+bAOTNKrm7LMpNf/PA6vnL2N/n1r24H4IDDn8bbP/NGai98RsnVSdKIMtiPJoO9pKp5+IFH6P74RVz8uUtYv654S+2sIw+g9qKjOe6Uo5kz7/DS16k/dP/DXPffN/CLS66j50fXc//vHgBgz6fszuv+6VWceubzaJncUmqNkjQKDPajyWAvqaruu+t+vvF//p1Lv3Iljz2y7k/9U3aawpEnHMozTprD0SfP4bDjDmbK1CmjWsujDz7KDT+5meVLV3L9lTdy2/I7NnvB1t5P25PT3vkS/vJtL2Snnatx868kDcFgP5oM9pKqbsP6Ddz001vp+dFyfvmj5dy2/I7Ntk+dNoXDn30Ic048nDnzDueI5xzKk578pIbOee9v/sCNV9/MymU3s/Lqm1m9Yg2bNj3xd9GUqZM58sTDOO6UY6i96GgOPGrGuL3ZV5KGwWA/mgz2kiaaB+97iBU/vonrrljJ9VeuZO2q3262PSKYecT+TNtlGns9dXg3rfZt6OO26+/g3jv/sFl/y+QWDjvuYI6eX/wvwRHPOdSZeUkTkcF+NBnsJU10D973EDf+9BZurM+u39pzGwCb+jZtNss+HK277cyRJxzGkfX/BTjsuKdX5vn6ktQAg/1oMthL0ubWr1vPb279Hb9ffQ99fZuGtW8E7H/oU5kx+2lMmuSLzyVpkIaD/eSRqEKSNDFMnTaVg+bO5KC5M8suRZI0iDP2WxERDwO3lF2HNIS9gfvKLkLaAv98aqzyz6bGsmmZOaeRAzhjv3W3ZGat7CKkwSKixz+bGqv886mxyj+bGssiouH13y5ylCRJkirAYC9JkiRVgMF+6zrLLkDaAv9saizzz6fGKv9saixr+M+nN89KkiRJFeCMvSRJklQBBntJkiSpAgz20jgSEQdExHci4sGIeCgiLoqIGWXXJUmSymewH0JEfDAibo2ITRHx8rLrkQAiohW4AjgceAPweuAQYGlEPKnM2jRxeH2UpOGLiG9GxIqIuC4ifhERzxuN8/iCqqFdBnQBF5ZdiDTAm4GDgMMy838AImIF8GvgLcCnS6xNE4fXR0kavrdk5h8BIuIY4PKI2DszN43kSZyxH0Jm/iwzby+7DmmQlwI/6w/1AJm5GrgaeFlpVWlC8fqo8SQivh4RuZWfS8uuURNDf6ive/JoncdgL40fRwIrh+i/ETiiybVI0njwcYrr5k+A59R/TqlvOwdoK6kuTUARcX5E3A58FzhtpGfrwWAvjSd7Ag8M0X8/sEeTa5GkMS8zVwL7A1fX/7fpZ8C6+uYfZOYdpRWnCScz35mZBwELgXMjYupIn8NgL0mSKikiZgK7A9cP6D4aSOCGEkqSyMxLKCbkjhrpYxvsd5Dr9lSCBxh6Zn5LM/lSKbw+agx5Rr1dMaDvGOCOzHyohHo0AUXEzhFx4IDPzwH2Am4fNK7ha6dPxRlCRHwIOBOYDsyJiM8Cx2fmbwYM+zjFv7QeAN5T73sycAnFur0LmlawJoobKdbZD3YEcFOTa9EE5fVR48xciqU3twzoO4bNZ/Cl0bYz8PWI2BXYCDxKscZ+8KRcw9dOZ+yHkJkfysz9M3OnzNy7/uvfDBrjuj012/eA4yPioP6OiJgFnFjfJo06r48aZ54B3JiZfQARMQWYjcFeTZSZ92fmczJzTmYenZknZuYVQ4xr+NpZ+WAfEc/fxn9r9P9cOczjum5PzfZvwB3AxRHxsoh4KXAxcCfOgGoM8fqoMeQZbP7n8CnAFIqHDkhjykhcOyfCUpyfUvzrfFt6h3lc1+2pqTLz0Yg4GTgf+BoQwOXAOzLzkVKLkzbn9VGlq7+R+2DgcwO676V4/OUnIuKBzPxaKcVJQ2v42ln5YJ+ZvcDN2zs+InI7jhm4bk8lyMy1wGll16GJyeujxpPMfBRoGdS3jlF4Eom0Nc28dlZ+Kc5wZWZs66c+1HV7kiYUr4+SNHzNvHYa7AeJiD0i4gcRcWtEXB8Rl0bE04cY6ro9SROK10dJGr5mXjsN9v9bAp/JzEMz8xnAD4AvDRwwYN3ewDVQA9ftvb5ZxUpSE3l9lKTha9q1MzK3uexnQouIGvCdzJxVdi2SNJZ4fZSk4RvNa6cz9tv2DopHCkqSNvcOvD5K0nC9g1G6dlb+qTiNiIizgYOAtrJrkaSxxOujJA3faF87DfZbEBEfAE4FXlh/ZKYkCa+PkrQjmnHtNNgPof6vqf7f+AfLrkeSxgqvj5I0fM26dnrz7CARcSTFHci3Af1v89yYmbXyqpKk8nl9lKTha+a102AvSZIkVYBPxZEkSZIqwGAvSZIkVYDBXpIkSaoAg70kSZJUAQZ7SZIkqQIM9pIkSVIFGOwlSZKkCjDYS5IkSRVgsJckSZIqwGAvSRoREfGXEXFNRBwcEd0R8UhE/DIinl52bZI0ERjsJUkNi4j3AN8Dfgq8FHgxsAg4BPhaiaVJ0oQRmVl2DZKkcSwiDgduAD6bme+MiKXAfZn5qoj4IvAWYJfMfLTUQiWp4pyxlyQ1qg3YCHwwIgI4FvhVfdtv621LGYVJ0kRisJckNeok4CeZ+QjwdGBX4Lr6tqcC92TmQyXVJkkThsFektSoGcDv6r9+Zr3tD/bzKdbeS5JG2eSyC5AkjXsPAgfWf/1M4HeZeU9EvAY4GHh9aZVJ0gTijL0kqVHfBuZFxAeA5wK/i4iPAF8F3p+Zvyy1OkmaIJyxlyQ16mxgT+A9FOvrN1L8/XJGZn69zMIkaSLxcZeSpBEREQcBtwEvy0zX1UtSk7kUR5I0Umr1tqfUKiRpgjLYS5JGSo3ixtm7yi5EkiYil+JIkiRJFeCMvSRJklQBBntJkiSpAgz2kiRJUgUY7CVJkqQKMNhLkiRJFWCwlyRJkirAYC9JkiRVgMFekiRJqoD/DxoPYAjHcP+xAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = 0.11   # arc-length between two successive points (resolution of solution points)\n",
    "\n",
    "z0 = np.array([0.5*np.pi, 0]) # unloaded state\n",
    "zp = z0\n",
    "dz= np.array([0, ds])  # assumed step to the next point (determines direction)\n",
    "\n",
    "stable_phi = []\n",
    "stable_mu =[]\n",
    "nonstable_phi = []\n",
    "nonstable_mu = []\n",
    "stable_prev = True   # we known by intution that the initial state is stable\n",
    "\n",
    "for k in range(17):\n",
    "    ZS = root(sys, zp+dz, args=(zp,ds), jac=dsys, method='hybr') \n",
    "    zs =ZS.x   # extract value from solution structure\n",
    "    stable = d2Epot_dphi2(zs[0], zs[1]) > 0\n",
    "    if stable:\n",
    "        stable_phi.append( zs[0] )\n",
    "        stable_mu.append( zs[1] )\n",
    "    else:        \n",
    "        nonstable_phi.append( zs[0] )\n",
    "        nonstable_mu.append( zs[1] )\n",
    "        \n",
    "    if not stable==stable_prev:\n",
    "        branchpoint = 0.5*(zp+zs) # stability change appeared somewhere in between\n",
    "    \n",
    "    # move on\n",
    "    dz = zs - zp \n",
    "    zp = zs\n",
    "    stable_prev = stable\n",
    "  \n",
    "# plot points of mu-phi curve\n",
    "plt.plot(stable_phi, stable_mu, 'bo')\n",
    "plt.plot(nonstable_phi, nonstable_mu, 'bo', fillstyle='none')\n",
    "plt.xlabel('$\\\\varphi$')\n",
    "plt.ylabel('$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = ['$-\\\\frac{1}{2}\\pi$', '$0$', '$\\\\frac{1}{2}\\pi$', '$\\pi$', '$\\\\frac{3}{2}\\pi$']\n",
    "plt.xticks(tick_pos, labels)\n",
    "\n",
    "# plot distance constraint for initial point\n",
    "plt.plot(z0[0], z0[1], 'bo')   # initial point\n",
    "phi_range = np.linspace(-np.pi, np.pi)\n",
    "mu_range =  np.linspace(-1.1*ds, 1.1*ds)\n",
    "PHI, MU = np.meshgrid(phi_range, mu_range)\n",
    "ZERO = np.zeros(PHI.shape)\n",
    "DS = np.sqrt(distance_squared(PHI, MU, ZERO, ZERO)) # refer to initial point (delta_phi=delta_mu=0)\n",
    "plt.contour(PHI+phi0, MU, DS, [ds])  # distance constraint\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A filled circle denotes a stable equilibrium and an empty circle an instable equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ellipse around the initial point shows where the next point, in this case the first point, has to be located, as dictated by the distance constraint. We will find intersections of the constraint-ellipse with the $\\phi$-$\\mu$-curve for all directions, even when the $\\phi$-$\\mu$-curve goes horizontal and the formulation for fixed $\\mu$ fails. For a quantitative analysis and to receive some inspiration for setting the pseudo arc-length metric it is advisable to calculate the determinant of the Jacobian (in function ``dsys``)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that the upright position is an equilibrium and at $\\mu=1$ it becomes unstable.\n",
    "Now we take a closer look at this point (stability loss) and try to follow the stable equilibrium branching off into the buckled state. Therefore we start the iterations from a state without load increase but rotated out of the vertical position ($\\varphi\\ne\\pi/2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvYAAAH3CAYAAADQVRsYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsfklEQVR4nO3dfZCmaV0f+u9vZheXNqACC4cU9vQCMY4Mi0hLUMZzGJVIsIIcSQJzWkrNYBNPsingRFBajppDUwsopDREbTMEoR4bI1KinIQDBUNwkBdbxGVlBRV2Rl7URZC3ZsO+XOeP+2m3p+mZnZ6np+/uuz+fqq5r7t/99mtq6uE71173/VRrLQAAwN52oO8GAACAyQn2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAV/TdwG72hCc8ob3pTW/quw0AAIavJr2AGfsL+NSnPtV3CwAAcFEEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHAIABEOwBAGAABHsAABgAwR4AAAZAsAcAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAdgWy8vJkSPJwYPduLzcd0cA+8sVfTcAwN63vJwsLCQnTyZHjyanTycnTnT7jh/vtzeA/cKMPQATW1zsQv2xY8mVV3bjyZNdHYCdIdgDMLGbbupm6tc7erSrw34wGiUzM8mBA904GvXdEfuRYA/AxA4f7pbfrHf6dFeHoRuNkvn55MyZpLVunJ8X7tl5gj0AE1tY6NbUnzqV3HZbN5440dVh6BYWktXVc2urq/7+s/M8PAvAxNYekL3uum75zeHD3fp6D86yH5w9u7U6XC6CPQDb4vhxQZ79aXq6W36zWR12kqU4AAATWFxMpqbOrU1NeSsUO0+wBwCYwNxcsrSUHDqUVHXj0lJXh51UrbW+e9i1Zmdn28rKSt9tAAAwfDXpBczYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAAMg2AMAwAAI9gAAMACCPQAADIBgDwAAA9BrsK+qB1bVm6qq9dkHAADsdVf0deOq+v4kL0ty2xbPe1yS1yc5u8nuRyR5Rmvt5PjYn07yL5N8esNx72it/dutdQwAALtXnzP2z0vy+CTvvIRzf7u19s3rf5L8aJJb04X+9f7vjccK9QDbb3k5OXIkOXiwG5eX++4IYH/pbcY+yWNba7dX1VbPuynJL29Sf0aS17XWPjNxZwBsyfJysrCQnDyZHD2anD6dnDjR7Tt+vN/eAPaL3mbsW2u3X+J5f9Vae9f6WlXdK8lTkyxtR28AbM3iYhfqjx1LrryyG0+e7OowBKNRMjOTHDjQjaNR3x3BV+pzxn47HU9ytrX2u5vse0JVPT3J/dOt539jkutba6s72SDAkN10UzdTv97Ro10d9rrRKJmfT1bHyeHMmW47Sebm+usLNhrK6y6fkc1n61eTfDHJU1prD09yIskPJHlrVV25g/0BDNrhw93ym/VOn+7qsNctLNwV6tesrnZ12E32fLCvqocneXiSV2/c11p7SWvtRGvtb8fbf5juod3HJPkX57nefFWtVNXKLbfccvkaBxiQhYVuTf2pU8ltt3XjiROCD8NwdrP38F2gDn0ZwlKcH0nym621ja+0PJ/3jMfHJPmKFXKttaWMZ/9nZ2e9Xx/gIqw9IHvddd3ym8OHu/X1HpxlCKanu+U3m9VhN9nTM/ZV9VVJ5rL5W3JSVVdvUr5jPB68XH0B7EfHjyc33pjccUc3CvUMxeJiMjV1bm1qysPh7D67PthX1X2r6h7n2f2UJH91nodmk+RMVW0M8I8aj+/blgYBgEGbm0uWlpJDh5Kqblxa8uAsu8+uDvZVdU2Sjyd5w3kOeUaSX7nAJe6Z5GfWwn1VHUpyfZIPJfm1bWwVABiwubnk5puTO+/sRqGe3ai3NfZV9dJ03zw7Pd5+/3jXo1trXx7/+UtJPp3kE5uc/5B06+T/2QVuM5fk/0jy/nG4n0rypiQv8LpLAACGpFrzfOj5zM7OtpWVlb7bAABg+GrSC+zqpTgAAMDFEewBAGAABHsAABgAwR4AAAZAsAcAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHYFssLydHjiQHD3bj8nLfHQHsL1f03QAAe9/ycrKwkJw8mRw9mpw+nZw40e07frzf3gD2CzP2AExscbEL9ceOJVde2Y0nT3Z12M1Go2RmJjlwoBtHo747gktXrbW+e9i1Zmdn28rKSt9tAOx6Bw8mt97ahfo1t92WXHVVcscd/fUFFzIaJfPzyerqXbWpqWRpKZmb668v9q2a9AJm7AGY2OHD3fKb9U6f7uqwWy0snBvqk257YaGffmBSgj0AE1tY6NbUnzrVzdSfOtVtC0jsZmfPbq0Ou52HZwGY2NoDstddl9x0UzdTv7jowVl2t+np5MyZzeuwF5mxB2BbHD+e3Hhjt6b+xhuFena/xcVuTf16U1Me+mbvEuwBgH1pbq57UPbQoaSqGz04y17mrTgX4K04AADsEG/FAQAABHsAABgEwR4AAAZAsAcAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHAIABEOwBAGAABHsAABgAwR4AAAZAsAdgWywvJ0eOJAcPduPyct8dAewvgj0AE1teThYWkl/4heTWW7txYUG4px+jUTIzkxw40I2jUd8dwc6o1lrfPexas7OzbWVlpe82AHa9I0e6MH/s2F21U6eS665Lbryxv77Yf0ajZH4+WV29qzY1lSwtJXNz/fUFF6EmvoBgf36CPcDFOXiwm6m/8sq7arfdllx1VXLHHf31xf4zM5OcOfOV9UOHkptv3uluYEsmDvaW4gAwscOHk9Onz62dPt3VYSedPbu1OgyJYA/AxBYWkhMnuuU3t93WjSdOdHXYSdPTW6vDkFzR582r6oFJ/kuS72mtTfyfHwDox/Hj3XjddclNN3Uz9YuLd9Vhpywubr7GfnGxv55gp/QW7Kvq+5O8LMltl3Du25PcP8mXN+x6WWvt1RuO/cdJXpjknkmuTPKrSV7cWrvzEtoG4DyOHxfk6d/aA7ILC93ym+npLtR7cJb9oM8Z++cleXyShSQPvYTzn9hau/lCB1TV0SRvTPIvWmu/VVVfn+Q9Se6d5Ccu4Z4AwC43NyfIsz/1ucb+sa21P73M93hxkne31n4rSVprf5Hk5Un+r6r6+5f53gAAsGN6C/attdsv5/XH6/e/PcmpDbvelm5JzpMu5/0BAGAn7eW34jynqt5ZVX9SVe+oqh/esP/h4/GjG+pr29de3vYAAGDn7NVg/7dJ/izJ45I8LMnPJ/nFqvrZdcfcbzx+fsO5nxuP993swlU1X1UrVbVyyy23bFvDAABwOe3JYN9ae3Jr7edba7e11u5orb0uyckkz66qid5U21pbaq3NttZmr7766u1pGAAALrM9GezP4z3pfp9vHW9/ajzea8Nx9x6Pf7MTTQEAwE7Yc8G+qu5RVV+zya47xuPB8fiB8Tiz4bhrxuMN29waAAD0ZtcH+6q6b1XdY13p25P8100OfdR4/MMkaa19Msm70q3DX+9Yui/F+p3t7RQAAPqzq4N9VV2T5ONJ3rBh13dV1feuO+5xSZ6Z5DUb3o3/3CTfXlVPGh/3oCTPTvJzrbWPX8bWAQBgR/X2zbNV9dJ03zw7Pd5+/3jXo1trXx7/+UtJPp3kE+tOfV+6wP78qnpRkq9O8uUkL0zy0vX3aK2drqp/muSFVfXCJPdI8ook11+O3wkAAPpSrbW+e9i1Zmdn28rKSt9tAAAwfDXpBXb1UhwAAODiCPYAADAAgj0AAAyAYA8AAAMg2AMAu9polMzMJAcOdONo1HdHsDsJ9gBsi+Xl5MiR5ODBblxe7rsjhmA0SubnkzNnkta6cX5euIfNCPYATGx5OVlYSH7hF5Jbb+3GhQXhnsktLCSrq+fWVle7OnAu77G/AO+xB7g4R450Yf7Ysbtqp04l112X3Hhjf32x9x040M3Ub1SV3HnnzvcDl5H32APQv5tuSo4ePbd29GhXh0lMT2+tDvuZYA/AxA4fTk6fPrd2+nRXh0ksLiZTU+fWpqa6OnAuwR6AiS0sJCdOdMtvbrutG0+csA6ayc3NJUtLyaFD3fKbQ4e67bm5vjuD3cca+wuwxh7g4i0vd7OoN93UzdQvLCTHj/fdFcCeMfEae8H+AgR7AAB2iIdnAQAAwR4AAAZBsAcAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHAIABEOwBAGAABHsAABgAwR4AAAZAsAcALrvRKJmZSQ4c6MbRqO+OYHiu6LsBAGDYRqNkfj5ZXe22z5zptpNkbq6/vmBozNgDsC2Wl5MjR5KDB7txebnvjtgtFhbuCvVrVle7OrB9zNgDMLHl5S6knTyZHD2anD6dnDjR7Tt+vN/e6N/Zs1urA5fGjD0AE1tc7EL9sWPJlVd248mTXR2mp7dWBy6NYA/AxG66qZupX+/o0a4Oi4vJ1NS5takp//CD7SbYAzCxw4e75TfrnT7d1WFuLllaSg4dSqq6cWnJg7Ow3QR7ACa2sNCtqT91Krnttm48ccLDkdxlbi65+ebkzju7UaiH7efhWQAmtvaA7HXXdctvDh/ulll4cBZg51Rrre8edq3Z2dm2srLSdxsAAAxfTXoBS3EAAGAABHsAABgAwR4AAAZAsAcAgAEQ7AEAYAB6DfZV9cCqelNVeTUPAABMoLf32FfV9yd5WZLbtnje1yU5keRpSe6Z7h8nf5Hk37fWTm849lVJjib5wobLvLa1dv2ldQ4AALtPn19Q9bwkj0+ykOShWzjvmUn+XZLHtdZurKqDSX4uyf+oqn/SWnvzhuOf0Vp7+3Y0DAAAu1WfS3Ee21r700s895daazcmSWvtjiQ/nuSOJNdtV3MAALCX9DZj31q7/RJPfekm17q1qj6T5Osm6woAAPamPfdWnNbaHeNZ+r8zXnd/vyRv3+SUp1fVO6rqg1X17qp6TlX1uQQJAAC23Z4L9ufxjCR/neTlG+qfT/JXSZ6Q5GFJnp9uTf9rz3ehqpqvqpWqWrnlllsuU7sAsPuNRsnMTHLgQDeORn13BFxItdbvmybHb675wdZaXeL5D0vyjiRPuZiHZKvqx5K8JMnR1to7L3Ts7OxsW1lZuZS2AGBPG42S+flkdfWu2tRUsrSUzM311xcM2CVl4fX29Ix9VU0n+Z0kP7KFN9+8Zzw+5rI0BQADsLBwbqhPuu2FhX76Ae7eng32VfX1Sd6S5Lmttddvsv9gVd1nk1PX1ucfvJz9AcBedvbs1upA/3Z9sK+q+1bVPTbU1kL9Qmvtdevqv73usK9P8q5NLvmo8fi+7e4VAIZienprdaB/uzrYV9U1ST6e5A3rag9KcirJjUmuqqofWPtJcu2GS3xDVc2vO/faJD+Rbk3+Wy93/wCwVy0udmvq15ua6urA7tRbsK+ql1bV+5M8abz9/vHP+tn5LyX5dJJPrKs9J8lDkjwlyWs2/Kz3iXTfUvvUqrqhqj6U5PVJXp3kia3vp4YBBmZ5OTlyJDl4sBuXl/vuiEnMzXUPyh46lFR1owdnYXfr/a04u5m34gBcnOXl7qHKkyeTo0eT06eTEye62d3jx/vuDmBP2N9vxQFgd1hc7EL9sWPJlVd248mTlm0A7CQz9hdgxh7g4hw8mNx6axfq19x2W3LVVckdd5z/PAD+jhl7APp3+HC3/Ga906e7OgA7Q7AHYGILC92a+lOnupn6U6e6bV9mBLBzrui7AQD2vrUHZK+7Lrnppm6m3oOzADvLGvsLsMYeAIAdYo09AAAg2AMAwCAI9gAAMACCPQAADIBgDwAAAyDYAwDAAAj2ALDHjUbJzExy4EA3jkZ9dwT0wRdUAcAeNhol8/PJ6mq3feZMt50kc3P99QXsPDP2ALCHLSzcFerXrK52dWB/EewBYA87e3ZrdWC4BHsA2MOmp7dWB4ZLsAeAPWxxMZmaOrc2NdXVgf1FsAeAPWxuLllaSg4dSqq6cWnJg7OwH1Vrre8edq3Z2dm2srLSdxsAAAxfTXoBM/YAADAAgj0AAAyAYA8AAAMg2AOwLZaXkyNHkoMHu3F5ue+OAPaXK/puAIC9b3m5+6bTkyeTo0eT06eTEye6fceP99sbwH5hxh6AiS0udqH+2LHkyiu78eRJ71IH2Eled3kBXncJcHEOHkxuvbUL9Wtuuy256qrkjjv66wtgD/G6SwD6d/hwt/xmvdOnuzoAO0OwB2BiCwvdmvpTp7qZ+lOnuu2Fhb47A9g/BHsAJnb8eLee/rrruuU3113Xbe/nB2dHo2RmJjlwoBtHo747AobOGvsLsMYegEsxGiXz88nq6l21qalkaSmZm+uvL2BXs8YeAHabhYVzQ33SbVuaBFxOgj0AbLOzZ7dWB9gOgj0AbLPp6a3VAbaDYA8A22xxsVtTv97UlC/sAi4vwR4AttncXPeg7KFDSVU3enAWuNy8FecCvBUHAIAd4q04AACAYA8AAIMg2AMAwAAI9gAAMACCPQCDNBolMzPJgQPdOBr13RHA5dVrsK+qB1bVm6rKq3kA2DajUTI/n5w5k7TWjfPzwj0wbL0F+6r6/iTvSvKQSzz/WVX1waq6oareV1VPPs9xc1X1R+PjbqyqH5mgbQD2gIWFZHX13NrqalcHGKorerz385I8PslCkodu5cSq+vEk/y7JP2qt/XlVPT7Jf6uqJ7XW/vu6456W5JVJvqO19t6qujbJu6rqQGvtl7ftNwFgVzl7dmt1gCHocynOY1trf7rVk6rqa5O8IMl/aq39eZK01t6S5M1JfnbdcQeSvDTJr7fW3js+7oYkr0ryoqr6qkl/AQB2p+nprdUBhuCign1VPb2qjozD8sZ9lxSQW2u3X8p5SZ6QZCrJqQ31tyX5pqr6xvH2o5M86DzH3SfJsUu8PwC73OJiMjV1bm1qqqsDDNXFztj/apI/SvKFqvr9qlqqqh+tqm9LMl9Vpy9fi1/h2vH40Q31j27Yf7HHATAwc3PJ0lJy6FBS1Y1LS10dYKguNtjfJ8l3J/nJJB9K8tgkv5DkdJL/kORhl6O587jfePz8hvrnxuN9t3jcOapqvqpWqmrllltumahRgP1kt71ecm4uufnm5M47u1GoB4buoh6eba39bbolLX+3rKWqrk7yb5M8O8nxy9FcH1prS0mWkmR2dtZrOAEuwtrrJdfeRLP2eslEoAbYKZf88Gxr7ZbW2gvSLdPZyWD/qfF4rw31e4/Hv9nicQBMyOslAfq3HW/FeXOSJ2/DdS7WDeNxZkP9mg37L/Y4ACbk9ZIA/bvYt+K8oKqeWFUP3GT3/ZN8cXvbOufe962qe6wrvSnJapLHbTj0WJIPttb+ZLz93iQfO89xn07y9u3uFWC/8npJgP5d7Iz9jyV5Y5KPVdUnq+r/raoXVtXPJPnpdF82te2q6pokH0/yhrXaeL3//5PkX1fVg8fHfXeS70n3pVVrx9057vtpVTU7Pu7hSX4oyfNba7dejp4B9iOvlwTo38V+8+zXJPmGJN8y/nlUkv8zydeO97+sqn4gyR8kWWmtvf7uLlhVL033zbPT4+33j3c9urX25fGfv5Rudv0T689trV1fVbcmeWNV3Z7kjiT/fP23zo6Pe21VHUzyyqpKun/IPGf8gCwA22TtAdmFhW75zfR0F+o9OAuwc6q1S3/xy3jG/FG5K/B/S5L7tNYObk97/ZqdnW0rKyt9twEAwPDVpBe42Bn7TbXWPpLkI0l+4+86qrKiEgAAdth2vBXnHK0170AAAIAdtu3BHgAA2HmCPQAADIBgDwAAAyDYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAAMg2AMAwAAI9gAAMACCPQAADIBgD8C2GI2SmZnkwIFuHI367ghgf7mi7wYA2PtGo2R+Plld7bbPnOm2k2Rurr++APYTM/YATGxh4a5Qv2Z1tasDsDMEewAmdvbs1uoAbD/BHoCJTU9vrQ7A9hPsAZjY4mIyNXVubWqqqwOwMwR7ACY2N5csLSWHDiVV3bi05MFZgJ1UrbW+e9i1Zmdn28rKSt9tAAAwfDXpBczYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAAMg2AMAwAAI9gAAMACCPQAADIBgDwAAAyDYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAbIvRKJmZSQ4c6MbRqO+OAPaXK/puAIC9bzRK5ueT1dVu+8yZbjtJ5ub66wtgP6nWWj83rrp/kpcnmR2XPpDkWa21j93NeT+U5Pokf7lh15VJvinJd7fW3jo+9lVJjib5woZjX9tau/7uepydnW0rKyt3dxjAvjcz04X5jQ4dSm6+eae7AdiTatIL9DJjX1X3SPKWJB9O8rAkLckrk5yqqke21jYG8Y1+qbX20xuueTzJzyZ5+4Zjn9Fa21gDYBudPbu1OgDbr6819j+Y5Nokz2ut3d5auyPJ85I8OMmP3s2570jy+k3qz0jyyvG1ANhB09NbqwOw/foK9k9Jcra19pG1QmvtL5N8cLzvvFprH2mt3bC+VlUPTvK/JfnPl6FXAO7G4mIyNXVubWqqqwOwM/oK9tcm+egm9Y8mefglXO8ZSd7cWttkhWeeXlXvqKoPVtW7q+o5VeWhYYBtNDeXLC11a+qrunFpyYOzADupr4B7vyR/sEn9c0mmquqerbUvXcyFqupguqU9/3qT3Z8f/1yX5EtJjiX5jSTfnuSfned680nmk2Taf0MGuGhzc4I8QJ+G8B777x2Pb9y4o7V2XWvt+a211dZ5W7o36jylqh672cVaa0uttdnW2uzVV199GdsGAIDt01ew/1SSe21Sv3eS1YudrR97RpL/0lq7/SKPf894fMwW7gEAALtaX8H+hiQzm9SvSfc++4tSVQ9M8j1JfmWTfQer6j6bnLb21pyDF3sfAADY7foK9q9PcqiqZtYKVfWAJIeT/Ob6A6vqAVV1vj5/OMnbzvPQ7Ncnedcm9UeNx/dttWkAANit+gr2r0o3M//iqrpiHNyvT/dWnF9cO2i8Dv4TSV6x8QJVVUn+ZZKlC9znG8YPw66dc22Sn0j3Lvy3Tv5rAADA7tBLsG+tfTnJ49Mti/lgkpvSra//zg3fOvuFJJ9N8slNLnMsyVSS3znPbT6R5JlJnlpVN1TVh9L9l4JXJ3lia61tx+8CAAC7Qcm35zc7O9tWVlb6bgMAgOGrSS8whNddAgDAvifYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAAMg2AMAwAAI9gBsi9EomZlJDhzoxtGo744A9pcr+m4AgL1vNErm55PV1W77zJluO0nm5vrrC2A/MWMPwMQWFu4K9WtWV7s6ADtDsAdgYmfPbq0OwPYT7AGY2PT01uoAbD/BHoCJLS4mU1Pn1qamujoAO0OwB2Bic3PJ0lJy6FBS1Y1LSx6cBdhJ1Vrru4dda3Z2tq2srPTdBgAAw1eTXsCMPQAADIBgDwAAAyDYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAAMg2AMAwAAI9gAAMACCPQAADIBgDwAAAyDYAwDAAAj2AAAwAII9ANtiNEpmZpIDB7pxNOq7I4D95Yq+GwBg7xuNkvn5ZHW12z5zpttOkrm5/voC2E/M2AMwsYWFu0L9mtXVrg7AzhDsAZjY2bNbqwOw/QR7ACY2Pb21OgDbT7AHYGKLi8nU1Lm1qamuDsDOEOwBmNjcXLK0lBw6lFR149KSB2cBdlK11vruYdeanZ1tKysrfbcBAMDw1aQXMGMPAAAD0Fuwr6r7V9Woqj40/nldVT3oIs+9uarev8nPd29y7FxV/VFV3VBVN1bVj2z/bwMAAP3q5QuqquoeSd6S5MNJHpakJXllklNV9cjW2hfu7hqttW++iPs8bXzd72itvbeqrk3yrqo60Fr75Ul+BwAA2E36mrH/wSTXJnlea+321todSZ6X5MFJfnQ7blBVB5K8NMmvt9bemySttRuSvCrJi6rqq7bjPgAAsBv0FeyfkuRsa+0ja4XW2l8m+eB433Z4dJIHJTm1of62JPdJcmyb7gMAAL3rK9hfm+Sjm9Q/muThF3OBqnpJVf1eVX24qt5cVU/a5B5r19x4j/X7AQBgz+sr2N8vyec3qX8uyVRV3fNuzv/rJO9LcjTdGv03JHlDVf2bDffIJvf53Hi875Y6BgCAXWxPvu6ytfbo1tprW2t3ttZua629Isl/S7d2/qpJrl1V81W1UlUrt9xyy/Y0DAAAl1lfwf5TSe61Sf3eSVZba1+6hGu+Z3zNh627Rza5z73H499sdpHW2lJrbba1Nnv11VdfQhsAALDz+gr2NySZ2aR+TZIPXOjEqrpnVf29TXbdMR4PrrtHNrnPNRv2AwDAntdXsH99kkNVNbNWqKoHJDmc5DfXH1hVDxi/unLNU5P83CbXfFSS/5nuzTpJ8t4kH0vyuA3HHUvy6SRvv+TuAQBgl+kr2L8q3cz8i6vqinFwvz7dG2t+ce2gqnpskk8kecWG849X1beuO+6pSZ6c5CVrX27VWrszyY8leVpVzY6Pe3iSH0ry/NbarZflNwMAgB708s2zrbUvV9Xjk7w83Qx7S3Jjku/c8K2zX0jy2SSfXFf77+m+eOo/VdWVSb42yWeS/KvW2tKG+7y2qg4meWVVJd0/ZJ6z8TgAANjrqrXWdw+71uzsbFtZWem7DQAAhq8mvcCefN0lAABwLsEeAAAGQLAHAIABEOwB2BajUTIzkxw40I2jUd8dAewvvbwVB4BhGY2S+flkdbXbPnOm206Subn++gLYT8zYAzCxhYW7Qv2a1dWuDsDOEOwBmNjZs1urA7D9BHsAJjY9vbU6ANtPsAdgYouLydTUubWpqa4OwM4Q7AGY2NxcsrSUHDqUVHXj0pIHZwF2UrXW+u5h15qdnW0rKyt9twEAwPDVpBcwYw8AAAMg2AMAwAAI9gAAMACCPQAADIBgDwAAAyDYAwDAAAj2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAAMg2AMAwAAI9gAAMACCPQAADIBgDwAAAyDYA7AtRqNkZiY5cKAbR6O+OwLYX67ouwEA9r7RKJmfT1ZXu+0zZ7rtJJmb668vgP3EjD0AE1tYuCvUr1ld7eoA7AzBHoCJnT27tToA20+wB2Bi09NbqwOw/QR7ACa2uJhMTZ1bm5rq6gDsDMEegInNzSVLS8mhQ0lVNy4teXAWYCdVa63vHnat2dnZtrKy0ncbAAAMX016ATP2AAAwAII9AAAMgGAPAAADINgDAMAACPYAADAAgj0AAAyAYA8AAANwRV83rqr7J3l5ktlx6QNJntVa+9jdnPfAJP8qyf+e7h8mVyT5YJKfaq19YMOxb09y/yRf3nCZl7XWXj3p7wAAALtFL8G+qu6R5C1JPpzkYUlaklcmOVVVj2ytfeECp/9Uku9M8l2ttb+oqquSvCbJe6rqH20M90me2Fq7edt/CQAA2EX6Worzg0muTfK81trtrbU7kjwvyYOT/OhFnP+S1tpfJElr7dYkP57knknmL1O/AACwq/UV7J+S5Gxr7SNrhdbaX6ZbUvOUuzn336Sb3V/vE+Px67atQwAA2EP6CvbXJvnoJvWPJnn4hU4cz/DfuaH8DePx7Zuc8pyqemdV/UlVvaOqfnjL3QIAwC7XV7C/X5LPb1L/XJKpqrrnFq83n+SP0621X+9vk/xZkselW8v/80l+sap+dovXBwCAXa23t+Jsl6r6riRPTfK/ttb+5/p9rbUnbzj8dVV1LMmzq+rnW2tnN7nefMZr9aenpy9P0wAAsM36mrH/VJJ7bVK/d5LV1tqXLuYiVfWIJL+a5EmttQ9e5L3fk+73/tbNdrbWllprs6212auvvvoiLwkAAP3qK9jfkGRmk/o16d5nf7eq6tokv5Xkaa2139tk/z2q6ms2OfWO8XjwojoFAIA9oK9g//okh6pqZq1QVQ9IcjjJb64/sKoeUFUHNtSuTfKGJE9vrZ0e1x5YVb+87rBvT/JfN7n3o8bjH076SwAAwG7RV7B/VbqZ+RdX1RXj4H59urfi/OLaQVX12HSvsnzFutrDk7w1yf+XZKaqfqCqfiDdOvt/uOE+31VV37vu3McleWaS17TW/nT7fy2A/Ws0SmZmkgMHunE06rsjgP2ll4dnW2tfrqrHJ3l5unfXtyQ3JvnODd86+4Ukn03yyXW1n0n3Vp1njn/W+x/r/vy+JM9N8vyqelGSr07y5SQvTPLS7fttABiNkvn5ZHW12z5zpttOkrm5/voC2E+qtdZ3D7vW7OxsW1lZ6bsNgF1vZqYL8xsdOpTcfPNOdwOwJ9WkF+hrKQ4AA3L2K14efOE6ANtPsAdgYuf72g9fBwKwcwR7ACa2uJhMTZ1bm5rq6gDsDMEegInNzSVLS92a+qpuXFry4CzATvLw7AV4eBYAgB3i4VkAAECwBwCAQRDsAQBgAAR7AAAYAMEeAAAGQLAHAIABEOwBAGAABHsAABgAwR4AAAZAsAcAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEegG0xGiUzM8mBA904GvXdEcD+ckXfDQCw941Gyfx8srrabZ85020nydxcf30B7Cdm7AGY2MLCXaF+zepqVwdgZwj2AEzs7Nmt1QHYfoI9ABObnt5aHYDtJ9gDMLHFxWRq6tza1FRXB2BnCPYATGxuLllaSg4dSqq6cWnJg7MAO6laa333sGvNzs62lZWVvtsAAGD4atILmLEHAIABEOwBAGAABHsAABgAwR4AAAZAsAcAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHAIABEOwBAGAABHsAABiAaq313cOuVVWfT/KhvvuATdwvyaf6bgLOw99Pdit/N9nNrmqtHZnkAldsVycD9aHW2mzfTcBGVbXi7ya7lb+f7Fb+brKbVdXKpNewFAcAAAZAsAcAgAEQ7C9sqe8G4Dz83WQ38/eT3crfTXazif9+engWAAAGwIw9AAAMgGAPAAADINjDHlJVX19Vr6uqz1bV56rq9VU13XdfAED/BPtNVNULqurDVXVnVT25734gSapqKsnbknxjkh9M8vQk/yDJqar66j57Y//w+QiwdVX161V1Q1X9YVW9t6q+63LcxxdUbe4tSUZJXtl3I7DOjyR5cJJ/2Fr7sySpqhuS/GmSZyZ5WY+9sX/4fATYume21v42SarqkUneWlX3a63duZ03MWO/idbau1trH+m7D9jgSUnevRbqk6S19tEk70zyfb11xb7i85G9pKp+raraBX7e3HeP7A9roX7say7XfQR72DseluTGTep/nOSbdrgXgL3gRek+N383ybeNf54w3nd9kvme+mIfqqqXV9VHkvxmkqds92x9ItjDXnKfJJ/ZpP7pJF+3w70A7HqttRuTPCjJO8f/tendSW4d735ja+3m3ppj32mtPbu19uAkc0leUlX32O57CPYAwCBV1aEkX5vkj9aVvzlJS/KBHlqCtNbelG5C7uHbfW3B/hJZt0cPPpPNZ+bPN5MPvfD5yC7yiPF4w7raI5Pc3Fr7XA/9sA9V1T2r6pp129+W5L5JPrLhuIk/O70VZxNV9dNJnpHk6iRHquo/JnlMa+1j6w57Ubp/aX0myXPHta9J8qZ06/Z+eccaZr/443Tr7Df6piQf3OFe2Kd8PrLHXJtu6c2H1tUemXNn8OFyu2eSX6uqeyW5PckX062x3zgpN/Fnpxn7TbTWfrq19qDW2le11u43/vPHNhxj3R477beTPKaqHrxWqKqZJI8d74PLzucje8wjkvxxa+2OJKmqK5McjmDPDmqtfbq19m2ttSOttW9urT22tfa2TY6b+LNz8MG+qr77bv6zxtrP27d4Xev22Gm/kuTmJG+oqu+rqicleUOSv4gZUHYRn4/sIo/IuX8P/5ckV6Z76QDsKtvx2bkfluL8Xrp/nd+d1S1e17o9dlRr7YtV9Z1JXp7kNUkqyVuTPKu19oVem4Nz+Xykd+Nv5H5IklesK9+S7vWXL66qz7TWXtNLc7C5iT87Bx/sW2urSf7kYo+vqnYR16xYt0cPWmtnkzyl7z7Yn3w+spe01r6Y5OCG2q25DG8igQvZyc/OwS/F2arWWt3dz/hQ6/aAfcXnI8DW7eRnp2C/QVV9XVW9sao+XFV/VFVvrqqHbnKodXvAvuLzEWDrdvKzU7D/Si3Jf2itfUNr7RFJ3pjkP68/YN26vfVroNav23v6TjULsIN8PgJs3Y59dlZrd7vsZ1+rqtkkr2utzfTdC8Bu4vMRYOsu52enGfu796x0rxQE4FzPis9HgK16Vi7TZ+fg34oziar6qSQPTjLfdy8Au4nPR4Ctu9yfnYL9eVTVTyZ5YpJ/PH5lJgDx+QhwKXbis1Ow38T4X1Nr/8N/tu9+AHYLn48AW7dTn50ent2gqh6W7gnkP0+y9m2et7fWZvvrCqB/Ph8Btm4nPzsFewAAGABvxQEAgAEQ7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHAIABEOwB2BZV9U+r6l1V9ZCqWq6qL1TV71fVQ/vuDWA/EOwBmFhVPTfJbyf5vSRPSvJPkpxI8g+SvKbH1gD2jWqt9d0DAHtYVX1jkg8k+Y+ttWdX1akkn2qt/fOq+qUkz0zy91prX+y1UYCBM2MPwKTmk9ye5AVVVUkeleQPxvs+Ph4P9tEYwH4i2AMwqccl+d3W2heSPDTJvZL84Xjf30/y1621z/XUG8C+IdgDMKnpJJ8c//lbxuNasD+Wbu09AJfZFX03AMCe99kk14z//C1JPtla++uqelqShyR5em+dAewjZuwBmNRvJDlaVT+Z5DuSfLKq/n2SVyd5fmvt93vtDmCfMGMPwKR+Ksl9kjw33fr629P9/8sPtdZ+rc/GAPYTr7sEYFtU1YOT/HmS72utWVcPsMMsxQFgu8yOx5VeuwDYpwR7ALbLbLoHZz/RdyMA+5GlOAAAMABm7AEAYAAEewAAGADBHgAABkCwBwCAARDsAQBgAAR7AAAYAMEeAAAGQLAHAIAB+P8B9mrvLKPqViEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "zp = branchpoint\n",
    "dz = np.array([0.5, 0.0])   # guide direction to find other branch\n",
    "for k in range(10):\n",
    "    ZS = root(sys, zp+dz, args=(zp,ds), jac=dsys, method='hybr') \n",
    "    zs =ZS.x\n",
    "    if d2Epot_dphi2(zs[0], zs[1]) > 0:\n",
    "        stable_phi.append( zs[0] )\n",
    "        stable_mu.append( zs[1] )\n",
    "    else:\n",
    "        nonstable_phi.append( zs[0] )\n",
    "        nonstable_mu.append( zs[1] )\n",
    "    dz = zs -zp\n",
    "    zp = zs\n",
    "\n",
    "plt.plot(stable_phi, stable_mu, 'bo')\n",
    "plt.plot(nonstable_phi, nonstable_mu, 'bo', fillstyle='none')\n",
    "plt.xlabel('$\\\\varphi$')\n",
    "plt.ylabel('$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = ['$-\\\\frac{1}{2}\\pi$', '$0$', '$\\\\frac{1}{2}\\pi$', '$\\pi$', '$\\\\frac{3}{2}\\pi$']\n",
    "plt.xticks(tick_pos, labels)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In accordance with our initial intuitive expectation we recognize the buckled state branching off when the critical load is exceeded.\n",
    "By symmetry there is also a mirrored branch (buckling to the other side). It is not shown here, but if you imagine it, then you may understand why this kind of bifurcation is called _pitchfork_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is much more to discover. In terms of physics it may be interesting to study the influence of an imperfection ($\\varphi_0\\ne \\pi/2$), whereas in terms of numerics it is desireable to automate the pseudo arc-length metric, to predict branched curves and to make sure all branches have been found. For further reading we recommend the textbooks _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_ by J. Guckenheimer and P. Holmes and _Elements of Applied Bifurcation Theory_ by Y. Kuznetsov."
   ]
  }
 ],
 "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": 4
}