{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 3.2: Newton incremental iteration scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from bmcs.api import PullOutModel\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate the idea of a non-linear solver even in simple terms, let us consider a known function defining the non-linear residuum. The principle assumption of the incremental time stepping scheme is that we know some initial value of the state variable that satisfies the governing equations, i.e. $R(u_0) = 0$.\n",
    "\n",
    "Then, we assume to know the derivatives of the sought function with respect to the state variable $u$. The time stepping algorithm provides us a means how to travel through the space of state variables along an admissible path satisfying the governing equations. The director of the travel is the pseudo-time variable $t$.\n",
    "\n",
    "In order to illustrate the concept, let us consider the simple function \n",
    "\\begin{equation}\n",
    "f(t)= \\mathrm{sin}(u)\n",
    "\\end{equation}\n",
    "Think of this equation as the equilibrium requirement of our mechanical model. On the left hand side, there is the prescribed history of loads $f(t)$. On the right hand side, the force of response of the structure. Now, for prescribed history on the right hand side $f(t)$, we want to find the corresponding history of displacements satisfying the equlibrium $u(t)$. \n",
    "\n",
    "Of course, we might solve this equation just by evaluating\n",
    "\\begin{align}\n",
    "u = \\arccos{f(t)}\n",
    "\\end{align}\n",
    "But this is not possible in a general case of a discretized structure. \n",
    "We thus pretend, that the inversion of the function $\\sin(u)$ is not available and define the residuum of the problem as\n",
    "\\begin{equation} \\label{eq:residuum}\n",
    "R = \\mathrm{sin}(u) - f(t) = 0.\n",
    "\\end{equation}\n",
    "In a numerical code, the values of the function $f(t)$ is given in incremental steps.\n",
    "For each increment of $f(t)$ an iteration loop must be performed to find the value of $u$ satisfying the residuum.\n",
    "\n",
    "To get the numerical algorithm, we first expand the residuum using the first two terms of the Taylor series as\n",
    "\\begin{align}\n",
    "R(u^{k+1},t) = R(u^{k},t) + \\left.\\frac{ \\partial R(u) }{ \\partial u }\\right|_{u^k} \\Delta u^{k+1} = 0.\n",
    "\\end{align}\n",
    "\n",
    "In the considered case of $\\sin{u}$, the derivative of the residuum with respect to $u$ is calculated as\n",
    "\\begin{equation}\n",
    "\\left.\\frac{ \\partial R(u) }{ \\partial u } \\right|_{u^k} = \\mathrm{cos}(u^k)\n",
    "\\end{equation}\n",
    "\n",
    "The iteration loop can then be obtained by solving the expanded residuum for the  increment of the displacement  $\\Delta u^{k+1}$ as\n",
    "\\begin{equation}\n",
    "\\Delta u^{k+1} = - \\left[ \\left.\\frac{ \\partial R(u) }{ \\partial u }\\right|_{u^k} \\right]^{-1} R(u^k)\n",
    "\\end{equation}\n",
    "The new value of the  variable $u$ can be given as\n",
    "\\begin{equation}\n",
    "u^{k+1} = u^k + \\Delta u^{k+1} \n",
    "\\end{equation}\n",
    "\n",
    "The last two equations are repeated until the residuum is (almost) zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sinus function as a pseudo model of a nonlinear structural response "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(u):\n",
    "    return np.sin(u)\n",
    "\n",
    "def df_du(u):\n",
    "    return np.cos(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithm implementation\n",
    "Note that the algorithm is implemented as a double loop. The outer loop defines the load levels. The inner loop iterates over the values of residuum until the residuum gets smaller than the required threshold value prescribing the accuracy of the calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of iterations = 0\n",
      "number of iterations = 3\n",
      "number of iterations = 3\n",
      "number of iterations = 3\n",
      "number of iterations = 4\n",
      "number of iterations = 10\n"
     ]
    }
   ],
   "source": [
    "K_max = 100 # max number of iterations\n",
    "u = 0.0 #  initial value of displacement\n",
    "f_levels = np.linspace(0, 0.99999, 6) # prescribed load levels\n",
    "u_levels = []  # calculated displacement values for each load level\n",
    "for f_level in f_levels: # loop over the load levels\n",
    "    for K in range(K_max): # iteration loop\n",
    "        R = f(u) - f_level # value of the residuum\n",
    "        if np.fabs(R) < 1e-8: # residuum equal to zero?\n",
    "            break # stop iteration\n",
    "        dR = df_du(u) # derivative of the residuum\n",
    "        d_u = - R / dR # increment of displacement\n",
    "        u += d_u # update the total displacement       \n",
    "        if K == K_max - 1:\n",
    "            raise ValueError('No convergence')\n",
    "    print('number of iterations =',K)\n",
    "    u_levels.append(u) # record the found solution for current load level f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Print the obtained values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u [0.0, 0.20135587954930018, 0.4115124817003838, 0.6434936070044599, 0.9272818847867171, 1.5663241841617594]\n",
      "f [0.       0.199998 0.399996 0.599994 0.799992 0.99999 ]\n"
     ]
    }
   ],
   "source": [
    "print('u', u_levels)\n",
    "print('f', f_levels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the obtained results with exact solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb9fbfca690>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3wUdf7H8deXkAqBAIEEQldq6ESxoShqQDxBxK4cCGJDvTvFdmf5nXeWwwIqiIhYTkVP5RAVCSoq2AGpoYOUJCaEEgiQsuX7+2NzGjHAEjaZ3c37+XjkYXZ32HkTMm/nMfPZGWOtRUREQl8tpwOIiEhgqNBFRMKECl1EJEyo0EVEwoQKXUQkTNR2asWJiYm2devWTq1eRCQkLVmyZKe1tnFFrzlW6K1bt2bx4sVOrV5EJCQZY7Ye7jUdchERCRMqdBGRMKFCFxEJEyp0EZEwoUIXEQkTRy10Y8x0Y8wOY8yqw7xujDHPGGM2GmNWGGN6BT6miEjom7U0m9Mfm0+bez7i9MfmM2tpdkDf35899FeAAUd4fSDQruxrDPD88ccSEQkvs5Zmc+/MlWQXFGGB7IIi7p25MqClftQ5dGvtAmNM6yMsMhh4zfquw/udMSbBGNPUWvtzgDKKiAQna6Gk8Nev0v1Qsg9K9vu+d5eA1w2eUrI/W8sobxGRtT0s9rZnobcbRS4P4zPWMaRnSkDiBOKDRSnA9nKPs8qe+12hG2PG4NuLp2XLlgFYtYhIFbAWDuTD3u2wN8v3tS8HDuz0PX8g/9fvvS6/3vIWgEjf95PdF7HQ2w2AnIKigMUORKGbCp6r8K4Z1tqpwFSAtLQ03VlDRJzj9foKe9fGX792boCCrbA3Gzwlv12+dizUbQx1GkO9ZtC0G8QlQlwjiKkH0fEQXQ9vZF3W7Pby+ZYiPl1fwLYCF25TG2pFctBjcBNB+dpslhAbsL9SIAo9C2hR7nFzICcA7ysiEhj7d0DuSshbBXmZvq+dG35b2lHxkHgiNO0OHS+E+i2gfnPfV0ILiEkAU9H+K5S6vXy3eRdzV+Xyyeo88gtLiIwwnH5iS644J5lzOyfx1Yad3DtzJW6X55c/FxsZwbj0DgH7awai0GcDY40xbwF9gL06fi4ijtm/A7IWQ/YSyPnRV+QH8n99Pb4ZJKVC236Q2A4anQiN2kHdJoct7IocLHWzYH0+GZl5fLomj8JiN3FREZzdoQnnpyZxdscm1IuJ/GX5/x0nH5+xjpyCIpolxDIuvUPAjp+DH4VujJkB9AMSjTFZwIOUHQmy1k4B5gAXABuBg8DIgKUTETkSjxt+Xg7bvoXsxZC1BPZu871mIqBJZ2h3PiR1geQuvv/GNaz06vYedPHZ2jzmrsplwYZ8il1eEuIiSU9NZkBqMme0SyQmMuKwf35Iz5SAFvih/JlyufIor1vKjveLiFQpdynkLIWtX8GWr2H7975pEoD6LaF5b+gzBlLSfIdOouKOe5U79hWTsTqPeZm5fLtpF26vJbleDJentSA9NZmT2zSkdkRwfEbTscvniogclbWQvxY2fgabPoOt34K7bCqkcSfofgW0Oh1anQbxyQFb7dZdB8jIzGXuqlyWbi/AWmiTWIfRfdsyoEsy3VLqU6uW/4dnqosKXUSCS9Ee2PS5r8A3fQ77yj54k9gBeg2H1mf4CrxOYsBWaa1lbW7hLyW+NrcQgNRm9fjLue1J75JMuyZ1McdwjN0JKnQRcd6eLbB2DqybA1u/AeuB6PrQ9iw462444RzfpEkAeb2Wpdv3kJHpOya+bfdBjIGTWjXkb4M6kZ6aTIuGx3/Ipjqp0EWk+lkLPy+DtR/5inxHpu/5xp3g9Nuh/QBI6Q0Rga0ol+e344U7ysYLTzshkZv6ncC5nZJoHB8d0HVWJxW6iFQPa30jhJkzIfO/vr1yUwtangrn/xM6DIRGJwR8tUWlHr5cn8+8zFw+XZPHvmI3sZERnN2xMempyb8bLwxlKnQRqVo71sCqmb4i37XRN07Y9izoewd0GAR1GgV8lXuLXMwvGy/8cr1vvLB+bCTndU5mQJdk+h5lvDBUqdBFJPD258PKd2D5m769cozvZOapt0CniwJ6QvN/dhQWMy8zj4xy44VJ9aK5rNx4YWSQjBdWFRW6iASGqxjWz4XlM2DDJ74Tm017wIDHIfViiE8K+Cq37Trom0zJzOXHbXuwFlo3imNU3zYMSE2me/OEoBwvrCoqdBE5PnmZsOQVWPE2FO+F+KZw2ljofiU06RTQVVlrWZdXyNxVuWRk5rHm530AdG5ajz/1b8+ALsm0Twr+8cKqokIXkWNXetB3YnPJy5C1CCKifIdSelzlu0ZKrcAdn/Z6LcuyCshYlUtGZi5bdvnGC3u3bBCy44VVRYUuIv7bsQYWT4flb0PJXt9Frc7/p29vPIAnN10eL99v3k1Gpq/EdxSWULuW4bQTE7n+zLac1zmJJvExAVtfuFChi8iReT2+Y+PfT4GfFvj2xjsPht4jfB+7D9DhjWKXhwXr85mbmctna3awt8hFbGQEZ7VvzIAuvvHC+rHhMV5YVVToIlKxoj2w9HX4YSoUbIN6zaH/g9DrjwHbG99b5OLztTvIyMzli3X5FLk81I+NpH+nJqSnJnNmu8bERoXfeGFVUaGLyG/t2gTfTvJNq7gO+vbCz/+Hb2Y8AJ/c3FFYzCer88jIzOPbTTtxeSxN4qMZ1rs56anJ9Gkb/uOFVUWFLiI+WYvh64mw5gOIiISul0GfG3y3WjtO23cf/OXCV0vKxgtbNYrjutPbkN4lmR41bLywqqjQRWoyrxc2fuIr8q1fQ0x96PsXOPmG45obt9ayPm9/2XhhLqvLxgs7lY0XpndJokNSfI0dL6wqKnSRmsjj9o0dLnwS8tf4jo+nP+K7PG10fKXe8pfxwsxcMlb9Ol7Yq2UD/nqBb7ywZSONF1YlFbpITeJx+T4AtPBJ2L3Zd4u2i6dCl6G+wyzHyOXx8sNPv44X5u3zjReeekIjRvdty/mdk2hST+OF1UWFLlITuEtg2Rvw1dO+iZXkbnD5674TnbWO7QRkscvDwg07mbsql8/W5lFw0EVMZK1fxgvP6ZBE/TiNFzpBhS4Sztyl8OOrsPApKMzx3Wvzgid8N04+huPX+4p/O154sNRDvZjanNspifNTkzmrvcYLg4EKXSQcedy+scMv/wV7t/muOT5ksu9j+X4WeX5hCZ+u8V2C9puy8cLG8dFc3DOFAV2SOaVtI40XBhkVukg48Xp91x3//BHYvQma9YQ/PA0n9PeryP83XjgvM49FW3djLbRsGMfI09uQnppEzxYNNF4YxFToIuHAWt/9OOf/A3ashiapcMWb0OGCIxa5tZYNO/aTscp3CdrMHN94YcfkeG47px0DuiTTMVnjhaFChS4S6rb/APPuh+3fQaMT4ZKXIHXoYU92er2WFdl7mbsql3mZuWzeeQCAXi0TuO+CjqSnJtOqUZ3q/BtIgKjQRULVzg3w6UOw9kOomwQXPg09h1f48Xz3b8YL88jdV/zLeOHIM9pwfuckkjReGPJU6CKhpjAPvnwMlrwKkbFw9l/hlJshuu5vFit2efhqw07mlt0c+X/jhWe2a8xdXTrQv6PGC8ONCl0kVLiK4btJvhFEdzGcNArOvAvqNv5lkcJiF/PX7mBeZh6fr9vBwVIP8WXjhempSZzZvjFxUdrsw5X+ZUWCnbWw+n345H7fh4I6Xgjn/R0anQDAzv0lfLo6j7mZuXyzcRelHi+JdaMZ0jOFAam+8cKo2hovrAlU6CLBLGcZZNznu3BWUhcYPhvankXWnoNkfPUTGZm5LN6yG6+FFg1j+eNprUhPTaZnywZEaLywxlGhiwSj/fnw2UOw9A2Ia4S9cCIbU4aQsSafuR8tZFX2r+OFY89px4DUZDo11XhhTadCFwkmHjcsfgnm/xPrOkh+1+t5I/pyPvhyP5vzvwagZ8sE7h3oGy9snajxQvmVX4VujBkATAQigGnW2scOeb0+8DrQsuw9n7DWvhzgrCJhZ9bSbMZnrCOnoIgB9X7i0ZhXSdi3nk3xJ3GfdzjfL2pERK0dnNq2ESNPa815nZNJrq/xQqnYUQvdGBMBTALOA7KARcaY2dba1eUWuwVYba39gzGmMbDOGPOGtba0SlKLhIFZS7O5d+ZK6rp28UTkm1xS+hXZJY242/UnvtjThzPbN+HJ1GT6d2pCQlyU03ElBPizh34ysNFauxnAGPMWMBgoX+gWiDe+A3h1gd2AO8BZRcLKk3NXc6l3DndG/4doXDzrHsJk90XExsWz9J5zNF4ox8yf35gUYHu5x1lAn0OWeQ6YDeQA8cDl1lrvoW9kjBkDjAFo2bJlZfKKhIX9W3/k2aK76BG5mQWerjzgHsEW2xSA4oMulblUij+/NRWdNreHPE4HlgHnACcAnxhjFlpr9/3mD1k7FZgKkJaWduh7iIS/kv3kzX6QRpnTSTHx3Fo6lg+8p1J+M2uWEOtcPglp/nzaIAtoUe5xc3x74uWNBGZan43AT0DHwEQUCQ+eNR9R+FQvkjKn8WHEubxzynt8GtGX8mUeGxnBuPQOzoWUkObPHvoioJ0xpg2QDVwBXHXIMtuA/sBCY0wS0AHYHMigIiGrMI+D7/+FuI0fku1twdzWz3HdVVdQLyaSZsm/Trk0S4hlXHoHhvRMcTqxhKijFrq11m2MGQtk4BtbnG6tzTTG3Fj2+hTgYeAVY8xKfLsbd1trd1ZhbpHgZy0sn0HpR/cQ4TrIRHsFrQbfw596t/llkSE9U1TgEjB+nXmx1s4B5hzy3JRy3+cA5wc2mkgI27MV9+zbqf3T5yz3tufVxDu46+qLaNkozulkEsZ0Kl0kkLxeWPQink8eotTt5WHXCOqfeSNPn9tB99+UKqdCFwmUXZuws27BbP+Wr73dmBBzM/f8MZ2T2zR0OpnUECp0keNVtlduP3mQg55aPFB6I8Wpl/Hyxd10AwmpVip0keOxZwu8Pxa2LOQrenK/dww3D+3Lpb2b68qHUu1U6CKVYS0sno6ddz8lHssDrutZmzyYl6/sRRtdAVEcokIXOVb7cmDWzbD5c5ZEdOf2olFcdFYf3j23ve4MJI5SoYsci1UzsR/+GXdpMf/wXEdG5CCeGt2D005IdDqZiApdxC/Fe2HOOFjxNpsiO3B90Rjad+7Bx0O70aCOLm0rwUGFLnI0W76C/96Id18OL5jLeK54MH8d0o0rT26hE58SVFToIofjLoX5D2O/eZZd0SmMLn4QV9NevH9FT05sUtfpdCK/o0IXqcjODfDeKPh5OR9FDmDc3su4pm8n7kzvQHTtCKfTiVRIhS5SnrWw9HXsx3dRQiR/dt/J4shTmTqqO33bNXY6ncgRqdBF/qdoD3zwJ1g9i9XRPbhu72i6durI3Eu60ahutNPpRI5KhS4CsPVbmHk93n0/M6nW1Uw6MIi/Dk7lmlNa6cSnhAwVutRsXg8sfBL7xaPsiWrKiOIHKU3qwewre9I+Kd7pdCLHRIUuNVdhHswcDT8t4PPIfty691ouPa0T9wzsSEykTnxK6FGhS820aT525hg8RYXc77mRTyL689yIHpzdsYnTyUQqTYUuNYvHDV88il34JDmRrRhRNI5m7Xry8aXdaRyvE58S2lToUnPszYb3RsO2b3jf9OeBg8P506DujDitNbVq6cSnhD4VutQMm+Zj3xuNq/gg40pvZnXiAN66oiedm9VzOplIwKjQJbx5vbBgPPaLR9laqwWjiu7h1D6nMvuCzsRG6cSnhBcVuoSvA7uwM6/HbPqM9719edyM4e/X9uG8zklOJxOpEip0CU/bF+H9z3A8hfk84BrF9taXMevyHiTVi3E6mUiVUaFLeLEWfngRb8Z9/OxtyM2u/+PC9IH884w2OvEpYU+FLuGj9CCe2bcRseod5nt68Uz9O3jkyr50SanvdDKRaqFCl/CwezMlb1xN5K41jHddxp5et/DWH7oQF6Vfcak59NsuIc+uz8D1zmiKSz38ydzH4CuvZUCXpk7HEql2KnQJXV4vxZ89StTX49ngbcmU5P/jgasH0LR+rNPJRByhQpfQVFRAwRsjSciaz0xPX/L7PcaEs1OJ0IlPqcFU6BJyXLlrKHzlMuKLspkQPYazr7mPoS0bOB1LxHF+FboxZgAwEYgApllrH6tgmX7ABCAS2GmtPSuAOaWGmrU0m/EZ68gpKKJZQiwPtt/KGSvvxeONZGqbCYy+8irqRmu/RAT8KHRjTAQwCTgPyAIWGWNmW2tXl1smAZgMDLDWbjPG6BqkctxmLc3m3pkrKXJ5MHi5uPBNzl3+Lqtpzc8DXuKWU3s7HVEkqPiza3MysNFauxnAGPMWMBhYXW6Zq4CZ1tptANbaHYEOKjXP+Ix1FLk8xFHMk5HPMzBiETM9ZzAh5mYWqMxFfsefQk8Btpd7nAX0OWSZ9kCkMeYLIB6YaK197dA3MsaMAcYAtGzZsjJ5pQbJKSiiudnBtMgnaWeyeNh1DS95BmJcTicTCU7+FHpFYwO2gvfpDfQHYoFvjTHfWWvX/+YPWTsVmAqQlpZ26HuI/Mb5dTbyqPtf1MIy3HUPX3u7AtAsQWOJIhXxp9CzgBblHjcHcipYZqe19gBwwBizAOgOrEekEn6a+yzPuf+PrTaJUa472WqTAYiNjGBcegeH04kEp1p+LLMIaGeMaWOMiQKuAGYfssz7QF9jTG1jTBy+QzJrAhtVagSPi59eu4k23/2NH2v34Jtz3sJdvw0GSEmI5dGhXRnSM8XplCJB6ah76NZatzFmLJCBb2xxurU20xhzY9nrU6y1a4wxc4EVgBffaOOqqgwu4cce3E32i5fTZs8PzK5zCWfdNJk+dWMY3s/pZCKhwVjrzKHstLQ0u3jxYkfWLcHHnbeWgpcuIb4kl/80vZPLRt9FdG3dUUjkUMaYJdbatIpe0ycyxHHF6z7F+/ZwrCeCt1Mnc82wy3TtcpFKUKGLowq/mkrcp3ez0ZvCqn5TGX7OaU5HEglZKnRxhtdDwft3kbB8Gl/aHngumc4l3U9wOpVISFOhS/UrKaTg38NJyJrPm2YQnUc+Q49WiU6nEgl5KnSpXgXb2ffyJdQt2MCEmBsZOuZBWjaKczqVSFhQoUu1sdlLKXp1GJQc4JEGf+fW62+gQZ0op2OJhA0VulQLz5oP8bwzit2eukxvNYm7rh1CTKTGEkUCSYUuVa70q0nU/vSvrPa24fOez/K3wadrLFGkCqjQpep4PRR9cBexS6cx13MSO857hj+f2dnpVCJhS4UuVaNkPwdn/JG4LZ/ykncQKZeNZ3hXXYNFpCqp0CXwCvM4+MpQonet5hEzmvRRf6N3K93zU6SqqdAlsPLXUfTyxXBgJ/fF3MeN199Mm8Q6TqcSqRFU6BI4W7+h5N+Xc8AFjzZ4nPtGX0WjutFOpxKpMVToEhDelTPxzhxDlieRaa3G849rLyA2SmOJItVJhS7Hx1pcXz1D5GcPsNjbgc97TOQfQ04hQmOJItVOhS6V5/VQ/OHdxPz4Ih96+pDffyJ3ndURY1TmIk5QoUvluIo5+PYo4jZ+yHTvIJKGjWdkd40lijhJhS7HrqiA/a9eRt3c73nC/JGzrnuIk1o3dDqVSI2nQpdjszebwulDiC7YxEPRf+Ha6+/ghMZ1nU4lIqjQ5VjsWMOB6YOhaB9/r/8wt18/msbxGksUCRYqdPGLd8s3lPz7Mva7I3g+ZQL3jRhGXJR+fUSCibZIOSpX5mx49zpyPIn8N/UZ/jbsXGpH1HI6logcQoUuR3Tw2+nEZNzBcm9blp05lTv699JYokiQUqFLxaylYN6jJHz7OF96e1B40TRGprVzOpWIHIEKXX7P62Xnu38mcfUrfEBfmgyfxlknJjudSkSOQoUuv+UuIe+1kSRt+4gZERdx0pjnODGpvtOpRMQPKnT5VUkhP08dRtNd3zE9biQX3vgYTerFOJ1KRPykQhcA7IGd5E7+A433r2Va4p1cOeZe6kTr10MklGiLFUp3b2f3lEE0LMnhzdb/ZMTwGzWWKBKCVOg1XGHWaopfHkycu5A5PSZx7ZDLNJYoEqJU6DVY/vrviZxxKbW8lsX9/s3FZ5/ndCQROQ5+FboxZgAwEYgApllrHzvMcicB3wGXW2vfDVhKCYhZS7MZn7GOnIIizq+zgSc9j7HP1iF38Fuc0+skp+OJyHE66oFSY0wEMAkYCHQGrjTGdD7Mco8DGYEOKcdv1tJs7p25kuyCIs6ttZhn3P8gx9uQT099jV4qc5Gw4M+Zr5OBjdbazdbaUuAtYHAFy90KvAfsCGA+CZDxGesocnm4uNZCno+cwBrbkstKH+CFZSVORxORAPGn0FOA7eUeZ5U99wtjTApwMTDlSG9kjBljjFlsjFmcn59/rFnlOOQUFHFtxDyejnqe772duLr0PgqIJ6egyOloIhIg/hR6RSMP9pDHE4C7rbWeI72RtXaqtTbNWpvWuHFjfzPK8bKWP8d8wMORrzDP05vrXOM4QCwAzRJinc0mIgHjz0nRLKBFucfNgZxDlkkD3iobd0sELjDGuK21swKSUirPWla8fBu3MYP/es5gnGsM7rJ/9tjICMald3A4oIgEij976IuAdsaYNsaYKOAKYHb5Bay1bay1ra21rYF3gZtV5kHA6yHzhRF02/YaX9QfjOeiySQlxGOAlIRYHh3alSE9dWNnkXBx1D10a63bGDMW3/RKBDDdWptpjLmx7PUjHjcXZ1h3CeumXE3qzk+Y1+gazrnpGWrXjmDYSa2cjiYiVcSvOXRr7RxgziHPVVjk1toRxx9LjoctPcimycPoWPA1HyXdxIAbHiWilj79KRLu9EnRMGNLCtny7EW0LVzK7BZ3cuF1f6OWylykRlChhxHvgT1kTRpEiwNrmN32AQYP/7OuyyJSg6jQw4SncAe5ky4guWgzH3Z4lMFX3qAyF6lhVOhhwLVnO7uev4CGJbnM6TqBIcOGOx1JRBygQg9xpfk/sfeFgdRxFTCv12SGDL7U6Ugi4hAVeggryV3LgRcHEeku4os+LzL4gj84HUlEHKRCD1HFWSspmX4hXo+X7854lT+cp2uZi9R0KvQQdHDLYjyvDqHIW5vlZ7/BBf3OdDqSiAQBFXqI2b/ha8ybw9jrjWPd+W+QfvopTkcSkSChQg8hhWvmU/vtK8mzCWwZNIP+J/dyOpKIBBEVeojYu/JjYt4bzjbbhLwhb9GvZ1enI4lIkFGhh4CCpe8T9/51bLTN2TfsP5zRVZe8FZHfU6EHud2L3qHeRzewxram5Ip3OKXTCU5HEpEgpUIPYru+fYP6GWNZadvBNe+Q1k6XvhWRw/PnBhfigB0LXqJBxi0spRNRI2bRU2UuIkehQg9CeZ8/T5P5f+EHuhJ/3X9JbdPM6UgiEgJ0yCXI/DzvGZp+cz9fmd40vf5tTmimm2mLiH9U6EEk5+MnaPb9w3xZ62Ra3vAf2iQ1cDqSiIQQFXqQyPrwEZovfpz5EafR7sa3aNG4vtORRCTEqNCDwPZZ/0eLZU/xWe0z6XTzmzRrGO90JBEJQSp0J1nL1pn302rls3wSeTbdb3mDJgl1nE4lIiFKhe4Ua/npP/fQZs0U5kWdS++x/6ZRvTinU4lICFOhO8FaNs+4g7brX2JuzABOufVVEurEOJ1KREKcCr26WcvG12/nxE2v8nHsIE6/7WXqxUY7nUpEwoAKvTpZy/rXxtL+p9f5OO4i+t42nboxkU6nEpEwoUKvLtay7pWb6LB1BnPrXky/W6cRG60fv4gEjhqlOni9rJk+hk5Z7zC33qX0GzuFmCj96EUksNQqVc3rJfPFUaT+PJOMhCs455bJREVGOJ1KRMKQCr0KWa+HVS+MpGve+2Q0vIb+Nz9D7doqcxGpGir0KmK9HlY+P4Ju+bP5JHE45940kYgIXdxSRKqOX4VujBkATAQigGnW2scOef1q4O6yh/uBm6y1ywMZNBTMWprN+Ix1/FxwgCeipzHUfMFnSSPpP+YpaqnMRaSKHbXQjTERwCTgPCALWGSMmW2tXV1usZ+As6y1e4wxA4GpQJ+qCBysZi3N5t6ZKylxuXi89lSGmgVM9FxCyz53qsxFpFr40zQnAxuttZuttaXAW8Dg8gtYa7+x1u4pe/gd0DywMYPf+Ix1lLhc/CtyKpfWXsAE91Cedl3CE/PWOx1NRGoIfwo9Bdhe7nFW2XOHMwr4uKIXjDFjjDGLjTGL8/Pz/U8ZAnILDjA+8gWGRSzgadclTHAPAyCnoMjhZCJSU/hT6KaC52yFCxpzNr5Cv7ui1621U621adbatMaNw+dOPF63m6eiX+CSiIU85RrGRM8lv7zWLCHWwWQiUpP4c1I0C2hR7nFzIOfQhYwx3YBpwEBr7a7AxAt+XrebZc9dxWCzkKfdl/KM5+JfXouNjGBcegcH04lITeLPHvoioJ0xpo0xJgq4AphdfgFjTEtgJnCttbbGHDT+X5n3KsjgqxY30HroQ6QkxGKAlIRYHh3alSE9j3R0SkQkcI66h26tdRtjxgIZ+MYWp1trM40xN5a9PgV4AGgETDbGALittWlVF9t51uNm+aSr6VWQwcLmYzjjuscxxnBxrxp3PlhEgoSxtsLD4VUuLS3NLl682JF1Hy/rcbNs0tX03D2XBc3H0HfUvyj7H5mISJUyxiw53A6zBqSPkW/P/FpfmadcrzIXkaChQj8G1uNm+eRr6bF7Dl80G03f0eNV5iISNFTofrIeNyueH06PXXP4vOlozrr+CZW5iAQVFbofrNfDiuf/SPedH/F58ij6jVGZi0jwUaEfha/MR9B954fMT7qOfjc8qTIXkaCkQj8C6/WwcspIuufPZn7SCPqNUZmLSPBSoR+O18vKF0bRbcf7zG8ynH5jntZVE0UkqKmhKuL1smLqKLrl/Zf5idfQ74aJKnMRCXpqqUNZy4oXx9AtdyafJV5Nv5ueVZmLSEhQU5VnLStevIFuP7/D/EZX0O+m51TmIhIy1Fb/Yy0rpt1Et5y3+bzh5Zx18/O6B6iIhBQ1FvjKfPotdMuewRcNhnHmLVNU5iISctRa1rLi5dvotmK0SS0AAAeSSURBVP0Nvki4hDNumaoyF5GQVLOby1pWvPInum17jS/rD+GMsS9Su3aE06lERCql5ha6tax49S902/oKX9a7iNNuna4yF5GQVmMLffm/76LbluksqDeI0257mUiVuYiEuBpZ6Mtfv5fum6eyMH4gp9z6GpG1/bm1qohIcKtxhb78jb/RfeNkvq6bTp/bXicqUmUuIuGhRhX68hkP0H3Ds3xT51zSVOYiEmZqTKEvf/vvdF83kW/jzqHXbTOIjopyOpKISEDViEJf9s4jdF/zJN/F9aPn7W8RE60yF5HwE/aFvuzdx+mR+Tg/xPalx+3/ISY62ulIIiJVIqwLffnMJ+ix6hEWxZxOt9vfVZmLSFgL20JfPmsC3Vc8zJLoU+hy+3vExMQ4HUlEpEqFZaEvn/0s3Zc9yI/RJ9Hp9pnExsY6HUlEpMqFXaEv/3AyXZfcz7Ko3nS8fRZxcXWcjiQiUi3CqtCXf/QCXRfdx8roHrS77X3i4uo6HUlEpNqETaEvnzudLj/cTWZUN064bTZ16sY7HUlEpFqFRaEvz3iV1G/vYG1UKq1vnU3duvWcjiQiUu1CvtBXfPJvOn/zZ9ZHdqTFrR8SXy/B6UgiIo4I6UJf8dkMOn51O5si25Ey9iPq1WvgdCQREcf4dXUqY8wAYCIQAUyz1j52yOum7PULgIPACGvtjwHOyqyl2YzPWEdOQRF/iFvFeM+/2FK7LU1vmUP9hIaBXp2ISEg56h66MSYCmAQMBDoDVxpjOh+y2ECgXdnXGOD5AOdk1tJs7p25kuyCIvrWWs54z3jW2RYsP/tl6jdoFOjViYiEHH8OuZwMbLTWbrbWlgJvAYMPWWYw8Jr1+Q5IMMY0DWTQ8RnrKHJ5OLVWJlMjn2Kjbca1pfcy4av8QK5GRCRk+VPoKcD2co+zyp471mUwxowxxiw2xizOzz+2Is4pKAJgh03ge28nrim9l73U/eV5EZGazp9CNxU8ZyuxDNbaqdbaNGttWuPGjf3J94tmCb6P72+yKfzRdQ97qPeb50VEajp/Cj0LaFHucXMgpxLLHJdx6R2IjfztjZxjIyMYl94hkKsREQlZ/hT6IqCdMaaNMSYKuAKYfcgys4HhxucUYK+19udABh3SM4VHh3YlJSEWA6QkxPLo0K4M6fm7IzsiIjXSUccWrbVuY8xYIAPf2OJ0a22mMebGstenAHPwjSxuxDe2OLIqwg7pmaICFxE5DL/m0K21c/CVdvnnppT73gK3BDaaiIgci5D+pKiIiPxKhS4iEiZU6CIiYUKFLiISJozvfKYDKzYmH9hayT+eCOwMYJxAUrbKUbbKUbZjF6y5wL9sray1FX4y07FCPx7GmMXW2jSnc1RE2SpH2SpH2Y5dsOaC48+mQy4iImFChS4iEiZCtdCnOh3gCJStcpStcpTt2AVrLjjObCF5DF1ERH4vVPfQRUTkECp0EZEwEdSFbowZYIxZZ4zZaIy5p4LXjTHmmbLXVxhjegVRtqvLMq0wxnxjjOkeLNnKLXeSMcZjjBkWTNmMMf2MMcuMMZnGmC+DIZcxpr4x5gNjzPKyXFVyRdHDZJtujNlhjFl1mNed3A6Ols3J7eCI2cotV63bgT+5Kr0NWGuD8gvfpXo3AW2BKGA50PmQZS4APsZ3x6RTgO+DKNtpQIOy7wcGU7Zyy83HdxXNYcGSDUgAVgMtyx43CZJc9wGPl33fGNgNRFXTz+1MoBew6jCvO7Id+JnNke3An2zl/u2rezs42s+s0ttAMO+hB8XNqSubzVr7jbV2T9nD7/Ddxak6+PNzA7gVeA/YUU25/M12FTDTWrsNwFpbHfn8yWWBeGOMAeriK3R3NWTDWrugbH2H49R2cNRsDm4H/vzcwIHtwI9cld4GgrnQA3Zz6ipwrOsdhW8PqjocNZsxJgW4GJhC9fLn59YeaGCM+cIYs8QYMzxIcj0HdMJ3a8WVwO3WWm81ZPOHU9vBsarO7eCoHNwOjqbS24BfN7hwSMBuTl0F/F6vMeZsfL/IZ1RponKrrOC5Q7NNAO621np8O5zVxp9stYHeQH8gFvjWGPOdtXa9w7nSgWXAOcAJwCfGmIXW2n1VmMtfTm0HfnNgO/CHU9vB0VR6GwjmQg+Km1Mfhl/rNcZ0A6YBA621u6ohl7/Z0oC3yn6JE4ELjDFua+2sIMiWBey01h4ADhhjFgDdgaosdH9yjQQes76DmhuNMT8BHYEfqjCXv5zaDvzi0HbgD6e2g6Op/DZQXScoKnHioDawGWjDryeqUg9ZZhC/PRn0QxBla4nvHqunBdvP7ZDlX6H6Tgb583PrBHxWtmwcsAroEgS5ngceKvs+CcgGEqvx37U1hz+J5sh24Gc2R7YDf7Idsly1bQd+/MwqvQ0E7R66DaKbU1cy2wNAI2By2R6A21bDFd78zOYIf7JZa9cYY+YCKwAvMM1ae8Sxs+rIBTwMvGKMWYmvOO+21lbLJViNMTOAfkCiMSYLeBCILJfNke3Az2yObAd+ZnPE0XIdzzagj/6LiISJYJ5yERGRY6BCFxEJEyp0EZEwoUIXEQkTKnQRkTChQhcRCRMqdBGRMPH/0t0iBh+fHSUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(u_levels,f_levels,'o-')\n",
    "u_analytical = np.linspace(0, u_levels[-1], 50)\n",
    "plt.plot(u_analytical, f(u_analytical))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Questions and tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 What if the derivative is not available?\n",
    "Is it possible to trace the $\\sin(u)$ function even without the derivative? Sometimes it is difficult to obtain the derivative of the function so that other strategies are necessary. Examine the algorithm behavior by setting a value of the derivative, i.e.  $\\partial f(u) / \\partial u = \\cos(0)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 Can the algorithm trace the descending branch of the function?\n",
    "Can the present algorithm trace the $\\sin(u)$ function also in the range $u > \\pi/2$. \n",
    "Why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4 Can the algorithm handle an infinite derivative?\n",
    "Use the algorithm above to trace the function $\\sqrt{u}$. Why doesn't it work? How to fix it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(u):\n",
    "    return np.sqrt(u)\n",
    "\n",
    "def df_du(u):\n",
    "    return 1./2./np.sqrt(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 Can the algorithm handle a concave function?\n",
    "Use the algorithm above to find the solution to the functions $u^2$. Why doesn't it work? How to fix it? Use the algorithm in combination with the function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(u):\n",
    "    return u**2 + .5*u\n",
    "\n",
    "def df_du(u):\n",
    "    return 2*u + .5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5 Visualization of the iteration process\n",
    "Apply an extnded version of the algorithm that records the iterations staps showing what the algorithm went through. Use the modified Newton method - using a constant derivative to see the all the loops that the algorithm went through. Don't forget to increase the number of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of iterations = 0\n",
      "number of iterations = 5\n",
      "number of iterations = 4\n",
      "number of iterations = 4\n",
      "number of iterations = 3\n",
      "number of iterations = 3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfZzNdf7/8cdrxjXJVaXIDLbkIiqyJGNOJbShFmEnKjERlVrkZtxqaXXl1q1WUY1WtTXq2279tkR0xbieZohEQjJoal2Vi7BjzPv3x4yzQzMM5nM+Z+Y877fbud3mnPfHOc9zzJzn5+q8jznnEBGRyBXldwAREfGXikBEJMKpCEREIpyKQEQkwqkIREQiXDm/A5yuOnXquNjYWL9jiIiUKitWrNjlnDuvsLFSVwSxsbFkZGT4HUNEpFQxs8yixrRrSEQkwqkIREQinIpARCTCqQhERCKcikBEJMKpCEREIpyKQEQkwqkIRETCnHMOL78ywLMiMLMZZrbDzL4uYtzMbIqZbTKzr8zsKq+yiIiUNikpKcTGxmJmVKlShXbt2nn2WF5uEbwGdD3JeDfgkvxLIvCih1lEREqNlJQUEhMTyczM+zDw4cOHSU9PJyUlxZPH86wInHMLgT0nWaQn8A+XZzlQw8wu9CqPiEhpkZSUxMGDB4+7zTlHUlKSNw94bN+TFxcgFvi6iLEPgWsLXP8MaFPEsolABpARHR3tgEIvzjkXExPjYmJinMv7h4VeCo7PnTvX3XbbbUUue2y8QoUKzjnnKlSoUOhyBcdvu+02N3fu3CLvs+D4scx6TnpOek56Tqe6mJk7U0CGK+K92pyHByDMLBb40DnXopCx2cATzrnF+dc/A8Y451ac7D7btGnjNOmciJRlsbGxwd1CBcXExLBly5Yzuk8zW+Gca1PYmJ9nDW0HLi5wvT6Q5VMWEZGw8dhjjxEdHX3cbVWqVGHSpEmePJ6fRfABMDD/7KF2wF7n3I8+5hER8V12djazZs3i6NGjnHvuuUDelkBycjIJCQmePKZn30dgZm8B8UAdM9sOPAqUB3DOvQTMAW4CNgEHgbu8yiIiUhocPHiQXr16MXfuXCZPnsyoUaNC8rieFYFzrv8pxh0w3KvHFxEpTfbu3cvNN9/MkiVLSE5OZsiQISF77FL3DWUiImXNzp076dKlC2vWrOGtt96ib9++IX18FYGIiI+2b99O586d2bJlC++//z433XRTyDOoCEREfLJp0yZuuOEG9uzZw7x584iLi/Mlh4pARMQHa9asoXPnzhw9epT58+fTunVr37Jo9lERkRBbvnw5nTp1oly5cixcuNDXEgAVgYhISH322WfccMMN1KpVi8WLF9O0aVO/I6kIRERC5djB4IYNG7Jo0SJiY2P9jgSoCEREQuKNN96gV69eXHnllaSmpnLhheEz2bKKQETEY1OnTmXgwIF06tSJTz/9lFq1avkd6TgqAhERjzjnePzxxxkxYgQ9e/Zk9uzZVKtWze9Yv6EiEBHxgHOOhx9+mKSkJG6//Xb++c9/UqlSJb9jFUqfIxARKWFHjx5l2LBhTJ8+neHDhzNlyhSiosJ3vTt8k4mIlELZ2dkkJCQwffp0kpKSeP7558O6BEBbBCIiJebgwYP06dOHOXPm8PTTTzN69Gi/IxWLikBEpATs27eP7t27s2jRopBPI322VAQiImdp165ddO3aldWrV/syjfTZUhGIiJyFH374gc6dO/P999/7No302VIRiIicoU2bNtG5c2d2797t6zTSZ0tFICJyBtasWcONN97IkSNHfJ9G+myF9zlNIiJhKC0tjU6dOhEVFRUW00ifLRWBiMhp+Pzzz7n++uuD00g3a9bM70hnTUUgIlJMJ04j3bBhQ78jlQgVgYhIMbz55pv06tWLVq1ahd000mdLRSAicgrTpk1jwIABYTuN9NlSEYiIFME5xxNPPMHw4cPp0aMHs2fP5pxzzvE7VolTEYiIFMI5x9ixYxk3bhwJCQn861//CttppM+WPkcgInKCo0ePcu+995KcnMy9995bKmYQPRtl95mJiJyBI0eOcPvtt5OcnMy4ceN44YUXynQJgLYIRESCDh06RJ8+fZg9ezZPPfUUY8aM8TtSSKgIREQ4fhrpl19+mcTERL8jhYyKQEQiXsFppGfOnEm/fv38jhRSKgIRiWhpaWm0a9cOgA8//JA//OEPPicKPU+PgJhZVzP71sw2mdnYQsbPNbNZZrbazNaa2V1e5hERSUlJITY2lqioKC644IJgCfzxj3+MyBIAD7cIzCwamAp0BrYD6Wb2gXNuXYHFhgPrnHPdzew84FszS3HOZXuVS0QiV0pKComJiRw8eBCAHTt2BMc6derkVyzfeblF0BbY5JzbnP/G/jbQ84RlHHCOmRlQDdgD5HiYSUQiWFJSUrAE5H+8LIJ6wLYC17fn31bQC0BTIAtYAzzgnMs98Y7MLNHMMswsY9WqVZhZoReA2NhYYmNjj/27Qi8Fx+fNm0ffvn2LXPbYeMWKFQGoWLFiocsVHO/bty/z5s0r8j4Ljh/LrOek56Tn5P1z2rp164lvL0Fr164tcqysM+ecN3ds1gfo4pwbnH99ANDWOXdfgWV6Ax2Ah4DGwCdAK+fcvqLut02bNi4jI8OTzCJStsXGxpKZmfmb2y+88EKysrJ8SBQ6ZrbCOdemsDEvtwi2AxcXuF6fvDX/gu4C3nN5NgHfA5d5mElEItikSZOoXLnycbdVrlyZyZMn+5QoPHhZBOnAJWbW0MwqAP2AD05YZitwPYCZXQA0ATZ7mElEIljPnj1p3rx58HqDBg2YPn06CQkJPqbyn2dnDTnncsxsBDAPiAZmOOfWmtnQ/PGXgMeA18xsDWDAw865XV5lEpHItWXLFnr06MHatWt57rnnuP/++4PHFSKdpx8oc87NAeaccNtLBX7OAm70MoOISGpqKr179yYnJ4ePPvqIG2/U205BZXtKPRGJeMnJydxwww3Url2btLQ0lUAhVAQiUiYdOXKEESNGcM8999C5c2fS0tK49NJL/Y4VllQEIlLm7N69m65duzJ16lRGjRrFrFmzOPfcc/2OFbY06ZyIlClr166lR48ebN++nddff52BAwf6HSnsqQhEpMyYNWsWCQkJVK1aldTU1OCEcnJy2jUkIqWec44nn3ySnj17cumll5Kenq4SOA3aIhCRUu3QoUMMHjw4+IUyM2bM+M2nh+XktEUgIqXWDz/8QFxcHG+99RaPP/44M2fOVAmcAW0RiEiplJaWxq233sr+/fv597//TY8ePfyOVGppi0BESp0333yTTp06UalSJZYtW6YSOEsqAhEpNY4ePcqYMWMYMGAA7du3Jz09nRYtWvgdq9TTriERKRX27t3Ln/70J+bMmcO9997Lc889R/ny5f2OVSaoCEQk7G3cuJEePXqwadMmXnzxRYYOHep3pDJFRSAiYe3TTz/ltttuIyoqik8++YT4+Hi/I5U5OkYgImHJOceUKVPo2rUr9erVIz09XSXgERWBiISd7OxshgwZwgMPPMDNN9/M0qVLadiwod+xyiwVgYiElR07dnD99dfz97//naSkJN577z3OOeccv2OVaTpGICJhY9WqVfTs2ZOdO3fy9ttv07dvX78jRQRtEYhIWHj33Xfp0KEDubm5LFq0SCUQQioCEfFVbm4uEyZMoHfv3rRs2ZL09HRat27td6yIol1DIuKbX3/9lTvuuIN3332XO+64g5dffpmKFSv6HSviqAhExBeZmZn07NmTNWvW8Mwzz/Dggw9iZn7HikgqAhEJucWLF/PHP/6R7OxsZs+eTdeuXf2OFNF0jEBEQurvf/871113HTVr1iQtLU0lEAZUBCISEjk5OTzwwAMMHjyYQCDA8uXLadKkid+xBBWBiITAnj176NatG1OmTOHBBx9k9uzZ1KxZ0+9Ykk/HCETEU9988w09evRg69atzJgxg7vuusvvSHICFYGIeGbOnDn079+fSpUqMX/+fK655hq/I0khtGtIREqcc47Jkydz880307hxY9LT01UCYUxbBCJSovbt28dtt93GvHnz6NOnD6+++ipVq1b1O5achIpAREpMeno6bdu2BeCxxx4jKSlJHxIrBbRrSERKxLRp04IlADB+/HiVQCnhaRGYWVcz+9bMNpnZ2CKWiTezVWa21sxSvcwjIiUnJSWF2NhYoqKiqFKlCsOHDwegQYMG2hVUynhWBGYWDUwFugHNgP5m1uyEZWoA04AezrnmQB+v8ohIyUlJSSExMZHMzEyccxw6dAiAp556ij599Gdc2ni5RdAW2OSc2+ycywbeBnqesMyfgPecc1sBnHM7PMwjIiUkKSmJgwcP/ub2qVOnsmXLltAHkrPiZRHUA7YVuL49/7aCLgVqmtkCM1thZgMLuyMzSzSzDDPLWLlyJWZW6AWgRo0a1KhR49i/K/RScHzmzJl07NixyGWPjUdF5b1UUVFRhS5XcLxjx47MnDmzyPssOH4ss56TnlNpek6ZmZmF/amydetW3n33XU0dUcqYc86bOzbrA3Rxzg3Ovz4AaOucu6/AMi8AbYDrgcrAMuAPzrkNRd1vmzZtXEZGhieZRaR4YmNjCy0DM+P111/n9ttvD5aJhAczW+Gca1PYmJdbBNuBiwtcrw9kFbLMXOfcr865XcBCoJWHmUSkBEyaNIkqVaocd1tUVBTPPfccAwYMUAmUMl4WQTpwiZk1NLMKQD/ggxOWeR/oaGblzKwK8HvgGw8ziUgJSEhIYNCgQcHdUjVr1uS1117j/vvv9zmZnAnPPlDmnMsxsxHAPCAamOGcW2tmQ/PHX3LOfWNmc4GvgFzgFefc115lEpGzt2fPHoYPH87bb79N27Zt+cc//qFjAqWcZ8cIvKJjBCL+mTt3LoMGDWLnzp08+uijjB07lnLlNEFBaeDXMQIRKSMOHDjAsGHD6NatG7Vq1SItLY3x48erBMoIFYGInNSSJUu44oorePnllxk1ahQZGRlcddVVfseSEqQiEJFC/fe//2Xs2LHExcWRm5vLggULmDx5MpUqVfI7mpSwU27XmVl74HagI3AhcAj4GpgNvOmc2+tpQhEJudWrVzNgwADWrFnDkCFDeOaZZzjnnHP8jiUeOekWgZl9BAwm78yfruQVQTNgPFAJeN/MengdUkRCIycnhyeeeIKrr76aHTt28OGHH5KcnKwSKONOtUUwIP+DXgUdAFbmX54xszqeJBORkNq0aRMDBw5k2bJl9O7dmxdffJE6dfTnHQlOukVQSAmc0TIiEr6cc7z44ou0atWKb775hpSUFN555x2VQAQp1rlfZrYfOPaBgwpAeeBX51x1r4KJiPd++OEH7r77bubNm8eNN97IjBkzqFfvxLkhpawrVhE4547bQWhmt5A3zbSIlELOOd566y2GDx9OdnY206ZNY+jQoZojKEKd0emjzrl/A9eVcBYRCYFdu3bRt29fEhISaNq0KatWrWLYsGEqgQhW3F1DfyxwNYq8qaNL19wUIsLs2bMZPHgwu3fv5oknnmD06NFER0f7HUt8VtzPh3cv8HMOsIXfftuYiISp/fv38+c//5np06dz+eWXM3fuXFq10ozvkqe4xwju8jqIiHhj4cKF3HnnnWzZsoWHH36YCRMmULFiRb9jSRg51QfKxptZrZOMX2dmN5d8LBE5W4cPH2b06NHEx8djZixcuJAnn3xSJSC/caotgjXALDM7TN4HyHaS94niS4ArgE+Bxz1NKCKnbeXKlQwcOJC1a9cydOhQJk+eTLVq1fyOJWHqVGcN9XbOdSBviom15H3BzD7gTfK+f/hB59xOjzOKSDHl5OTw17/+ld///vfs2bOHOXPm8OKLL6oE5KROtUXQ2sxigAQgcMJYZfImoBORMPDtt99yxx13kJaWRr9+/Zg6dSq1ahW5Z1ck6FRF8BIwF2gEFPxaMCPv9NFGHuUSkWLKzc1l2rRpjBkzhkqVKvHWW2/Rr18/v2NJKXLSInDOTQGmmNmLzrlhIcokIsW0bds27rrrLj777DO6devGK6+8wkUXXeR3LCllivXJYpWASHhxzvHGG29w+eWXs3z5cl5++WVmz56tEpAzom8oEylldu7cSe/evRk4cCAtWrRg9erVJCYmaooIOWMqApFS5IMPPqBFixZ8+OGHPPXUU6SmptK4cWO/Y0kpV9wpJkTER/v27WPkyJG8+uqrtGrVik8//ZTLL7/c71hSRmiLQCTMLViwgJYtW/L6668zbtw4vvjiC5WAlCgVgUiYOnToEA8++CCBQIDy5cuzePFiJk2aRIUKFfyOJmWMdg2JhKGMjAwGDBjA+vXruffee3n66aepWrWq37GkjFIRiISRXbt2cd555wFQr1694FdIinhJu4ZEwsT06dODJQCwZs0alYCEhLYIRHyWnZ1N7dq1OXDgQPC2J598kpo1a/qYSiKJtghEfPT5559TsWLFYAmsXr3a50QSiVQEIj5wzhEIBLj++usB6NKlC7m5uVxyySU+J5NIpF1DIiG2YcMGmjRpEryemppKXFycj4kk0nm6RWBmXc3sWzPbZGZjT7Lc1WZ21Mx6e5lHxA8pKSnExsYSFRVF9erVgyVQu3ZtsrOzVQLiO8+KwMyigalAN6AZ0N/MmhWx3FPkfQuaSJmSkpJCYmIimZmZOOfYv38/AEOGDGHXrl2UL1/e54Qi3m4RtAU2Oec2O+eygbeBnoUsdx/wLrDDwywivkhKSuLgwYO/uf3jjz/+zW3OOT766KNQxBI5jpdFUA/YVuD69vzbgsysHnAred+EViQzSzSzDDPLWLlyJWZW6AWgRo0a1KhR49i/K/RScHzmzJl07NixyGWPjUdF5b1UUVFRhS5XcLxjx47MnDmzyPssOH4ss55T2XxOmZmZhf5OZ2ZmkpWVxaxZs47L1qtXLwC2bNlysj8JkZLlnPPkAvQBXilwfQDw/AnL/BNol//za0DvU91v69atnUhpERMT48j7WtfjLjExMc455xYuXOji4uIc4Bo0aOBeeeUVl52d7W9oKZOADFfE+6qXWwTbgYsLXK8PZJ2wTBvgbTPbAvQGppnZLR5mEgmpSZMmUaVKleNuq1KlCoMGDaJr167ExcWxYcMGnn/+eTZs2MDdd9+t4wYScl6ePpoOXGJmDYEfgH7Anwou4JxreOxnM3sN+NA5928PM4mEVEJCApB3rGDr1q3UrVuXiy66iEcffZTatWvz9NNPM3z48N+UhUgoebZF4JzLAUaQdzbQN8A7zrm1ZjbUzIZ69bgi4SYhIYGPP/6Yfv368dNPP7Fx40YmTpzI5s2bGT16tEpAfOfpB8qcc3OAOSfcVuiBYefcnV5mEfHDli1beOyxx3j99depWLEiY8eOZdSoUdSqVcvvaCJB+mSxiAeysrKYNGkS06dPJyoqivvuu4+xY8dywQUX+B1N5Dc015DIaVq2bBkTJ04kNzf3N2M7d+5k1KhRNG7cmOTkZAYNGsTGjRt59tlnVQIStrRFIFJMv/76K+PHj+dvf/sbzjkSEhJo3LgxAL/88gvPPPMMzz33HAcPHuT222/nkUceCY6LhDMVgUgxzJ8/n8GDB7N582auuOIKVq1aBcCBAweYMmUKkydP5pdffqFPnz5MmDCBpk2b+pxYpPi0a0jkJPbt28ewYcO47rrriIqKIjU1lYceegiAZ599lkaNGpGUlMS1117Ll19+yTvvvKMSkFJHRSBShLlz59KiRQuSk5P585//zOrVq2nXrh0vvZR34tvUqVNp1aoVy5YtY9asWVxxxRU+JxY5MyoCkRP8/PPP3HXXXXTr1o1q1aqxZMkSnnzySf75z39y2WWXsXTpUgBeeeUVPvnkE9q1a+dzYpGzoyIQKeD999+nWbNmvPHGGyQlJbFixQoyMzNp0aIFd955JzVr1qR///4AxMfH+xtWpISoCETIO+2zf//+3HLLLVxwwQV88cUXtG3blnbt2tGvXz+io6N59913ycjIoFu3bn7HFSlROmtIIppzjnfeeYcRI0awd+9eJk6cSJs2bbj33ntJS0ujcePGvPnmm8EyECmLtEUgEevHH3+kV69e9OvXj4YNG/L888/z2WefcdNNN5GVlcX06dP55ptvSEhIUAlImaYtAok4zjneeOMNRo4cyaFDh+jbty+7d+9m6NCh1K1bl+eff54hQ4ZQsWJFv6OKhISKQCLKtm3buOeee/joo4+oUaMGjRs35v/+7/+oVauWpoSWiKUikIjgnGP69OmMGjUq+AXyv/zyC7m5uUyYMIGRI0dSvXp1n1OK+ENFIGXe5s2bGTJkCJ9//nnwtipVqnD//fczevRoTQktEU9FIGVWbm4uL7zwAg888EDwtgoVKjBs2DDGjh1L3bp1fUwnEj5UBFImffvtt/Ts2ZNvv/02eFtiYiLjx4/n4osvPsm/FIk8On1UypScnBzGjRvHZZddFiyBAQMGsGnTJl5++WWVgEghtEUgZcby5ctp37598HrHjh156aWXaNasmY+pRMKfikBKvf379xMTE8PPP/8cvG3lypVceeWVPqYSKT20a0hKrezsbBITE6levXqwBGbNmoVzTiUgchpUBFLqHD16lOTkZCpWrMj06dMBGDlyJM45br75Zp/TiZQ+2jUkpUZubi7/+te/6Nu3b/C2unXrsm7dOmrWrOljMpHSTVsEEvacc8yaNYvLLrvsuBKYO3cuP/74o0pA5CypCCRsOef49NNPad++PT169GDjxo0A3HPPPezbt48uXbr4nFCkbNCuIQlLS5YsYfz48SxYsCB4W0xMDK+//jqdOnXyL5hIGaQtAgkrK1as4KabbuLaa689rgQeeugh1q1bpxIQ8YC2CCQsfP311zz66KO89957REX9b/3ksssuY8aMGcd9UExESpa2CMRXGzduJCEhgZYtW/LJJ59w1VVXUblyZaKjoxk3bhxffvmlSkDEY9oiEF9s3bqViRMn8tprr1GhQgUGDRrEtm3b+Pjjj2nZsiWvvvoqV111ld8xRSKCikBC6scff+Txxx8nOTkZgOHDh9O4cWP++te/8ssvvzBhwgTGjh1LhQoVfE4qEjlUBBIS69ato3nz5sHriYmJDB48mCeeeIIpU6bQpk0bPvvsMy6//HIfU4pEJk+PEZhZVzP71sw2mdnYQsYTzOyr/MtSM2vlZR4JjZSUFGJjY4mKiqJ+/fqY2XEl8Oqrr9KhQwe6dOnCnDlzeOqpp1i2bJlKQMQnnhWBmUUDU4FuQDOgv5mdOB/w90An51xL4DEg2as8EhopKSkkJiaSmZmJc44ffvghODZixAgAnn76ae644w6aNWvG6tWrGTNmDOXKaeNUxC9ebhG0BTY55zY757KBt4GeBRdwzi11zh2bO3g5UN/DPBICSUlJHDx48De3x8TEcMsttwCQmZnJlClTWLhwIU2aNAl1xDOSm5vLl19+ybPPPsu0adP8jiNSorwsgnrAtgLXt+ffVpS7gY8KGzCzRDPLMLOMFStWYGaFXgAqVapEpUqVjv27Qi8Fx6dOnUqzZs2KXLbg+Mnus+B4s2bNmDp1apHLFRw/lrmsPKfMzMxC/3MzMzO54YYbqFChAmvWrGH37t1MnDgRgIsuuqjQ+/zLX/4SHM/KymLWrFlFPv6x8e7duwPQvXv3QpcrOD5r1iyysrKKvM8pU6YwceJEKlWqRJ06dbjqqqt46KGHWL58OW3btiU2NvYkv84ipYc557y5Y7M+QBfn3OD86wOAts65+wpZNgBMA651zu0+2f22adPGZWRkeBFZSkBsbGyhZRATE8OWLVtCH6iYnHOsXbuWBQsWMH/+fFJTU9m9O+9XsVGjRsTHxxMIBIiPj6d+fW24SuljZiucc20KG/Nyx+x2oOAXxNYHsk5cyMxaAq8A3U5VAhL+Jk2aRGJi4nG7h6pUqcKkSZN8TPVbzjnWr1/P/Pnzg2/8O3fuBPJKq3v37sE3/gYNGvicVsRbXhZBOnCJmTUEfgD6AX8quICZNQDeAwY45zZ4mEVCJCEhAcg7VrB161YaNGjApEmTgrf7xTnHhg0bmD9/PgsWLGDBggX85z//AaB+/fp07dqVQCBAIBDQLh+JOJ7tGgIws5uA54BoYIZzbpKZDQVwzr1kZq8AvYBj+xJyitp0OUa7hqQ4nHN89913wTX+BQsW8OOPPwJ5xxyOre0HAgEaNWoUPB4iUladbNeQp0XgBRWBFMY5x/fffx98058/f37w1NW6desG3/QDgQC/+93v9MYvEcevYwQinsrMzDxujX/r1q0AnH/++cTHxwff/Js0aaI3fpGTUBFIqbFt27bg2v78+fODZyHVrl2b+Ph4xowZQyAQoGnTpnrjFzkNKgIJW1lZWcet8X/33XcA1KpVi06dOvHggw8SCARo3rz5cd9hICKnR0UgYeOnn346bo3/2HcU16hRg7i4OEaMGEEgEODyyy/XG79ICVIRiG927NgRfONfsGAB69evB6B69erExcUxdOhQ4uPjadWqFdHR0T6nFSm7VAQSMrt27SI1NTW4xr9u3ToAqlWrRseOHRk0aBCBQIArrrhCk9CJhJD+2sQze/bsITU1NbjWv2bNGgCqVq3Ktddey8CBA4mPj6d169Z64xfxkf76pMT8/PPPLFq0KLjG/9VXX+Gco3LlynTo0IF+/foRCARo06YN5cuX9zuuiORTEcgZ27t3L4sWLQqu8X/55Zc456hUqRLXXHMNEyZMIBAI0LZtW331pEgYUxFIse3fv5/FixcH1/hXrlxJbm4uFSpUoH379jz66KPBN/5j02KLSPhTEUiRDhw4wJIlS4Jr/BkZGRw9epTy5cvTrl07kpKSCAQCtGvXjsqVK/sdV0TOkIpAgg4ePMjSpUuDa/zp6enk5ORQrlw52rZty9ixYwkEArRv354qVar4HVdESoiKIIIdOnSIZcuWBc/jT0tL48iRI0RHR3P11VczevRo4uPj6dChA1WrVvU7roh4REUQQQ4fPkxaWlpwjX/58uVkZ2cTFRVF69atg1M2dOjQgXPOOcfvuCISIiqCMuy///0vX3zxRXCNf9myZRw+fPK+8/AAAAlFSURBVJioqCiuvPJK7r//fuLj4+nYsSPVq1f3O66I+ERFUIZkZ2eTkZERXONfunQphw4dwsxo1aoVw4YNIxAI0LFjR2rUqOF3XBEJEyqCUuzIkSOsWLEieFbP4sWLg98V3LJlSxITE4mPjycuLo5atWr5nFZEwpWKoBTJycnhyy+/DK7xL168mAMHDgDQvHnz4Fw9cXFx1KlTx+e0IlJaqAjC2NGjR1m1alVwjX/RokXs27cPgKZNmzJgwAACgQCdOnXi/PPP9zmtiJRWKoIwkpuby1dffRVc41+4cCF79+4FoEmTJvTv3z/4xl+3bl2f04pIWaEi8FFubi5ff/11cI0/NTWVn3/+GYDf/e533HbbbcHv3r3ooot8TisiZZWKIIScc6xbty54OueCBQvYvXs3AI0aNeLWW28lEAgQHx9P/fr1fU4rIpFCReAh5xzr168/7lu4du7cCUBMTAzdu3cPrvHHxMT4nFZEIpWKoAQ559i4ceNxa/w//fQTAPXr16dr167BNf6GDRv6nFZEJI+K4Cw45/juu++Oe+PPysoC4KKLLuL6668nPj6eQCBAo0aNMDOfE4uI/JaK4DQ459iyZUvwrJ4FCxawfft2AC644AICgUBwjf+SSy7RG7+IlAoqglPIzMwM7uOfP38+W7duBeC8884LvukHAgGaNGmiN34RKZVUBCfYvn37cWv833//PQC1a9cmPj6eMWPGEB8fT7NmzfTGLyJlQsQXQVZW1nFr/N999x0AtWrVolOnTowcOZJAIEDz5s2JioryOa2ISMmLuCL46aefjjudc8OGDQCce+65dOrUiREjRhAfH0/Lli31xi8iEaHMF8GOHTtITU0NrvGvX78egOrVqxMXF0diYiKBQIBWrVoRHR3tc1oRkdArc0Wwa9cuUlNTg2v9a9euBaBatWp07NiRQYMGER8fz5VXXkm5cmXu6YuInDZzznl352Zdgb8B0cArzrknTxi3/PGbgIPAnc65lae4TxcTE8OkSZNISEhgz549LFy4MLjGv2bNGgCqVq3KtddeGzyrp3Xr1nrjF5GIZWYrnHNtCh3zqgjMLBrYAHQGtgPpQH/n3LoCy9wE3EdeEfwe+Jtz7venuF8HUK5cOerVq8fWrVtxzlG5cmU6dOgQPKXz6quvpnz58p48NxGR0uZkReDlKnJbYJNzbnN+iLeBnsC6Asv0BP7h8tpouZnVMLMLnXM/nurOc3Jy+M9//sOECRMIBAJcffXVVKxY0YvnISJSpnlZBPWAbQWubydvrf9Uy9QDjisCM0sEEk98gMOHD/PII4+sKJG0JasOsMvvECcR7vkg/DOGez4I/4zhng/KVsYiZ7b0sggK+7TVifuhirMMzrlkIBnAzDKK2rwJF+GeMdzzQfhnDPd8EP4Zwz0fRE5GL0+U3w5cXOB6fSDrDJYREREPeVkE6cAlZtbQzCoA/YAPTljmA2Cg5WkH7C3O8QERESk5nu0acs7lmNkIYB55p4/OcM6tNbOh+eMvAXPIO2NoE3mnj95VjLtO9ihySQr3jOGeD8I/Y7jng/DPGO75IEIyevo5AhERCX+aTEdEJMKpCEREIlzYFoGZdTWzb81sk5mNLWTczGxK/vhXZnZVmOW7zMyWmdl/zWxUKLOdRsaE/NfuKzNbamatwixfz/xsq8wsw8yuDWW+4mQssNzVZnbUzHqHMl/+Y5/qdYw3s735r+MqM3sknPIVyLjKzNaaWWoo8xUno5mNLvD6fZ3/f10rjPKda2azzGx1/mtYnOOt/+OcC7sLeQeXvwMaARWA1UCzE5a5CfiIvM8itAPSwizf+cDVwCRgVJi+htcANfN/7haGr2E1/nccqyWwPtxewwLLfU7eyQ+9wy0jEA98GOrfwdPIV4O8GQca5F8/P9wynrB8d+DzcMoHjAOeyv/5PGAPUKG4jxGuWwTB6Smcc9nAsekpCgpOT+GcWw7UMLMLwyWfc26Hcy4dOBKiTCcqTsalzrmf868uJ+9zHOGU74DL/80GqlLIhw39zpjvPuBdYEcow+Urbka/FCffn4D3nHNbIe9vJwwzFtQfeCskyfIUJ58DzsmfyLMaeUWQU9wHCNciKGrqidNdxit+PnZxnW7Gu8nbwgqVYuUzs1vNbD0wGxgUomzHnDKjmdUDbgVeCmGugor7/9w+f7fBR2bWPDTRgOLluxSoaWYLzGyFmQ0MWbo8xf5bMbMqQFfyij9UipPvBaApeR/IXQM84JzLLe4DhOu8zCU2PYVH/Hzs4ip2RjMLkFcEodwHX9zpRf4f8P/MLA54DLjB62AFFCfjc8DDzrmj5s93WBcn40ogxjl3IH/G338Dl3ieLE9x8pUDWgPXA5WBZWa23Dm3wetw+U7n77k7sMQ5t8fDPCcqTr4uwCrgOqAx8ImZLXLO7SvOA4TrFkG4T09RGqbGKFZGM2sJvAL0dM7tDlE2OM3X0Dm3EGhsZnW8DlZAcTK2Ad42sy1Ab2Camd0SmnhAMTI65/Y55w7k/zwHKB/C17G4f8tznXO/Oud2AQuBUJ64cDq/i/0I7W4hKF6+u8jbveacc5uA74HLiv0IoTwocxoHR8oBm4GG/O/gSPMTlvkDxx8s/iKc8hVY9i/4c7C4OK9hA/I+1X1NmOb7Hf87WHwV8MOx6+GS8YTlXyP0B4uL8zrWLfA6tgW2hup1LGa+psBn+ctWAb4GWoTTa5i/3Lnk7XuvGob/xy8Cf8n/+YL8v5U6xX2MsNw15LybniJk+cysLpABVAdyzWwkeUf6i7WpFoqMwCNAbfLWYgFyXIhmWixmvl7kzUV1BDgE9HX5v+lhlNFXxczYGxhmZjnkvY79QvU6Fiefc+4bM5sLfAXkkvdthl+HIl9xM+YveivwsXPu11BlO418jwGvmdka8laOH3Z5W1fFoikmREQiXLgeIxARkRBREYiIRDgVgYhIhFMRiIhEOBWBiEiEUxGIiEQ4FYGISIRTEYicJTOLNbOvC1wfZWZ/8TGSyGlREYiIRDgVgYhIhFMRiJy9HI7/W6rkVxCRM6EiEDl7/wHON7PaZlYRuNnvQCKnIyxnHxUpTZxzR8xsIpBG3jzw632OJHJaNPuoiEiE064hEZEIpyIQEYlwKgIRkQinIhARiXAqAhGRCKciEBGJcCoCEZEI9/8Bic09XU+WG10AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "K_max = 100 # max number of iterations\n",
    "u = 0.0 #  initial value of displacement\n",
    "f_t = np.linspace(0, 0.99999, 6) # prescribed load levels\n",
    "u_t = []  # calculated displacement values for each load level\n",
    "f_t_iter = [0.0] \n",
    "u_t_iter = [0.0] # interim values of displacement\n",
    "f_t_levels = [] # used for visualization of load levels\n",
    "u_t_levels = [] # used for visualization of load levels\n",
    "for f_level in f_t: # loop over the load levels\n",
    "    for K in range(K_max): # iteration loop\n",
    "        R = f(u) - f_level # value of the residuum\n",
    "        if np.fabs(R) < 1e-8: # residuum equal to zero?\n",
    "            break # stop iteration\n",
    "        dR = df_du(u) # derivative of the residuum\n",
    "        d_u = - R / dR # increment of displacement\n",
    "        u += d_u # update the total displacement\n",
    "        u_t_iter.append(u)# record data for visualization of iterations\n",
    "        f_t_iter.append(f_level)\n",
    "        u_t_iter.append(u) \n",
    "        f_t_iter.append(f(u))\n",
    "        f_t_levels.append([f_level,f_level]) \n",
    "        u_t_levels.append([0,u]) \n",
    "        if K == K_max - 1:\n",
    "            raise ValueError('No convergence')\n",
    "    print('number of iterations =',K)\n",
    "    u_t.append(u) # resord the found solution for current load level f\n",
    "plt.plot(u_t, f_t, 'o', color='black', lw=2) # plot solved increments as bullets\n",
    "plt.plot(u_t_iter, f_t_iter, color='black') # plot iterations as solid liness\n",
    "plt.plot(np.array(u_t_levels).T, np.array(f_t_levels).T, '-.',  \n",
    "                color='black', lw=1) # plot load levels as dash-dotted lines\n",
    "plt.xlabel('u')\n",
    "plt.ylabel('f(u)')\n",
    "plt.xlim(0)\n",
    "plt.ylim(0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6 Solve the problem using the newton method provided in scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import newton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pdoc newton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}