{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solving equations with SageMath" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%display latex" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exact solutions : `solve`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us start with a simple example: solving $x^2-x-1=0$ for $x$:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]\$" ], "text/latex": [ "$\\displaystyle \\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]$" ], "text/plain": [ "[x == -1/2*sqrt(5) + 1/2, x == 1/2*sqrt(5) + 1/2]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(x^2-x-1 == 0, x)" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [], "source": [ "eq = x^2-x-1 == 0" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle x^{2} - x - 1 = 0\$" ], "text/latex": [ "$\\displaystyle x^{2} - x - 1 = 0$" ], "text/plain": [ "x^2 - x - 1 == 0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle x^{2} - x - 1\$" ], "text/latex": [ "$\\displaystyle x^{2} - x - 1$" ], "text/plain": [ "x^2 - x - 1" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq.lhs()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 0\$" ], "text/latex": [ "$\\displaystyle 0$" ], "text/plain": [ "0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq.rhs()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]\$" ], "text/latex": [ "$\\displaystyle \\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]$" ], "text/plain": [ "[x == -1/2*sqrt(5) + 1/2, x == 1/2*sqrt(5) + 1/2]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(eq, x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solutions are returned as a list:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]\$" ], "text/latex": [ "$\\displaystyle \\left[x = -\\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}, x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\\right]$" ], "text/plain": [ "[x == -1/2*sqrt(5) + 1/2, x == 1/2*sqrt(5) + 1/2]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol = solve(eq, x)\n", "sol" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\$" ], "text/latex": [ "$\\displaystyle x = \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}$" ], "text/plain": [ "x == 1/2*sqrt(5) + 1/2" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[1]" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x == 1/2*sqrt(5) + 1/2\n" ] } ], "source": [ "print(sol[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The access to the value of the solution is via the `rhs()` method:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}\$" ], "text/latex": [ "$\\displaystyle \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2}$" ], "text/plain": [ "1/2*sqrt(5) + 1/2" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[1].rhs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A numerical approximation (recall that `n` is a shortcut alias for `numerical_approx`):" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1.61803398874989\$" ], "text/latex": [ "$\\displaystyle 1.61803398874989$" ], "text/plain": [ "1.61803398874989" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n(sol[1].rhs())" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1.6180339887498948482045868343656381177203091798058\$" ], "text/latex": [ "$\\displaystyle 1.6180339887498948482045868343656381177203091798058$" ], "text/plain": [ "1.6180339887498948482045868343656381177203091798058" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n(sol[1].rhs(), digits=50)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A new equation involving the above solution:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2} = \\phi\$" ], "text/latex": [ "$\\displaystyle \\frac{1}{2} \\, \\sqrt{5} + \\frac{1}{2} = \\phi$" ], "text/plain": [ "1/2*sqrt(5) + 1/2 == golden_ratio" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[1].rhs() == golden_ratio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Asking Sage whether this equation always holds:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\mathrm{True}\$" ], "text/latex": [ "$\\displaystyle \\mathrm{True}$" ], "text/plain": [ "True" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bool(_)" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [], "source": [ "y = var('y')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we may form the system:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(x^{2} + x y + 2 = 0, y^{2} = {\\left(x + y\\right)} x\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(x^{2} + x y + 2 = 0, y^{2} = {\\left(x + y\\right)} x\\right)$" ], "text/plain": [ "(x^2 + x*y + 2 == 0, y^2 == (x + y)*x)" ] }, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\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]\$" ], "text/latex": [ "$\\displaystyle \\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]$" ], "text/plain": [ "[[x == -1/2*I*sqrt(10) + 1/2*I*sqrt(2), y == -I*sqrt(2)], [x == 1/2*I*sqrt(10) + 1/2*I*sqrt(2), y == -I*sqrt(2)], [x == -1/2*I*sqrt(10) - 1/2*I*sqrt(2), y == I*sqrt(2)], [x == 1/2*I*sqrt(10) - 1/2*I*sqrt(2), y == I*sqrt(2)]]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol = solve((eq1, eq2), x, y)\n", "sol" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left[x = -\\frac{1}{2} i \\, \\sqrt{10} + \\frac{1}{2} i \\, \\sqrt{2}, y = -i \\, \\sqrt{2}\\right]\$" ], "text/latex": [ "$\\displaystyle \\left[x = -\\frac{1}{2} i \\, \\sqrt{10} + \\frac{1}{2} i \\, \\sqrt{2}, y = -i \\, \\sqrt{2}\\right]$" ], "text/plain": [ "[x == -1/2*I*sqrt(10) + 1/2*I*sqrt(2), y == -I*sqrt(2)]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[0]" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle y = -i \\, \\sqrt{2}\$" ], "text/latex": [ "$\\displaystyle y = -i \\, \\sqrt{2}$" ], "text/plain": [ "y == -I*sqrt(2)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[0][1]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle -i \\, \\sqrt{2}\$" ], "text/latex": [ "$\\displaystyle -i \\, \\sqrt{2}$" ], "text/plain": [ "-I*sqrt(2)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[0][1].rhs()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle -1.41421356237310i\$" ], "text/latex": [ "$\\displaystyle -1.41421356237310i$" ], "text/plain": [ "-1.41421356237310*I" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n(sol[0][1].rhs())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximate solution: `find_root`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us consider the following transcendental equation:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle x \\sin\\left(x\\right) = 4\$" ], "text/latex": [ "$\\displaystyle x \\sin\\left(x\\right) = 4$" ], "text/plain": [ "x*sin(x) == 4" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eq = x*sin(x) == 4\n", "eq" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left[x = \\frac{4}{\\sin\\left(x\\right)}\\right]\$" ], "text/latex": [ "$\\displaystyle \\left[x = \\frac{4}{\\sin\\left(x\\right)}\\right]$" ], "text/plain": [ "[x == 4/sin(x)]" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(eq, x)" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 8.962126200824814\$" ], "text/latex": [ "$\\displaystyle 8.962126200824814$" ], "text/plain": [ "8.962126200824814" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "find_root(eq, 0, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "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": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGFCAYAAAAPa6wiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABC50lEQVR4nO3deZyNdf/H8ddlhrFkRritjSKSUSHOmGxpkbqlUqm0qLTpR0X7erfcRft9V5R0t6pulcioKEmKkJ0ItzYkt6hmLBkx1++PT9wIszjX+V7nnPfz8TiPkzHnXG+dWT7nu3y+nu/7iIiIiMjelXEdQERERCTsVDCJiIiIFEEFk4iIiEgRVDCJiIiIFEEFk4iIiEgRVDCJiIiIFEEFk4iIiEgRVDCJiIiIFEEFk4iIiEgRQl0weSbd8zzPdRYRERFJXqmxupDnebcBA4AnfN/vBxR5JkteXh4ZGRnk5eUFHU9ERESSU7EGZWIywuR5XgS4Epgfi+uJiIiIRFPgBZPneQcArwFXAL8EfT0RERGRaIvFCNNg4D3f9z8q6hMLCgrIz8/f5SYiIiLiWqAFk+d55wFHA7cV5/MHDhxIRkbGjltmZmaQ8URERESKxfP9Itdel+6JPS8TmAmc5Pv+vD8+9gkwd2+LvgsKCigoKNjx5/z8fDIzM8nLyyM9PT2QnCIiIsWxZYvdly0L2rudUIr1aga5S64lUAOYtVNXgBSgg+d5fbdu3UpKSsouD0hLSyMtLS3ASCIiIkUrLIQvvoCRI+Hzz2HhQvj1V/u7cuXgsMOgRQvo2hVOOQUOOMBpXImBIEeYKgMH7/bhF4HFwEO+7y8o6jny8/N3tBXQCJOIiARt61Z48UV47DFYsgRq1IDjjoNmzaB2bShTBvLyYNEimDoV5s2DypXhiivgxhvtcyTuuB1h8n1/PfDlzh/zPG8jsM73/S/3/CgRERE3JkyA666zYujMM+GZZ6BDB9htMmQX33wDzz8PTz8Nzz0H990HfftCasy6HEqshLrTt4iISNA2bYIrr4QTT4QDD4SZM2HECBtZ2lexBNCgATzwgBVOF14I118P7dvD8uWxyS6xE9iUXDEUeWFNyYmISJC+/dbWIX3zDTzxBFx++f4t6J46Fc47DzZs+F/RJaEXnk7fIiIiYTNtGmRnw+bNNqp0xRX7v/vtmGNgzhxo2RI6d4Z//zs6WcU9FUwiIpJ0pkyBTp2gcWMrnLKyovfcVavCe+9Bjx5wwQUwbFj0nlvc0bI0ERFJKlOmwMkn2yjQe+9BpUrRv0bZsrbbrmxZuOQSKF8euneP/nUkdlQwiYhI0pg5M/hiabsyZWDoUPjtN7joIqhTB9q2De56Eiwt+hYRkaSwfDm0bg0HH2wtBIIslnZWUAAnnWTNL6dNg4YNY3NdKTYt+hYREQHIz4cuXWxqLDc3dsUSQFoajBpla5u6d7dF5hJ/VDCJiEhC83249FIbYXr/feveHWtVq8Jbb8FXX0H//rG/vuy/UBZMgwcPJisri0gk4jqKiIjEuX/+086Ee/llaNLEXY5mzeDJJ2HIELUbiEdawyQiIglr+nRo1w769YNHHnGdxka7LrwQxoyxNU2Zma4TCcVcw6SCSUREEtLGjdC8OVSrBp99Zlv8w+DXX6FpUzjqKJsi3N9mmbLftOhbRESS1403wqpV1jgyLMUSQJUq1m5g3DibJpT4oIJJREQSztixtlbo0UehUSPXaf6sSxfrzdS/P6xZ4zqNFIem5EREJKFs2GBHnTRpYqM4YZ3yWrcODjsMTj8dXnjBdZqkpik5ERFJPvfcAz/9ZCNMYS2WwNZWPfCAHaEybZrrNFIUjTCJiEjCmDfPjj35+9/htttcpynatm2QnW3//cUXkJLiNk+S0giTiIgkj8JC6N0bGjeGG25wnaZ4UlJg0CCYPRtee811GtkXFUwiIpIQ/vUvm9oaMgTKlXOdpviOOQbOOgvuukvHpoSZCiYREYl7+flw553Qsye0b+86Tck98AD88AM8/bTrJLI3KphERCTuDRxou+MGDHCdpHQaN4bLLrPCKS/PdRrZExVMIiIS177/Hv7xD2tUWbeu6zSld/fd8Ntv4TjCRf4slAWTDt8VEZHiuv12OPBAuPlm10n2T5060LevHdD788+u08ju1FZARETi1owZti3/uefg8stdp9l///0vHHII3HqrjThJTOjwXRERSWwnnwzLl8OCBYnTw+i66+CVV2yqUb/6YkJ9mEREJHF9/jl88IF19k6UYglsanHTJhg82HUS2ZlGmEREJC516gSrV1t37zIJ9vb/6qthxAgbZapY0XWahKcRJhERSUyffQYffWSjS4lWLIGNMv38Mwwb5jqJbKcRJhERiTvHH28FxezZiVkwAZx9Nnz5JSxalLj/xpDQCJOIiCSeTz6BiRPh3nsTu5Do3x+WLIFx41wnEdAIk4iIxJmOHWH9epg5E7xijQ3EJ9+HnByoXNmmHyUwGmESEZHEMnUqTJpkB9UmcrEE9u/r3x8mTID5812nEY0wiYhI3OjWDb76KnnW9fz+Oxx6KJxwArz4ous0CUsjTCIikjgWL4bRo+Gmm5KjWAIoWxauuQZefx1++sl1muSWJF9yIiIS7x55BGrXhgsvdJ0ktnr1sum5l15ynSS5hbJg0uG7IiKys1WrrCdRv36QluY6TWxVq2YtBoYOhcJC12mSl9YwiYhI6N18Mzz7rJ0bl5HhOk3sTZ4M7dvD+PFw4omu0yQcrWESEZH4t369FUu9eydnsQTQti1kZdn/B3FDBZOIiITaK6/Axo3Qt6/rJO54Hlx1Fbzzjp2fJ7GngklERELL92HQIGsnkJnpOo1bF11ku+ZeeMF1kuSkgklERELro4+sncA117hO4t6BB8K552rxtysqmEREJLSefBKOOsoWPAtcdhl8/z18+qnrJMlHBZOIiITS11/De+/Btdcm/jEoxdW2LTRoYOu6JLZUMImISCgNHmzTUOef7zpJeHge9OwJb71lC+EldgItmDzPu83zvBme5633PG+N53nveJ7XOMhriohI/NuwwRY3X3EFVKjgOk24XHSR/f955x3XSZJL0CNMxwKDgRygE5AKfOh5XqWArysiInFs2DDrv3T11a6ThE+DBramS9NysZUa5JP7vn/yzn/2PO9SYA3QMsjriohI/NreSuCMM+Dgg12nCaeePa0v0w8/QN26rtMkh1ivYdreo/XnPf1lQUEB+fn5u9xERCS5fPYZLFoEffq4ThJe3btDuXLw2muukySPmBVMnud5wOPAZN/3v9zT5wwcOJCMjIwdt8xk71ImIpKEnnsOGjaE445znSS8MjJsBO7ll21EToIXyxGmQcBRQI+9fcJtt91GXl7ejtuKFStil05ERJz7+WfbAXbFFWolUJSePW0kbvZs10mSQ6BrmLbzPO8p4DSgg+/7K/f2eWlpaaSlpcUikoiIhNCwYdbF+pJLXCcJv06doFYtePVVaKmVwYELuq2A53neIOBM4Hjf978N8noiIhK/fN+m4844A2rUcJ0m/FJTbS3Tm2/qqJRYCHpKbjBwIXA+sN7zvFp/3NRVQ0REdjF1KixcaNNxUjznnQerVsHkya6TJL6gC6arsZ1xnwA/7nQ7N+DriohInBk6FOrXhxNOcJ0kfuTkQGYmvPGG6ySJL9CCyfd9by+3l4K8roiIxJdff7WppSuugDI6tKvYypSBc8+1hfJbt7pOk9j0ZSkiIs69/jps2aLF3qVx7rnw00/wySeukyQ2FUwiIuLcSy9Bly5Qu7brJPGnZUs7LmX4cNdJEpsKJhERcWrhQpgxQ6NLpeV5tvh75EgbpZNgqGASERGnXn4ZqlWzESYpne7d4ZdfNC0XJBVMIiLizNat1qzy/PPtbDQpnWbNbIfhqFGukyQuFUwiIuLMhx/C6tWajttfngfdusE776iJZVBUMImIiDMvvwxHHAEtWrhOEv/OPNOKz2nTXCdJTKEsmAYPHkxWVhaRSMR1FBERCcgvv9iIyCWX6KDdaDjmGDtbbuRI10kSk+f7vqtrF3nh/Px8MjIyyMvLIz09PRaZREQkRp55Bq65BlautF/0sv9694bx42HZMhWhJVCs/1OhHGESEZHE99JLcPLJKpai6cwz4ZtvYP5810kSjwomERGJucWL4YsvtNg72jp2hIwM7ZYLggomERGJuddes1/sp57qOkliKVcOunbVOqYgqGASEZGY8n07O+6ss6B8eddpEk+3brBgga1jkuhRwSQiIjH1xRe2zub8810nSUydO0OFCpqWizYVTCIiElOvv24LvTt2dJ0kMVWqZEXTO++4TpJYVDCJiEjMbN0Kb7xhh8WmpLhOk7i6drUGlmvXuk6SOFQwiYhIzEycCP/9r6bjgvbXv9oRKePGuU6SOFQwiYhIzLz+OjRsCK1auU6S2GrVgkgE3n3XdZLEoYJJRERi4rff4O234YIL1IU6Fk491UaYfv/ddZLEoIJJRERi4v33Yf166NHDdZLkcOqpkJcHU6a4TpIYQlkw6fBdEZHE8/rr0LIlNG7sOklyaNEC6tTRtFy06PBdEREJ3K+/2rqaAQPg+utdp0keV14Jn35qR9HIXunwXRERCYdRo2DLFjj3XNdJkkuXLrBkCfznP66TxD8VTCIiErjXX7dGlXXruk6SXE44AdLS4L33XCeJfyqYREQkUD/+CB9/rN5LLhxwABx3nNYxRYMKJhERCdSbb1pX77POcp0kOZ16KkyaZDvmpPRUMImISKDefNPONjvwQNdJklOXLnYkzYQJrpPENxVMIiISmJUr4fPPoXt310mS1yGHwGGHwQcfuE4S31QwiYhIYN5+G8qWhdNOc50kuXXubAWTu05C8U8Fk4iIBOatt+Ckk6BKFddJklvnzvD997B0qesk8UsFk4iIBGLlSjuWQ9Nx7nXsCOXKaVpuf6hgEhGRQGyfjjv9dNdJpFIlaNdOBdP+UMEkIiKB0HRcuHTuDJ98Aps3u04Sn0JZMOnwXRGR+PbDD5qOC5uTT4ZNm2DyZNdJ4pMO3xURkah74gm46SZYs0YjTGHh+3Y0zQUXwCOPuE4TKjp8V0RE3HjrLejUScVSmHieTZFqHVPpqGASEZGo2j4dd845rpPI7jp3hgULYNUq10nijwomERGJKu2OC69OnWyk6cMPXSeJPyqYREQkqjQdF17Vq0PLlpqWKw0VTCIiEjU//GC7sLQ7Lrw6d4bx42HbNtdJ4osKJhERiZqRIyE1VdNxYda5M6xbB3PmuE4SX1QwiYhI1IwaBSecAAce6DqJ7E3r1tb5e8IE10niiwomERGJirVrYdIk6NbNdRLZl3LloEMHFUwlpYJJRESiYswYa46o6bjwO/54W2tWUOA6SfxQwSQiIlExahS0aQO1arlOIkU54QT47TeYNs11kvihgklERPbb+vXW20fTcfGhWTOoWlXTciWR6uKinud5eXl5f/p4QUEBBTuND65fvx6wM+VERCS8Ro2y6Z0TTwT9yI4P7dpZP6Ybb3SdxK2MjIx0YL1fxOG6Tg7f9TwvHfhzxSQiIiISexm+7++z1HdVMHl5eXmFu3989xGmH3/8kezsbBYtWkTdunUDzxWJRJgxY4auU0L5+flkZmayYsUK0tPTA71WrP5NsbxWLK4Ty9cIEuv/XSyvE6/fSwUF0KAB9O0Lt90W7LX2Rd9LJbNsmXX9fustuOOOxHmNoGSvU0ZGRgbFGGFyMiVXVKjdVa5cOSZfmCkpKbrOfkhPTw/8erH8NyXi6xSL1wgS7/+dvpf2bexY2LABzj8f9vSUifg6JcL3UosWULcuTJ2amK8RFO91KmpkaTst+t5Jnz59dJ2Qi+W/Sa9T6SXa/zu9Rvs2ahQceigccUTw19oXvU4l43m2W27CBL1GxeFkSu4PRV545cqVO4bUDjrooFhkklLIz88nIyODvLy8mL5zkOLTaxQf4vF12rYNateGiy+GRx5xnSZ48fga7csrr9hr99NPdjBvoijh6+QV5zlDPcKUlpa2y72EU1paGnfffbdepxDTaxQf4vF1+vxz+2V75pmuk8RGPL5G+3L88XY/caLbHNEWxOsU6hGmRKvkRUQSzfXXw/DhsHIllAn1W3DZm8aN4bjjYMgQ10mcif8RJhERCS/fh5Ej7SgUFUvxa/s6Jtk3fYmLiEipzJ0L33+fPNNxieqEE6zFwPLlrpOEmwomEREplVGjoEoV6NjRdRLZHx072o65jz92nSTcVDBJqQwcOJBIJELlypWpUaMGZ5xxBkuWLHEdS4owcOBAPM+jX79+rqPITn744QcuvPBCqlWrRsWKFWnevDmzZs1yHatII0fCqadC2bKukwRv69at3HnnndSvX58KFSrQoEED7rvvPgoL/9SDOe5Uq2Zny8Xjwu9PP/2Url27UqdOHTzP45133tnl733f55577qFOnTpUqFCBjh07snDhwlJdSwWTlMqkSZPo06cP06ZNY/z48WzdupWTTjqJjRs3uo4mezFjxgyGDh3KUUcd5TqK7OSXX36hbdu2lC1blrFjx7Jo0SIee+wxqlSp4jraPv3nP7BwYfJMxz300EMMGTKEQYMG8dVXX/Hwww/zyCOP8NRTT7mOFhXHHguTJrlOUXIbN26kWbNmDBo0aI9///DDD/P4448zaNAgZsyYQa1atejUqdOOs2pLxPd9V7ci5eXl+YCfl5dXnE8Xh9asWeMD/qRJk1xHkT1Yv36936hRI3/8+PH+scce61933XWuI8kfbrnlFr9du3auY5TYQw/5foUKvr9xo+sksdGlSxe/V69eu3zszDPP9C+88EJHiaJr1CjfB9//7jvXSUoP8EeNGrXjz4WFhX6tWrX8Bx98cMfHNm/e7GdkZPhDhgzZ5aHFuYVyhGnw4MFkZWURiURcR5Fiysuzs5SrVq3qOInsSZ8+fejSpQsnnnii6yiym9zcXFq1akX37t2pUaMGLVq04LnnnnMdq0gjR0LnzlCxousksdGuXTsmTJjA0qVLAZg3bx6TJ0/mr3/9q+Nk0dG+vd3H4yjT3nz77besXr2ak046acfH0tLSOPbYY/n8889L/HyhLJj69OnDokWLYnbIquwf3/e5/vrradeuHUfs7WwEcWb48OHMnj2bgQMHuo4ie/DNN9/wzDPP0KhRIz744AN69+7NtddeyyuvvOI62l6tWgXTp0O3bq6TxM4tt9xCjx49OPzwwylbtiwtWrSgX79+9OjRw3W0qKhWDY48MrEKptWrVwNQs2bNXT5es2bNHX9XEk4O35XE0rdvX+bPn8/kyZNdR5HdrFixguuuu44PP/yQ8uXLu44je1BYWEirVq0YMGAAAC1atGDhwoU888wz9OzZ03G6PRszBlJSbMF3snjjjTd49dVXef3112natClz586lX79+1KlTh4svvth1vKg49lg7SDnReN6ufSl93//Tx4ojlCNMEj+uueYacnNzmThxos77C6FZs2axZs0aWrZsSWpqKqmpqUyaNIknn3yS1NRUtm3b5jpi0qtduzZZWVm7fKxJkyYsD3FTnNxcaNcOkmkG/qabbuLWW2/lvPPO48gjj+Siiy6if//+CTVye+yx8PXX1rU9EdSqVQvgT6NJa9as+dOoU3GoYJJS8X2fvn37MnLkSD7++GPq16/vOpLswQknnMCCBQuYO3fujlurVq244IILmDt3LikpKa4jJr22bdv+qSXH0qVLOfjggx0l2rcNG6wr9GmnuU4SW5s2baLMbu3MU1JSEqKtwHYdOth9okzL1a9fn1q1ajF+/PgdH9uyZQuTJk2iTZs2JX4+TclJqfTp04fXX3+d0aNHU7ly5R0VfEZGBhUqVHCcTrarXLnyn9aVVapUiWrVqmm9WUj079+fNm3aMGDAAM455xy++OILhg4dytChQ11H26Px46GgALp2dZ0ktrp27coDDzxAvXr1aNq0KXPmzOHxxx+nV69erqNFTY0akJVlBdMFF7hOUzwbNmxg2bJlO/787bffMnfuXKpWrUq9evXo168fAwYMoFGjRjRq1IgBAwZQsWJFzj///JJfrLjb6QK4FUltBcILOzz5T7cXX3zRdTQpgtoKhM+YMWP8I444wk9LS/MPP/xwf+jQoa4j7dUll/h+kyauU8Refn6+f9111/n16tXzy5cv7zdo0MC/4447/IKCAtfRourqq33/sMNcpyi+iRMn7vF30cUXX+z7vrUWuPvuu/1atWr5aWlpfocOHfwFCxbs/jTFqls83/f3s74rtSIvnJ+fT0ZGBnl5eaSnp8cik4iI7MW2bVCrFlx2GTz4oOs0EoQ33oDzzrOdkLVru04TM8VaAa41TCIiUizTpsHatcm3fimZHHus3X/6qdscYaSCSUREiiU3F/7yF2jd2nUSCUqtWtC4MXzyiesk4aOCSUREiiU313ovaXNlYovXc+WCpoJJRESKtHQpLF6s6bhkcOyx8NVXsGaN6yThooJJRESKNGYMpKVBp06uk0jQtI5pz0JZMOnwXRGRcBkzBk48ESpVcp1Egla3Lhx6qAqm3amtgIiI7NO6dVCzJjz9NFx5pes0EguXXgrz5sHs2a6TxITaCoiIyP4bO9Z6MCXTYbvJrl07K5jy810nCQ8VTCIisk+5uRCJQJ06rpNIrLRvD4WFMHWq6yThoYJJRET2qqAAxo3T7rhk06iR9dyaPNl1kvBQwSQiIns1aRKsX6+CKdl4nk3LffaZ6yThoYJJRET2KjcXDj4YjjzSdRKJtXbtYPp02LLFdZJwUMEkIiJ75PtWMJ12mo04SHJp3x42b06anXJFUsEkIiJ7NG8erFih6bhk1bw5VKyoabntVDCJiMge5eZCejp06OA6ibhQtizk5Gjh93YqmEREZI9yc+GUU6BcOddJxJX27WHKFGsxkOxUMImIyJ+sXAmzZmk6Ltm1a2ed3hcvdp3EPRVMIiLyJ+++CykpNsIkyat1a/s60LRcSAsmHb4rIuJWbq6tXTrwQNdJxKXKlW3xtwomHb4rIiK72bABqlWDhx6Cfv1cpxHX+veHd96Bb791nSQwOnxXRERK7sMPrVlh166uk0gYtGsH331n69qSmQomERHZRW4uNG0Khx7qOomEQdu2dj9litscrqlgEhGRHbZtswXf2h0n29WqBQ0bqoGlCiYREdlh6lTbRq6CSXbWrp0WfqtgEhGRHXJzoUYNyM52nUTCpE0bWLAA1q93ncQdFUwiIrJDbq4t9i6j3w6ykzZtrNv3jBmuk7ijbwkREQFgyRK7aTpOdtekCWRkwOefu07iTmAFk+d5h3ie97zned96nveb53lfe553r+d5OpVIRCSExoyB8uXhxBNdJ5GwKVMGjjlGBVNQDv/j+a8CmgL9gd7AgACvKSIipZSbC506QcWKrpNIGLVpY5sCkvUg3sAKJt/3x/m+f6nv+x/6vv+N7/u5wKPAmUFdU0RESmftWuuzo+k42Zs2beDXX5P3IN5Yr2HKAH6O8TVFRKQI779vIwennuo6iYRVdrZNzSXrtFzMCibP8w4FrgGG7O1zCgoKyM/P3+UmIiLBy821k+lr1XKdRMKqcmU48kiblktGJS6YPM+7x/M8v4hbq90eUwcYB7zl+/6/9vbcAwcOJCMjY8ctMzOz5P8iEREpkc2bYdw4TcdJ0dq0Sd4RJs/3/ZI9wPOqA9WL+LTvfN/f/Mfn1wEmAtOBS3zf375c7E8XLigooKCgYMef8/PzyczMJC8vj/T09BLlFBGR4hk3Dk45xRoTHnGE6zQSZq++ChddZGveqlVznSZqvOJ8UmpJn9X3/bXA2mIl8Ly6WLE0C7h0p2Jpj9LS0khLSytpJBER2Q+5uVC/vh24K7IvbdrY/bRp0KWL2yyxVuKCqbj+GFn6BFgO3Aj8xfOsiPN9f3VQ15XYKCiAlSvhxx9tOH/LFruBNTfLyIA6daBmTfCKVbuLiAu+bwXT2Wfre1WKVr++/Vz//HMVTNF0EtDwj9vK3f5O35ZxorAQFi60dvhz5sDcudYJ+Kefivf4Aw6wU66bNLHDGzt0gKwsHbsgEhZz5sAPP2j9khSP5/2vH1OyCaxg8n3/JeCloJ5fgrNuHYweDWPHwief2Fy150HjxtC8uTW2q1cPDjrIRpEqVoRy5exWWAj5+fDLLzYCtWyZ3ebNg7fegq1boWpVOOEEOOcc+Otf1SRPxKXcXBsRbt/edRKJF8ccA/fcYz/PU4McdgmZJPqnyr6sWwcjRtht4kQrfHJyoHdvOP54679RqVLxnqtGDbuPRHb9+MaNNu/96afw3nvQvbs9Z7du0LevbWkWkdjKzbU3LmXLuk4i8aJNG9i0CebPh6OPdp0mdkq8Sy6Kirxwfn4+GRkZ2iUXEN+3AuaZZ+DNN+3dwnHHwVlnwRlnBN+PZdkyeOMNeOEF+OYbK8r69bORp5SUYK8tIrBihY0W//vfcN55rtNIvNi8GdLT4fHH7c1uAijWMiGtJElCGzfCs89Cixb2TmHKFLjvPlvAPX68jSrFonldw4Zwxx2wdKkd+pmeDuefb9uahw9P3vOKRGJlzBibUjn5ZNdJJJ6ULw8tWyZfPyYVTEkkLw8GDIBDDoH/+z+7HzcO/vMfuPlm+Mtf3ORKSbHjGMaPt8XlDRpAjx42Tz5zpptMIskgNxeOPRaqVHGdROJNMjawVMGUBH79Fe66Cw4+2EaSuneHr7+Gd96Bzp3DtWOtVStb3zRpkg37ZmfDlVfawnMRiZ716229onbHSWm0aQPffw+rVrlOEjsh+lUp0VZQYHPMhx4Kjz0Gl11ma4WeftpGl8KsQweYNQuefNJ212VlWYEnItHxwQfWO61rV9dJJB4dc4zdJ1N7ARVMCaiwEIYNszYAN99sI0rLllnRVKeO63TFl5pqCwoXL7Zvzm7drOhbv951MpH4l5trB6nWr+86icSjOnVs1iKZpuVCWTANHjyYrKwsIrvvS5ciTZ9u2/N79rRFeV9+CUOGxFehtLuaNW106fnnbTdfs2Ywe7brVCLxa+tWm/rWdJzsj2RbxxTKgqlPnz4sWrSIGTNmuI4SN376CS6/3HonbdsGkyfD22/D4Ye7ThYdnge9elkDzKpVoW1beOUV16lE4tOUKfDzz3D66a6TSDzLybFO8QUFrpPERigLJim+bdtsTdJhh1mBNHiw7TRr29Z1smA0aGDFYI8ecPHFNmW3/Qw7ESme3FyoXdtGoUVKKyfHiqV581wniQ0VTHFs4UIrjPr2tYMzly61dgGJ3vSxfHmbnnvmGRg61FoS5Oe7TiUSH3zfjj7q2jVcO2Ql/jRrZkdiTZ/uOkls6NslDhUUwN13W+PJvDz47DN47jl3fZRc8DxrsPnBB/DFF7arLpm2t4qU1uLF1lZE03Gyv9LS7GiUadNcJ4kNFUxxZupU+wIdMABuvRXmzk3c6bfiOO44KxjXrbOddF995TqRSLiNHm0HXh9/vOskkghat1bBJCGzaZOds9a2LRxwgO0Su+8+q/CT3ZFHWiGZnm5di5NlPl2kNHJzrWFt+fKuk0giyMmx/n4//eQ6SfBUMMWBGTNsVOnZZ+HRR20b55FHuk4VLgcdBJ98ApmZNuqkDZYif/bf/9pogNoJSLTk5Nh9MqxjUsEUYr//DvfcY1NNlSvb9s3rr0/8Rd2lVa0aTJhgrRQ6dVKvJpHdvfeerf/r0sV1EkkUBx8MNWqoYBKHFi+2pmD332/nwH3+eeL0VApSlSq2ELxxY5t2WLTIdSKR8Bg92n6uJNMGEQmW59koUzKsY1LBFDKFhXZ+WosWdgTI1Km2I65sWdfJ4kflyjB2rPWZ6dTJ5tdFkt2mTTB+vKbjJPpat7bdyoWFrpMESwVTiKxYASedBNddB1dcYVNKOh2mdKpWtV8OlSrBCSfAypWuE4m4NWEC/Pab2glI9OXkWC+8xYtdJwmWCqYQ8H149VVbyL14MXz4oY0yVazoOll8q1kTPvrI3vWcdBL88ovrRCLujB5tU9WHHeY6iSSaSMSm5hJ9Wi6UBVMyHb67bh2ccw5cdJF1rF6wwKaRJDrq1bMCdPVqOOssHaMiyamwEMaM0XScBKNyZWjaNPELJs/3fVfXLvLC+fn5ZGRkkJeXR3p6eiwyxdT778Nll9kv8SFDoHt314kS16efWiHaowe8+KK9GxJJFtOm2W7bzz6Ddu1cp5FEdMUVto4pTvvgFes3QihHmBLdhg1w1VW2tbdFC/jySxVLQevQwQqll1+2nYciySQ3F6pXt6JJJAg5Ofa7bMMG10mCk+o6QLKZPBkuuQR+/NFGla68UqMdsXL++fDtt3DnndCwoY02iSSD0aNtyl893CQorVvb1O/MmdCxo+s0wdAIU4xs3gw332wjHTVr2rDlVVepWIq122+39WKXXRa3Q8ciJbJsmfUj0/olCVKTJraWKZHXMalgioHZs6FVK3jiCXjwQVtP07Ch61TJyfPsiJnDD4du3eDnn10nEgnWmDF25qQ2k0iQUlJst5wKJimV33+He++1ocpy5WDWLBtl0rC4WxUqwMiR1jekRw/Yts11IpHg5ObCiSfaod0iQcrJsSNS3O0lC5YKpoAsWmQLLP/+d7jtNqu6jzjCdSrZ7pBDYPhw69N0112u04gE4+efbWecpuMkFnJyrIXL8uWukwRDBVOUbdsGjz0GRx8NGzfa0Sb33WcjTBIuJ55oU6QDB8Lbb7tOIxJ9779vP5NOPdV1EkkGrVvbfaJOy6lgiqKFC6F9e7jpJujTR0ebxIMbb4Szz4ZevWwHnUgiyc2F7GyoU8d1EkkGNWpA/fo2LZeIVDBFQUGBHZDbooUNgU+aZKNMFSq4TiZF8Tz417/s7LkePWzdmUgiKCiAceM0HSex1bq1RphkLyZPhubNbVrn1lth7lwbZZL4kZFh65lmzdJ6Jkkcn3wC69erYJLYysmx2ZVEPIZKBVMp5eXB1VdbcVSlin2B3HcflC/vOpmURuvW8MAD8NBDdvacSLzLzbXNDdpsIrGUk2Ojm3Pnuk4SfaEsmMJ8+G5hIbzwgp34/eqr8NRTNsqkH0rx78Yb4aSTrLHl6tWu04iUnu9bwXT66WqOK7HVvDmULWvnyiUaHb5bAtOmwTXXWOv3Cy6w0Yi6dV2nkmj673+hWTO7jR0LZUL5lkJk3+bMsZ26EybA8ce7TiPJJjvbmgO/8orrJMWmw3ej5ccf4eKLra9SYaGNKL36qoqlRFSzph3Q++GHdtafSDwaPdrW5mk9pbiQnQ0zZrhOEX0qmPbh11/toNZGjayfydChNszYtq3rZBKkzp2hd29rD7Fsmes0IiU3apT1Xipb1nUSSUbZ2bB4sa31TSQqmPZg40ZraFi/Pjz+OPTtC0uXwhVX6FiTZPHIIzbadOmlOjpF4ss338D8+XZWoogL25cfz5zpNke0qWDayebNMGgQHHoo/O1vcOGF8PXXVjwdeKDrdBJLBxwAL70EU6bAP//pOo1I8Y0aZbt1Tz7ZdRJJVo0bQ+XKibfwWwUT8MsvMGAAHHwwXHcdnHKKjSg99RTUru06nbjSoQP07w933GFnA4rEg5EjbbdnpUquk0iyKlPGRplUMCWQRYusQKpXz3oodesGS5bAiy9a/xKR+++3qdmLL4atW12nEdm31avt/EpNx4lr2dkqmOLe5s22w619e2ja1Do8X3cdfP+97Ypq2NB1QgmTChVs19zs2TY1KxJmo0fbu/uuXV0nkWSXnQ2rVsEPP7hOEj0xKZg8z0vzPG+u53m+53nNY3HNnW3bBhMn2s6nunWtMWG5cvDmm7BihY0i1KwZ61QSL7Kz4bbb4N57YcEC12lE9m7UKJtKrlbNdRJJdtnZdp9I7QViNcL0MLAqRtcCYN06K4guvxwOOsiat33wge10W7LEGrp1726Fk0hR7rrL2ktccYV2zUk45eXBxx9rOk7CoW5dWwOcSNNyqUFfwPO8U4CTgLOAU4K6zurVMH263SZMsKrW9yErC84/H845xypeHRMgpZGWZn242reHZ56xVhMiYfLee/D773DGGa6TiJhEW8cUaMHkeV5N4DngDGBTSR+/ebPdr1xpI0H5+XZ0xfbb11/DV1/Zbd06+9yaNaFjR5t+69TJRpdEoqFdO/u6uu02O6MrM9N1IpH/GTUKWrXS16WER3a2HSFWWJgYx0wFVjB5nucBLwFDfN+f6XneIUU9pqCggIKCgh1/fvRRq5iaNv3z51apYruXDj/cttBmZdk2xsxMjSJJcB580BbW9ulj9/pakzD47Tc7+/D2210nEfmf7Gwb6Fi61H5Xx7sSF0ye590D3F3Ep0WANkA6MLC4zz1w4EDuvffenT5yKAAjRkDVqtYIq2ZNqFHDpkhEYi0jw5qbnnUWvP02nH2260Qi8NFHdkKB1i9JmLRqZfczZiRGweT5vl+yB3hedaB6EZ/2HTAc6ArsfIEUYBvwmu/7PXd/0O4jTPn5+WRmZpKXl0d6enqJcooEqVs3mDbNenmpC7y41qsXfP65nd8lEiaHH27LY556ynWSfSrWXEGJC6ZiX93z6mEjTNvVAT4Azgam+76/oqjnyM/PJyMjQwWThM4PP0CTJtCjBzz7rOs0ksy2boVatWwH58Bij+eLxEbPnrYzffp010n2qVgFU2DLsHzfX+77/pfbb8DSP/7qa9/3VwZ1XZFYqFvXfjk991zofxBIgps82Ta9aDpOwigSgblzYafJo7iVAOvWRdzo3RtatID/+z/1ZhJ3Ro2y3cDb14uIhEl2NmzZAvPnu06y/2JWMPm+/53v+57v+3NjdU2RIKWkwNNP27EpmpYTF3wf3nnHei8lwrZtSTzNmkHZsonRj0nfYiL7oXVr6yZ/xx2wZo3rNJJsZs+G5cs1HSfhVb68FU2JcESKCiaR/TRwoL27v+UW10kk2YwaZS1XOnRwnURk7xKl47cKJpH9VL26FU0vvQRTprhOI8lk1Cjo2hVSAz/kSqT0srOt5UVenusk+0cFk0gUXHaZ7Qb5v/+zbd4iQVuyxPqAaTpOwi4SsfV2s2a5TrJ/VDCJRMH2BeALFtjhvCJBGzECDjjAjoYSCbPGje2kjnifllPBJBIlrVrZAvC//Q3WrnWdRhLdiBFw6qlQoYLrJCL7lpJiPx9VMInIDvffb0PPf/ub6ySSyJYts2aAOstQ4kV2dvzvlAtlwTR48GCysrKIRCKuo4iUSI0acPfd1pdp3jzXaSRRjRgBFSvCKae4TiJSPNnZsHIlrFrlOknpBXaWXDEUeWGdJSfx6Pff4aijoGZNmDgRvGKdUiRSfC1bwqGHwptvuk4iUjwrV0JmpjVaPf1012n+xO1ZciLJqmxZ+Oc/YdIkGwkQiaZvvrGGlZqOk3hSty7Urh3f65hUMIkEoHNn649z442waZPrNJJI3n7bFnr/9a+uk4gUn+dZewEVTCLyJ48/Dj/+CI895jqJJJK33rK1Swcc4DqJSMlsX/hdWOg6SemoYBIJSMOGcO218NBDVjiJ7K/vv7dfOJqOk3iUnW3dvpctc52kdFQwiQTozjvt8Mm77nKdRBLBiBGQlmb9l0TiTatWdh+v03IqmEQCVKWKtRl44QW1GZD9N2IEnHyydU0WiTcHHgiHHaaCSUT2ondvaNQIbrjBmlqKlMaKFTBtGnTv7jqJSOllZ6tgEpG9KFsWHn0UJkyA9993nUbi1dtvQ7lymo6T+BaJwJw5sGWL6yQlp4JJJAZOPRWOP97aDPz+u+s0Eo9GjLCDdjMyXCcRKb3sbCuW5s93naTkVDCJxIDnWXuBJUtg6FDXaSTerFgBU6ZoOk7iX/PmkJoan+fKqWASiZHmzeHSS20R+K+/uk4j8eTNN2133BlnuE4isn/Kl4dmzWD6dNdJSi6UBZMO35VE9fe/w2+/wYABrpNIPBk+HLp0AR2pKYkgEtEIU9T06dOHRYsWMSMe/4+K7EOdOnDLLfDEE3YmmEhR/vMfmDkTzjvPdRKR6IhE4KuvYP1610lKJpQFk0giu/FGqF4d7rjDdRKJB2+8YcegdOniOolIdGRnW4uV2bNdJykZFUwiMVaxItx3n02zzJzpOo2E3fDhcPrp9nUjkgiaNIFKleKvH5MKJhEHLr4YsrLg5pvVzFL27ssvYeFCTcdJYklJgaOPjr91TCqYRBxITbVDeSdOhHHjXKeRsBo+3I7XOekk10lEoiseF36rYBJxpEsX6NDBFoFv2+Y6jYSN71vBdNZZ1uFbJJFkZ8N338FPP7lOUnwqmEQc8Tx45BFYsACGDXOdRsJm1iz4+mtNx0li2t41KJ5GmVQwiTiUnW3dm++6y/oziWw3fDjUqAEdO7pOIhJ99etDtWoqmESkBAYMgNWr4cknXSeRsCgstHYC55xj691EEo3nxd86JhVMIo41bAi9e8PAgbBunes0EgZTpsDKlZqOk8QWiVhrgXjZKayCSSQE7rrLFn7ryBQBm47LzIRjjnGdRCQ4kYgt+l6+3HWS4lHBJBICNWrYbrlBg2zniCSvLVusYOrRA8roJ7QksHhb+B3Kb0cdvivJqH9/WwR5552uk4hLY8fCzz/DRRe5TiISrFq1bCQ1Xjp+e767ycMiL5yfn09GRgZ5eXmk65huSQJDh8JVV9mW8qOPdp1GXDj7bGsnMGeO6yQiwTvrLHuDMHGi0xhecT4plCNMIsmqVy84/HCbnpPk8+uvMGaMRpckeWRn2xvEeGjeq4JJJERSU+HBB+Gjj+DDD12nkVh76y3YutXWL4kkg0gE1q+HJUtcJymaCiaRkDntNGjb1g7mLSx0nUZiadgwOPFEqF3bdRKR2GjZ0u7jYeG3CiaRkNl+ZMq8efDaa67TSKx89x189pmm4yS5ZGRA48YqmESklI45Bs4803bMbd7sOo3EwquvQqVK0K2b6yQisZWdHR875VQwiYTUwIHwww/Wm0kSm+/bdNyZZ1rRJJJMIhEbUd+yxXWSfVPBJBJShx0GV14JDzxg224lcc2YAUuXajpOklMkYsXS/Pmuk+xb4AWT53ldPM+b7nneb57nrfU8b2TQ1xRJFHffDb//bqNNkriGDYM6deD4410nEYm95s1th3DY1zEFWjB5nncWMAx4EWgGtAVeD/KaIomkZk246SZ46in4/nvXaSQIv/9uR6Gcfz6kpLhOIxJ75cvDUUeFfx1TYAWT53mpwBPATb7vD/F9f6nv+0t83x8R1DVFEtENN0CVKnZArySe99+HtWuhZ0/XSUTciUSSe4TpaKAuUOh53hzP8370PG+s53lNA7ymSMI54AC4917bRaXjMhLP889Dq1Zw5JGuk4i4k50NixZZE8uwCrJgavDH/T3A/cCpwC/AJM/zqu7pAQUFBeTn5+9yExG47DJbBK4jUxLLjz/aCFOvXq6TiLgVidhu0dmzXSfZuxIXTJ7n3eN5nl/ErdVOz/2A7/tv+74/C7gUO3S3+56ee+DAgWRkZOy4ZWZmlvofJpJIUlPhoYdg/HgdmZJIhg2DsmV1FIpIkyZQsWK4p+U83/dL9gDPqw5UL+LTvgOOAT4G2vu+P3mnx08HPvJ9//bdH1RQUEBBQcGOP+fn55OZmUleXh7p6eklyimSaHwfOnSwIetZs7RAON75vv2SaNXKpltFkl2HDnYs0BtvxPzSXnE+KbWkz+r7/lpgbZFX97xZQAHQGJj8x8fKAocAe9zvk5aWRlpaWkkjiSSF7UemHHOMHZmiRcLxbepUO3D06addJxEJh0gERoa48VBga5h8388HhgD3ep53kud5jYFn/vjrt4K6rkgiy8mBs8+2I1N++811Gtkfzz8PhxwCHTu6TiISDpGInan400+uk+xZ0I0rbwKGY72YZgAHA8f7vv9LwNcVSVgDBthi4aeecp1ESmvDBpt2uPRSKKPzFkQAK5gAZs50m2NvAv1W9X3/d9/3b/R9v6bv++m+73fyfX9hkNcUSXSNGkHv3lY4rVvnOo2UxltvwaZNcMklrpOIhEeDBlC1angbWOq9jUgcuusu2LbNzpmT+PP889CpE9Sr5zqJSHh4XrgbWKpgEolDNWpYT6ZBg+Dbb12nkZJYuBCmTLHeWiKyq+0FUwk38MeECiaRONW/P1SvDnfc4TqJlMSzz9oZgWec4TqJSPhkZ8OaNbB8ueskf6aCSSROVaoE990H//53eBdJyq42bYJXXrHO3uXKuU4jEj7bF36HcVpOBZNIHLvkEsjKgptvDucQtuzqjTcgPx+uuMJ1EpFwqlULDjpIBZOIRNn2I1MmToSxY12nkaIMGQKdO0P9+q6TiIRXdrYKJhEJQJcu1vzwxhth61bXaWRv5syx7dK9e7tOIhJukYgtMygsdJ1kV6EsmAYPHkxWVhaR7ZOZIrJXngePPw6LF8PQoa7TyN48+yzUrWsFrojsXSRiZ2YuWeI6ya5KfPhuFBV54fz8fDIyMnT4rkgx9OoFubmwbBlUqeI6jexs/XqoUwduuAHuucd1GpFwy8uzn2EvvxyzMzOLdfhuKEeYRKTk7r8fNm+2ewmXl1+2s/8uv9x1EpHwy8iAxo3D1/FbBZNIgqhTx5pZPvmkjTJJOBQW2rl/Z55pu39EpGhh7PitgkkkgdxwgzVFvPlm10lkuw8+gKVL4brrXCcRiR+RCMydC1u2uE7yPyqYRBJIxYrw4IMwahR88onrNAI24nf00dCmjeskIvEjO9uKpfnzXSf5HxVMIgmmRw/7YXP99XZAr7izeDGMG2ejS16xlpWKCEDz5tZnLkzTciqYRBJMmTLwj39Y359XXnGdJrk99ZQdlHzuua6TiMSX8uXhyCNVMIlIwNq0sV/St98OGza4TpOcfv3Vdsf17g1paa7TiMSf7Oxw7ZRTwSSSoB58EH75xY5Okdh74QVbg6HO3iKlE4nAV1+F502fCiaRBHXIIbaO6dFH4dtvXadJLtu2waBBcM45ULu26zQi8SkSsbYcs2e7TmJUMIkksNtvh2rVrHCS2Bk1yopUtRIQKb2sLNv5G5ZpORVMIgnsgAPgscfgnXdst5YEz/fh4YfhuOPsHbKIlE5qqrXkCMvC71AWTDp8VyR6zjkHOnaEa6+FggLXaRLfJ5/YD3g1DxXZf2Hq+K3Dd0WSwJdfWl+T+++HW291nSaxnXwyrF5tbR3Ue0lk/wwfbr3lfvoJqlcP7DI6fFdEzBFHwDXXwN//DitXuk6TuObNs6NQbr5ZxZJINGyfaArDKJMKJpEkcc89ULky3Hij6ySJ6+GHbXfiOee4TiKSGBo0gKpVVTCJSAxlZNgv9DfegIkTXadJPN99Z/9vb7jBFquKyP7zPGjVSgWTiMTYhRdaF/BrroHff3edJrE88ghUqQK9erlOIpJYtnf8drfk2qhgEkkiZcpYQ8VFi+ycM4mOlSvhX/+y0aWKFV2nEUkskQisWQMrVrjNoYJJJMm0aAF9+sDf/gbLl7tOkxgeesh6XvXt6zqJSOIJy8JvFUwiSeiBB2xNU58+7oe5492qVfDcc9ZNvXJl12lEEk/t2nDQQe47fqtgEklC6ek2JffuuzBypOs08e2hh6BCBVsXJiLBCEMDSxVMIkmqWzc47TT7RZ+X5zpNfPrxRxg6FPr3tyJURIIRicDMmXYYrysqmESSlOfZKFN+Ptxxh+s08envf7fRJR2yKxKs7GxYvx6WLHGXQQWTSBKrV8+OS3n6aZg82XWa+LJsma1duv12Ww8mIsFp2dLuXU7LhbJg0uG7IrFzzTWQk2P9gzZtcp0mftx1F9SsaQvnRSRYVarAYYe5LZh0+K6IsGSJHc579dXw+OOu04TfnDlw9NE2wnT55a7TiCSHiy6CpUth+vSoP7UO3xWR4mnc2Kbm/vlPmDLFdZrwu/12e7d7ySWuk4gkj0gE5s6FLVvcXF8Fk4gA0K+fTc1deqmm5vZlwgQYN856WenMOJHYiUSsWFqwwM31VTCJCAApKfDii9b9+667XKcJp61bbUdcu3Zw1lmu04gkl+bN7U2KqwaWKphEZIftU3P/+Iem5vZkyBA7h++JJ6wtg4jEToUKcOSR7hZ+q2ASkV30729TcxdfbH1PxKxbZ+fvXXaZLfgWkdhz2fFbBZOI7CIlBV55Bf77X7j2WtdpwuNvf4Nt22ztkoi4EYnYKO+GDbG/tgomEfmThg1h0CB46SV44w3XadxbsMCm4+6+G2rUcJ1GJHllZ9vxKLNnx/7agRZMnucd5nneaM/z1nqel+953hTP844L8poiEh09e8J558FVV8H337tO405hoTWnbNQI+vZ1nUYkuWVl2VomF9NyQY8wvQekAscDLYG5wLue59UK+Loisp88D555xjrsXnCB7RBLRs89B599ZiNM5cq5TiOS3FJTbQ1hQhVMnudVBxoCD/q+P9/3/f8AtwIVgaZBXVdEoqdKFXj1VZg2LTlbDfzwA9x8s3Xz7tjRdRoRAZuWc9FaIMgRpnXAV0BPz/MqeZ6XClwF/BeYFeB1RSSK2rWDgQPhwQchN9d1mti65hob/n/4YddJRGS7SAS+/RbWro3tdQPrU+v7vu95XidgNLAeKMSKpZN93/91T48pKCigoKBgx5/z8/ODiiciJXDjjfD557auadYsOPRQ14mC9/bbMGoUvPkmHHig6zQisl0kYvczZ8LJJ8fuuiUeYfI87x7P8/wibq08z/OAp4E1QHsgGyue3vU8r/aennvgwIFkZGTsuGVmZu7Pv01EosTzbMfcX/5iHa5/+811omCtXm2L3bt1g7PPdp1GRHZ26KH2JibW03Ke7/sle4CtTapexKd9B7QFPgQO9H1/x1CR53n/AZ73fX/g7g/a0whTZmYmeXl5pKenlyiniETf/PnW1PLMM2HYsMTsdu370KWLbVtesMCKRBEJl86dbRPGmDFRebpi/SQr8ZSc7/trgSJnDj3Pq/jHfxbu9leF7GVkKy0tjbS0tJJGEpEYOeooO2/uvPPsGJVEXAg+ZAiMHQvvvqtiSSSsIhH417/sDU6s3rgFueh7KvAL8LLnec3+6Mn0CFAfazcgInHo3HPhvvus8/Xw4a7TRNecOXD99dC7t40yiUg4RSJ2GsGKFbG7ZmAF0x8jUScDBwAfAzOBdsDpvu/PC+q6IhK8O++ECy+ESy6BqVNdp4mOX36x9UpZWXb4sIiEV3a23ceyH1OgjSt935/p+35n3/er+b6f7vv+Mb7vjw3ymiISPM+z4fBIBE4/Hb75xnWi/VNYaDsAf/4ZRoyA8uVdJxKRfaldG+rWTaCCSUQSV1qabbvPyIBOnWDVKteJSu/BB23N0quvQv36rtOISHFEIrHdKaeCSURKrXp1GD8eCgqsaIp1I7loGDnSphjvukvrlkTiSXa29YUr3H1rWUBUMInIfjnkEPjoI/jpJ2sil5fnOlHxTZ9u5+R17w733OM6jYiURCQC+fmwZElsrqeCSUT22+GHw4cfwtdfwwknxMdI09Kl0LUrtGwJL78MZfTTUCSuRCK2nnL69NhcTz8iRCQqmjeHiRNh+XI49thwr2n6+ms4/nibUnznHS3yFolHGRm2q3XatNhcTwWTiERN8+bw2Wc2TN6+vR2QGTbff2/FUsWKMGGCFU0iEp9at1bBJCJxqnFjmDzZprhat7YCKiwWLoR27SAlBT7+2LYmi0j8ysmxI4w2bgz+WqEsmAYPHkxWVhaR7UcSi0hcOfhga2jZtKmN5jz7rOtEVsS1a2eHdk6eDAcd5DqRiOyvnBzbJTdzZvDXCmXB1KdPHxYtWsSMWHakEpGoql7dFoJfdZUdNXLVVbBpk5ssr7wCJ54IzZrBp59CnTpucohIdGVlwQEHxGZaLpQFk4gkhrJlYdAgeO45K1patIhtZ95Nm+Dqq+Hii+H882HcOKhSJXbXF5FgpaTYbrlY7JRTwSQigbv8cjvYtnJlaNMGbrrJFoYHaepUW4T+0kswZAg8/7x2w4kkopwc+373/WCvo4JJRGLi8MPth9o998DgwbY4/KWXYOvW6F7nxx+hVy9o2xaqVoW5c2060POiex0RCYecHFi9GlasCPY6KphEJGbKloU77oDFi63twKWXwmGH2aLw/d3lsmoV3HILNGoEo0fDU0/Z4u7GjaOTXUTCqXVruw96HZMKJhGJuXr14M03bZouErF1RjVrQs+e8P77sGFD8Z5n3Tr497+tY3e9evDMM9C3LyxbBn36QGpqsP8OEXGvZk07oinodUyeH/Sk394VeeH8/HwyMjLIy8sjPT09FplExIHvvoPXXoNhw+xcqNRUW3/UtCk0aGALtStUsEXcP/9sI1QLF8KiRbZuIRKxabgePaz7r4gklx497JSBKVNK9fBiTdirYBKR0PB9K4YmTbLh9cWLrZjKz4fffrPu3BkZNo2XlQWtWtmBv2oTIJLcnnjCpuTz86FcuRI/XAWTiCQO39fCbRHZs2nT4JhjrG1Jq1YlfnixfrJoDZOIxAUVSyKyNy1a2MhSkAu/VTCJiIhIXEtLs6JJBZOIiIjIPrRunYQFkw7fFRERkZLIyYGvv4a1a4N5fi36FhERkbj37bfWhuTdd6FLlxI9VIu+RUREJDkccgjUqBHctJwKJhEREYl7nmfrmILq+K2CSURERBJCTo4VTIWF0X9uFUwiIiKSEHJyrNv34sXRf24VTCIiIpIQIhGbmgtiWk4Fk4iIiCSEypXt0O4gFn6rYBIREZGEkZOjgklERERkn3Jy4MsvYcOG6D6vCiYRERFJGK1b2y65mTOj+7wqmERERCRhNGlia5miPS2ngklEREQSRkoKXHstNGwY3efVWXIiIiKSzOL3LLnBgweTlZVFJBJxHUVEREREI0wiIiKS1OJ3hElEREQkTFQwiYiIiBRBBZOIiIhIEVyuYSqS53npQB6Q4ft+vus8IiIikpzCXjB5QGVgvR/moCIiIpLQQl0wiYiIiISB1jCJiIiIFEEFk4iIiEgRVDCJiIiIFEEFk4iIiEgRVDCJiIiIFEEFk4iIiEgRVDCJiIiIFOH/AZL2PVczqkWbAAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(x*sin(x)-4, (x, 0, 10))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can get it by forcing the search in the interval $[6,8]$:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 6.901412609021304\$" ], "text/latex": [ "$\\displaystyle 6.901412609021304$" ], "text/plain": [ "6.901412609021304" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "find_root(eq, 6, 8)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(6.901412609021304, \\begin{array}{l}\n", "\\verb| |\\verb|converged:|\\verb| |\\verb|True|\\\\\n", "\\verb| |\\verb|flag:|\\verb| |\\verb|'converged'|\\\\\n", "\\verb| |\\verb|function_calls:|\\verb| |\\verb|10|\\\\\n", "\\verb| |\\verb|iterations:|\\verb| |\\verb|9|\\\\\n", "\\verb| |\\verb|root:|\\verb| |\\verb|6.901412609021304|\n", "\\end{array}\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(6.901412609021304, \\begin{array}{l}\n", "\\verb| |\\verb|converged:|\\verb| |\\verb|True|\\\\\n", "\\verb| |\\verb|flag:|\\verb| |\\verb|'converged'|\\\\\n", "\\verb| |\\verb|function_calls:|\\verb| |\\verb|10|\\\\\n", "\\verb| |\\verb|iterations:|\\verb| |\\verb|9|\\\\\n", "\\verb| |\\verb|root:|\\verb| |\\verb|6.901412609021304|\n", "\\end{array}\\right)$" ], "text/plain": [ "(6.901412609021304,\n", " converged: True\n", " flag: 'converged'\n", " function_calls: 10\n", " iterations: 9\n", " root: 6.901412609021304)" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "find_root(eq, 6, 8, full_output=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximate solution to a system: `mpmath.findroot`" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle x \\sin\\left(x\\right) - y = 0\$" ], "text/latex": [ "$\\displaystyle x \\sin\\left(x\\right) - y = 0$" ], "text/plain": [ "x*sin(x) - y == 0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\$\\displaystyle y \\cos\\left(y\\right) - x - 1 = 0\$" ], "text/latex": [ "$\\displaystyle y \\cos\\left(y\\right) - x - 1 = 0$" ], "text/plain": [ "y*cos(y) - x - 1 == 0" ] }, "metadata": {}, "output_type": "display_data" } ], "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left( x, y \\right) \\ {\\mapsto} \\ x \\sin\\left(x\\right) - y\$" ], "text/latex": [ "$\\displaystyle \\left( x, y \\right) \\ {\\mapsto} \\ x \\sin\\left(x\\right) - y$" ], "text/plain": [ "(x, y) |--> x*sin(x) - y" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f1(x,y) = eq1.lhs()\n", "f1" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left( x, y \\right) \\ {\\mapsto} \\ y \\cos\\left(y\\right) - x - 1\$" ], "text/latex": [ "$\\displaystyle \\left( x, y \\right) \\ {\\mapsto} \\ y \\cos\\left(y\\right) - x - 1$" ], "text/plain": [ "(x, y) |--> y*cos(y) - x - 1" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f2(x,y) = eq2.lhs()\n", "f2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For solving a system, we use the [mpmath](http://mpmath.org/) toolbox:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "import mpmath" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\begin{array}{l}\n", "\\verb|[-0.642871006593746]|\\\\\n", "\\verb|[|\\verb| |\\verb|0.385398467857987]|\n", "\\end{array}\$" ], "text/latex": [ "$\\displaystyle \\begin{array}{l}\n", "\\verb|[-0.642871006593746]|\\\\\n", "\\verb|[|\\verb| |\\verb|0.385398467857987]|\n", "\\end{array}$" ], "text/plain": [ "matrix(\n", "[['-0.642871006593746'],\n", " ['0.385398467857987']])" ] }, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left[\\left[-0.642871006593746\\right], \\left[0.385398467857987\\right]\\right]\$" ], "text/latex": [ "$\\displaystyle \\left[\\left[-0.642871006593746\\right], \\left[0.385398467857987\\right]\\right]$" ], "text/plain": [ "[[mpf('-0.64287100659374617')], [mpf('0.38539846785798665')]]" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol.tolist()" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(-0.642871006593746, 0.385398467857987\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(-0.642871006593746, 0.385398467857987\\right)$" ], "text/plain": [ "(-0.642871006593746, 0.385398467857987)" ] }, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle -5.55111512312578 \\times 10^{-16}\$" ], "text/latex": [ "$\\displaystyle -5.55111512312578 \\times 10^{-16}$" ], "text/plain": [ "-5.55111512312578e-16" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f1(x1,y1)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle 1.11022302462516 \\times 10^{-16}\$" ], "text/latex": [ "$\\displaystyle 1.11022302462516 \\times 10^{-16}$" ], "text/plain": [ "1.11022302462516e-16" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f2(x1,y1)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(-0.6428708247411281,\\,0.38539856760460756\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(-0.6428708247411281,\\,0.38539856760460756\\right)$" ], "text/plain": [ "(-0.6428708247411281, 0.38539856760460756)" ] }, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(-0.6428708247411281, 0.38539856760460756\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(-0.6428708247411281, 0.38539856760460756\\right)$" ], "text/plain": [ "(-0.6428708247411281, 0.38539856760460756)" ] }, "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": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(-3.0233688164127415 \\times 10^{-07}, -1.0387405890988077 \\times 10^{-07}\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(-3.0233688164127415 \\times 10^{-07}, -1.0387405890988077 \\times 10^{-07}\\right)$" ], "text/plain": [ "(-3.0233688164127415e-07, -1.0387405890988077e-07)" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f1(x2,y2), f2(x2,y2)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\displaystyle \\left(-1.81852617853195 \\times 10^{-7}, -9.97466205743258 \\times 10^{-8}\\right)\$" ], "text/latex": [ "$\\displaystyle \\left(-1.81852617853195 \\times 10^{-7}, -9.97466205743258 \\times 10^{-8}\\right)$" ], "text/plain": [ "(-1.81852617853195e-7, -9.97466205743258e-8)" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x1 - x2, y1 - y2" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.6.beta4", "language": "sage", "name": "sagemath" }, "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.10" } }, "nbformat": 4, "nbformat_minor": 4 }