{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Root Finding \n", "\n", "## Overview, Objectives, and Key Terms\n", " \n", "In [Lecture 19](ME400_Lecture_19.ipynb), SymPy was applied to perform symbolic differentuation, and in [Lecture 21]( ME400_Lecture_21.ipynb), finite-difference approximations were developed to perform numerical differentiation. We'll need both techniques as we begin to solve **nonlinear equations** and **optimization problems**. Both subjects are rich, so we'll touch on only the basics, but you'll have tools at your disposal to tackle such problems in varied applications.\n", " \n", "### Objectives\n", "\n", "By the end of this lesson, you should be able to\n", "\n", "- Find one or more roots of a one-dimensional, nonlinear equation $f(x) = 0$ using the bisection and Newton methods.\n", "- Find local extrema of a function $f(x)$ using the bisection and Newton methods.\n", "- Use `fsolve` to solve nonlinear systems.\n", "- Use `minimize` to solve nonlinear optimization problems.\n", "\n", "### Key Terms\n", "\n", "- nonlinear equation\n", "- transcendental equation\n", "- graphical solution\n", "- bisection method\n", "- Newton's method\n", "- quadratic convergence\n", "- second-order convergence\n", "- order of convergence\n", "- secant method\n", "- Steffensen's method\n", "- nonlinear systems\n", "- Jacobian matrix\n", "- `scipy.optimize.fsolve`\n", "- optimization\n", "- extremum\n", "- critical point\n", "- objective function\n", "- `scipy.optimize.minimize`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Nonlinear Equations\n", "\n", "Linear equations are characterized by *linear combinations* or the unknowns. For example, the system of equations for $x$ and $y$\n", "\n", "$$\n", "\\begin{split}\n", " ax + by &= 1 \\\\\n", " cx + dy &= 2 \n", "\\end{split}\n", "$$\n", "\n", "is linear because $x$ and $y$ appear with only constant coefficients. There are no $x^2$ terms, or $\\sin(y)$ terms. Just multiples of $x$ and $y$. Any deviation from this pattern results in **nonlinear equations**. For example, quadratic equation $ax^2 + bx + c = 0$ is nonlinear, and we know how to find the **roots** of this equations. A root is the solution $x$ to any equation (generally, nonlinear) $f(x) = 0$. Our goal now is to solve such equations using the tools at our disposal." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analytic Solutions\n", "\n", "In rare cases, nonlinear equations have closed-form solutions (i.e., you can write it down explicitly). The quadratic equation is one example. In fact, solutions exist for the general cubic (add a $x^3$ term) and quartic (add both $x^3$ and $x^4$). The roots are hideous (compared to the quadratic equation roots), but they are known. \n", "\n", "***\n", "\n", "**Exercise**: Use SymPy to find the roots of the cubic equation $ax^3 + bx^2 + cx + d = 0$. Can it also find roots to *quintic* equation $ax^5 + bx^4 + cx^3 + dx^2 + ex + f = 0$?\n", "\n", "*Solution*: The cubic is easy, but SymPy can't do the quintic. In fact, the roots to general polynomials of degree 5 and higher cannot be determined explicity, a fact [proven by Abel](https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem).\n", "\n", "***\n", "\n", "For other problems, we often have to have a bit of luck.\n", "\n", "***\n", "**Exercise**: Find all values of $x$ such that $\\sin(ax) - b = 0$.\n", "\n", "*Solution*: The equation has no solution if $|b| > 1$ (because $-1 \\leq \\sin(ax) \\leq 1$). If $|b| \\leq 1$, then $x = \\frac{\\sin^{-1}b}{a} + n\\pi$ for $n = 0, \\pm 1, \\pm 2, \\ldots$. Such an equation is called **transcendental** because it involves a transcendental function (i.e., $\\sin(x)$). Transcendental functions cannot be defined in terms of polynomial of finite degree; in other words, if one needs an infinitely long [Taylor series](ME400_Lecture_20.ipynb) to represent the function, it is transcendental. Transcendental equations almost always require some form of numerical evaluation.\n", "\n", "**Exercise**: Solve the nonlinear system of equations $y - x = 4$ and $x^2 + y = 3$.\n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Graphical Solutions\n", "\n", "When analytic methods fail, solutions (albeit, approximate) to nonlinear equations may be found **graphically**. For single equation in $x$, it suffices to plot $f(x)$ over the range of interest and identify the roots.\n", "\n", "***\n", "\n", "**Exercise**: Revisit $\\sin(ax) - b$ for $a=1$ and $b=1/4$ and find all roots between $x = -5$ and $x = 5$.\n", "\n", "*Solution*:\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "x = np.linspace(-5, 5, 100)\n", "ax = plt.figure(1, (8,3)).gca()\n", "plt.plot(x, np.sin(x)-0.25, x, 0*x,'r--')\n", "plt.xlabel('x')\n", "plt.ylabel('f(x)')\n", "ax.set_xticks(np.arange(-5, 6, 1))\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because `plt.plot` does not by default draw a line for $y=0$, a dashed, red line is superimposed to provide more obvious intersections between the function $f(x)$ and the horizontal axis $y = 0$. It appears (from inspections) that the roots of interest are near -3.5, 0.3, and 2.8. One could certainly zoom in somewhat to improve the precision of the estimates, but arbitrary precision will require other techniques.\n", "\n", "***\n", "\n", "**Exercise**: A simple, somewhat realistic model of a spherical nuclear reactor consists of two regions: (1) an inner spherical volume of *fuel* of radius $R$, in which neutrons induce fission, leading to energy release and additional neutrons, and (2) an outer, infinitely large *reflector* surrounding the fuel, which minimizes the number of neutrons that escape the fuel. The reactor is *critical* if the number of neutrons produced and lost is constant over time. For this simple, two region model, the condition for criticality is\n", "\n", "$$\n", " BR \\cot(BR) = 1 - \\frac{D_r}{D_f} \\left ( \\frac{R}{L_r} + 1 \\right ) \\, ,\n", "$$\n", "\n", "where $D_f$, $D_r$, and $L_r$ are fixed, material properties of the fuel and reflector. Hence, in order for this criticality condition to be satisfied, an appropriate value of $B$ must be determined. Here, $B$ is called the *buckling* and is a function of several fuel material properties. For $R=20$, $D_f = 0.9$, $D_r = 1.1$, and $L_r = 1.7$, estimate the first positive value of $B$ that satisfies the condition.\n", "\n", "***\n", "\n", "**Exercise**: The solution of transient heat-conduction problems often involves a technique called *separation of variables*. Application of this technique to an infinite slab of thickness $2L$ subject to certain initial and boundary conditions leads to the transcendental equation\n", "\n", "$$\n", " \\cot \\lambda L = \\frac{\\lambda L}{hL/k} = \\frac{\\lambda L}{\\text{Bi}} \\, ,\n", "$$\n", "\n", "where $\\lambda$ is a dimensionless, undetermined parameter, $h$ is the heat transfer coefficient, $k$ is the thermal conductivity, and $\\text{Bi}$ is the *Biot number*. The Biot number quantifies how hard it is for heat to flow *within* a body relative to how hard it is to flow *through the outer surface* of the body. Only for certain values of $\\lambda$ can the equation be solved. Determine these values graphically by plotting the left-hand and right-hand sides of the equation as functions of $\\lambda L$.\n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Determining Roots Numerically\n", "\n", "The numerical solution of $f(x) = 0$ can be challenging since it requires we know at least a little about the desired solution, particularly the range in which the root is contained. We'll explore three basic schemes for solving $f(x)$: the [bisection method](http://mathworld.wolfram.com/Bisection.html), [Newton's method](http://mathworld.wolfram.com/NewtonsMethod.html), and the secant method.\n", "\n", "### Bisection\n", "\n", "Application of bisection to finding a root of $f(x)$ requires a key piece of information: a range $x \\in [L, R]$ in which a *single* root lives. The algorithm is actually quite similar to [binary search](ME400_Lecture_14). A central point $C$ is selected, and the function is evaluated at that point. Then, if $f(C)$ has the same *sign* (i.e., $+$ or $-$) as $f(A)$, we know that the root must be between $C$ and $B$. Why? If there is only one root $x^*$ between $L$ and $R$, then $f(x)$ must have a constant sign for $x \\in [A, x^*)$ and the opposite sign for $x \\in (x^*, B]$. If, instead, $f(x)$ changes signs between $x=A$ and $x=x^*$ then it *must* cross the $x$ axis, at which point there *must* be a second root. \n", "\n", "***\n", "\n", "**Exercise**: Write a `lambda` function for $x$ that returns `True` if $x \\geq 0$ and `False` otherwise.\n", "\n", "***\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An important difference between bisection and binary search is that now the search can require an unlimited number of comparisons because the search range $[A, B]$ is split in half repeatedly until a root is *converged* to within a desired tolerance $\\tau$. In other words, roots will very rarely be found *exactly*; rather, we often must accept *approximate* values, as was first introduced [earlier](Loops_in_Python.ipynb). \n", "\n", "In pseudocode, the bisection method can be written as follows (where the notation from binary search is adopted as closely as is reasonable):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "'''Bisection method for finding an isolated root of f between a and b'''\n", "Input: f, a, b, tau\n", "\n", "# Define the left and right boundaries\n", "Set L = a\n", "Set R = b\n", "\n", "# Define the central point\n", "Set C = (a+b)/2\n", "\n", "# Go until f(C) is close enough to zero (i.e., that\n", "# C is close enough to the root between a and b\n", "While |f(C)| > tau \n", " If sign(f(C)) == sign(f(L)) then\n", " # the root must be between C and R\n", " L = C\n", " Otherwise\n", " # the root must be between L and C\n", " R = C\n", " C = (L+R)/2\n", "Output: C\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Application of this algorithm to find the first positive root of $\\sin(x) - 2/5$ is shown graphically below for the first two iterations. At each step, the search space is reduced by half, illustrated by the shrinking highlighted rectangles. By the second iteration, $C = 0.375$, which is visually quite close to the root of interest $x = \\sin^{-1}(2/5) \\approx 0.4115$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEPCAYAAACKplkeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xdc1dUfx/HXYSpoDhDcmSs1cyKKlqXpL9PSSkvNkSM1K8WRIrn33pm5lXJk5qocZalZWa7UtKE4cKOiIkNlnd8fF02WXOByv/fC5/l48PDe+z3f7/18geub8x3nKK01QgghhCU4GF2AEEKInENCRQghhMVIqAghhLAYCRUhhBAWI6EihBDCYiRUhBBCWIyEihBCCIuRUBFCCGExEipCCCEsxsnoAqzN09NTlylTxugyhBDCIg4ezL5t16798PscvK61LpLeOrkuVMqUKcOBAweMLkMIISxCqezb9sP/VSqlQsxZRw5/CSGEsBgJFSGEEBYjoSKEEMJiJFSEEEJYTK47Uf8ot2/f5urVq8TGxj547dlnKxAWZvlvk4dHHHv2nLT4dkXu4e7uTsmSJXFwkL8Nhe2QUEl0+/ZtQkNDKVGiBHnz5kUlXlIRFpY97xcW5kTlypWzZ+Mix0tISODixYtcv34dLy8vo8sR4gH5EyfR1atXKVGiBG5ubg8CRQhb5eDggLe3N+Hh4UaXIkQSEiqJYmNjyZs3r9FlCGE2Z2dn4uLijC5DiCTk8NdDrN1DeejUjRCZoIiPh+vXja5DiP9IqBgoI+dXq1V7ijlz5vH8889nqV2LFi/Rtm07Ond+O933LFeuDAsWLKZJkybmFyqsSilwdTW6CiH+I6FiJ44ePZ7hdqNHj+LUqWCCgj5/8Nq33261eG1CCHGfnFMRQghhMRIqdqJcuTLs2LGD0aNH0a7dm3Tp0pmCBfNTrdpTSQbIvN9u27ZtTJo0gbVrv6BAgXzUqlUdgMaNn2fJksUAnDp1iiZNGuPl5YG3tyedOnXg1q1bBuydECKjjh49TK9e6R/GtjYJFTv09debefPNdoSF3eLll1vSt+8HKdo0a9aMIUM+4s032xIeHsmhQ0dStNFaExAQyPnzlzh27G/Onz/PmDGjrLAHQojMSEhIYOvWb3j55cY880xNNm/+yuiSUrDpUFFKNVNK/auUClZKDXlEuzZKKa2U8rFmfUZp0OAZmjdvjqOjIx07duLo0ZSBYY7y5cvTtGlTXF1dKVKkCP37D+Cnn3ZbuFohRFbduXOHZcsWUqdOFdq2fYVTp04yduwU/vrrvNGlpWCzJ+qVUo7APKApcAHYr5TarLX+K1m7/EBf4HfrV2mMokWLPnjs5ubG3bt3iYuLw8kpYz/Oq1ev0q9fX37+eQ8REREkJCRQqFAhS5crhMiksLDrLFr0CQsXfsz169eoWbM2S5as4tVX2+Ds7Gx0eamy5Z6KLxCstT6ttY4B1gCtUmk3FpgC3LVmcfYgvftuhg4NRCnFH38c5ebN2wQFfY7W2krVCSHScubMaQYO/IAqVUozYcJIfHzqsmXLLnbt2s8bb7S32UAB2w6VEsDDfbsLia89oJSqCZTSWn9jzcLshbe3N2fPniUhISHV5REREeTLl4+CBQty8eJFpk2bauUKhRAPO3z4EF26tKNmzQosX76QNm3as2/fcdau/ZpnnnnOLoaQsuVQSe279+DPaKWUAzATGJjuhpTqqZQ6oJQ6cO3aNQuWaNvatHkDAC8vD+rUqZVi+fDhI/njj0MULlyAli1b8Nprr1u7RCFyPa01u3f/SKtW/6Nhw9rs2LGVPn0GcuzYWebNW0KlSlWMLjFDlK0e7lBK+QGjtNYvJj4PBNBaT0x8XgA4BUQmrlIUuAG01FqnOQm9j4+PTm2O+r///jvVUYOz8w8DGbZJZNU///xN6dIy2rU9SkhIYMuWzUyfPpGDB/fh7V2U997rR7du71KgQAGzt1O+PFy9avn6vL3hypX/niulDmqt070YymZP1AP7gQpKqSeAi0A74K37C7XW4YDn/edKqV3Ah48KFCGEMFpcXBzr1q1hxoyJ/PPPX5QpU5ZZsz7lrbfeJk+ePBneXnCw+W3v3QNPz/TbZYXNhorWOk4p9QGwHXAElmqtjyulxgAHtNabja1QCCHMd+/ePVauXM6sWZM5e/YMVapUZfHilbz++psZvnLTltn0nmittwBbkr02Io22z1ujJiGEyIjo6GiWLVvInDlTuXz5Ej4+dZk8eTYvvtgiR87aadOhIoQQ9ioiIoLFiz9h7tzpXL9+jWeffZ4FC4J47rnGdnEVV2ZJqAghhAWFh4ezYMFc5s2byc2bN3jhhRcZPHgYfn7PGF2aVUioCCGEBYSHhzN//mw++WQmt27d4sUXWxAQMAIfH1+jS7MqCZV0eHtDaOijWmj82Isv+8hPBBHkZx++7MWP1G+1+W+7Qgj7lzxMWrRoRUDACGrUSHlvWG4goZKOh6/TBtMUwA4Opgdq6RLUtCmmi8RjY01fzs6mLy8v9IeD0d26m54LIXKU27dvM3/+bObNm8GtW7d4+eVXCQgYQfXqNY0uzVA579IDa4iMxKFpY9SggagzZ1BRUaiYGJTWpn+jokyvDxqIQ9MXIDIy/W1m0fDhw/D29qREiaLpN86kjRs3UKZMKQoUyMcff/xBtWpPsWvXrmx5r4kTJ9Cz5zvZsm1ztWjxEkFBKwytQdieiIgIpk2bQNWqZRg/fgT16zdkz55DrFq1IdcHCtjwHfXZJaN31CcXGx2L04uNYf9+1L176bbXrq5Qx5eE73/Ith7L+fPnqVy5IqdPh+Dl5cXZs2cpX/4J7t6Ntej17xUrlmPatBm0bJnauJ7my2h92bU/D0tt6mV7IHfUW090dDSLFs1j5szJ3LgRRrNmLxMYOIqaNWsbXZrZsnLzo7l31EtPJYMcli2BQ4fMChTA1O7QQdSypdlWU0hICB4eHnh5eVlke3FpjB8TEhJClSpPZWkbRrClWoT9uXfvHp9+Opdq1coyfPhgatb04ccff2ft2q/tKlCsRUIlI7TGYfoUVHR0hlZT0dGoqVMgC73CyZMnUbFiOQoWzM/TT1dh48YNAOzYsYNmzZpy6dIlChTIR7duXWjUqCEAHh4FKVAgH3v37gVg2bKlVK1aGU/PQrz00ouEhIQ82L6Tk+KTT+ZRqVIFKlWqkOS97927R4EC+YiPj6dWrepUrFgO+G/qYjD9pf/mm23o3LkjhQo9xooVy9m3bx916/pQqNBjFC/uzcCBAwDSrO9ho0ePonPnjo9sn9H96d/fnzJlSlGo0GP4+tZmz549AGZNvZyQkMD48eMoW/ZxihXzokuXzoSHhwOmnpSTkyIoaAVPPFEab29PJkwYn4GfrrBFcXFxrFixmJo1KzB4cF8qVqzEtm0/sWHDtlx3RVdGSKhkxN69mR+57Wqoaf1MKleuHLt27eHGjXCGDx9J584duXz5Mk2aNOGbb7ZSvHhxwsMjWbp0OTt3/gRAWNgtwsMj8fPzY9OmjUyaNIEvv1zPlSvXeOaZZ+nQoX2S99i0aSO//vo7f/6ZZB40XF1dCQ83nRc6dOgIJ06cSrXGzZs30bp1G8LCbvHWWx3o39+fPn38uXnzNidOnOKNN94ESLW+R7HU/vj41OHgwcNcu3aDdu3eol27N7h7965ZUy+vWLGcoKDl7Nixk5MnTxMZGZliGudffvmZv/76l++++4Fx48bw999/P3K/hG1KSEhg3bo11KlThT59euDtXYxNm77n2293Ur/+s0aXZ/MkVDJi3z7TFV6ZEReHOrA/02/dps0bFC9eHAcHB958sy0VKlRg//59Zq+/cOECAgICqVy5Mk5OTgQGfsSRI4eT/HUfEBBI4cKFyZs3b6ZqrFfPj1atXsXBwYG8efPi7OxMcHAw169fJ1++fNSrVy9T201NZvanQ4eOeHh44OTkxIABA7l37x7//vuvWe+3atVK+vUbQNmyZcmXLx/jx0/kiy/WJDm0Nnz4SPLmzUv16tWpVq16pqd5FsbQWrN16zc880xNunVrT548eVizZhM//vgbjRo1ydF3wVuShEpGRERkPlRiYkzrZ9JnnwVRu3YNPDwK4uFRkGPHjnH9+nWz1z93LoT+/f0frF+kSGG01ly8ePFBm1KlSmW6vtTWX7RoCSdPnuCppypRr14dvvnGcnOpZWZ/ZsyYTtWqlSlcuAAeHgUJDw83+3t4+fIlSpd+/MHzxx9/nLi4OEIfuokp+TTPkVa46k9Yxq+/7uHFF5+lbdtXiIqKYvHilfzyy2GaN28pYZJBcp9KRuTPb7qCKyYm4+u6uJjWz4SQkBB69erBd9/9gJ+fH46OjtSuXSPNqX9T+xCULFmKwMChvPVWhzTfJ6sfnuTrV6hQgZUrV5OQkMCGDetp27YNV6+GZfh9LLE/e/bsYerUyXz33Q889dRTODg44OlZ6MH3ML2aihUrzrlz//WCzp07h5OTE97e3ly4cCFD+yNsx7FjRxk1KpDvvttC0aLFmDXrUzp16mbT0/XaOumpZISvb+YvC3ZyQvvUydSqUVFRKKUoUqQIAMuXL+PYsWNpti9SpAgODg6cPn36wWu9er3L5MkTOX78OGC6C3jdui8zVY+5Vq78nGvXruHg4ECBAgUBcHR0TLW+R7HE/kRERODk5ESRIkWIi4tj7Ngx3L59+8Hy9KZebteuPbNnz+TMmTNERkYybJjpHExOGrI8Nzl79gw9enSiQYMa7Nv3K6NHT+Lw4WC6deslgZJFEioZ4ecHmb1s19vbtH4mVKlShf79B/LMM34UL+7Nn3/+Sf36DdJs7+bmRmDgUBo2bICHR0F+++03Xn31NQYNCqBDh3YUKvQY1atXZdu2rZnbFzNt376NatWeokCBfAwY4M+qVWvIkydPqvU9iiX258UXX6RZs5eoXLkiZcs+Tp48eZIcHktv6uWuXbvRoUMnGjVqSPnyT5AnTx5mz56bye+MMEpY2HUCAvpRu/aTbNq0jn79BnPkyGn69w/Azc3N6PJyBLn5MZG5Nz/Gz/sUh8EDM3RZsXZzQ0+bge7ZK0O1CpEeufnRPNHR0XzyySxmzZpMZGQknTp1IzBwFMWLlzC6NKuSmx9tUELX7lCrlulOeTNoV1eoXRvdtVs2VyaESC4+Pp7PPltKzZoVGDNmKM8+24jffvuTuXMX5bpAsRYJlYxydibhm61QxxedTndZu7mBry8JX2+RQSWFsCKtNdu3b6FBgxq8/353SpYszfbte1i9eiOVKlUxurwcTUIlM/LlI+H7H0yHtJ4oi3Z3R7u6opUy/evuji5bFj1tBgnf/QD58hldsRC5xpEjf9CqVVPeeKMFd+/eJSjoS3bs+DXXTJJlNLl0JbOcndE9e6F79IS9e003NkZEQP786Dq+UK8eyPXtQljNxYsXGDNmKGvWfEahQoWZMmUO3br1wsXFxejSchUJlaxSCurXR9evb3QlQuRKkZGRzJw5mY8/nk58fDx9+37IwIEfUbBgQaNLy5UkVIQQdik+Pp6VK5czduwwQkOv0Lp1O0aNmsjjj5cxuLLcTUJFCGF3du/+kcDA/hw7dhRfXz9WrtyAr6/lxpYTmScn6oUQdiM4+CTt2rXilVde4PbtcJYtW8P33/8igWJDJFQyofah/bwf/C+XY8ybqMsarDGdsC1r3Pi/uU+Sk6mJ7d+tW7f46KOB1K37FD/99COjRk3kwIF/aN26rQz4aGPk8FcmHImK4u/oaFaEhvK2tzfDSpehmIt5N0Nmh/PnzzNz5vRsn07YXgUGfvTgsVFTE3/7bfYOiZNTxcfHs2LFYsaOHcaNG2F06tSN4cPH4e2dO/94sgfSU8mkGK25m5DAsitXqLD/d0N7LtaaTljI98aa9uzZxbPP1qJfv3epVKkKP/10kI8/XiyBYuMkVLLIWuFi5HTCAHfv3qVz5454eXng4VGQevXqPJhLJDw8nB49ulOyZDFKly7B8OHDiI+Pf7Du4sWLqFq18oPaDx06BJjGW2vc+Hk8PApSrdpTfP315gfrdOvWhT593ueVV1pQsGB+/PzqcurUfzNOfv/99zz1VCUKFy5A374fpDkNAMjUxPYmJOQsnTu/QYsWjbh9O5ygoC/ZsmUX1avXNLo0YQYJFQvJ7nAxcjphgKCgFYSHh3P27HmuXg1j3rxPH8yo2LXr2zg5OfHvv8EcOPAH33//3YP/QNet+5IxY0axbFkQN2/eZsOGzXh4eBAbG8urr75C06b/4/Llq8yaNZdOnTokmYlxzZrVDB8+kuvXb1K+fHmGDx8KwPXr13nzzdaMGTOO0NDrlC1bjl9//cWs76NMTWy7oqOjGTduBHXqVGb79m8ZOnQM+/f/zauvtpHzJnZEQsXCsitcjJ5O2NnZmbCwMIKDgxMnCavNY489RmhoKNu2bWXGjFm4u7vj5eVFv379+eKLNQAsWbKYDz8cTJ06dVBKUb58eR5//HF+++03IiMjCQgYgouLC40bN6ZFi5dZs2b1g/d87bXX8fX1xcnJifbtO3DkyGEAtm7dQuXKVWjdug3Ozs74+/dLMutiRsnUxMbSWrN+/Vp8fCoxZcpYXn75NQ4e/JeAgOGZntpaGEdCJZvcD5dFly/z1t8p//LPKKOnE+7YsRP/+9+LdOjQjlKlihMQMJjY2FhCQkKIjY2lZMliD7bdu3cvrl27CsCFC+cpW7Zciu1dvnyJUqVK4eDw369g6dKPc+nSf/U8fOz84el5L126lKRWpRQlS2Z+KmSZmtg4x4//ycsvN6ZLl7YULuzBtm0/sXTpqiz9PIWx5NKgbOKiFI5K8bZ3UYY99B9IZtjCdMLOzs6MGDGSESNGcvbsWV55pTlPPvkkL73UHFdXV0JDr6d6NVXJkqU4ffpUiteLFSvO+fPnSUhIeBAs58+fo0KFimnW8N+6xTh//vyD51prLlw4/4g1/iNTE9uGW7duMX78CBYtmkeBAgWZOXM+Xbr0wNHR0ejSRBZJT8XCXJQir4MD3YoW42SdunxcviJFs3i5sS1MJ7xz507+/PNP4uPjeeyxx3B2dsbR0ZFixYrRtOn/GDRoILdv3yYhIYFTp06xe/duALp3f4cZM6Zx8OBBtNYEBwcTEhJC3bp1cXd3Z+rUKcTGxrJr1y6++eZr2rZtl24tzZu34K+/jrNhw3ri4uKYO3cOV65cMWs/ZGpiYyUkJPDZZ0upVasiixbNo2vXXhw6dILu3d+VQMkhJFQsJDvC5D5bmE44NPQKbdu2oVChx6hatTINGz5Hhw6mK6qWLw8iJiaGp5+ugqdnIdq2bcOVK5cB07mgwMChdOr0FgUL5qd161e5ceMGLi4ubNiwmW3btuLt7UmfPu+xfHkQlSpVSrcWT09P1qz5ko8+GoKXlwfBwScf+f2w9PdGpibOnEOHDtCkSX3ef7875cpVYPfuA8yc+QkeHh5GlyYsSKYTTmTudMKxseDyy64Hz5Mf5rJUkAhhDnuYTvjGjRuMHTuUpUsXUKSIF2PGTKF9+05yRZcBrDGdcM7vb2cTCRMhHu3+oa6RI4cQHn6L3r39CQwcRYECBYwuTWQjCZVMqO7ujt9jBSRMhEjD4cOHGDDgPQ4c+B0/v2eYPn0eVatWM7osYQUSKplwsFYdo0sQwiaFh4czbtxwFi2ah4eHJwsWBNGuXUc51JWL2PSJeqVUM6XUv0qpYKXUkFSWD1BK/aWUOqqU+kEplbVrd4UQmaK1Zu3aVdSu/SSLFs2je/feHDz4r5w7yYVsNlSUUo7APOAloArQXilVJVmzPwAfrXU1YB0wxbpVCiFOnPiHli2b8M47HShV6nF27drP9Okfy3S+uZTNhgrgCwRrrU9rrWOANUCrhxtorXdqraMTn/4GlMzKG6Z1X4EQtsjoKzfv3LnDmDHD8POrxpEjh5g5cz47dvxKjRopL6MWuYcth0oJ4OHbpC8kvpaW7kCmJ61wd3fn4sWLxMTEGP5hFSI9Wmtu3AjDxSWPIe///ffbqFu3KtOmjad163YcOPCP3MAoANs+UZ/agdhU/7dXSnUEfIDn0ljeE+gJULp06VTfrGTJkly/fp2QkJBHzpkRHw9yiFjYAheXPHh6ZqlznmGXL19iyJB+bNjwJRUqPMk33/xIw4aNrFqDsG22HCoXgIdHlSsJXEreSCnVBBgKPKe1TnU4YK31QmAhmG5+TK2Ng4MDXl5e6U50df06uMpVxCKXiY+PZ8mSTxk9OpCYmBiGDRuLv/8gXOXDIJKx5VDZD1RQSj0BXATaAW893EApVRNYADTTWl+1folC5HxHjx7G378XBw/uo1GjpsyY8QnlypU3uixho2z2nIrWOg74ANgO/A2s1VofV0qNUUq1TGw2FcgHfKmUOqyU2pzG5oQQGRQVFcWwYYN47jkfzp07y5Ilq9i4cbsEingkW+6poLXeAmxJ9tqIhx43sXpRQuQC3323lQEDenPuXAhduvRg9OjJFCpUyOiyhB2w6VARQljX1auhBAT046uv1vDkk5XZvn0Pfn7PGF2WsCM2e/hLCGE9WmuCgpbg41OJr79ez9ChY/j55z8kUESGSU9FiFzu5MkT+Pv35Oefd9OgQUPmzFlIhQpPGl2WsFPSUxEil4qJiWHq1PHUr1+NY8eOMHfuIr79dqcEisgS6akIkQsdOLCPPn3e4fjxP3n11TZMnToXb++iRpclcgDpqQiRi0RFRREYOIAmTfy4efMGq1dvJCjoSwkUYTHSUxEil9i5cwf+/j05e/YM3bu/y6hRk2QWRmFx0lMRIoe7efMm773XjVatmuLk5MzWrbuZOXO+BIrIFhIqQuRgmzevx9e3CqtXBzFgwBB+/fUIDRo0NLoskYPJ4S8hcqDQ0CsMGtSHjRvXUb16Tdat20L16jWNLkvkAtJTESIH0VqzalUQdepUYevWrxk5cgI//vi7BIqwGumpCJFDXLhwHn//Xnz//Vbq1q3PvHlLqFixktFliVxGQkWIbFa+PFzNhokZvLwgONjUO1m+fBHDhn1IfHw8kybNolevD2QWRmEICRUhsll2BMr97YaEnKVPn3fYtesHGjZsxNy5i3niibLZ84ZCmEFCRQg7Vq9eVRwcHJg5cz5du/bEwUFOkwpjSagIYTiNH3vxZR/5iSCC/OzDl734AeqRa9at24C5cxdRqlRp65QqRDokVIQwiBOxdGMJg5mCN1dxIhYXYonBmTicCcWLKQxmKd2JwznVbWzYsA2lHh08QliT9JWFMIBa+Dst/GcwqPA4ynGGfESRhxgc0OQhhnxEUY4zzGAgP/AC7kSmvh0JFGFjJFSEsDInYtEV7rC1eU2eXrWE9/z9uVy4cKpt3YnGl31soTlOxFq5UiEyTkJFCCvrxhIAYlxcuOvqypLmzSm7alWa4ZKHe9TmIF1Zau1ShcgwCRUhrEozmClJXjEnXNyJJoApgLZirUJknISKEFbkx168Sf3GlfTCxZtQ/NhrrVKFyBQJFSGsyJd96Z4bSStcnIijDvutVKkQmSOXFGeQOh+C0veMLsNu/BUfQ1SuP2JT9cGj/ETgYuYJ9xgXFwAWtGzJ8TJl+LH/APITkS0VCmEpEioZpGLuQaF8RpdhN6Ji7lDQQcagui+C/MSkcc9Jci4xMTgmJNB12zaGBwURiwsR5M/mCoXIGgkVIaxoH75p3sh4X/IwKXrzJgCRuLOfOtYoU4hMy1CoKKXqAc2AekBxIC9wHfgX2A1s1FrftHSRQuQUe/EjFK9Ul6UVJveF4p04dIsQtsusUFFKvQ18CDwF3AaOAieBO0BhoC7QCZinlFoLjNZan8mWioWwa4opDE7ySnphAhCFG5MZTHpjgQlhtHRDRSl1BPACgoDOwGGtdYpTr0qpAsDLQAfguFKqq9b6CwvXK4TdW0p34BezwgTgDq4coDbL6GbdQoXIBHN6KsuAT7XWdx/VSGsdDqwEViqlqgNFLVCfEHbt3r2UVwrG4Yw6mZeXjv/OjKC5lL15Mc31o3DjALVpwZZ0z8UIYQvSDRWt9ayMblRrfQQ4kqmKhMgh/vn7OMOHDMB0yjEp3bMu31ILLwoSwBS8CcWJOJyJIRYX4nAiFG8mM5hldJNAEXYjoyfqa2mtD2VXMULkBPHx8axYuoBPP55FwUKF0mwXhzOL6MUieuLHXuqwP8l8Kr9RDzmHIuxNRi8p3qmUelVrvTNbqhHCzp0/F8KIwIEcOXyQps1aEDh8LI0bpLeWYi/12Ut9a5QoRLbK6DAtq4AtSqnWyRcopZ5RSv1smbKEsC9aa9Z/uZp2rzfn1KkTjJs8k0nT5lKwYNo9FSFyogyFita6NzARWKOUehdAKfW0Uupr4CdAPkEi1wm7fo3+H/Rg3KiPqFqtBms3bKP5y68+mEDLwyM+W97XK/XbXYQwVIbvqNdaj1FKXQTmK6XaAw2A80A3TJcdC5Fr7N65gzEjhhAVGcHAgOG079gFB4ekf6t9/9N5s7d3KzISn8pV028ohI3KcKgopQoDFYF44FngV+B5rXWchWsTwmZFR0Uxfco4Nqxbw5OVqjBu2SrKla9odFlCGC6jV3+NBPonrjcdCAY+BWYAfS1enRA26OiRPxg+ZAAXzofwdrde9O7THxcXV6PLEsImZLSnMhRYjGkYllAApdQ5YINSyhvoqLW22ETaSqlmwGzAEVistZ6UbLkrpkNutYEwoK3W+qyl3l+Ih8XFxbFk4TwWfzqXIl7eLFy+mto+dY0uSwibktFQqay1PvXwC1rrH5VSjYAtwDbgBUsUppRyBOYBTYELwH6l1Gat9V8PNesO3NRal1dKtQMmA20t8f5CPOxcyFmGDenPsaOHeenlVgwZNob8+R8zuiwhbE5Gr/46lcbrh4BngDIWqOk+XyBYa31aax0DrAFaJWvTCliR+Hgd8IK6f8mNEBagtWbj+rW0b92CkLOnmTh1DuMnz5JAESINFptPRWsdrJSy5N1bJTBdVXbfBUyjIafaRmsdp5QKBzwwDcefLVr27oRyTjrp1OvNXqZHhy5E37lD656dUqzT4bU36Ph6W67fuEEn/54plr/TvhOtm7fiwuWL9Bjsn2J5n65bigLgAAAYPElEQVQ9ad74f5w4HYz/yCEplg/u3ZdG9Rty9O9jBEwYlWL5yP4B1KtVh98O7Wf0zMkplk/+aBTVKldl568/MWX+nBTLZ4+eRMWy5dny43fMXbYwxfJFU2ZTslgJvtqyicWrP0uyLCIhnhmz5lOoUGE2b1jH15vWpVh/zvxl5M2bl7WrP+P77d+m3P7yNQAELVvInt0/Jlnm6pqHjxcsN7WbP4d9v/+aZHmBAoWYNns+AHNnTuHokaQDQnh5F2X8ZNNIRFMnjuHEv/91hOPi4rh2NZRLFy/gU6cenkW8WLd2JevWrnzQpuKTVRgUOAKAoQH9uBp6Jcn2q1WvRZ/+plGJP/TvTXh40kEjfevWp0dv0+nID3p1ISo6ivxu7g+WN2v2Mn37fghA8+bPp/jevPbam/To8R7R0dG0adM8xfIOHbrQoUMXwsKu06lTmxTLu3fvTevWbblw4Tw9U/nd7dNnIC+99AonT/6Lv3+vFMsHDRpGo0ZNOHr0MEOG9EuxfOTICdStW5/ff/+V0aM/SrF80qRZVKtWg507dzB16rgUy2fPXkCFCk+ydevXzJ07PcXyhQs/o2TJUnz11RcsWTI/xfLPPluHh4cnK1cuZ+XK5SmWr1u3BTc3NxYt+oQNG9amWL5lyy4A5syZxrZt3yRZlidPXtav3wrA5Mlj2b37hyTLCxf24PPPvwJg1KhA9u3bm2R58eIlWbz4cwACAvrx55+HkywvX74ic+aYPm99+/YkOPhEkuVPP12DyYm/u++805FLly4kWe7r68eoURMB6NixNTduhLFhw64U+2hp6fZUlFKblFI1zdmY1jpUKZVHKTXg/n0sWZBajyP56MjmtEEp1VMpdUApdeDatWtZLEvkBrdvh/PXsaNcvnSRvgMCmL/kc/K6uRldlhA2T6Uyin3SBkrNAXoBhzHdUb8HOPrwJcRKqeKYDle9ArwOXAS6aa33ZbowpfyAUVrrFxOfBwJorSc+1GZ7Ypu9Sikn4ApQJLWh+e/z8fHRBw4cyGxZ3Pj9BC4ynbDZ9tvZdMIxMfeYN3sany1fTJmy5Rg/eRaVq1jvvhG5T0Vkp3v3wNMzc+sqpQ5qrX3Sa2fO4a8YoBHQHhgJFAC0Uuo2cA/TXfTOmHoN+4B+wGda64TMlf7AfqCCUuoJTCHVDngrWZvNwNvAXqAN8OOjAkXkPu3DLlLN2ZV33AtSxPHRv+6ng0/y0WB/Tvz7N2+07Ui/QR+RN29eK1UqRM5gTqj0A9ZqrfsopSIwXeHlBxQD8mC6lPcf4CetdYilCks8R/IBsB3TJcVLtdbHlVJjgANa683AEuAzpVQwcANT8AjxwL/xsZyOj2Xz3Sha5nFPNVy01qz7YiUzpozDzT0fs+YtpuHzFrmIUYhcx5xQucF/Y3oFYJqHPuXZ3mygtd6C6VLlh18b8dDju8Ab1qhF2C/TjVOajXcjU4TLzRthjB4ewE+7fqD+M88xatxUPIsUMbhiIeyXOaHyMzBNKVUE0yEuObwk7FLycPG9dZvjgwYSeS6EQYEjaPvW2ynG7RJCZIw5ofIBpntBVmAKlB1KqaPAHw99HbfknfRCZKf74bLHzRU172OaJDjQpEQZCRQhLCDdT5HW+pLWuimme0IU8AVwGWiGaciWg0CEUuqQUmpJdhYrhEW5uKBdXNiVx4mWYReZePs61+JlXFQhssLsmx+11leUUhuAmVrrvwGUUvmAGkDNxK9a2VKlENnofs9l3d1ITsXHsrhQMYMrEsJ+ZeiOeq1162TPIzGdc5EZH4XdcgYcULTM404P94JGlyOEXbPYMC1C2JvkYeKZzn0sQoj0yZlJkSPFxcXxyZyUY0WBKUxcUbyaJx9fe5Qg8DFPCRQhLEQ+SSLHuXTxAh8N9ufo4UPQ7r8jttIzESL7yadK5CjfbfuGcaM+Ag0Tp84hEAkTIaxJPl0iR7gTHc3USaPZ+NVanq5ekwlTZlOiZCmWJ479JWEihHXIp0zYvRP//MWQQX0JOXOarj3e4933++Hs7AzAao8SBlcnRO4ioSLsltaaL1YFMXPqBAoULMj8xZ/jW8+S88QJITJKQkXYpVu3bjJ62GB279zBMw0bMXr8VAoV9jC6LCFyPQkVYXcOHvidoYP7cSMsjIEBw3mrU1eUSm0SUCGEtUmoCLsRHx/Pok/nsvjTuZQsVZoVq9dbdVZGIUT6JFSEXQi9cpmhAf04dGAfLVq+xpBhY3B3l2mdhbA1EirC5u3euYNRQwcRExPDmAnTeLlV6/RXEkIYQkJF2KyYmHvMnj6Z1Z8v48lKVZg0fS6PlylrdFlCiEeQUBE26VzIGYZ82Jd//jpG+45d8B84BBcXV6PLEkKkQ0JF2Jyt32xi/OihODk5MWPuAp5v/D+jSxJCmElCRdiMO9HRTJk4mk3r11Kjpg8Tps6maLHiRpclhMgACRVhE4JP/suQgX04czqY7j3fp9f7/XBykl9PIeyNfGqFobTWbPzqC6ZMGEW+fPn5ZFEQdf2eMbosIUQmSagIy9Mahz8O4Hj0MI/fvoVLvvxEP12dqBq14KE73yMjI5gwehjbtmymrl8Dxk6ciWeRIgYWLoTIKgkVYTmxsTitW4PL4k9QYdchPo4nYmPRzs5oRyfiPDwI7daLsNff5O/Ew10XL5zjff8P6fpObxwcZCJSIeydhIqwiAZXTuO3Zw8jFn5M8UuXkixTMTFADI4XoikxZTxxyxby/uXLuHh4sHD5amrV9jWmaCGExcmfhiLrYmM5mhDHsjq1Kbd0Ke/5+3O5cOFUmzrevUPx8+fY5ebGmi82SqAIkcNIqIgsc1q3BoAYFxfuurqypHlzyq5alWa45AGqxsRQbucPVq5UCJHdJFRE1miNy+JPkrxkTrg43r2D95JPQWtrViuEyGYSKiJLHP44YDopn4r0wsUpLAz3w4esVaoQwgokVESWOB49DPFxj2yTVriouDjcjh21UqVCCGuQq79E1kRFQmysWU1jXFwAWNCyJcfLlGHXgAHEREZwKyE+Oyu0K+4ygaWwcxIqImvc84Gzs1lNXWJicExIoOu2bQwPCgJnZ0o8VhAvl7zZXKQdiYlHzjIJeyahIrIkvloNcHz0r1HyMCl68yYAOq8b8U9Xt0aZQggrkVARWRJfoza3Eg9rJZdWmNynPYuQUNPHGmUKIaxEQkVk2o1bN3l3SH9Khd9K8np6YQKg8+Yl5p3eScYCE0LYPwkVkSn7Dh+kS//eXLl2lUmBowDzwgRAu7gS/1Q14lq3s2LFQghrkFARGaK15uPlCxkxbQIlvIvx/aqN1K5WgxVXTlN//6+M+GQexZKN/ZVk/bx5iX+qGncXfmb2CX4hhP2wyVBRShUGvgDKAGeBN7XWN5O1qQHMBx4D4oHxWusvrFtp7nIz/Ba9Awfw7Q/bafHCi8yfOINCBQoC8EvRsvBqKZziHUlY9N8oxcTGmsLD0QntWYSYd3qbeigSKELkSDYZKsAQ4Aet9SSl1JDE5wHJ2kQDnbXWJ5VSxYGDSqntWutbyTcmsu7g0cO83e9dLoZeZmLgSN5/uwcq+fkQZ2fi2nUirm1H03wqfx4x3cfino/4ajVIqFFbzqEIkcPZaqi0Ap5PfLwC2EWyUNFan3jo8SWl1FWgCCChYkFaaxauXE7gpNEULeLF9pXr8a1R+9ErKUVCrTok1KpjlRqFELbDVkPFW2t9GUBrfVkp5fWoxkopX8AFOGWN4nKL25ERfDD0QzZs+4ZmjZqwYNIsChcsZHRZQggbZlioKKV2AEVTWTQ0g9spBnwGvK21TkijTU+gJ0Dp0qUzWGnudPTvY3T278XZC+cZPfAj+snMjEIIMxgWKlrrJmktU0qFKqWKJfZSigFX02j3GPAtMExr/dsj3mshsBDAx8dHRsF4BK01y79cxaCxwylcsBBbgtZR30cm0hJCmMdW//TcDLyd+PhtYFPyBkopF2ADEKS1/tKKteVYkVFR9Bzcl77DB9OgTl1+2bhdAkUIkSG2GiqTgKZKqZNA08TnKKV8lFKLE9u8CTQEuiilDid+1TCmXPv3T/AJnn+jBV98vYFhfT9k/aLPKeLhaXRZQgg7Y5Mn6rXWYcALqbx+AHgn8fHnwOdWLi1HWrPpK/xHBuDu5s7mZat53u9Zo0sSQtgpmwwVYR13791l8PgRLPtiJQ3q1GPZ9HkU807t2gkhhDCPhEoudSrkDJ39e3H07+MM7PUBw/oOwslJfh2EEFkj/4vkQpu2f8t7Hw3E0dGRLxesoNnzaV6IJ4QQGWKrJ+pFNoiJiWHIhJF07NuTimXL8/OG7RIoQgiLkp5KLnHh8kU693uX/YcP0btTd8YNHoZLGpNrCSFEZkmo5ALf/7STdwb1ITY2lqBZn/LaS68YXZIQIoeSw185WHx8PONmT6V1z04U8yrK7q+2SKAIIbKV9FRyqGth1+k28H127f2Zjq+3ZfqI8bjlzWt0WUKIHE5CJQf69cA+uvTvzc3wW8wbP53ObWTaXiGEdcjhrxxEa83sJZ/SvHMb3PLm4ce1myVQhBBWJT2VHOLW7XB6Bw7gmx3baPVic+aNn06B/I8ZXZYQIpeRUMkBjvx1jE59e3L+8sW0p/oVQggrkMNfdkxrzYovV/NC25bExMaw9bN1fNClpwSKEMIw0lOxU9F37jBgzEesXL+Wxg0asnjaxxQp7GF0WUKIXE5CxQ6dPHOKTn178dfJfwj8YAAB7/XD0dHR6LKEEEJCxd5s2v4tvQMH4OzszPpFn9Pk2eeNLkkIIR6Qcyp2IjY2liETR9Gxb0+eLF+Bnzdsl0ARQtgc6anYgUuhl+ns/y6//3FABoMUQtg0CRUbt/PXn+g28APu3L3D8pmf0Lp5K6NLEkKINMnhLxuVkJDAlE9m0arbWxTx8GT3uq0SKEIImyc9FRsUdvMGPQf7891PP9L2ldeZPWYy7m5uRpclhBDpklCxMQePHqaTf09Cr11j5sgJdG/fWW5mFELYDTn8ZSO01ixetYL/vfUaSim+W7WBd956WwJFCGFXpKdiA6Kio/EfEcAXX6/nfw0bs3DKbDwKFTa6LCGEyDAJFYOdOB1Mx749+Sf4BMP9B/Hhu31xcJAOpBDCPkmoGGjD1q9576OBuLq6smnpKhrVb2h0SUIIkSXyJ7EBYmNjGTJhJJ37vUuVJyvxy8btEihCiBxBeipW9vDd8e917s7YQXJ3vBAi55BQsaJde/fQdcD73L17V+6OF0LkSHL4ywoSEhKY+ukcWnV7C8/CHuxat0UCRQiRI0lPJZvdDL9FzwB/tu3cwRsvv8qcMVPI5+5udFlCCJEtJFSAokUhNNTc1hXN3m7hgvfI716PS1evMH3EeHrIzYxCiBxOQoWMBErG3LjlilveeLavXE+d6rWy502EEMKGSKhksz3rt+NZWO6OF0LkDhIqZtH4sRdf9pGfCCLIzz582Ysf8OjDWRIoQojcRELlEZyIpRtLGMwUvLmKE7G4EEsMzsThTCheTGEwS+lOHM5GlyuEEIaTUEmDO5Fs4SVqcYh8RCdZlocYIIZ8nGEGA+nAKpqzhSjyGVOsEELYCLlPJRVOxLKFl/Blf4pASc6daHzZxxaa40SslSoUQgjbZJOhopQqrJT6Xil1MvHfQo9o+5hS6qJS6mNLvX83llCLQ+Thnlnt83CP2hykK0stVYIQQtglmwwVYAjwg9a6AvBD4vO0jAV2W+6tNYOZkm4PJTl3oglgCqAtV4oQQtgZWw2VVsCKxMcrgFdTa6SUqg14A99Z6o392Is3VzO1rjeh+LHXUqUIIYTdsdVQ8dZaXwZI/NcreQOllAMwHRhkyTf2ZV+mz404EUcd9luyHCGEsCuGXf2llNoBFE1l0VAzN/EesEVrfT69oU+UUj2BngClS5d+ZNv8ROCSyVBxJob8RCR5LeZmZKa2JXIn7eKKNu9UnhAZ5myFOx8MCxWtdZO0limlQpVSxbTWl5VSxSDV41F+wLNKqfeAfICLUipSa53i/IvWeiGwEMDHx+eRJz0iyE8MzomXDWdMLC5EkD/Ja4Xrmj9WmBBC2DtbvU9lM/A2MCnx303JG2itO9x/rJTqAvikFigZtQ/fxBsZMx4qcTixnzpZLUEIIeyWrZ5TmQQ0VUqdBJomPkcp5aOUWpydb7wXP0JTnsIxSyjeiUO3CCFE7mSToaK1DtNav6C1rpD4743E1w9ord9Jpf1yrfUHlnl3xRQGE4lbhtaKwo3JDCa9scCEECIns8lQMdpSunOIWtzF1az2d3DlALVZRrdsrkwIIWybhEoq4nCmOVvZhy9R6fRYonBjH760YIsMKimEyPUkVNIQRT5e4Af6M4NTlCUSd+7iSjyKu7gSiTunKEt/ZtCEH2QwSSGEwHav/rIJcTiziF4soid+7KUO+5PMp/Ib9ZBzKEII8R8JFbMo9lKfvdQ3uhAhhLBpSuvcNQCiUuoaEJL01dq1s+8dDx7Mvm1bjSdw3egirCi37S/IPucWWdnnx7XWRdJrlOtCJauUUge01j5G12FNuW2fc9v+guxzbmGNfZYT9UIIISxGQkUIIYTFSKhk3EKjCzBAbtvn3La/IPucW2T7Pss5FSGEEBYjPRUhhBAWI6GSCqVUM6XUv0qpYKVUiuH0lVKuSqkvEpf/rpQqY/0qLcuMfR6glPpLKXVUKfWDUupxI+q0pPT2+aF2bZRSWill91cKmbPPSqk3E3/Wx5VSq6xdo6WZ8btdWim1Uyn1R+Lvd3Mj6rQUpdRSpdRVpdSxNJYrpdScxO/HUaVULYsWoLWWr4e+AEfgFFAWcAGOAFWStXkP+DTxcTvgC6PrtsI+NwLcEh/3zg37nNguP/AT8BumOXsMrz2bf84VgD+AQonPvYyu2wr7vBDonfi4CnDW6LqzuM8NgVrAsTSWNwe2YhoOpB7wuyXfX3oqKfkCwVrr01rrGGAN0CpZm1bAisTH64AXVHpzGtu2dPdZa71Tax2d+PQ3oKSVa7Q0c37OAGOBKcBdaxaXTczZ5x7APK31TQCtdWqzrtoTc/ZZA48lPi4AXLJifRantf4JuPGIJq2AIG3yG1AwcYZdi5BQSakEcP6h5xcSX0u1jdY6DggHPKxSXfYwZ58f1h3TXzr2LN19VkrVBEpprb+xZmHZyJyfc0WgolLqF6XUb0qpZlarLnuYs8+jgI5KqQvAFqCPdUozTEY/7xkiY3+llFqPI/klcua0sSdm749SqiPgAzyXrRVlv0fus1LKAZgJdLFWQVZgzs/ZCdMhsOcx9Ub3KKWqaq1vZXNt2cWcfW4PLNdaT1dK+QGfJe5zQvaXZ4hs/f9LeiopXQBKPfS8JCm7ww/aKKWcMHWZH9XdtHXm7DNKqSbAUKCl1vqelWrLLuntc36gKrBLKXUW07HnzXZ+st7c3+1NWutYrfUZ4F9MIWOvzNnn7sBaAK31XiAPpjGyciqzPu+ZJaGS0n6gglLqCaWUC6YT8ZuTtdkMvJ34uA3wo048A2an0t3nxENBCzAFir0fZ4d09llrHa619tRal9Fal8F0Hqml1vqAMeVahDm/2xsxXZSBUsoT0+Gw01at0rLM2edzwAsASqnKmELlmlWrtK7NQOfEq8DqAeFa68uW2rgc/kpGax2nlPoA2I7pypGlWuvjSqkxwAGt9WZgCaYucjCmHko74yrOOjP3eSqQD/gy8ZqEc1rrloYVnUVm7nOOYuY+bwf+p5T6C4gHBmmtw4yrOmvM3OeBwCKlVH9Mh4G62PMfiUqp1ZgOX3omnicaCaZpabXWn2I6b9QcCAaiga4WfX87/t4JIYSwMXL4SwghhMVIqAghhLAYCRUhhBAWI6EihBDCYiRUhBBCWIyEihBCCIuRUBFCCGExEipCCCEsRkJFCAMppcorpWKVUqOTvT5fKRVh52ONiVxIQkUIA2mtg4HFQP/EsbZQSo0AugGv2flYYyIXkmFahDCYUqooptkJPwH+wTQTYXut9VpDCxMiE2RASSEMprW+opSahWlgQyegrwSKsFdy+EsI23AScAX2aq3nGV2MEJkloSKEwZRSjTHNVbMXaKCUqm5wSUJkmoSKEAZSStXCNDHWYkxzYJwDJhhZkxBZIaEihEGUUuWBrcB3QB+tdQwwGmiulGpoaHFCZJJc/SWEARKv+PoVU8/kRa31vcTXHYFjwE2tdX0DSxQiUyRUhBBCWIwc/hJCCGExEipCCCEsRkJFCCGExUioCCGEsBgJFSGEEBYjoSKEEMJiJFSEEEJYjISKEEIIi5FQEUIIYTH/ByA72Ioh3QslAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import nonlinear_plots\n", "nonlinear_plots.bisection_root()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Apply the bisection algorithm *by hand* to solve $f(x) = x^3 + 2x^2 - 5x - 6 = 0$ for $x\\in [1, 4]$. Find the root to within $\\tau = 10^{-1}$. Write out all values of $L$, $R$, $C$, and $f(C)$ for each step.\n", "\n", "***\n", "\n", "**Exercise**: Implement the bisection algorithm in Python and use it to solve $\\sin(x) - 2/5$ for $x \\in [0, 1]$.\n", "\n", "***\n", "\n", "**Exercise**: Modify the bisection algorithm to accept a desired number of iterations $n$ instead of the tolerance $\\tau$. Implement this modified algorithm in Python. Then, show the error in the approximate root of $f(x) = \\sin(x) - 2/5$ for $x\\in [0, 1]$ as a function of $n$.\n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Newton's Method\n", "\n", "A more power way to find roots of $f(x) = 0$ is [Newton's method](http://mathworld.wolfram.com/NewtonsMethod.html), sometimes called the Newton-Raphson method. Like bisection, Newton's method produces a sequence of approximations for a root. The values of the sequence are increasingly close to the root. Unlike bisection, Newton's method requires not a range in which a single root lives but an initial guess for what the root is. \n", "\n", "Let the root of interest be $x_r$. Assume we have a good way to guess the root (intuition, a graph, etc.), and call that initial approximation $x_0$. Then, for some $\\Delta$, $x_r = x_0 + \\Delta$,Taylor series provides the following relationship:\n", "\n", "$$\n", " f(x_0) + \\Delta \\frac{df}{dx}\\Bigg |_{x=x_0} + \\mathcal{O}(\\Delta^2) = 0 \\, .\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Formally, this equation is still exact because $\\mathcal{O}(\\Delta^2)$ captures *all* terms proportional to $\\Delta^2$ and higher powers of $\\Delta$. However, if $\\Delta$ were small enough that the $\\Delta^2$ terms can be ignored (i.e., our guess $x_0$ were close enough to $x_r$), then we're left with the *approximate* relationship\n", "\n", "$$\n", " f(x_r) = 0 \\approx f(x_0) + \\Delta f'(x_0) \\, ,\n", "$$\n", "\n", "from which it follows that\n", "\n", "$$\n", " \\Delta = -\\frac{f(x_0)}{f'(x_0)} \\, .\n", "$$\n", "\n", "We expect $x_1 = x_0 + \\Delta = x_0 - f(x_0)/f'(x_0)$ to be a *better* approximation than $x_0$. Moreover, the same process can be applied to determine a sequence of estimated roots $x_2$, $x_3$, and so on. This process is Newton's method, and is presented in the following pseudocode:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "'''Newton's method for finding the root f given x_0'''\n", "Input: f, fp, x_0, tau\n", "Set x = x_0\n", "While |f(x)| > tau do\n", " # Compute the Newton \"step\" \n", " Set Delta = -f(x)/fp(x)\n", " # Compute the next value of x\n", " Set x = x + Delta\n", "Output: x\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Consider $f(x) = \\cos(x) = 0$. With $x_0 = 1$, apply *one* step of Newton's method. What is the result (i.e., $x_1$)?\n", "\n", "*Solution*: $x_1 = x_0 - \\frac{\\cos(x_0)}{-\\sin(x_0)} \\approx 1.642$\n", "\n", "\n", "***\n", "\n", "**Exercise**: Apply Newton's method *manually* to solve $f(x) = x^3 + 2x^2 - 5x - 6 = 0$ to within $\\tau = 10^{-5}$. Use $x_0 = 3/2$. Write out $n$, $x_n$, $f(x_n)$, $f'(x_n)$ at each step $n = 0, 1, \\ldots$.\n", "\n", "*Solution*:\n", "\n", "```\n", "n x f(x) f'(x)\n", "0 1.5000000000000000e+00 -5.625e+00 7.750e+00\n", "1 2.2258064516129030e+00 3.807e+00 1.877e+01\n", "2 2.0229637924064581e+00 3.487e-01 1.537e+01\n", "3 2.0002760690096917e+00 4.142e-03 1.500e+01\n", "4 2.0000000406383567e+00 6.096e-07 1.500e+01\n", "5 2.0000000000000009e+00 1.421e-14 1.500e+01\n", "```\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The exercise above illustrates an important feature of Newton's method. Once Newton's method \"get's going,\" the number of correct digits usually *doubles* at each step. For example, step 3 got 3 digits right, while step 4 got 7, and step 5 got 15. That's pretty quick convergence. \n", "\n", "Such convergence is called **quadratic convergence** or **second-order convergence**. In other words, if the error is $\\epsilon_n$ at step $n$, the error at step $n+1$ will be $\\epsilon_{n+1} = a \\epsilon_n^2$ for some constant $a$. More generally, the **order of convergence** of a method is the value of $p$ for which $\\epsilon_{n+1} = c\\epsilon_{n}^p$ for some constant $c$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Use the last three values of $x$ in the last exercise to estimate the order of convergence of Newton's method.\n", "\n", "*Solution*: Using logarithms, we have $\\log \\epsilon_4 = p \\log \\epsilon_3 + \\log c$ and $\\log \\epsilon_5 = p \\log \\epsilon_4 + \\log c$. Those are two equations for the two unknowns $p$ and $\\log c$. We can set up and solve the $2\\times 2$ system via:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9990273906811293\n" ] } ], "source": [ "e_3 = 2.0002760690096917e+00-2\n", "e_4 = 2.0000000406383567e+00-2\n", "e_5 = 2.0000000000000009e+00-2\n", "import numpy as np\n", "A = np.array([[np.log(e_3), 1], [np.log(e_4), 1]])\n", "b = np.array([np.log(e_4), np.log(e_5)])\n", "p = np.linalg.solve(A, b)[0]\n", "print(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "That's pretty darn close to 2, and Newton's method is therefore confirmed to exhibit second-order convergence (for this problem).\n", "\n", "***\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Knowledge of a method's order of convergence is useful when verifying its implementation. If we saw less than second-order convergence for the exercise above, then the implementation might be suspect.\n", "\n", "### Secant Method\n", "\n", "If we don't know the derivative of $f(x)$, then Newton's method can't be used. Of course, the bisection method is one option, but we can adapt Newton by leveraging the finite-difference methods from [Lecture 21](ME400_Lecture_21.ipynb). The result is the **secant method**. The algorithm is nearly identical to Newton's method, but $f'(x)$ is approximated when needed. The classical secant method begins with two values of $x$ and the corresponding function values. The finit\n", "\n", "```\n", "'''Secant method for finding the root f given x_0'''\n", "Input: f, x_0, tau\n", "Set x = x_0\n", "Set x_0 = x - tau\n", "While |f(x)| > tau do\n", " # Approximate the derivative\n", " Set fp = (f(x)-f(x_0))/(x-x_0)\n", " # Compute the Newton \"step\"\n", " Set Delta = -f(x)/fp\n", " # Store old x, and compute next value\n", " Set x_0 = x\n", " Set x = x + Delta\n", "Output: x\n", "```\n", "\n", "***\n", "\n", "**Exercise**: Apply the secant method *manually* to solve $f(x) = x^3 + 2x^2 - 5x - 6 = 0$ to within $\\tau = 10^{-5}$. Use $x_0 = 3/2$. Write out $n$, $x_n$, $f(x_n)$, $f'(x_n)$ at each step $n = 0, 1, \\ldots$. How does it compare to Newton's method?\n", "\n", "***\n", "\n", "**Exercise**: Implement the secant method as a Python function.\n", "\n", "***\n", "\n", "**Exercise**: Implement a function `newton(f, x_0, fp=None, tau=1e-8, delta=1e-8)` that applies Newton's method when `fp` is provided and the secant method if `fp` is not provided. Here, `delta` is to be the step used in the finite-difference approximation to $f'(x)$.\n", "\n", "***\n", "\n", "**Exercise**: Use the secant method *without* defining $f'(x)$ to solve the following equations:\n", "\n", " 1. $\\tan(x)/x^2$ for $x_0 = 1$.\n", " 2. $s^2 = 3$ for $s_0=1.5$. (What is this doing?)\n", " 3. $2\\cosh(x/4)-x = 0$ with $x_0 = 4$.\n", "\n", "***\n", "\n", "**Exercise**: Another method that can be used when the derivative is not available is [Steffensen's method](https://en.wikipedia.org/wiki/Steffensen%27s_method). Steffensen's method produces a sequence of approximate roots according to $x_{n+1} = \\frac{f(x_n + f(x_n)) - f(x_n)}{f(x_n)}$. Try this method on nonlinear equations from the last exercise. Does it converge quadratically like Newton's method?\n", "\n", "***\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving Systems of Nonlinear Equations Numerically\n", "\n", "The solution of nonlinear systems is considerably more challenging than single-variable equations. However, the Newton (and secant) methods described so far are up to the challenge. A system of nonlinear equations can be written generically as $\\mathbf{f}(\\mathbf{x}) = \\mathbf{0}$, where the bolded names indicate vector quantities. Whereas Newton's method for single-variables equations requires the derivative $f'(x)$, Newton's method for systems of equations requires the **jacobian matrix**. For a system of $n$ unknows, the jacobian $\\mathbf{J}$ is a square matrix defined as\n", "\n", "$$\n", " \\mathbf{J}(\\mathbf{x}) = \\begin{bmatrix}\n", " \\frac{\\partial f_0}{\\partial x_0} & \\frac{\\partial f_0}{\\partial x_1} & \\ldots & \\frac{\\partial f_{0}}{\\partial x_{n-1}} \\\\\n", " \\frac{\\partial f_1}{\\partial x_0} & \\frac{\\partial f_1}{\\partial x_1} & \\ldots & \\frac{\\partial f_{1}}{\\partial x_{n-1}} \\\\\n", " & & \\ddots & \\\\\n", " \\frac{\\partial f_{n-1}}{\\partial x_{0}} & \\frac{\\partial f_{n-1}}{\\partial x_{1}} & \\ldots & \\frac{\\partial f_{n-1}}{\\partial x_{n-1}} \n", "\\end{bmatrix}\n", "$$\n", "\n", "where $x_i$ indicates the $i$th unknown (or element $i$ of $\\mathbf{x}$).\n", "\n", "Then, Newton's method leads to the sequence\n", "\n", "$$ \n", " \\mathbf{x}_{n+1} = \\mathbf{x}_{n} - \\mathbf{J}^{-1}(\\mathbf{x}_n) \\mathbf{f}(\\mathbf{x}_n) \\, .\n", "$$\n", "\n", "Here, $\\mathbf{f}(\\mathbf{x}_n)$ and $\\mathbf{J}(\\mathbf{x}_n)$ indicate the function and its jacobian are evalaluated at $\\mathbf{x}_n$. \n", "\n", "Note that if we let $\\mathbf{\\Delta}_n = \\mathbf{J}^{-1}(\\mathbf{x}_n)\\mathbf{f}(\\mathbf{x}_n)$, then it becomes clear that $\\mathbf{J}(\\mathbf{x}_n)\\mathbf{\\Delta}_n = \\mathbf{f}(\\mathbf{x}_n)$. In other words, each Newton step for a system of *nonlinear* equations requires the solution of a *linear* system of equations where the Newton step $\\mathbf{\\Delta}_n$ is the unknown.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Consider the nonlinear system of equations $x-y = 4$ and $x^2 + y = -3$. Derive the jacobian matrix, and apply one step of Newton's method using $x_0=1/2$ and $y_0=-3$. \n", "\n", "*Partial Solution*: Here, we have \n", "\n", "$$\n", "\\mathbf{f}(\\mathbf{x}) = \\left[ \\begin{array}{c}\n", " x-y - 4\\\\\n", " x^2 + y + 3\n", "\\end{array}\\right ] = \n", "\\left[ \\begin{array}{c}\n", " 0\\\\\n", " 0\n", "\\end{array}\\right ] \\, ,\n", "$$\n", "\n", "where $\\mathbf{x} = [x, y]^T$.\n", "\n", "The jacobian is, therefore, defined as\n", "\n", "$$\n", "\\mathbf{J}(\\mathbf{x}) =\n", "\\left[ \\begin{array}{cc}\n", " 1 & -1 \\\\\n", " 2x & 1 \\\\\n", "\\end{array}\\right ] \\, .\n", "$$\n", "\n", "With $\\mathbf{x}_0 = [1/2, -3]^T$, Newton gives $\\mathbf{x}_1 = \\mathbf{x}_0 - \\mathbf{J}^{-1}(\\mathbf{x})\\mathbf{f}(\\mathbf{x}_0)$. \n", "\n", "\n", "***\n", "\n", "**Exercise**: What's the Jacobian matrix for the *linear* system $ax+by = 1$ and $cx+dy=2$? How many Newton iterations are required to solve the system for *any* initial values $x_0$ and $y_0$? \n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### `scipy.optimize.fsolve`\n", "\n", "Although implementation of a Newton (or secant) method for systems would be possible using NumPy, it is more practical to use the tools included in SciPy. In particular, the `scipy.optimize.fsolve` function implements a method that combines Newton's method with some other techniques applied when Newton's method cannot reliably converge (e.g., what happens when $f'(x)$ vanishes near the root of interest?). \n", "\n", "The signature (from `help(fsolve)`) is\n", "```python\n", "fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)\n", "```\n", "The arguments of interest here are `func`, `x0`, `args`, and `fprime`; refer to the documentation for more on the the other arguments. From `help(fsolve)`:\n", "```\n", "func : callable ``f(x, *args)``\n", " A function that takes at least one (possibly vector) argument.\n", "x0 : ndarray\n", " The starting estimate for the roots of ``func(x) = 0``.\n", "args : tuple, optional\n", " Any extra arguments to `func`.\n", "fprime : callable(x), optional\n", " A function to compute the Jacobian of `func` with derivatives\n", "```\n", "Hence, `func` represents our nonlinear system function $\\mathbf{f}(\\mathbf{x})$, and `args` represents any values needed to define what $\\mathbf{f}(\\mathbf{x})$ does. The initial guess is `x0`. If available, `fprime` should compute $\\mathbf{J}(\\mathbf{x})$; if not provided, `fsolve` approximates the Jacobian using a forward-difference approximation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "**Exercise**: Solve the nonlinear system of equations $x-y = 4$ and $x^2 + y = -3$ using `fsolve` $x_0=1/2$ and $y_0=-3$.\n", "\n", "*Solution*:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.6180339887498933 -3.381966011250107\n" ] } ], "source": [ "from scipy.optimize import fsolve\n", "\n", "def f(z):\n", " x, y = z # unpack z\n", " return [x-y-4, x**2+y+3]\n", "\n", "z0 = [1/2, -3]\n", "x, y = fsolve(f, z0)\n", "print(x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Further Reading\n", "\n", "Checkout out the [SciPy documentation](https://docs.scipy.org/doc/scipy/reference/optimize.html) on optimization and root finding." ] } ], "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }