oving02_81.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Possible solution to exercise 8.1 using Python.\n",
    "\n",
    "One way is to use Python's symbolic algebra package SymPy.\n",
    "We begin by importing SymPy and defining the symbols - I will use $x$,$y$,$z$, and $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "x,y,z,alpha = sp.symbols('x y z alpha')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define f symbolically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 x^{2} + x y - 6 x + y^{2} + y z - 7 y + z^{2} - 8 z + 9$"
      ],
      "text/plain": [
       "2*x**2 + x*y - 6*x + y**2 + y*z - 7*y + z**2 - 8*z + 9"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = 2*x**2 + x*y + y**2 + y*z + z**2 - 6*x - 7*y - 8*z + 9\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us calculate the gradient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}4 x + y - 6\\\\x + 2 y + z - 7\\\\y + 2 z - 8\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[    4*x + y - 6],\n",
       "[x + 2*y + z - 7],\n",
       "[    y + 2*z - 8]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = sp.Matrix([f.diff(x),f.diff(y),f.diff(z)])\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Search direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = sp.Matrix([1,2,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now verify that $\\nabla f(0,0,0)^T p < 0$, that is that $p$ is a descent direction at $0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle -20$"
      ],
      "text/plain": [
       "-20"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df0 = df.subs([(x,0),(y,0),(z,0)])\n",
    "df0.dot(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now solve for $\\alpha$ satisfying the sufficient decrease condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0 \\leq \\alpha \\wedge \\alpha \\leq 0.5$"
      ],
      "text/plain": [
       "(0 <= alpha) & (alpha <= 0.5)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c1 = 4/5\n",
    "f_alpha = f.subs([(x,0+alpha*p[0]),(y,0+alpha*p[1]),(z,0+alpha*p[2])])\n",
    "f0      = f_alpha.subs(alpha,0)\n",
    "sp.solve(f_alpha <= f0 + c1*alpha*df0.dot(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore the sufficient decrease condition is satisfied, if $\\alpha \\in [0,0.5]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course not everything can be done using symbolic computations.  One could e.g. plot a picture and see approximately where the lines intersect: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc40892ad30>"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0mElEQVR4nO3dd3QUVRvH8e8NCQRS6AYl9B5ICCQ0qRGQKkVAQECqoYmogKDAS1dBFBAJvUkLHWnSJDSpKQRSpIpSBekBacl9/5iogAFD2uxmn885c8hmZ2d+1/XMkzvlXqW1RgghhO2xMzuAEEIIc0gBEEIIGyUFQAghbJQUACGEsFFSAIQQwkbZmx3gReTKlUsXLFgwSZ+9c+cOTk5OKRvIwkmbbYO02TYkp80hISF/aK1zP/17qyoABQsWJDg4OEmf3bFjB7Vq1UrZQBZO2mwbpM22ITltVkr9mtDv5RSQEELYKCkAQghho6QACCGEjbKqawBCCJGWHj58yLlz57h3757ZUciaNSvR0dHPXcfR0RF3d3ccHBwStU0pAEII8Qznzp3DxcWFggULopQyNcvt27dxcXF55vtaa65evcq5c+coVKhQorYpp4CEEOIZ7t27R86cOU0/+CeGUoqcOXO+UG9FCoAQQjyHNRz8//KiWW2iAIStnML16YOJjX1kdhQhhLAYNlEAbs2dRvPAvWyrmofoX0PMjiOEEBbBJgpAjbXhbGldk7oHr3KvSgUmLevHg9gHZscSQghT2UQBUHZ2ZOwxnFvLF1Lihj1tO39NjwGlCL6QtGElhBAiLd25c4d+/fqxf//+FN2uTRSAv2Rr0Y4sYUfJ7JaX6ZNOM//diny8ZQB3H941O5oQQjzTtGnTePDgAXv27EnR7dpUAQCgRAlcwiKhQX0mb9SUGjieCpM92Xlmp9nJhBAiQZs2baJo0aJ4e3un6HZt80GwrFlxWLsBRoyg88iR+Fw9T6NztWj8Wg/G1h2LayZXsxMKISzMB5s+4PClwym6Te883kysP/G569y7d4/Y2FgOHz7MwIEDU3T/ttcD+IudHYwYAWvW4HkjI1FzMnNs5XRKB5Rmw/ENZqcTQggATpw4QWxsLMWLF0/0EA+JZZs9gMc1bYo6eBCX5s35ccEJxjV7SOObjXnb620m1Z9Eriy5zE4ohLAA//WXemq5cuUKx48fZ+HChSm+bdvtATyuZEk4cADVpAkDV/7Okd1l2BC2jFJTShEYEYjW2uyEQggbdeHCBVq0aEFcXBzXr19P0W2bWgCUUh8qpSKVUhFKqSVKKUfTwri6wsqV8PnneO6I4tLKQtR4+AptV7alaWBTzt86b1o0IYRtevToEaGhoVy6dIkPP/yQDBkypOj2TSsASqm8wPuAr9a6DJABaGNWnvhQMGgQ/PADjr9fZcW4X1mdpQvbTm/DI8CDGSEziNNxpkYUQtgOe3t7vv76a1asWMH8+fNxdU3ZG1TMPgVkD2RWStkDWYALJucxvP46hISgihal2cdzOHe9CxXcytN9fXdqf1ebk9dOmp1QCCGSzbQCoLU+D4wHfgMuAje11lvMyvMvBQvCnj3QpQs5vprC1iUOzK/2FaEXQ/Ga6sVXe78iNi7W7JRCCJFkyqwLnEqp7MBKoDVwA1gOrNBaL3xqPX/AH8DNzc0nMDAwSfuLiYnB2dk5SZ99ecMGik2axIPs2dkz+AOG2a1n79W9lHQpycclPqaQU+ImX0hryWmztZI224a0anPWrFkpWrRoqu8nMWJjYxN1DeDkyZPcvHnzid/5+fmFaK19/7Wy1tqUBWgFzH7s9TtAwPM+4+Pjo5MqKCgoyZ/VWmsdHKx1wYJaOzjouG+/1UuOLNa5xuXSDiMd9PCg4fr+o/vJ234qSHabrZC02TakVZujoqLSZD+JcevWrUStl1BmIFgncEw18xrAb0BlpVQWZcxiUBt4/oSXZvLxgZAQqFMH9d57tBm7gehOwbQq3YrhO4fjM8OHg+cPmp1SCCESzcxrAAeAFUAocDQ+ywyz8iRKjhywfj2MGgWLF5OrVkMWeQxlXdt1XP/zOlVmV6H/lv4yuJwQwiqYeheQ1nqY1rqk1rqM1rqD1vq+mXkSxc4OhgyBLVvgyhWoUIHGoTFE9Y7Cv7w/X+37Cs+pngT9EmR2UiGEeC6zbwO1XnXqQGgoeHlB27a49vuUqXUnsqPjDuyUHa999xr+6/y5ee/mf29LCCFMIAUgOdzdYccO6NcPpkyBatWoqfMT3iOcAa8OYHbYbDwCPFh7bK3ZSYUQ6cCqVat4//33U2x7UgCSy8EBxo+H1avhxAkoX54sG7cyru44DnQ7QM7MOWka2JS2K9ty+c5ls9MKIaxYWFgY5cuXT7HtSQFIKc2aGaeEihQxfv7oI3xzeRHsH8wov1Gsil6FxxQPFh1ZJIPLCSFeyPHjx2natCkTJ05k+PDhTJw4MUW2K8NBp6TCheGnn6B/f5gwAfbuJePSpQypMYQ3S71J17Vdab+6PUsiljC10VTyZc1ndmIhRGJ98AEcPpyy2/T2hv84mN+/f5+33nqLadOm8fbbb7N37148PDzo0aMHjo7JGz9TegApLVMmmDwZli+H6GjjC/7+ezxye7Cn8x4m1ptI0JkgSgeUZlrwNBlcTgjxXFu3bqVs2bLkyZMHV1dX8uTJg6OjI7GxyR+KRnoAqaVlSyhXDlq3Nk4J9e1LhrFj6Vu5L01KNMF/vT89N/RkScQSZr0xi2I5i5mdWAjxPCl02uVFHT58GE9PTyIiIvDy8uLy5cu4uLjg5OSU7G1LDyA1FSlinBJ6/32YNAmqVoVTpyiUvRBb2m9hdpPZhF8Kx2uaF+N+GsejuEdmJxZCWBhXV1d+/vnnvwvAsGHD6N27d4psWwpAasuUyTj4r14Np04ZvYKlS1FK0aVcF6J6R1G/aH0GbhtI5VmVOfL7EbMTCyEsSPv27Tlx4gRjx45l6tSp5MiRgz59+qTItqUApJVmzYwLSGXKQJs24O8Pd+/yissrrHprFctaLuPsrbP4zPDhf0H/4/4jy38oWgiR+nLkyMHOnTtxd3fn4MGDjBkzBmP4tOSTApCWChSAnTvhk09g1iyoUAGOHkUpRavSrYjqFUXbMm0ZtWsU5WeUZ/+5/WYnFkJYgPv37xMTE0POnDlTdLtSANKagwN89hls3gxXr0LFijBtGmhNziw5+a75d2x8eyO379/m1dmv8uGmD7nz4I7ZqYUQJsqUKRNHjqT86WEpAGapWxfCw6FmTejZ07hr6No1ABoUa0Bkr0h6+vZk4oGJeE71ZNvpbSYHFkKkN1IAzOTmBhs3wpdfwrp1ULYs7NoFgEsmF6Y0msKuTruwt7On7oK6dP2+Kzfu3TA3sxAi3ZACYDY7O+PJ4b17wdER/Pzgf/+DR8YtodULVCe8RzgDqw5kfvh8PKZ4sObnNeZmFsKGWNPQLS+aVQqApfD1NcYS6tDBmHCmRg345RcAMjtk5os6X3Dw3YO4ObvRfGlz3lr+Fr/H/G5yaCHSN0dHR65evWoVRUBrzdWrV19oeAhTnwRWSmUDZgFlAA100VrvMzOTqVxcYN48qF8fevQwTgkFBED79gCUf7k8B7sd5Mu9XzJi5wi2nd7GxPoT6eDVIcVuCxNC/MPd3Z1z585x5coVs6Nw7969/zy4Ozo64u7unuhtmj0UxCRgk9a6pVIqI5DF5DyWoU0bqFIF2rUzegQ//GDMN5AtGw4ZHPi0+qd/Dy7XcU1HlkQsYXrj6eTPmt/s5EKkKw4ODhQqVMjsGADs2LGDcuXKpeg2TTsFpJRyBWoAswG01g+01jfMymNxChQwJpsZNQqWLjV6A7t3//12yVwl2d15N5MbTGb3r7spHVCaKQenyOByQohEU2ad21JKeWNMAh8FlAVCgL5a6ztPrecP+AO4ubn5BAYGJml/MTExODs7JyeyaVyiovAYMwbHS5f4rW1bznTsiHZw+Pv9S/cu8dXxrwi+Hoynqyf9S/Qnf5b8Vt3mpJI22wZp84vx8/ML0Vr7/usNrbUpC+ALPAIqxb+eBIx63md8fHx0UgUFBSX5sxbh1i2tu3bVGrT28dE6OvqJt+Pi4vS8sHk6+xfZdaZRmfTnuz/XW3/calJY81j995wE0mbbkJw2A8E6gWOqmXcBnQPOaa0PxL9eAaTcXGfpjYuLMXzEypXG3UHlyxsXiON7cEopOnp3JKp3FI2LN+aTHz+hV1gvDl86bG5uIYTFMq0AaK0vAWeVUiXif1Ub43SQeJ4334SjR43bRHv3hoYN4eLFv9/O45yHFW+tYEWrFfxx/w98Z/gy+MfB3Ht0z8TQQghLZPZzAH2ARUqpI4A38Jm5cazEK68YdwZ9+61xodjTE1aseGKVFh4tmFdhHh3KduCzPZ9Rbno59p7da05eIYRFMrUAaK0Pa619tdZeWutmWuvrZuaxKkoZPYCwMChYEFq1Mm4ZvXHj71VcHVyZ23Qum9tv5s+Hf1JtTjXe/+F9Yh7EmBZbCGE5zO4BiOQqWRL27TOGj1iyxOgNbHty4LjXi7xORK8I3qv4Ht8e/JYyAWXYfHKzSYGFEJZCCkB64OAAI0YY4wk5ORkjjfbpg929f877O2d05psG37C7824c7R2pv6g+ndZ04tqf10wMLoQwkxSA9KRiRWM8ob594dtv8X33XaN38Jiq+atyuMdhPq32KQuPLMRjigcro1aaFFgIYSYpAOlNliwwcSJs347dw4dQrRoMGgT3/5li0tHekTG1xxDsH8wrLq/QcnlLWixrwcXbF5+9XSFEuiMFIL3y8+PQ7NnQpQuMHQs+PhAS8sQq3nm8OfjuQb6o/QUbjm/AI8CDuWFzrWLkQyFE8kkBSMdinZxg5kxj0pnr16FSJeNi8YMHf69jb2fPwGoDCe8RTpmXytBlbRfqLazHmRtnzAsuhEgTUgBsQYMGEBFhjC46apQxGX1Y2BOrlMhVgp2ddjKl4RT2ndtHmYAyfHPgG2LjYk0KLYRIbVIAbEX27DB/PqxdC5cvGxeMhw17ojdgp+zoVaEXkb0iqVGgBn039aX63OpEX4k2MbgQIrVIAbA1b7wBkZHGnAMjRxq9gdDQJ1bJnzU/G97ewILmCzh29Rje070ZvWs0D2MfmhRaCJEapADYohw5YMECozdw5YrRGxg8+Ik7hZRStPdqT3TvaJqVbMbQoKFUmFmB0Iuhz9mwEMKaSAGwZX/1Bjp0gM8+g3LlYP/+J1Z5yekllrZcyurWq7l85zIVZ1Zk0LZB/PnwT5NCCyFSihQAW5c9O8ydawwuFxMDr74KH30Ed56Yl4dmJZsR1TuKTt6dGPvTWLyne7P7193P2KgQwhpIARCG+vWNO4V69IAJE8DLC7Zvf2KVbI7ZmNVkFls7bOVB7ANqzKtB7w29uX3/tkmhhRDJIQVA/MPV1ZhkZscOsLOD2rWhW7cnRhgFqFO4DhE9I/ig0gdMDZ5K6YDS/HDiB1MiCyGSTgqA+LeaNeHIEfj4Y5g3D0qVglWrnljFKaMTE+pPYG/XvbhkcqHh4oa8s/odrt69ak5mIcQLkwIgEpY5szGExMGDkCcPtGhhzEZ24cITq1V2r0yofyhDawxlScQSSk0pxbLIZTKchBBWwPQCoJTKoJQKU0qtNzuLSED58kYR+OIL40JxqVIwbRrExf29Sib7TIz0G0mIfwgFshWg9YrWNF/anAu3Lzxnw0IIs5leAIC+gDxqaskcHGDgQGMuYl9f6NnTmJM4MvKJ1bzcvNjXdR/j6oxj86nNeEzxYHbobOkNCGGhTC0ASil3oBEwy8wcIpGKFjVmG5s7F6KjjecGhg6FxyaesbezZ0DVARzteRTvPN50W9eNugvqcvr6aRODCyESosz860wptQL4HHAB+mutGyewjj/gD+Dm5uYTGBiYpH3FxMTg7OycjLTWJzXb7HDjBkUCAsizdSt33d05/sEH3PDxeWKdOB3H+ovrmX56OnE6jq6FutI8b3MyqAypkgnke7YV0uYX4+fnF6K19v3XG1prUxagMRAQ/3MtYP1/fcbHx0cnVVBQUJI/a63SpM1bt2pdtKjWoHWHDlr//vu/Vjl786xutKiRZji60sxKOuL3iFSLI9+zbZA2vxggWCdwTDXzFFBVoIlS6gwQCLymlFpoYh6RFHXqGLeMDh4MgYHGJPUzZz5xkdjd1Z11bdex+M3FnLx2knLTyzFy50gexD54zoaFEKnNtAKgtf5Ea+2utS4ItAG2a63bm5VHJEPmzDB6NISHG08Q+/tD9epGYYinlKKtZ1uie0fT0qMlw3YMw2eGD4fOHzIxuBC2zRLuAhLpRalSEBRkPDx2/LhxC2m/fnD7n6EicjvlZnGLxaxts5brf16n8uzKDNgygLsP75qXWwgbZREFQGu9QydwAVhYIaWgY0c4dgy6doWvvzYKw/Ll8NgNB2+UeIPIXpF0K9eN8fvGU3ZaWXae2WlicCFsj0UUAJEO5cgB06fDvn2QOze89RbUq2cUhnhZHbMy/Y3pbH9nO3E6jlrza9FzfU9u3b9lXm4hbIgUAJG6KleG4GCYPNl4otjTEz799Inhpv0K+XG051H6VenHjNAZlA4ozYbjG0wMLYRtkAIgUl+GDPDee8Zf/23awOefG6eFVqz4+7RQFocsjH99PPu67iObYzYaL2lMu1XtuHLnisnhhUi/pACItOPmBt99B7t3GxPRtGoFr79uPFUcr2LeioT4hzC85nCWRy7HI8CDJUeXyHASQqQCKQAi7VWrBiEh8M03cOiQcevogAFwyzj3nzFDRobVGkZo91AKZy/M26vepmlgU87fOm9ycCHSFykAwhz29tCnj3G7aMeOMH48lChh9BDiHyIr81IZ9nbZy9evf82209vwCPBgRsgM4nTcf2xcCJEYUgCEuV56CWbNggMHoEABoxi8+qrRMwAy2GXgwyofcrTnUXxe9qH7+u7U/q42J6+dNDm4ENZPCoCwDBUrwt69xkNkZ84Yrzt3hosXASiSowg/vvMjM9+YSejFULymevHV3q+IjYs1NbYQ1kwKgLAcdnZGD+D4cWP+gcWLoXhxYzKae/dQStGtfDeiekVRt0hd+m/tT5XZVTj6+1GzkwthlaQACMvj6moc9CMj4bXX4JNPwMPDmJdYa/K65mVN6zUEtgjkzI0zlJ9RnmFBw7j/6L7ZyYWwKlIAhOUqWhS+/x62bgUnJ2NeYj8/CA1FKUXrMq2J6h1F69KtGblrJD4zfIi6FWV2aiGshhQAYfnq1IGwMAgIMHoFvr7G9YELF8iVJRcL31zI+rbruXn/Ju+FvcdHmz/izoM7/71dIWycFABhHeztjbmIT5wwRhhdtAiKFYMRI+DOHRoVb0Rkr0gav9yYCfsn4DXNi+2/bDc7tRAWTQqAsC7ZssGXXxpPDzdsCMOHGxeK587F1d6Jj4p/xI6OO7BTdtT+rjb+6/y5ee+m2amFsEhSAIR1KlLEGGJ6zx7Ilw+6dIHy5ckeHEzNgjUJ7xHOgFcHMDtsNh4BHqw9ttbsxEJYHCkAwrpVrWoMOR0YCLduUXbAAKhfnyxRJxhXdxwHuh0gZ+acNA1sSpsVbbh857LZiYWwGKYVAKVUPqVUkFIqWikVqZTqa1YWYeWUgtat4eefOdmzpzHsdLly0KkTvrFuBPsHM7LWSFb/vBqPKR4sOrJIBpcTAnN7AI+AflrrUkBloLdSysPEPMLaZcrEubfeglOnoH9/o1dQvDgZPx3KUK/3COseRrGcxWi/uj2NlzTm7M2zZicWwlRmTgp/UWsdGv/zbSAayGtWHpGOZM8O48YZ8w+0bm1cNC5cGI+569nTZisT6k1gx5kdlA4ozbTgaTK4nLBZyhK6wkqpgsAuoIzW+tZT7/kD/gBubm4+gYGBSdpHTEwMzs7OyUxqXaTNBqfTpyk8YwY5Dxzgfq5cnOnYkbBa3ow/NZGQGyF4ZfViQPEBuGdxNyl18sj3bBuS02Y/P78QrbXvv97QWpu6AM5ACPDmf63r4+OjkyooKCjJn7VW0uan7NypdZUqWoPWxYvruKVL9ZzgWTrr51m142hHPXbPWP0w9mGaZU0p8j3bhuS0GQjWCRxTTb0LSCnlAKwEFmmtV5mZRdiAGjXgp5+M4SUcHFCtW9PZP4BTxadQv0g9Bm4bSOVZlQm/FG52UiHShJl3ASlgNhCttf7arBzCxigFTZpAeLgx+cy1a+Rs0Z5VAdcIKjySs7fO4jvTl6Hbh8rgciLdM7MHUBXoALymlDocvzQ0MY+wJRkyQIcOxoXiKVNQJ09S653/cXZzGQY7vs7o3aPxnu7N3rN7zU4qRKox8y6gPVprpbX20lp7xy8bzcojbFTGjNCrF5w8CV9+ScawcIZ/vJFLuytT4MwNqs2pRt8f+hLzIMbspEKkOHkSWAiALFmMZwdOn4ZRo3ALjmbT+EsEby3MtvXf4DnVk22nt5mdUogUJQVAiMe5usKQIfDLLzBkCOWPXCZiqmLyd1fp83Vdun7flet/Xjc7pRApQgqAEAnJnh1GjYJffkF98gmNjsURGQB1hs6l6bDirI5ebXZCIZJNCoAQz5MzJ4wZgzpzBruBg3jrlCM7xv/Bn2+9yYcT6nEp5pLZCYVIMikAQiRGrlzw+edkOPMrun9/Wp7KyFf9trD/1fysXT5aBpcTVkkKgBAvInduMoz7koy/neda3+7UPRFHk7eGsrdCHi7uWG92OiFeyDMLgFJqQfy/MkyzEE/LlYtcE6aR+dwlDrzbAI/Iy7zs9wa/vupB3O5dZqcTIlGe1wPwUUoVALoopbIrpXI8vqRVQCEsmV3OXFSasZHbx44yt1UxModHY1ejJndfrQCbN4OcGhIW7HkFYBqwCSiJMVjb40tw6kcTwnrkz1+GTkuPsWXbdD5pnJlrkcFQvz7a1wdWroQ4GXJaWJ5nFgCt9TfamKxljta6sNa60GNL4TTMKIRVUErRvoo/fZecZsDU5nR7A86ejYKWLcHDA2bPhvsyvpCwHP95EVhr3TMtggiRXuRxzsOSt1fRcOxKqgzIRptWivOxN6BbNyhcGMaPh1u3/nM7QqQ2uQtIiFTyZqk3iXgvmiztOuHe7ne69MrLjQJ5YMAAyJ8fPvkELl40O6awYVIAhEhF2TNnZ07TOWzpsIXtxezJXi+UcRNa8rDOa8a0lQULQteuEBVldlRhg6QACJEG6hapS0SvCPpU7MOgmyspWi2E3VtmGaeFliyB0qWhUSPYvl3uHBJpRgqAEGnEOaMz3zT4ht2dd5PFIQs19nShU907XD8WDiNGwKFDULs2lC8PCxfCgwdmRxbpnNlTQtZXSh1TSp1USg0yM4sQaaVq/qqEdQ9jcPXBLDq6iFJLq7OyZWn47TeYOdO4U6hDB+P00GefwdWrZkcW6ZSZU0JmAKYADQAPoK1SysOsPEKkJUd7R0a/NppD7x4ir2teWi5vSYu17bjYuhFERMAPP4CnJwweDPnygb8/REaaHVukM2b2ACoCJ7XWp7XWD4BAoKmJeYRIc955vDnQ7QBf1P6CDcc34BHgwdzw+eh69YwniSMioH17WLAAypSBOnVg3TqIjTU7ukgHlFmjGCqlWgL1tdbd4l93ACpprd97aj1/wB/Azc3NJzAwMEn7i4mJwdnZOXmhrYy02bqcvXuW8cfHc+TmEXyz+9KveD/yOOYBwOHmTV5ev568a9aQ6Y8/+PPllznfrBmXGjTghlJW2+aksubvOamS02Y/P78QrbXvv97QWpuyAK2AWY+97gBMft5nfHx8dFIFBQUl+bPWStpsfWLjYvWUg1O082fO2mmMk560f5J+FPvonxUePNB62TKtq1XTGrTOnFmfb9RI68OHzQttAmv/npMiOW0GgnUCx1QzTwGdA/I99toduGBSFiEsgp2yo1eFXkT2iqRGgRr03dSX6nOrE30l2ljBwQFatYLduyEsDNq1w23bNvD2hurVjVtK5e4hkUhmFoBDQDGlVCGlVEagDbDWxDxCWIz8WfOz4e0NLGi+gGNXj+E93ZvRu0bzMPbhPyt5e8PMmexbtswYXuLiRXj7beOi8aefwpkzZsUXVsK0AqC1fgS8B2wGooFlWmu5zUGIeEop2nu1J7p3NM1KNmNo0FB8Z/oSciHkifUeubpCv35w/Dhs2gSVKsHYsca4Q40by0Vj8UymPgegtd6otS6utS6itR5jZhYhLNVLTi+xtOVSVrdezZU7V6g4qyIDtw7kz4d/PrminR3Uqwdr18Ivv8CQIRASAk2aGM8UjBgB586Z0gZhmeRJYCGsRLOSzYjqHUVn786M2zuOstPKsuvXZ8w+lj8/jBxpPFy2cqUxHPXw4VCgALzxhlEkHj1K0/zC8kgBEMKKZHPMxqwms9jWYRuP4h5Rc15NJpyYwK37zxhe2sEB3nzTeKbg1CkYNMjoFTRtahSJTz+FkyfTthHCYkgBEMIK1S5cm6M9j/Jh5Q9Zd2EdZQLKsPHExud/qHBhGDPG6BWsWQM+Psa1gmLFoFYt42Gzu3fTIr6wEFIAhLBSThmd+Lre13xb7ltcMrnQaHEjOqzuwB93/3j+B+3tjR7AunVGMRg9Gs6ehXfegTx5jGEn9u2TUUltgBQAIaych6sHof6h/K/G/wiMCMRjigfLIpf99YDl8+XNa4w3dOIE7NgBzZvDokXw6qtQsqQxGN3Zs6neBmEOKQBCpAOZ7DMxwm8EIf4hFMhWgNYrWtN8aXMu3E7ks5V2dlCzJsyfD5cuwZw54OZmFIcCBYwxiObPh9u3U7chIk1JARAiHfFy82Jf1318WfdLNp/ajMcUD2aFzkpcb+AvLi7QuTPs2mVcOB42zLittFMn4xRRu3bG8wZyF5HVkwIgRDpjb2dP/1f7c7TnUbzzePPuuneps6AOp6+ffvGNFS5sFICTJ2HPHmOegh9+gAYNjNNH778PBw7I9QIrJQVAiHSqaI6ibO+4nWmNpnHo/CE8p3oyYd8EYuOS8FSwUlC1KkybZgw5sWqVMfbQjBlQuTIULWo8eCZzG1sVKQBCpGN2yo7uvt2J6h2FX0E/PtryEVXnVCXycjJGXcmUybhYvGKFcb1g9myjp/D558bcxl5exu2m8nyBxZMCIIQNcHd1Z13bdSx+czGnrp+i3PRyjNw5kgexyRw5NFs26NIFtm6F8+fhm2+MawhDhhjPF/j6wrhxxjUEYXGkAAhhI5RStPVsS1SvKFp6tGTYjmH4zvDl0PlDKbODPHmgTx/46Sf49Vf48kvj7qKBA40eQoUKRjE4nYRrESJVSAEQwsbkdsrN4haLWdtmLdf+vEbl2ZXpv6U/dx+m4FPA+fND//5w8KBxwB871riOMHAgFCkC5coZD6BFR6fcPsULkwIghI16o8QbRPaKpFu5bny17yu8pnqx48yOlN9RoULw8cdGMfjlF2PugixZYOhQY5C6UqWMMYkOHZK7idKYFAAhbFhWx6xMf2M629/ZDoDffD+6r+vOzXs3U2eHBQsacxf89JNxzeDbb43bSceNg4oVjclsevUyBq+Tmc1SnRQAIQR+hfw40vMI/av0Z1bYLEoHlGbD8Q2pu9NXXoHevWHbNrh82XjSuFIl49/69SFXLnjrLVi4EK5dS90sNsqUAqCU+lIp9bNS6ohSarVSKpsZOYQQ/8jikIUvX/+SfV33kT1zdhovaUy7Ve24cudK6u88Rw5jMLqVK+GPP4yB6tq0+efhs9y5jecOxo6FiAg5VZRCzOoBbAXKaK29gOPAJyblEEI8pWLeioT4hzC85nCWRy7HI8CDwIjAFxtOIjkyZzamspwxw5jB7OBB47bSu3eN+Qw8Pancti307GkUijt30iZXOmRKAdBab4mfExhgP+BuRg4hRMIyZsjIsFrDCO0eSuHshWm7si1NA5ty7lYaTylpZ2fcPjpihDGRzblzMHMmt4sVM+YvaNLE6D3UqQNffSW9gxek0qyqPyuAUuuApVrrhc943x/wB3Bzc/MJDAxM0n5iYmJwdnZOck5rJG22Dand5lgdy6rzq5j9y2zslT3dC3en0cuNsFPmXUKMiYnBJWNGsh49So5Dh8h54ABOZ84AcD9XLq75+nLd15frPj48zJbNtJwpKTnfs5+fX4jW2vdfb2itU2UBtgERCSxNH1tnMLCa+EL0X4uPj49OqqCgoCR/1lpJm21DWrX55NWT2m+en2Y4uta8WvrE1RNpst+EJNjm337TetYsrVu10jp7dq2NvoDW5cpp3b+/1ps2aR0Tk+ZZU0pyvmcgWCdwTE21Eq61rqO1LpPA8j2AUqoj0BhoFx9QCGHBiuQowo/v/MiMxjMIvRiK51RPxu8dz6M4CxkWOl8+6NoVli2DK1eMawejR4OrK0yaZNxZlD071KgBw4cbw13fv292alOZdRdQfWAg0ERrLZOQCmEllFK86/MuUb2ieL3I6wzYOoAqs6tw5PcjZkd7UoYMxrWDwYONmc6uXzfmMPjgA+Ni8siRxgQ42bJB7dowapRREO7dMzl42jLrJN63gAuwVSl1WCk1zaQcQogkyOualzWt1xDYIpBfb/yKzwwfhgUN4/4jC/2L2skJ6tUzHjgLDoarV2H1aujRw3jGYNiwfwpCzZrGU8pbtqT7GdDszdip1rqoGfsVQqQcpRSty7SmduHafLj5Q0buGsnK6JXMbjKbSu6VzI73fNmzQ7NmxgJGEdizx+gF7NxpDG09erRxF1LZssZcCH8t+fKZmTxFyZPAQohkyZUlFwuaL2DD2xu4ef8mVWZX4aPNH3HngRXdn58jh3FL6fjxxphE168bPYDBg41iMWcOtG1rDHKXL5/xhPKECbBvn1WfNjKlByCESH8aFmtIZK9IBm0bxIT9E/j+2PfMaDyD2oVrmx3txbm4QN26xgLw8CGEhxsH/L17jWX5cuM9Bwejl1CpkjGeUcWKULy40XuwcJafUAhhNVwzuRLQKICdnXaSQWWgzoI6dFvbjRv3bpgdLXkcHIzJbfr0gSVLjPkOLlwwriN8+CE4OxtjGHXsaIxumi0b+PnBgAHG+sePQ1yc2a34F+kBCCFSXI0CNQjvEc7wHcMZv288G09sZGqjqTQt2dTsaCnn5ZefvI4QGws//2ycQgoONv6dPPmfW01dXIyeQrlyxuLtbQyHnSmTSQ2QAiCESCWZHTIztu5YWpVuRde1XWm2tBmtS7fmmwbf8JLTS2bHS3kZMhhzIpcuDZ06Gb97+BCiooxhLEJDISzMmEP5bvzd7/b2Ro/By8tYPD2NJW9eYwKdVCYFQAiRqnxf8eXQu4cY99M4Ru0axdbTW5lUfxLtPNuh0uAgZ6q/rg+ULWvMnQxGT+HkSeOawuHDxr+7dsGiRf98LmvWf4qJhweULo19Kjy0JtcAhBCpLmOGjAypMYSw7mGUyFmCDqs70HhJY87ePGt2tLSXIQOUKGHcSfTZZ7BhA/z2m3Er6s6dEBAA7doZvYNVq4xrDK+/jmtERIpHkR6AECLNeOT2YHfn3Uw5NIVPfvwEjwAPxtUZR3ff7qYOLmcR/hqmokaNJ39/+TJERXErJibFd2nj/8WFEGktg10G3q/0PhE9I6jsXpleG3tRa14tjl89bnY0y/TSS1CrFo9SYcRXKQBCCFMUyl6ILe23MKfJHI5ePorXVC/G7hlrOYPL2QApAEII0yil6FyuM1G9omhYrCGDfhxEpVmVCL8UbnY0myAFQAhhupddXmZV61WsaLWC87fO4zvTlyHbh3DvkfUOs2ANpAAIISxGC48WRPWOop1nO8bsHkO56eXYe3av2bHSLSkAQgiLkiNzDuY1m8emdpu4+/Au1eZUo+8PfYl5kPJ3wdg6KQBCCItUr2g9InpG0LtCbyYfnEyZgDJsObXF7FjpihQAIYTFcsnkwuSGk9nVeReO9o7UW1iPzt935vbD9D1RS1oxtQAopforpbRSKpeZOYQQlq1a/moc7nGYQVUHsSB8AZ2CO7EqepXZsayeaQVAKZUPqAv8ZlYGIYT1cLR35PM6n3Pw3YPkyJiDFsta0Gp5Ky7FXDI7mtUyswcwAfgY0CZmEEJYmfIvl2dqual89tpnrDu2Do8pHnwX/h1ay6HkRSkz/qMppZoAtbXWfZVSZwBfrfUfz1jXH/AHcHNz8wkMDEzSPmNiYnBOhUepLZm02TbYcpt/u/sbXx77kohbEVTIXoGPin9EHsc8ZsdLFcn5nv38/EK01r7/ekNrnSoLsA2ISGBpChwAssavdwbIlZht+vj46KQKCgpK8metlbTZNth6m2PjYvXkA5O10xgn7TTGSU8+MFnHxsWaFy6VJOd7BoJ1AsfUVDsFpLWuo7Uu8/QCnAYKAeHxf/27A6FKqfRZtoUQqcpO2fFexfeI7BVJtfzV6PNDH2rMrcGxP46ZHc3ipfk1AK31Ua31S1rrglrrgsA5oLzWWq7kCCGSrEC2AvzQ7gfmN5tP1JUoyk4ry+e7P+dh7EOzo1kseQ5ACJFuKKV4p+w7RPeO5o0Sb/Dp9k+pOKsiYRfDzI5mkUwvAPE9gQQvAAshRFK4ObuxvNVyVr61kksxl6gwswKfbPtEBpd7iukFQAghUsubpd4kqlcU75R9hy9++oKy08qy57c9ZseyGFIAhBDpWvbM2ZnTdA5b2m/hQewDqs+tznsb3+P2fRlOQgqAEMIm1C1Sl6M9j9K3Ul8CDgVQZmoZNp3cZHYsU0kBEELYDOeMzkysP5GfuvyEk4MTDRY1oOOajlz785rZ0UwhBUAIYXOq5KtCWPcwhlQfwuKji/GY4sGKqBVmx0pzUgCEEDYpk30mRr02iuB3g3F3dafV8la0WNaCi7cvmh0tzUgBEELYtLJ5yrK/237G1hnLxhMb8QjwYG7YXJsYXE4KgBDC5tnb2fNx1Y8J7xGO50uedFnbhXoL6/HL9V/MjpaqpAAIIUS84jmLs6PTDgIaBrD/3H7KTC3DpP2TiI2LNTtaqpACIIQQj7FTdvSs0JPIXpHUKliLDzZ/QPW51Ym6EmV2tBQnBUAIIRKQL2s+1rddz8LmCzl+9Tjlppdj9K7R6WpwOSkAQgjxDEop2nm1I6p3FM1LNmdo0FB8Z/oSciHE7GgpQgqAEEL8h5ecXiKwZSBrWq/hyp0rVJxVkYFbB/Lnwz/NjpYsUgCEECKRmpZsSlTvKLp4d2Hc3nGUnVaWXb/uMjtWkkkBEEKIF5DNMRszm8xkW4dtPIp7RM15Nem5vie37t8yO9oLM60AKKX6KKWOKaUilVLjzMohhBBJUbtwbY72PMqHlT9kesh0SgeUZuOJjWbHeiGmFACllB/G5PBeWuvSwHgzcgghRHI4ZXTi63pfs7frXlwzudJocSPar2rPH3etY44rs3oAPYEvtNb3AbTWl03KIYQQyVbZvTKh/qEMrTGUpZFL8ZjiwdKIpRY/nIRZBaA4UF0pdUAptVMpVcGkHEIIkSIy2WdipN9IQvxDKJCtAG1WtqH50uZcuH3B7GjPpFKrQimltgF5EnhrMDAG2A70BSoAS4HCOoEwSil/wB/Azc3NJzAwMEl5YmJicHZ2TtJnrZW02TZImy1PrI5lxbkVzDkzBwflQM8iPWmYpyFKqSRvMzlt9vPzC9Fa+/7rDa11mi/AJqDWY69PAbn/63M+Pj46qYKCgpL8WWslbbYN0mbLdeLqCV1zbk3NcLTfPD998urJJG8rOW0GgnUCx1SzTgGtAV4DUEoVBzIC1nHVRAghEqlojqJs77id6Y2nE3whGM+pnny972uLGVzOrAIwByislIoAAoGO8VVKCCHSFTtlh7+PP1G9o6hduDb9tvTj1TmvEnE5wuxo5hQArfUDrXV7rXUZrXV5rfV2M3IIIURacXd1Z22btSx+czGnr5+m/PTyjNgxggexD0zLJE8CCyFEGlFK0dazLVG9omhVuhXDdw7HZ4YPh84fMiWPFAAhhEhjuZ1ys+jNRaxru47rf16n8uzK9N/Sn7sP76ZpDikAQghhksbFGxPZK5J3y7/LV/u+wmuqFzvO7Eiz/UsBEEIIE2V1zMq0xtPY0XEHSin85vvRfV13bt67mer7lgIghBAWoGbBmoT3CGfAqwOYFTYLjwAP1h1bl6r7lAIghBAWIotDFsbVHceBbgfImTknTQKb0HZlW67cuZIq+5MCIIQQFsb3FV+C/YMZWWskK6NWUmpKKQ7fOJzi+5ECIIQQFihjhowMrTmUsO5h+LziQ97MeVN8H1IAhBDCgpV+qTSb228md6bcKb5tKQBCCGGjpAAIIYSNkgIghBA2SgqAEELYKCkAQghho6QACCGEjZICIIQQNkoKgBBC2ChlTTMxKqWuAL8m8eO5sL15h6XNtkHabBuS0+YCWut/PUlmVQUgOZRSwVprX7NzpCVps22QNtuG1GiznAISQggbJQVACCFslC0VgBlmBzCBtNk2SJttQ4q32WauAQghhHiSLfUAhBBCPEYKgBBC2Kh0VwCUUvWVUseUUieVUoMSeF8ppb6Jf/+IUqq8GTlTUiLa3C6+rUeUUnuVUmXNyJmS/qvNj61XQSkVq5RqmZb5Ulpi2quUqqWUOqyUilRK7UzrjCktEf9fZ1VKrVNKhce3ubMZOVOSUmqOUuqyUiriGe+n7PFLa51uFiADcAooDGQEwgGPp9ZpCPwAKKAycMDs3GnQ5leB7PE/N7CFNj+23nZgI9DS7Nyp/B1nA6KA/PGvXzI7dxq0+VNgbPzPuYFrQEazsyez3TWA8kDEM95P0eNXeusBVAROaq1Pa60fAIFA06fWaQp8pw37gWxKqZfTOmgK+s82a633aq2vx7/cD7inccaUlpjvGaAPsBK4nJbhUkFi2vs2sEpr/RuA1toW2qwBF6WUApwxCsCjtI2ZsrTWuzDa8SwpevxKbwUgL3D2sdfn4n/3outYkxdtT1eMvyCs2X+2WSmVF2gOTEvDXKklMd9xcSC7UmqHUipEKfVOmqVLHYlp87dAKeACcBToq7WOS5t4pknR45d9suNYFpXA756+zzUx61iTRLdHKeWHUQCqpWqi1JeYNk8EBmqtY40/EK1aYtprD/gAtYHMwD6l1H6t9fHUDpdKEtPmesBh4DWgCLBVKbVba30rlbOZKUWPX+mtAJwD8j322h3jr4MXXceaJKo9SikvYBbQQGt9NY2ypZbEtNkXCIw/+OcCGiqlHmmt16RJwpSV2P+v/9Ba3wHuKKV2AWUBay0AiWlzZ+ALbZwcP6mU+gUoCRxMm4imSNHjV3o7BXQIKKaUKqSUygi0AdY+tc5a4J34q+mVgZta64tpHTQF/WeblVL5gVVAByv+i/Bx/9lmrXUhrXVBrXVBYAXQy0oP/pC4/6+/B6orpeyVUlmASkB0GudMSYlp828YPR6UUm5ACeB0mqZMeyl6/EpXPQCt9SOl1HvAZoy7COZorSOVUj3i35+GcUdIQ+AkcBfjrwirlcg2/w/ICQTE/0X8SFvxSIqJbHO6kZj2aq2jlVKbgCNAHDBLa53grYTWIJHf8ShgnlLqKMapkYFaa6seIloptQSoBeRSSp0DhgEOkDrHLxkKQgghbFR6OwUkhBAikaQACCGEjZICIIQQNkoKgBBC2CgpAEIIYaOkAAiRSEqpM0qpXMldRwhLIQVACCFslBQAIRKglFoTP6hapFLK/6n3CiqlflZKzY8fk31F/NO3f+mjlApVSh1VSpWM/0zF+LkYwuL/LZGmDRIiAVIAhEhYF621D8aYQu8rpXI+9X4JYIbW2gu4BfR67L0/tNblgalA//jf/QzU0FqXw3gy+7NUTS9EIkgBECJh7yulwjHmT8gHFHvq/bNa65/if17IkyOsror/NwQoGP9zVmB5/ExPE4DSqRFaiBchBUCIpyilagF1gCpa67JAGOD41GpPj6Hy+Ov78f/G8s94W6OAIK11GeCNBLYnRJqTAiDEv2UFrmut78afw6+cwDr5lVJV4n9uC+xJxDbPx//cKUVSCpFMUgCE+LdNgL1S6gjGX+77E1gnGugYv04OjPP9zzMO+Fwp9RPG6JZCmE5GAxXiBSmlCgLr40/nCGG1pAcghBA2SnoAQghho6QHIIQQNkoKgBBC2CgpAEIIYaOkAAghhI2SAiCEEDbq/4EmVzUvE4JEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return 2*x[0]**2 + x[0]*x[1] + x[1]**2 + x[1]*x[2] + x[2]**2 - 6*x[0] - 7*x[1] - 8*x[2] + 9\n",
    "\n",
    "def df(x):\n",
    "    return np.array([4*x[0]+x[1]-6, x[0]+2*x[1]+x[2]-7, x[1]+2*x[2]-8])\n",
    "\n",
    "x0= np.zeros(3)\n",
    "p = np.array([1.0, 2.0, 0.0])\n",
    "c1= 4/5.\n",
    "\n",
    "def phi(alpha):\n",
    "    return f(x0 + alpha*p)\n",
    "\n",
    "def ell(alpha):\n",
    "    return f(x0) + c1*alpha*df(x0).dot(p)\n",
    "\n",
    "alpha = np.linspace(0,1.0,100)\n",
    "plt.plot(alpha, ell(alpha), 'g', label='$\\ell$')\n",
    "plt.plot(alpha, np.vectorize(phi)(alpha), 'r', label='$\\phi$')\n",
    "plt.xlabel('alpha')\n",
    "plt.ylabel('f')\n",
    "plt.grid(True)\n",
    "plt.legend(loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After this, we could run for example the classical bisection algorithm to find out more precisely, where the lines intersect.  Indeed, we can see that for $\\alpha = 0.2$ the green line is above the red (sufficient decrease is satisfied) whereas for $\\alpha = 1.0$ is below.  This could be our initial bracket."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current bracket: [2.000000e-01, 1.000000e+00]\n",
      "Current bracket: [2.000000e-01, 6.000000e-01]\n",
      "Current bracket: [4.000000e-01, 6.000000e-01]\n",
      "Current bracket: [4.000000e-01, 5.000000e-01]\n",
      "Current bracket: [4.500000e-01, 5.000000e-01]\n",
      "Current bracket: [4.750000e-01, 5.000000e-01]\n",
      "Current bracket: [4.875000e-01, 5.000000e-01]\n",
      "Current bracket: [4.937500e-01, 5.000000e-01]\n",
      "Current bracket: [4.968750e-01, 5.000000e-01]\n",
      "Current bracket: [4.984375e-01, 5.000000e-01]\n",
      "Current bracket: [4.992187e-01, 5.000000e-01]\n",
      "Current bracket: [4.996094e-01, 5.000000e-01]\n",
      "Current bracket: [4.998047e-01, 5.000000e-01]\n",
      "Current bracket: [4.999023e-01, 5.000000e-01]\n",
      "Current bracket: [4.999512e-01, 5.000000e-01]\n",
      "Current bracket: [4.999756e-01, 5.000000e-01]\n",
      "Current bracket: [4.999878e-01, 5.000000e-01]\n",
      "Current bracket: [4.999939e-01, 5.000000e-01]\n",
      "Current bracket: [4.999969e-01, 5.000000e-01]\n",
      "Current bracket: [4.999985e-01, 5.000000e-01]\n",
      "Current bracket: [4.999992e-01, 5.000000e-01]\n",
      "Current bracket: [4.999996e-01, 5.000000e-01]\n",
      "Current bracket: [4.999998e-01, 5.000000e-01]\n",
      "Current bracket: [4.999999e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Current bracket: [5.000000e-01, 5.000000e-01]\n",
      "Final alpha = 5.000000e-01\n"
     ]
    }
   ],
   "source": [
    "alpha_min = 0.2\n",
    "alpha_max = 1.0\n",
    "\n",
    "while alpha_max - alpha_min > 1.0E-10:\n",
    "    alpha = 0.5*(alpha_min + alpha_max)\n",
    "    print('Current bracket: [%e, %e]' % (alpha_min, alpha_max))\n",
    "    if ell(alpha)>phi(alpha):\n",
    "        alpha_min = alpha\n",
    "    else:\n",
    "        alpha_max = alpha\n",
    "\n",
    "print('Final alpha = %e' % alpha)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}