{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Solving equations with SageMath" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "%display latex" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Exact solutions : `solve`" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Let us start with a simple example: solving $x^2-x-1=0$ for $x$:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]$ " ] }, "execution_count": 2, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(x^2-x-1 == 0, x)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Note that the equation is formed with the operator `==`. \n", "\n", "The equation itself can be stored in a Python variable:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "eq = x^2-x-1 == 0" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}x^{2} - x - 1 = 0$ " ] }, "execution_count": 4, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eq" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}x^{2} - x - 1$ " ] }, "execution_count": 5, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eq.lhs()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}0$ " ] }, "execution_count": 6, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eq.rhs()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]$ " ] }, "execution_count": 7, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(eq, x)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "The solutions are returned as a list:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]$ " ] }, "execution_count": 8, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol = solve(eq, x)\n", "sol" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}$ " ] }, "execution_count": 9, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol[1]" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Each element of the solution list is itself an equation, albeit a trivial one. This can be seen by using the `print` function instead of the default LaTeX display:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x == 1/2*sqrt(5) + 1/2\n" ] } ], "source": [ "print(sol[1])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "The access to the value of the solution is via the `rhs()` method:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}$ " ] }, "execution_count": 11, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol[1].rhs()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "A numerical approximation, `n` being a shortcut for `numerical_approx`:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}1.61803398874989$ " ] }, "execution_count": 12, "metadata": { }, "output_type": "execute_result" } ], "source": [ "n(sol[1].rhs())" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}1.6180339887498948482045868343656381177203091798058$ " ] }, "execution_count": 13, "metadata": { }, "output_type": "execute_result" } ], "source": [ "n(sol[1].rhs(), digits=50)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "A new equation involving the above solution:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2} = \\phi$ " ] }, "execution_count": 14, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol[1].rhs() == golden_ratio()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Asking Sage whether this equation always holds:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}$ " ] }, "execution_count": 15, "metadata": { }, "output_type": "execute_result" } ], "source": [ "bool(_)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "### A system with various unknowns\n", "\n", "Since `x` is the only predefined symbolic variable in Sage, we need to declare the other symbolic variables denoting the unknowns, here `y`:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "y = var('y')" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Then we may form the system:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(x^{2} + x y + 2 = 0, y^{2} = {\\left(x + y\\right)} x\\right)$ " ] }, "execution_count": 17, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eq1 = x^2 + x*y + 2 == 0\n", "eq2 = y^2 == x*(x+y)\n", "(eq1, eq2)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left[x = -\\frac{1}{2} i \\, \\sqrt{10} + \\frac{1}{2} i \\, \\sqrt{2}, y = -i \\, \\sqrt{2}\\right], \\left[x = \\frac{1}{2} i \\, \\sqrt{10} + \\frac{1}{2} i \\, \\sqrt{2}, y = -i \\, \\sqrt{2}\\right], \\left[x = -\\frac{1}{2} i \\, \\sqrt{10} - \\frac{1}{2} i \\, \\sqrt{2}, y = i \\, \\sqrt{2}\\right], \\left[x = \\frac{1}{2} i \\, \\sqrt{10} - \\frac{1}{2} i \\, \\sqrt{2}, y = i \\, \\sqrt{2}\\right]\\right]$ " ] }, "execution_count": 18, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol = solve((eq1, eq2), x, y)\n", "sol" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Here again the solutions are returned as a list; but now, each item of the list is itself a list, containing the values for `x` and `y`:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{1}{2} i \\, \\sqrt{10} + \\frac{1}{2} i \\, \\sqrt{2}, y = -i \\, \\sqrt{2}\\right]$ " ] }, "execution_count": 19, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol[0]" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}y = -i \\, \\sqrt{2}$ " ] }, "execution_count": 20, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol[0][1]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-i \\, \\sqrt{2}$ " ] }, "execution_count": 21, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol[0][1].rhs()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-1.41421356237310i$ " ] }, "execution_count": 22, "metadata": { }, "output_type": "execute_result" } ], "source": [ "n(sol[0][1].rhs())" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Approximate solution: `find_root`" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Let us consider the following transcendental equation:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}x \\sin\\left(x\\right) = 4$ " ] }, "execution_count": 23, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eq = x*sin(x) == 4\n", "eq" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "`solve` returns a non-trivial equation equivalent to its input, meaning it could not find a solution:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = \\frac{4}{\\sin\\left(x\\right)}\\right]$ " ] }, "execution_count": 24, "metadata": { }, "output_type": "execute_result" } ], "source": [ "solve(eq, x)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "We use then `find_root` to get an approximate solution in a given interval, here $[0, 10]$:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}8.96212620082$ " ] }, "execution_count": 25, "metadata": { }, "output_type": "execute_result" } ], "source": [ "find_root(eq, 0, 10)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Actually `find_root` always return a single solution, even if there are more than one in the prescribed interval. In the current case, there is a second solution, as we can see graphically:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGFCAYAAAAPa6wiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmcDWT7x/HPGUYYjKXsS2m1pSwztoSyxtieaqZSCSlJVEob7aqn8vza5KFQMbRoxpIKCY/sW1G2UJGtMGM3zPn9cb1UCrOdc+6zfN+v17zO05jnnC/DzDX3fd3X7fF6vV5ERERE5IyiXAcQERERCXYqmERERESyoIJJREREJAsqmERERESyoIJJREREJAsqmERERESyoIJJREREJAsqmERERESyoIJJREREJAsqmERERESy4POCaejQoURFRXH//fef9te9Xi/p6enoRhYREREJFT4tmJYsWcLIkSOpXbv2GT9m//79xMbGsn//fl++tIiIiIjf+KxgOnDgALfccgujRo2iePHivnpaEREREed8VjDdc889dOjQgRYtWvjqKUVERESCQn5fPMmECRNYuXIlS5cu9cXTiYiIiASVPBdMW7dupX///syYMYPo6GhfZBIREREJKh5vHo+rpaam0qVLF/Lly/fHybcTJ07g8XjIly8fR48exePx/PHx6enpxMbG0rZtW/LnP7VeS0pKIikpKS9xRERERHwuzwXTwYMH+emnn0553+233061atUYNGgQ1apVO+XXThZMaWlpFCtWLC8vLSIi4lOZmbBiBcyeDT/8AFu2wJEjUKAAVK4M1apBy5Zw5ZUQpUmGESXPW3IxMTFUr179H+8rVarUP4olERGRYLR7N/z3v/DOO7B5MxQuDDVrQpUqEBMDR4/Chg0waRI88giULw933QW9e0Pp0q7TSyD4pOn77/66BSciIhKsDh6El1+2t8xMuPFGePddaNTIVpX+7tgxWLAAxo+HoUPtbdAgGDgQChUKfH4JnDxvyeWUtuRERCQYzJ0L3bvDtm3Qt6+tHJUqlf3//5498MIL8J//2HbdBx9Agwb+yytuaQdWREQiyokT8Pjj0KwZlCsHq1fbClNOiiWAkiXhpZfs/3/uudCkCQwbBrr5KzypYBIRkYhx8CAkJNhW2rPPwpw5cNFFeXvOSy6BefNgwAC4/364+27IyPBNXgkefulhEhERCTZ79sB119mK0GefQevWvnvu6Gj497/hssusGfy33yA52d4v4UErTCIiEvZ27YKrr4aNG21kgC+Lpb/q0cNO0k2ebA3kWmkKHyqYREQkrKWnQ9u2tuozbx7Uq+ff1+vQwYqmqVPhzjvV0xQuVDCJiEjYOnIEOnWCH3+EL76wLbNAaN8eRo+GMWPgmWcC85riX+phEhGRsJSZCd262dykL7+Eyy8P7OvffLNNCn/8cWsMT0wM7OuLb6lgEhGRsPT88/Dxx7Y9dtVVbjI8+qhdsdKzJ9SubVerSGhyNrjy5OW7unBXRER87fPPoV07GDwYnnzSbZYDByA+3v734sV21YqEHk36FhGRsLJ5M9StCw0bwpQpwXFJ7g8/QP36kJQEI0e6TiO5EQR/jURERHzj+HErSkqUsKtKgqFYAtuKe/VVGDUKpk93nUZyI0j+KomIiOTdCy/AkiUwbpwVTcGkVy+b/9SzJ+zd6zqN5JQKJhERCQvLlsFTT9klusF4Ca7HA++8A4cOQf/+rtNITqmHSUREQt7hw1CnDhQqBAsXQoECrhOd2bvv2kTwOXOgaVPXaSS7tMIkIiIhb/Bga/b+4IPgLpYAbr/dGtLvuUdXp4QSFUwiIhLSVq6EYcNgyBCoXt11mqxFRcGbb8L338Prr7tOI9mlLTkREQlZJ05Ao0Zw8CAsXx78q0t/de+9dnXK2rVQoYLrNJIVrTCJiEjIGjHChkGOGBFaxRLYHXOFCsETT7hOItmhgklERELSb7/Z1SM9ekDjxq7T5Fzx4lYsjR0Lq1e7TiNZUcEkIiIhacgQ8Hph6FDXSXKvd284/3wr/CS4qWASEZGQs3o1vP22nY477zzXaXKvQAF47jm7wmXePNdp5Gx0+a6IiIQUr9cmZm/eDGvWhF7v0t9lZkJcHJxzDvzvfzbgUoKPTsmJiEhImTYN2reHlBTo2NF1Gt+YPh3atYMZM+Daa12nkdNRwSQiIiEjMxOuuAJKloTZs8NnNcbrtVWmwoVtArgEH/UwiYhIyJg4Eb77Dp5/PnyKJbDfy+OPw9y5KpiClVaYREQkJBw/bpO8L7kEpk51ncb3vF648ko491yYOdN1Gvk7rTCJiEhIGDsWNmyAZ591ncQ/Tq4yzZoF33zjOo38nVaYREQk6B09aitL8fHw4Yeu0/hPZibUrAkXXQSTJ7tOI3+lFSYREQl6//0vbN0KTz/tOol/RUXBAw/YluOGDa7TyF+pYBIRkaB25Ig1eXfrBpdd5jqN/918s/Ux/d//uU4if6WCSUREgtqYMbBrFzz2mOskgVGwIPTpA6NHw549rtPISSqYREQkaB0/Di+9BP/6F1x8ses0gXP33XDiBIwc6TqJnKSCSUREgtbEiXYFyqBBrpMEVpkytjX3+uuQkeE6jYAKJhERCVKZmTB0KLRta/OJIs2AAbBtW3ifCgwlzgqmxMREEhISSE5OdhVBRESC2NSpdrnuI4+4TuJGzZp2r9ybb7pOIqA5TCIiEoS8XmjYEKKjYd4812ncmTQJunaFlSuhdm3XaSKbtuRERCTofP01LFoUuatLJ3XoAOXKwYgRrpOICiYREQk6L78Ml19u/UuRLDoaevaEDz6AAwdcp4lsKphERCSorFsHn31mTc8ej+s07vXsCQcPglp+3VLBJCIiQeW116B0aUhMdJ0kOFSuDO3awfDh1tslbqhgEhGRoLF3r032vvtum3gt5q67YMUKWLrUdZLIpYJJRESCxqhRNt37rrtcJwkubdrYSpOav93Jc8E0dOhQ4uLiKFasGGXKlKFz586sX7/eF9lERCSCHD9uk62TkqBsWddpgku+fNC9uw2xPHjQdZrIlOeCad68edx7770sWrSImTNnkpGRQatWrTh8+LAv8omISIT49FP45Re47z7XSYLTrbfC/v325ySB5/PBlb/99hulS5dm7ty5NGnS5B+/rsGVIiJyOo0b2zH6r792nSR4NWsG+fPDzJmuk0Qen/cw7du3D4/HQ8mSJX391CIiEqaWL4dvvoF+/VwnCW633w5ffQU//+w6SeTxacHk9Xrp378/TZo0oXr16r58ahERCWMjRkD58pCQ4DpJcPvXv6BwYXj/fddJIo9PC6Y+ffrw/fffM2HCBF8+rYiIhLH0dBg3Dnr1su0mObMiRaxoGjNGM5kCzWd/Nfv27ctnn33GvHnzKFeuXJYfn5iYSP6//ctISkoiKSnJV5FERCQEjBsHR45YwSRZu/12GDsW5s+H07QKi5/4pOm7b9++pKamMmfOHKpWrXrWj1XTt4iInOT1Qu3acOGFOv2VXZmZcMEFf07/lsDI85Zcnz59GDduHOPHjycmJoadO3eyc+dOjhw54ot8IiISxhYuhO++06DKnIiKsmtjPvoIMjJcp4kceV5hioqKwnOa2xFHjx7Nrbfe+o/3a4VJREROuu02mDcPNm60QkCyZ+VKuPJKu6S4bVvXaSJDnnuYMjMzfZFDREQizJ49MHEiPPWUiqWcql0bLrsMkpNVMAWK/oqKiIgTY8daP0737q6ThB6Px66Q+fRT0MUagaGCSUREAs7rhbffhq5doXRp12lCU1ISHDgA06a5ThIZVDCJiEjAzZkD69dD796uk4Suiy+GunVtW078TwWTiIgE3OjRcNFFcPXVrpOEtqQkW2FKS3OdJPypYBIRkYDavx8+/tgGMJ7mkLXkwI03wrFjmmEVCCqYREQkoD780BqVTzN5RnKoYkW46io7bSj+pYJJREQCavRoaNkSKlVynSQ8XH89zJoF+/a5ThLeVDCJiEjArF9vd6BplIDvdO5sE7+nTHGdJLw5K5gSExNJSEggWe39IiIRY8wYKF4cOnVynSR8VKgADRrAJ5+4ThLe8jzpO7cmTJigq1FERCLIiRPw3nt2sqtgQddpwkvXrvDEEzaXqUgR12nCk7bkREQkIGbMgG3btB3nD126wJEjMH266yThSwWTiIgExOjRUKMG1KvnOkn4qVoVrrhC23L+pIJJRET8bs8eSEmx1SXNXvKPrl1tiOWRI66ThCcVTCIi4nfJydbDdMstrpOEr65drYdpxgzXScKTCiYREfG70aPhuuugTBnXScJXtWr2pm05/1DBJCIifrV2LSxbpsnegdClC0yebHOZxLdUMImIiF+NGwexsbbCJP7VtSvs3Qtz5rhOEn5UMImIiN94vfDBB3Z9h2Yv+d8VV9iVM5r67XsqmERExG+++Qa2bIGbb3adJDJ4PNC+vRVMXq/rNOFFBZOIiPjNuHFQsSI0beo6SeRISIDNm2HNGtdJwosKJhER8Ytjx2DiRFtditJ3m4Bp1gxiYrQt52u6fFdERPzi889tYKW24wKrYEFo1cpOy4nveLzewO5ypqenExsbS1pami7fFREJYzfeaCMFVq1ynSTyjBkDd9wBO3ZA6dKu04QHLZKKiIjPpaXZCocme7vRrp09TpvmNkc4UcEkIiI+N2kSHD0KSUmuk0Sm0qWhQQNty/mSCiYREfG5ceOs+bhiRddJIldCAnz5pS7j9RUVTCIi4lM7d8Ls2Vpdcq1DBzh0yD4XkncqmERExKc+/tjGCHTp4jpJZKteHS64QOMFfEUFk4iI+NTEiXDttVCqlOskkc3jsebv6dM19dsXVDCJiIjPbNsG//ufjRQQ99q2tatp1q1znST0qWASERGf+egjiI6GTp1cJxGA5s3hnHNslUnyRgWTiIj4zMSJ0Lo1FC/uOokAFC4MV1+tgskXVDCJiIhP/PQTLFyo7bhg064dzJkDBw+6ThLaVDCJiIhPfPih3WOWkOA6ifxV27Z2EbLGC+SNLt8VERGfmDjRVjOKFnWdRP7q4ouhalVty+VVflcvPGHCBF2+KyISJn78EZYtg4cecp1E/s7jsVWmzz6z8QIej+tEoUlbciIikmcffmgNxtdd5zqJnE7btrB5M6xf7zpJ6FLBJCIieTZxol3FERPjOomcjsYL5J0KJhERyZN162DVKp2OC2Ynxwt89pnrJKFLBZOIiOTJxInW6N22reskcjZt2sDcuXD4sOskoUkFk4iI5MnEidCxo40UkODVsiUcPWpX10jOqWCSLGn0Q+jQ5yq0hMPna/Vq+P77yNiOC/XPV40aUK4czJjhOon/+eNzpYJJshTqXyQiiT5XoSUcPl8TJ9o1KK1auU7if6H++fJ4bJVJBVPuqGASEZFc8XptnEDnzlCggOs0kh0tW8LKlbBrl+skoSfsC6ZA/kQQrq8VSOH6ZxiOn69w/fMLx88V+Of39d13Ntfnhhv8/1pnos9Xzlx7rT3OnOnf1zmdUP9cqWDSawWVcP0zDMfPV7j++YXj5wr88/v65BOIjYUWLfz/Wmeiz1fOlC0LtWr9uS2nz1X2+exqFK/Xy/79+7P8uPT09FMe/e348eN6Lb1WxLxWOP6e9FrB+1offgitW8ORI/bmz9c6E71WzjVtCpMmQVpa+Pye8vpaRYsWxZPFnTEer9frzWswsAIoNjbWF08lIiIiEjBpaWlZ3m/rs4IpJytMlSpV4pdfftHluyIiIeqVV+Dll2HTJihUyHUayYlDh6ByZXj2WbjrLtdpgkNAV5iy6+RKVHaqORERCU716kHVqrYtJ6GnRQu792/KFNdJQkfYN32LiIhvbdkCy5ZB166uk0hutWwJX38Nx465ThI6VDCJiEiOTJpkN9+3a+c6ieRWy5Zw4AAsXOg6SehQwSQiIjkyaZJN9i5a1HUSya0rr4SSJSNj6revqGASEZFs274dvvlG23GhLl8+uOYaFUw5oYJJTmvo0KHExcVRrFgxypQpQ+fOnVm/fr3rWJINQ4cOJSoqivvvv991FDmDX3/9lW7dunHuuedSuHBhateuzfLly13HypZPP7VvtgkJrpP4X2ZmJk888QRVq1alcOHCXHTRRTz77LOuY/nMNdfA0qWQjQPuQWnevHkkJCRQoUIFoqKimDx58j8+ZvDgwZQvX57ChQvTsmVLNm7cmOvXU8EkpzVv3jzuvfdeFi1axMyZM8nIyKBVq1YcPnzYdTQ5iyVLljBy5Ehq167tOoqcwb59+2jcuDHnnHMOX3zxBT/88AOvvPIKJUqUcB0tWz75xE5YhUjcPHnhhRcYMWIEb731FmvXruWll17ipZde4o033nAdzSeaN4cTJ2DePNdJcufgwYNcccUVvPnmm6cdCfDiiy/yxhtvMGLECBYvXkxMTAytW7fmWG473b0BlpaW5gW8aWlpgX5pyYPdu3d7PR6Pd968ea6jyBns37/fe8kll3hnzZrlbdasmXfAgAGuI8lpPPzww96mTZu6jpEru3d7vfnyeb1vv+06SWC0b9/e27Nnz1Pe17VrV2+3bt0cJfKtzEyvt0IFr/fBB10nyTuPx+NNTU095X3lypXzvvrqq3/8d1pamrdgwYLeiRMn5uo1tMIk2bJv3z48Hg8lS5Z0HUXO4J577qFDhw60+PvFXhJUpkyZQr169bjhhhsoU6YMderUYdSoUa5jZUtqKmRmQqdOrpMERqNGjZg1axYbNmwAYNWqVcyfP592YXI80OOxVaavvnKdxPc2b97Mjh07uOaaa/54X7FixYiPj2fBggW5ek6f3SUn4cvr9dK/f3+aNGlC9erVXceR05gwYQIrV65k6dKlrqNIFjZt2sTw4cN54IEHeOyxx1i0aBH9+vWjYMGC3HLLLa7jndWkSXDVVVCmjOskgTFo0CDS09O57LLLyJcvH5mZmTz33HMkJia6juYzLVrAuHGwd294bbPu2LEDj8dDmb/9ZS1Tpgw7duzI1XOqYJIs9enTh++//5758+e7jiKnsXXrVvr378+MGTOIjo52HUeykJmZSVxcHM888wwAtWvXZs2aNQwfPjyoC6a0NDtR9fLLrpMEzsSJExk/fjwTJkygevXqrFy5kvvuu4/y5cvTrVs31/F8onlz8HphzpzIWDn0er1ZXoFyJtqSk7Pq27cvn332GV9//TXlypVzHUdOY9myZezevZu6desSHR1NdHQ0c+bM4f/+7/8oUKAA3sDefiRZKFeuHNWqVTvlfdWqVePnn392lCh7pk6FjAzo3Nl1ksB56KGHeOSRR7j++uupUaMGN998MwMGDGDo0KGuo/nM+efb2+zZrpP4VtmyZfF6vezcufOU9+/atesfq07Z5axgSkxMJCEhgeTkZFcRJAt9+/YlNTWV2bNnU7lyZddx5AyuvfZavvvuO1auXMmqVatYtWoV9erV45ZbbmHVqlW5/mlK/KNx48asW7fulPetW7eOKlWqOEqUPZ98AnFxUKmS6ySBc+jQoX/8+4mKiiIzM9NRIv9o0SL8+pguuOACypYty6xZs/54X3p6OosWLaJRo0a5ek5nW3ITJkzQ5btBrE+fPiQnJzN58mRiYmL+qNJjY2MpWLCg43TyVzExMf/oLYuJiaFUqVL/WMkQ9wYMGEDjxo0ZOnQoN9xwA4sWLWLUqFGMHDnSdbQzOngQPv8cnnzSdZLA6tChA8899xyVKlWiRo0aLF++nGHDhtGzZ0/X0XyqeXN4913YvRvOO891muw7ePAgGzdu/GMVfdOmTaxatYqSJUtSqVIl+vfvz7PPPstFF13E+eefzxNPPEHFihXp2LFj7l4wt0f4cktjBUKDx+PxRkVF/eNt7NixrqNJNjRv3lxjBYLYtGnTvLVq1fIWKlTIW716de8777zjOtJZffyx1wte78aNrpME1oEDB7wDBgzwnn/++d7ChQt7L7roIu/gwYO9GRkZrqP51Nat9vn98EPXSXLm66+/Pu33qu7du//xMUOGDPGWK1fOW6hQIW+rVq28GzZsyPXrebzewDY4pKenExsbS1pamlaYRERCwE03wfffw8qVrpOIv1x6qW3NDR/uOknwUtO3iIic0dGj1vDdpYvrJOJPLVqEX+O3r6lgEhGRM/r6a7trLBKOnEey5s1h3Tr49VfXSYKXCiYRETmjlBS44AKoVct1EvGnZs3sUatMZ6aCSURETiszEyZPttUlTacIb6VLQ82aKpjORgWTiIic1tKltkWj7bjIEK73yvmKCiYRETmtlBQoVQpyOedPQkyLFrB5M/z0k+skwUkFk4iInFZKCnToAPl162hEuOoqe5wzx22OYKWCSURE/mHdOvjhB23HRZJSpay5f+5c10mCkwomERH5h9RUKFQIWrZ0nUQCqWlTFUxnost3RUTkH1JToXVrKFzYdRIJpKZNYcMG2L7ddZLgo8t3RUTkFDt2wIIFdiGrRJamTe1x7ly48Ua3WYKNtuREROQUU6bY3KX27V0nkUArWxYuuUTbcqejgklERE6RkmInps4913UScaFpU52UOx0VTCIi8of9+2HmTJ2Oi2RNm8KaNfDbb66TBBcVTCIi8ofPP4djx6BjR9dJxJWrr7bH//3PbY5go4JJRET+kJoKtWvbhbsSmSpXhipVtC33dyqYREQEgIwMmDpVq0uieUyno4JJREQAW1FIS1P/kti23MqV9vdBjAomEREB7HRc5cpwxRWuk4hrTZtCZiZ8843rJMFDBZOIiOD1Wv9Sp042g0ki20UX2Uwm9TH9SQWTiIiwfDls3artODEej23LqY/pTyqYRESElBQoUcIGVoqAbcstWQKHDrlOEhx0+a6IiJCSYleh5Hd2w6gEm6ZN4fhxu1dQdPmuiEjE27gRVq+Gp55ynUSCSfXqULKkDbC85hrXadzTlpyISIRLTYWCBaF1a9dJJJhERUGjRpr4fZIKJhGRCJeSAi1bQkyM6yQSbJo0gYULbWsu0qlgEhGJYLt22awdTfeW02ncGA4cgG+/dZ3EPRVMIiIRbOpUm8HUoYPrJBKM6tWDAgVg/nzXSdxTwSQiEsFSUmwVoXRp10kkGBUsaEWT+phUMImIRKwDB+DLLzWsUs6ucWMrmLxe10ncylPBdPz4cR5++GEuv/xyihQpQoUKFbjtttvYvn27r/KJiIiffPklHD2q/iU5uyZN4Ndf4aefXCdxK08F06FDh1i5ciVDhgxhxYoVfPrpp6xbt46O+tcnIhL0UlKgRg27N0zkTBo1ssdI72PyeL2+XWRbunQp8fHx/PTTT1SsWPEfv56enk5sbCxpaWkaXCki4sjx49a31KcPPPus6zQS7KpVg2bNYPhw10nc8XkP0759+/B4PBQvXtzXTy0iIj4ybx7s3av+Jcmexo21wuTTguno0aMMGjSIm266iSJFivjyqUVExIdSUqBCBahb13USCQVNmtj1Ofv2uU7iTo7ukhs/fjy9e/cGwOPxMH36dBo3bgxYA/j111+Px+PhrbfeyvK5EhMTyf+3Wx6TkpJISkrKSSQREckhr9cKpk6dwONxnUZCQePG9vdmwQJo29Z1Gjdy1MN08OBBdu7c+cd/V6hQgXPOOeePYmnLli189dVXlChR4ozPoR4mERG3VqyAOnVgxgy49lrXaSQUeL1Qtiz06hW5PW85WmGKiYmhatWqp7zvZLG0adMmZs+efdZiSURE3EtJgdhYuPpq10kkVHg8f85jilR56mE6ceIEXbt2Zfny5XzwwQdkZGSwc+dOdu7cSUZGhq8yioiID6WmQrt2EB3tOomEkiZNYPFiOHbMdRI3crTC9Hdbt25l6tSpAFxxxRUAeL1ePB4Ps2fPpmnTpnlPKAGTlgabN9vbli3w229w6JC9HT4M55xjt5kXKWLHkc8/394uvdR+TUSC3+bNsGoVPPaY6yQSaho3tu8FK1ZAfLzrNIGXp4KpSpUqnDhxwldZJICOHIFly2Dhwj/ftm7989cLFbKiKCYGChe2+4SOHbOrFPbvtxvOjx61j42Ohlq17LRNs2bQujWUKuXktyUiWUhNtctU27RxnURCzZVX2veC+fMjs2Dy+eDKrKjp253t2+1m8ilTYOZM+0mhUCGoX9/+8l9xBVStChdcYMXS2U7PZGbCzp2waROsXGnF1+LFsGYNREVBgwZ2+/lNN0HlyoH7PYrI2TVvbv/uP/vMdRIJRc2aQcmSMGmS6ySBp4IpzKWlwYcfwtix9lNBVJQtq3boYKdjatWC/HlaZzzVtm0wfbp9Mf7iCyvKWrWCO+6Azp3VMyHi0u+/2w9Dw4fDnXe6TiOh6LHHYNQo2LEj8kZS+HzSt7jn9cKcOZCUZMdA77rL+o7ee8+20ubOhYEDbXnVl8US2CC8nj3tp48dO2DkSNvCu/FGuPhieO01OHjQt68pItkzdap9fUhIcJ1EQlWTJvZ9ZONG10kCTwVTGDl2DD74AOrVs2XTlSvhqafgl1/g88+hW7fA9hYVLQo9etjK1sqV9g/t/vtt2+/NNyP3pIWIKykptv1etqzrJBKqGja0laVIvCZFBVMY2L8fXnzReo+6dYPzzrPtsO+/h4cegvLlXSeE2rWtmNu4Ea67Du69F6pXty/gIuJ/hw7Z1wXdHSd5Ubw41KwZmfOYVDCFsEOH4OWXbcVm8GCbq7J6ta0mtWoVnPvL558P774L334Ll1xifU0dO8LPP7tOJhLeTh70UMEkeRWpF/GqYApBR4/CG2/AhRfCI49Aly62cjNyJNSo4Tpd9tSsCdOmwccfw9Klttr02mt2+k5EfC8lBS67zOamieRFkyawdq3N6oskKphCiNcL48fbysx999m8o3XrYMQIqFTJdbqc83iga1f44Qfo3t1+T23b2vgDEfGdEydsnEjHjq6TSDho2NAeFy50myPQnBVMiYmJJCQkkJyc7CpCSFmyxKr6m2+2AZFr1sCYMbYdF+qKFYPXX7f+im+/hcsvty/uIuIb33xjqwHajhNfODmrb8EC10kCy8eHyrNvwoQJmsOUDdu3w6OPWnFUqxZ89ZUNngtHrVpZwdSjhx17fuIJePJJmx0lIrmXkmIn4+LiXCevcEVdAAAgAElEQVSRcODx2CpTpBVM+lYUpI4fh1dfte23KVNs0Nzy5eFbLJ103nl2dcPzz8Ozz1pTeHq661QiocvrtYIpIUE/fIjvNGxotzscP+46SeDon08QWrLErit58EG4/XbYsMGGT/p6yGSw8nismX3yZJg9265ZicQhaSK+sGaNXWGk7TjxpYYNbQjx6tWukwSOCqYgkp4O/frZYDmPBxYtst6eEiVcJ3OjfXv7CebECWjUyO6rE5GcSUmxSf8tWrhOIuGkXj37IT6StuVUMAWJadPsaP2778Irr1ihUL++61TuXXaZNaxWrWrTy2fOdJ1IJLSkpNiMtnPOcZ1EwknhwnZhuwomCZi0NLuYtn17Ox32/fcwYEDkbL9lR6lSMGuWnRJs184uExaRrP3yi63MapyA+EOkNX6rYHLoiy9sgOMnn9jK0rRpULmy61TBKSbGeppuuAESE2HsWNeJRILf5Mn2w1e7dq6TSDhq2ND6S3fvdp0kMFQwObB/P9x5J7RpA9WqwXff2eDGYLzKJJhER8N770HPnvbnpaJJ5OxSUuxkbfHirpNIOIq0AZba+AmwWbNsztDvv9uE7l69VCjlRFQUvP22HZU+WWTeeqvrVCLBZ98++Ppr+M9/XCeRcFWlis33WrAAOnRwncb/VDAFyIED8PDD8NZb9hPf11/bRbSSc1FRVmyCjV2IjoakJKeRRILOZ5/ZjJyEBNdJJFydHGD5zTeukwSGCqYAmDvXVkN27LBLc+++WwPk8upk0ZSRYStMJUrYFqeImJQUO/odivdMSuho2NBuZDh+PPwPK+nbth8dOgT9+9tx+AoV7NqPe+5RseQrUVEwapQVSl27Rs4+ukhWjh6F6dM1rFL8r2FD+1737beuk/ifLt/1k2++sRkVI0bYFSdffw0XXug6VfjJnx8mToQ6deC662wsg0ik++orawPQOAHxt7p1rS0iEsYLeLxerzeQL5ienk5sbCxpaWlhefnu4cMweLANn2zQwC7NveQS16nC3969cPXV9rh4MZQr5zqRiDu9e9sBkw0bdKhE/C8+Hi6+GD74wHUS/9LmkA8tWmQrHa+/Di+9BPPmqVgKlBIlbAsiM9O2IQ4fdp1IxI3MTJu/1KmTiiUJjEgZYKmCyQeOHoVHH7X7zooWheXL7eLcfPlcJ4ssFSrYN4qTc60Cu3YqEhwWL7YDJtqOk0Bp2NAueN61y3US/1LBlEeLF9uq0iuvwDPPWO9S9equU0WuunVtWXjiRHjqKddpRAIvJQXOPdd+gBMJhJMDLMN9lUkFUy4dOWJzlRo2tEsIly2zVaZwP1YZCrp0geees4IpTM8UiJxRSorNXtIKtwRKpUpQvrwKJjmNhQvhyittgu5zz9lfkpo1XaeSv3rkEejWzS42XrHCdRqRwFi7Ftat0zgBCayTAyxVMMkfDh+GgQOhcWPrVVqxAgYN0qpSMPJ4bKRD9eo2o2nPHteJRPwvNdVWvK+91nUSiTQNG8KSJTZMOFypYMqmk3OVXn8dhg5Vr1IoKFQIPvkE0tLgllvs9JBIOEtJgdat7e++SCA1amSLCqtWuU7iPyqYsrB3L9x1l60qlShhq0oPPaRVpVBx/vkwfjx8/rmawCW8bd9u7QLajhMX6tSBAgXCe1tOBdMZeL32jfayy+zx9ddh/nyoVs11Msmp1q3h6aftFOOXX7pOI+IfqanW6H3dda6TSCQ65xzr7V20yHUS/1HBdBobN9o32ZtvhqZN4YcfoG9fnToJZY8+Ci1bWiP4jh2u04j43qRJdm9lqVKuk0ikio8P7zs9VTD9xeHDthJRsyasXw/TpsFHH9lARAltUVHw/vv2ePPNcOKE60QivrNnD8yebSM1RFxp0AB+/BF++811Ev/Q5bvY9tvEibb99uyzcN99sGYNtGvnOpn4UunSMG6cfWMZOtR1GhHfmToVjh9X/5K4FR9vj4sXu83hLxF/+e7y5VYg/e9/Nuzt5ZftEkEJX4MH2/ysuXOtmV8k1HXqBLt3W5+liCteL5QpYwelnn7adRrfi9gtuS1b4LbboF49Own35ZfWNKliKfwNHmxLx926wf79rtOI5M2BA/DFF9qOE/c8HltlCtfG74grmHbuhHvvhUsusS8yb7wBK1daQ7BEhvz5rZ9p924YMMB1GpG8+fxzu6qpc2fXSUTsh9FFi8Jz7l3EFEzbt9vdb1Wr2uWsTz1lzWl9+mimUiSqWhWGDYN33rGVRZFQ9cknNlS3alXXSURshSktzQ5OhZuwL5g2boTevW2A4fDh1q+0aZPdNRYT4zqduNSjB3ToAL162cqjSKg5csQavrUdJ8Gifn3bmgvH8QJhWzAtXQo33giXXmorCE89BT//DM8/bxO7RTweGDnS/nevXtawKBJKZs2yHiYVTBIsYmNtwHM49jH5tGDq3bs3UVFRvPbaa7582mxLT7cLV+vWtSp36VJ46y1r8B40CIoXdxJLgliZMlY0TZli23MioWTSJOvH1L2WEkzCtfHbZwVTSkoKixcvpkKApzxmZtrx8J49oXx560mqUMG+Aa5fb9txBQsGNJKEmI4dbXuuf3/YvNl1GpHsOX7cVs+7dLHVUpFg0aABfPstHDrkOolv+aRg2rZtG/369WP8+PHkD0AHdUYGzJlj3+AqVYKrr4YZM6yp+6efYPJkaN9eV5lI9g0bZldK3HmntuYkNMybB7//Dl27uk4icqr4eLtNYdky10l8K8/Vjdfr5dZbb+Whhx6imp9ups3MtPvc5s+3UQAzZ9r2W9mycP311qvUsKFdeyGSG0WL2nZu27YwZgx07+46kcjZTZpkPzDWres6icipatSwQ1ULF8JVV7lO4zt5LpheeOEFChQoQN++fX2Rh4MHYd06WLvWiqQlS+wPPS3NCqL4eHjwQfvGVqeOiiTxnTZt4NZb4f777X+XK+c6kcjpZWbCp5/Cv/6l7TgJPvnz21DocOtjylHBNH78eHr37g2Ax+Nh6tSpvPbaa6xYsSLbz7F0qT32729LdocP22rR9u32tnv3nx9brpwVRQMH2gpS/fq2EiDiL6++aoMA+/a1+TYiwWjJEti2TafjJHjFx8P48a5T+FaO7pI7ePAgO/8ysObDDz/k8ccfx/OXH3FOnDhBVFQUlStXZtOmTf94jvfeS+e222IpVqwt+fLlJ18+q0YLFoT69ZO47rokqlWzcQCxsXn83YnkwkcfwQ03wMcfqz9EgtPDD8Po0fZDpno1JRh9+qkV9Fu32kGscJCny3f37t3L9u3bT3lfq1atuPXWW+nevTsXn+ZitmC7fFfk77xe+4e+YAF8/z2ULOk6kcifvF6787JFC/jvf12nETm9X3+1QumTT8JnJTRPHUAlSpSgevXqp7xFR0dTtmzZ0xZLIqHA47H5XUeOWD+TSDD57ju71ilcvglJeCpf3g4lhNPEb5+3THvUgShhoFw562caO9ZOZooEi0mToFgxW2ESCWbhNsAyT1tyuaEtOQkVXi+0bGl3D65eDYULu04kArVqQe3adom4SDB7+WUYMsROuYfDJfc6lC9yBie35rZtg6FDXacRsVErq1fb/DmRYNeggU37Xr3adRLfUMEkchaXXGL3EL70ks0HE3Hpo49stErr1q6TiGStTh07xRku23IqmESyMGgQVKxos5l0bYq49NFHkJCg+zElNBQubNvH4dL4rYJJJAuFCsEbb9iVPB9+6DqNRCptx0koCqfGbxVMItnQtq0d4x4wwCbTiwSatuMkFDVoYMX+vn2uk+SdCiaRbPrPf6xYGjzYdRKJRNqOk1AUH2+PS5a4zeELKphEsqlSJXjySXj9dVi50nUaiSTajpNQdfHFULx4eGzLqWASyYH77oNq1eDuu+3GeJFA0HachKqoKFtlCofGb2cFU2JiIgkJCSQnJ7uKIJJj0dEwfLj943/nHddpJFJ89BF06KDtOAlNJxu/Q/2UsSZ9i+TC7bfD1Kmwfr0u5xX/+uEHqF4dUlKgY0fXaURybto0aN/ebk244ALXaXJPW3IiufDCC3DsmPU0ifiTtuMk1MXF2WOo9zGpYBLJhbJl4Ykn7OqUcBn7L8FJ23ES6s47D6pWVcEkErHuu8++CPTvH/p78xKcTp6Ou+EG10lE8iYcBliqYBLJpQIFYNgwmDULUlNdp5FwpO04CRfx8bB8OWRkuE6SeyqYRPKgXTto0wbuvx+OHHGdRsKNtuMkXMTFwdGj8O23rpPkngomkTzweGyV6Zdf7FHEV9au1XachI8rr7SxLKG8LaeCSSSPLrsM+vWD556Dbdtcp5FwMWECFCum7TgJDwULQu3aKphEIt7gwVC4MAwa5DqJhAOvF8aPh86dtR0n4SPUG79VMIn4QGwsPP88fPABLFjgOo2EuuXLYcMGuOkm10lEfCc+Htatg337XCfJHRVMIj7SvTvUqWPbc7pnTvIiORlKl4YWLVwnEfGdkwMslyxxmyO3VDCJ+Ei+fPDaa7B0KYwd6zqNhKrMTOtfuv56yJ/fdRoR37n4YihePHS35VQwifhQ48aQlASPPgoHDrhOI6Fo3jw7PKDtOAk3UVG2yqSCKYcSExNJSEggOTnZVQQRvxg6FPbuhX//23USCUXJyVClCjRs6DqJiO+dbPwOxdsRPF5vYGOnp6cTGxtLWloaxYoVC+RLiwTMo4/Cf/4D69dDxYqu00ioOHYMypWDXr3sgmeRcDN1qg1j3bwZzj/fdZqc0ZaciB8MGmRXWjz2mOskEkpmzIA9e7QdJ+ErPt4eQ3FbTgWTiB8UKwZPPw3vvWdN4CLZkZwM1atDrVquk4j4x3nnwQUXqGASkb/o0QNq1IAHHgjN/XoJrEOHICXFDg14PK7TiPhPqA6wVMEk4if588Mrr8DcuZCa6jqNBLspU+DgQSuYRMJZfLwNZ83IcJ0kZ1QwifhR69bQpg0MHGgNvSJnkpxsR64vvNB1EhH/iouDI0fgu+9cJ8kZFUwifvbyy7BpE7z1luskEqz27oXp09XsLZHhyittBT7UtuVUMIn4WY0adkz86aftBJTI302aBMePww03uE4i4n+FCkHt2iqYROQ0nnrKviE+/bTrJBKM3n8fmje3GUwikSAUG79VMIkEQJkyNszyzTdtmKXISVu2wJw5cNttrpOIBE58PKxdC2lprpNknwomkQDp3x/Kl4eHHnKdRILJBx9ATAx07uw6iUjgxMXZ45IlbnPkhAomkQApWNCuu0hNtRUFEa/Xhpv+619QpIjrNCKBc8klEBsbWttyunxXJIASE+0nq/vvh8xM12nEtYULYcMGuPVW10lEAisqyr4WhlLBlN/VC0+YMEGX70rE8XhsmOVVV8H48XDLLa4TiUvvvQeVKkGzZq6TiARefDz897+20hoK0+21JScSYE2aQJcu1gR++LDrNOLKkSMwYYIVzVH6SiwRKC4Odu2Cn392nSR79M9UxIEXX4QdO2DYMNdJxJWpU2HfPm3HSeSKj7fHUNmWU8Ek4sBFF8E998DQobBzp+s04sLYsfYT9mWXuU4i4kbp0nD++SqYRCQLTzwB0dEwZIjrJBJou3bZVShaXZJIF0oDLFUwiThSsqQVTSNHwpo1rtNIICUnW99SYqLrJCJuxcfDsmWQkeE6SdZ8UjD98MMPdOzYkeLFi1OkSBHi4+PZunWrL55aJKzdcw9ccAEMHOg6iQTS2LHQvj2UKuU6iYhbcXF2AGL1atdJspbngunHH3/kqquuonr16sydO5fvvvuOJ554goIFC/oin0hYK1DAGsCnT4cZM1ynkUBYvRpWrNB2nAhAnTqQP39obMt5vF6vNy9PkJSURIECBRg7dmy2Pj49PZ3Y2FjS0tI0h0kEm0Fy1VWwfz8sXw758rlOJP50//122e62bVYwi0S6unXh8sth9GjXSc4uTytMXq+XadOmcfHFF9OmTRvKlClDgwYNSE1N9VU+kbB3cpjlt9/aVo2Er6NHbVjlbbepWBI5KVQav/NUMO3atYsDBw7w4osv0q5dO2bMmEHnzp3p0qUL8+bN81VGkbAXH28NwI8/DgcOuE4j/pKaCr//Dj16uE4iEjzi42HtWkhLc53k7HJUMI0fP56iRYtStGhRihUrxrp16wDo1KkT/fr14/LLL+fhhx+mffv2vP32234JLBKuhg6FPXvg5ZddJxF/GTUKGjWCatVcJxEJHnFx1pqwdKnrJGeXo7vkOnbsSIMGDf7473PPPZf8+fNT7W//+qtVq8b8+fPP+lyJiYnkz3/qyyclJZGUlJSTSCJh4/zz4b774N//hjvvhPLlXScSX9qyBWbOhHfecZ1EJLhceinExtq23DXXuE5zZjkqmGJiYqhateop76tfv/4fK00nrV+/nipVqpz1uXT5rsg/PfoovPuubc29+67rNOJLo0dDkSJw/fWuk4gEl6goqF8/+PuY8jxWYODAgUycOJFRo0bx448/8sYbbzB16lTuueceX+QTiSixsfDkkzBmDKxc6TqN+MqJE1YAJyVZ0SQipzrZ+J23c/v+leexAgBjxozh+eefZ9u2bVx66aU8/fTTtG/f/rQfq7ECImeXkQG1akHFijabyeNxnUjyavp0aNcOFi+2n6RF5FSTJ0PHjvDTT1C5sus0p+eTgiknVDCJZO3kF49p0+wbrYS2rl1h40ZbNVQBLPJPO3dC2bLw4YfBu22tu+REglCHDtCsGTz4IBw/7jqN5MWuXVYA9+ihYknkTMqUgSpVgruPSQWTSBA6Ocxy7Vo7ii6ha8wYm95+882uk4gEt2AfYKmCSSRI1akD3brB4MGQnu46jeTGiRPw9ttw4426aFckK/HxsGyZ9XEGIxVMIkHsueds8vcLL7hOIrnxxReweTPcfbfrJCLBLy4ODh+GNWtcJzk9FUwiQaxiRXjgARg2DH7+2XUayanhw+HKK+0nZxE5uzp1bPs6WLflVDCJBLmHHrL5TI895jqJ5MSWLXbKsU8fNXuLZEfhwnD55SqYRCSXihaFp5+GDz4I/ruW5E8jRkCxYjasUkSyJ5gbv1UwiYSAO+6AGjVsey6YJ+GKOXrUTjfefjvExLhOIxI64uPhhx+C86CLs4IpMTGRhIQEkpOTXUUQCRn588PLL8PcuZCa6jqNZOXjj+G339TsLZJTcXH2Q2EwrqZr0rdICGnTxiZGr1kD55zjOo2cSZMm9vmZNct1EpHQkpkJJUrAoEHwyCOu05xKW3IiIWTYMLtradgw10nkTFasgPnzrdlbRHImKsruWwzGPiYVTCIhpFo16NsXnn0Wfv3VdRo5nWHD7IqHjh1dJxEJTScbv4OtX1MFk0iIGTLEjt8OGuQ6ifzd9u0wYQL062d9ZyKSc/HxsGMHbN3qOsmpVDCJhJjixW0C+Pvvw8KFrtPIX735pvUu9ejhOolI6IqLs8dg25ZTwSQSgu64wyZI9+tnTZLi3uHDdm9cjx42aFREcqdsWahcWQWTiPhAvnzw2muwZAmMHes6jYCt+O3ZY0WsiORNMA6wVMEkEqKaNLEp0o88EpxD3iKJ1wv/+Q906gRVq7pOIxL64uNh2TI4ftx1kj+pYBIJYS+9BPv3wzPPuE4S2b74wqYTDxjgOolIeIiLg0OHbOZcsFDBJBLCKla0Fab/+z9Yt851msg1bBjUrWurfiKSd3XrWutBMG3LqWASCXEPPAAVKsB99wXf3JJI8O238OWXtrrk8bhOIxIeCheGWrVUMImIDxUqZCtMX3wBn3ziOk3keeEFG1R5ww2uk4iEl2Br/NbluyJhICHB3vr3t54mCYwff4SJE2HgQIiOdp1GJLzEx8P33wfPoRZdvisSJn76CapXh9694dVXXaeJDHfdBZ9+Clu22EqfiPjOmjVQsyZ89RU0b+46jbbkRMJGlSoweLDNZ1q1ynWa8Ld9O4webat6KpZEfO+yy6Bo0eDZllPBJBJGBgyASy+Fu+/WBHB/GzYMChaEPn1cJxEJT/nyQf36KphExA8KFIDhw2HBAnjnHddpwtfevfbn3KePrkER8aeTjd/BcAJYBZNImGnaFG69FR5+GHbvdp0mPL36qq3g9e/vOolIeIuPt+3vbdtcJ1HBJBKW/v1ve3zoIbc5wtHvv9s1KH37QpkyrtOIhLe4OHsMhm05FUwiYah0aRg6FMaMgblzXacJLy+/bI8DB7rNIRIJypWDSpVUMImIH/XqBQ0a2JiBI0dcpwkPu3bB66/DvffCuee6TiMSGYJlgKUKJpEwFRUFI0facMXnnnOdJjz8+9/25/rAA66TiESO+HhYuhSOH3ebQwWTSBirWRMee8yu7/j2W9dpQtvOnfDmm9boXaqU6zQikSMuDg4dsqnfLqlgEglzjzxiA+B69HD/E1ooGzrUxjYMGOA6iUhkqVvXZjK53pZTwSQS5goUsJlMy5fb6S7JuR9/hLfeskbvEiVcpxGJLDExtlqugklE/C4uzraSnngC1q93nSb0PPaYnTzU6pKIG8HQ+O2sYEpMTCQhIYHk5GRXEUQiytNP2/Hcbt20NZcTixfDxInwzDNQuLDrNCKRKT7eLuPdv99dBo/XG9iB4+np6cTGxpKWlkaxYsUC+dIiEW/hQmjcGJ580lab5Oy8XmjWzK5CWbHC+ihEJPBWr4ZatWD2bPs36YK25EQiSIMGtr301FN2TFfObsoUG/z50ksqlkRcqlYNihRxuy2nFSaRCJORAQ0bwsGD1gheqJDrRMHp6FG4/HKoXBm+/BI8HteJRCJbixZQvDhMmuTm9bXCJBJhoqPh/fdhyxYYNMh1muD16quwaZOdLFSxJOKe68ZvFUwiEahaNXjxRXjtNZgxw3Wa4PPzz/Dss9CvH9So4TqNiIAVTL/+Clu3unl9FUwiEapvX7jmGuje3Zqa5U8PPADFisGQIa6TiMhJcXH2uHixm9fPc8F08OBB+vbtS6VKlShcuDA1atRgxIgRvsgmIn4UFQVjxtiVA92724kwgZkz4eOP4ZVXrGgSkeBQvjxUrOhuWy7PBdOAAQP48ssvGT9+PGvXrqV///707duXqVOn+iKfiPhRxYpWNKWmago4WKN3375w9dWQlOQ6jYj8ncs+pjwXTAsWLOC2227jqquuonLlyvTq1YvatWuz2NWamYjkSEICPPggPPSQzWmKZM88Y43eb76pRm+RYBQfbyNRTpwI/GvnuWBq1KgRkydP5tdffwVg9uzZbNiwgdatW+c5nIgExvPPW3/ADTfA77+7TuPGihXwwgs20FON3iLBKS7ORqJ8/33gXzvPBdPrr79OtWrVqFixIgUKFKBdu3a8+eabNG7c2Bf5RCQAoqNhwgTrZ7rtNsjMdJ0osDIy4I47rFDSqAWR4FW3rvVfutiWy1HBNH78eIoWLUrRokUpVqwY8+fP57XXXmPRokVMnTqV5cuX88orr9CnTx+++uorf2UWET+oVMnmM02bZkfqI8mLL8J338G771rxKCLBqUgRqFnTTcGUo0nfBw8eZOfOnX/8d/ny5YmNjSU1NZU2bdr88f5evXqxbds2Pvvss388x8lJ323btiV//vyn/FpSUhJJ6rQUceqZZ2DwYDsp1rWr6zT+t2iR3a/38MPw3HOu04hIVu680/otv/02sK+bP+sP+VNMTAxVq1b947/3799PRkYGnr91R+bLl4/MLNb0J0yYoKtRRILQ44/breC33goXXghXXOE6kf/s3w833WTL/E8+6TqNiGRHfDy88w4cOGArToGSpx6mokWLcvXVVzNw4EDmzJnDli1bGDNmDO+99x5dunTxVUYRCSCPx7amqlWzE3R/WVQOO337wq5dMH68tuJEQkVcnPVZLlsW2NfNc9P3xIkTqV+/Prfccgs1atTgpZdeYujQodx5552+yCciDhQuDCkp1gzdpQscOeI6ke+NGwfvvWcjBC680HUaEcmu6tVtZSnQfUw56mHyhZM9TGlpadqSEwlyixZB8+bQujV89BHkz9EmfvBauRIaNbIerffe08wlkVDTvDmULAmffBK419RdciJyRvHxVihNmQJ33x0e16f89ht06mRbjv/9r4olkVDkYuK3CiYROavrrrMGy1GjbKhjKDt+HG680eZNffopFCrkOpGI5EZ8PGzbZm+BEiYL7CLiT7fdBrt3w8CBULy4XaUSarxeGDAA5syBWbOgcmXXiUQkt+Li7HHxYujcOTCvqYJJRLLlwQdh3z4rmo4fD72J2C++CG+8AW+/bZfrikjoqlDB3hYtUsEkIkHomWes8fuRR6xoevxx14myZ/RoyzxkCPTu7TqNiPhCoPuYVDCJSLZ5PDbgMX9+62c6csSKqGBunJ4wAXr1sunAQ4a4TiMivhIfb19/TpyAfPn8/3oqmEQkxx5/HAoWtO25bdvstFkwDn5MToZbboGbb4a33gruwk5EciYuzqZ9//CD3S/nbzolJyK58uCDNvxx3Dg7Sff7764Tneqtt6xY6tbNtuQC8ROoiAROvXoQFRW4bTlnBVNiYiIJCQkkJye7iiAieXTTTfDFF7B8ud3HFuirCk4nM9Mu0r3nHujXz0YiqFgSCT9FikCNGoErmDTpW0Ty7OefbWr2d9/ZVSN33OFm+2vfPujRw2Ysvfoq9O8f+AwiEji9etlogVWr/P9a2pITkTyrXBnmzbN5TT172jHf7dsDm2HJEqhTx2YsTZqkYkkkEsTHw+rV1svkbyqYRMQnChaEESPsbqeFC+2CzHfftS0yfzpyBJ56Cho3hvPOgxUr7OoTEQl/cXH2NWb5cv+/lgomEfGpLl1gzRpo3962x+rVg5kzff86Xi9Mnw6XXw7PPWcn9ubNgwsu8P1riUhwqlEDYmIC08ekgklEfK5UKXj/ffjf/+y+tpYt4aqrbKvsxIm8PXdmphVKjRpBu3ZQvjysXGlFU4ECvskvIqEhXz77oUwFk4iEtMaNrQoSaIoAAAJ/SURBVGhKSbH/7toVqlSB+++HBQuyXzxlZlpT5+DBtoLUrp29/4svYPZs2/4TkcgUqInfOiUnIgGzdCmMHQsffww7dkCxYtCggRU8F14IsbF2VPjECdizBzZtgnXrbKvt99+haFFITITu3e3/p0GUIjJpkv0wtm2brTj7iwomEQm4EyesMXzOHFtpWrcOtmyBjIw/PyYqyk7fXXQRNGwIzZvbY8GCzmKLSBDauhUqVbJxIv488KGrUUQk4PLls+26xo3/fJ/XC8eOwf79dlddkSL2KCJyNhUr2srSokUqmEQkAng8cM459iYikhPx8bZq7U9q+hYREZGQ1qiRTfw+ftx/r6GCSUREREJaw4Zw6BB8+63/XkOX74qIiEhIq1sXoqPhm2/89xo6JSciIiIhr0EDG08ybpx/nl9bciIiIhLyGjXy7wqTCiYREREJeY0a2Ty37dv98/wqmERERCTkNWxojwsW+Of5VTCJiIhIyKtQwW4H8Ne2nAomERERCQsNG2qFSUREROSsGjWyS76PHvX9c6tgEhERkbDQsKHdSblihe+fWwWTiIiIhIUrroBChfzTx6SCSURERMJCdDTUr6+CSUREROSsTjZ++/oek/y+fToRERERd9q0gT17rPG7YEHfPa/ukhMRERHJgrMtucTERBISEkhOTnYVQURERCRbtMIkIiIikgU1fYuIiIhkQQWTiIiISBZUMImIiIhkQQWTiIiISBYC3vTt9XrZv38/RYsWxePxBPKlRURERHIl4AWTiIiISKjRlpyIiIhIFlQwiYiIiGRBBZOIiIhIFlQwiYiIiGRBBZOIiIhIFlQwiYiIiGRBBZOIiIhIFv4f3b7fx9GnQgYAAAAASUVORK5CYII=" }, "execution_count": 26, "metadata": { }, "output_type": "execute_result" } ], "source": [ "plot(x*sin(x)-4, (x, 0, 10))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "We can get it by forcing the search in the interval $[6,8]$:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}6.90141260902$ " ] }, "execution_count": 27, "metadata": { }, "output_type": "execute_result" } ], "source": [ "find_root(eq, 6, 8)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(6.90141260902, \\begin{array}{l}\n", "\\phantom{\\verb!xxxxxx!}\\verb|converged:|\\phantom{\\verb!x!}\\verb|True|\\\\\n", "\\phantom{\\verb!xxxxxxxxxxx!}\\verb|flag:|\\phantom{\\verb!x!}\\verb|'converged'|\\\\\n", "\\phantom{\\verb!x!}\\verb|function_calls:|\\phantom{\\verb!x!}\\verb|10|\\\\\n", "\\phantom{\\verb!xxxxx!}\\verb|iterations:|\\phantom{\\verb!x!}\\verb|9|\\\\\n", "\\phantom{\\verb!xxxxxxxxxxx!}\\verb|root:|\\phantom{\\verb!x!}\\verb|6.901412609021304|\n", "\\end{array}\\right)$ " ] }, "execution_count": 28, "metadata": { }, "output_type": "execute_result" } ], "source": [ "find_root(eq, 6, 8, full_output=True)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Approximate solution to a system: `mpmath.findroot`" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}x \\sin\\left(x\\right) - y = 0$ " ] }, "execution_count": 29, "metadata": { }, "output_type": "execute_result" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}y \\cos\\left(y\\right) - x - 1 = 0$ " ] }, "execution_count": 29, "metadata": { }, "output_type": "execute_result" } ], "source": [ "eq1 = x*sin(x) - y == 0\n", "eq2 = y*cos(y) - x - 1 == 0\n", "for eq in [eq1, eq2]:\n", " show(eq)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( x, y \\right) \\ {\\mapsto} \\ x \\sin\\left(x\\right) - y$ " ] }, "execution_count": 30, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f1(x,y) = eq1.lhs()\n", "f1" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( x, y \\right) \\ {\\mapsto} \\ y \\cos\\left(y\\right) - x - 1$ " ] }, "execution_count": 31, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f2(x,y) = eq2.lhs()\n", "f2" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "For solving a system, we use the [mpmath](http://mpmath.org/) toolbox:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "import mpmath" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{l}\n", "\\verb|[-0.642871006593746]|\\\\\n", "\\verb|[|\\phantom{\\verb!x!}\\verb|0.385398467857987]|\n", "\\end{array}$ " ] }, "execution_count": 33, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f = [lambda a,b: f1(RR(a), RR(b)), lambda a,b: f2(RR(a), RR(b))]\n", "sol = mpmath.findroot(f, (0, 0))\n", "sol" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left[-0.642871006593746\\right], \\left[0.385398467857987\\right]\\right]$ " ] }, "execution_count": 34, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol.tolist()" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(-0.642871006593746, 0.385398467857987\\right)$ " ] }, "execution_count": 35, "metadata": { }, "output_type": "execute_result" } ], "source": [ "x1 = RR(sol.tolist()[0][0])\n", "y1 = RR(sol.tolist()[1][0])\n", "x1, y1" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}-5.55111512312578 \\times 10^{-16}$ " ] }, "execution_count": 36, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f1(x1,y1)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}1.11022302462516 \\times 10^{-16}$ " ] }, "execution_count": 37, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f2(x1,y1)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(-0.642870824741,\\,0.385398567605\\right)$ " ] }, "execution_count": 38, "metadata": { }, "output_type": "execute_result" } ], "source": [ "sol2 = minimize(eq1.lhs()^2 + eq2.lhs()^2, (-0.6,0.4))\n", "sol2" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(-0.642870824741, 0.385398567605\\right)$ " ] }, "execution_count": 39, "metadata": { }, "output_type": "execute_result" } ], "source": [ "x2, y2 = sol2[0], sol2[1]\n", "x2, y2" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(-3.02336881641 \\times 10^{-07}, -1.0387405891 \\times 10^{-07}\\right)$ " ] }, "execution_count": 40, "metadata": { }, "output_type": "execute_result" } ], "source": [ "f1(x2,y2), f2(x2,y2)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(-1.81852617853195 \\times 10^{-7}, -9.97466205743258 \\times 10^{-8}\\right)$ " ] }, "execution_count": 41, "metadata": { }, "output_type": "execute_result" } ], "source": [ "x1 - x2, y1 - y2" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.1", "name": "sage-8.1" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.14" } }, "nbformat": 4, "nbformat_minor": 0 }