{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "We review the second half of the course." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 3. Calculus" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write a function to make polynomials in a system.\n", " \n", "The *k*-th polynomial in the system is\n", "\n", "$$\n", " f_k(x_1,x_2, \\ldots, x_n) = x_k + \\sum_{i=1}^{n-k} x_i x_{k+i},\n", " \\quad k=1,2,\\ldots,n.\n", "$$\n", "\n", "Test your function on 8 variables." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 1" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[x1, x2, x3, x4, x5, x6, x7, x8]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X = [var('x' + '%d' % k) for k in range(1,9)]\n", "X" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "f = lambda X,k,n: X[k-1] + sum([X[i]*X[k+i] for i in range(0,n-k)])" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 + x6*x7 + x7*x8 + x1,\n", " x1*x3 + x2*x4 + x3*x5 + x4*x6 + x5*x7 + x6*x8 + x2,\n", " x1*x4 + x2*x5 + x3*x6 + x4*x7 + x5*x8 + x3,\n", " x1*x5 + x2*x6 + x3*x7 + x4*x8 + x4,\n", " x1*x6 + x2*x7 + x3*x8 + x5,\n", " x1*x7 + x2*x8 + x6,\n", " x1*x8 + x7,\n", " x8]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[f(X, k, 8) for k in range(1, 9)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a piecewise function ``int_inv_cub`` which\n", "as function of the end point $b$ is defined as\n", "\n", "$$\n", " \\int_{-1}^b \\frac{1}{x^3} dx\n", "$$\n", "\n", "and which must be correct for *any* value of the parameter $b$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 2" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "b = var('b')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle \\int_{-1}^{b} \\frac{1}{x^{3}}\\,{d x}\\)" ], "text/latex": [ "$\\displaystyle \\int_{-1}^{b} \\frac{1}{x^{3}}\\,{d x}$" ], "text/plain": [ "integrate(x^(-3), x, -1, b)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "i = integral(1/x^3, (x,-1,b), hold=True)\n", "i.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Without ``hold=True`` an error message prompts ``Is b+1 positive, negative or zero?``\n", "\n", "Let us assume that $b+1<0$ or $b<-1$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle -\\frac{1}{2 \\, b^{2}} + \\frac{1}{2}\\)" ], "text/latex": [ "$\\displaystyle -\\frac{1}{2 \\, b^{2}} + \\frac{1}{2}$" ], "text/plain": [ "-1/2/b^2 + 1/2" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "assume(b+1<0)\n", "negval = integral(1/x^3, (x,-1,b)).unhold()\n", "show(negval)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We needed another assumption $b > 0$.\n", "We reset to test other cases." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle \\int_{-1}^{b} \\frac{1}{x^{3}}\\,{d x}\\)" ], "text/latex": [ "$\\displaystyle \\int_{-1}^{b} \\frac{1}{x^{3}}\\,{d x}$" ], "text/plain": [ "integrate(x^(-3), x, -1, b)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "reset()\n", "b = var('b')\n", "assume(b+1>0, b>0)\n", "infval = integral(1/x^3, (x,-1,b), hold=True)\n", "infval.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get another error message for the above case.\n", "The value is infinity if 0 is in $[-1,b]$.\n", "\n", "We needed another assumption $b > -1$ and $b < 0$.\n", "We reset to test other cases." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle -\\frac{1}{2 \\, b^{2}} + \\frac{1}{2}\\)" ], "text/latex": [ "$\\displaystyle -\\frac{1}{2 \\, b^{2}} + \\frac{1}{2}$" ], "text/plain": [ "-1/2/b^2 + 1/2" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "reset()\n", "b = var('b')\n", "assume(b+1>0,b<0)\n", "negval = integral(1/x^3, (x,-1,b))\n", "negval.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far, we have the following cases for $b$:\n", " \n", "1. $-\\infty < b < -1$ : $-1/2/b^2 + 1/2$\n", "2. $-1 < b < 0$ : $-1/2/b^2 + 1/2$\n", "3. $b > 0 : \\infty$ \n", " \n", "Two cases are left to consider: $b = -1$ and $b = 0$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integral(1/x^3, (x,-1,-1))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle \\int_{-1}^{0} \\frac{1}{x^{3}}\\,{d x}\\)" ], "text/latex": [ "$\\displaystyle \\int_{-1}^{0} \\frac{1}{x^{3}}\\,{d x}$" ], "text/plain": [ "integrate(x^(-3), x, -1, 0)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "integral(1/x^3, (x,-1,0), hold=True).show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above integral is divergent and we will assign it the value $\\infty$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, this is the complete case analysis:\n", "\n", "1. $-\\infty < b < -1$ : $-1/2/b^2 + 1/2$\n", "\n", "2. $b = -1$ : 0\n", "\n", "3. $-1 < b < 0$ : $-1/2/b^2 + 1/2$\n", "\n", "4. $b \\geq 0$ : $\\infty$ " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "reset()\n", "b = var('b')\n", "int_inv_cub = piecewise([((-oo,-1),-1/2/b^2 + 1/2), \\\n", " ([-1,-1],0),((-1,0),-1/2/b^2 + 1/2),\\\n", " ([0,0],oo),((0,oo),oo)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us check all cases." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3/8" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "int_inv_cub(-2)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "int_inv_cub(-1)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-3/2" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "int_inv_cub(-1/2)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "+Infinity" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "int_inv_cub(0)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "+Infinity" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "int_inv_cub(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The arc length of continuous function $f(x)$ over an interval $[a,b]$ is\n", "\n", "$$\n", " \\int_a^t \\sqrt{1+[f'(x)]^2}.\n", "$$\n", " \n", "1. Compute the arc length of the positive half\n", " of the unit circle, i.e.: $f(x) = \\sqrt{1-x^2}$\n", " \n", " (*answer:* $\\pi$).\n", "\n", "\n", "2. Create a function (call it ``arc_length``) in $t$\n", " which returns a 10-digit floating-point approximation\n", " of the arc length of the positive half of the circle,\n", " for $x \\in [0,t]$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 3" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(\\displaystyle -\\frac{x}{\\sqrt{-x^{2} + 1}}\\)" ], "text/latex": [ "$\\displaystyle -\\frac{x}{\\sqrt{-x^{2} + 1}}$" ], "text/plain": [ "-x/sqrt(-x^2 + 1)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = var('x')\n", "f = sqrt(1-x**2)\n", "df = diff(f,x)\n", "df.show()" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "pi" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integral(sqrt(1+df^2),(x,-1,1))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def arc_length(t):\n", " \"\"\"\n", " Returns the arc length of the circle for x in [0, t].\n", " \"\"\"\n", " i = integral(sqrt(1+df^2),(x,0,t), hold=True)\n", " return i.n(digits=10)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.570796313279659" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arc_length(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the recurrence relation\n", "\n", "$$\n", " h(n) = 5 h(n-1) - 6 h(n-2), \\mbox{ for } n \\geq 2,\n", " \\mbox{ with } h(0) = 1 \\mbox{ and } h(1) = -2.\n", "$$\n", "\n", "Answer the following the questions.\n", "\n", "1. The generating function\n", "\n", " $$\n", " g(x) = \\frac{1-7x}{1-5x + 6x^2}\n", " $$\n", " \n", " defines $h(n)$ as the coefficient with $x^n$ in the Taylor\n", " expansion of $g(x)$. \n", " \n", " Use $g(x)$ to define $h$ as a function (call it ``t``) of $n$ \n", " which gives the value of $h(n)$.\n", "\n", "\n", "2. Write a function to compute $h(n)$, directly using the\n", " recurrence relation from above. \n", " \n", " Make sure your function can compute $h(120)$. \n", " \n", " Compare with the result of the first part of the question." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 4" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-76172*x^9 - 24964*x^8 - 8108*x^7 - 2596*x^6 - 812*x^5 - 244*x^4 - 68*x^3 - 16*x^2 - 2*x + 1" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = var('x')\n", "g = (1-7*x)/(1-5*x+6*x^2)\n", "taylor(g,x,0,9)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "def t(n):\n", " \"\"\"\n", " Returns the value of h(n) using the\n", " Taylor series of the generating function.\n", " \"\"\"\n", " tn = taylor(g,x,0,n)\n", " return tn.coefficient(x,n)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, -2, -16, -68, -244, -812, -2596, -8108, -24964, -76172]" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[t(k) for k in range(10)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we solve the second part of the question\n", "and use the recursion relation." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "def h(n, D={}):\n", " \"\"\"\n", " Applies memoization to compute h(n).\n", " \"\"\"\n", " if n in D:\n", " return D[n]\n", " else:\n", " if(n == 0):\n", " result = 1\n", " elif(n == 1):\n", " result = -2\n", " else:\n", " result = 5*h(n-1) - 6*h(n-2)\n", " D[n] = result\n", " return result" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, -2, -16, -68, -244, -812, -2596, -8108, -24964, -76172]" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[h(k) for k in range(10)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We observe that we get the same sequence of numbers as before." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-7188041199657724841646073178059495579561383474850002702724" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "h(120)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thanks to the memoization, the function ``h`` can efficiently compute `$h(120)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Legendre polynomials are defined by\n", "\n", "$$\n", " P_0(x) = 1, \\quad P_1(x) = x,\n", " \\quad P_n(x) = \\frac{2n-1}{n} x P_{n-1}(x)\n", " - \\frac{n-1}{n} P_{n-2}(x), \\ {\\rm for} \\ n \\geq 2.\n", "$$\n", "\n", "Write an *efficient, recursive* function ``legendre`` to compute\n", "$P_n(x)$. The function ``legendre`` takes on input\n", " the degree $n$ and the variable $x$.\n", "\n", " \n", "Compare the output of your ``legendre`` (50, $x$ )\n", "with the ``legendre_P`` (50, $x$ )." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 5" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "def legendre(n, x, D={}):\n", " \"\"\"\n", " Returns the n-th Legendre polynomial in x,\n", " with memoization.\n", " \"\"\"\n", " if (n, x) in D:\n", " return D[(n, x)]\n", " else:\n", " if(n == 0):\n", " result = 1\n", " elif(n == 1):\n", " result = x\n", " else:\n", " result = expand((2*n-1)/n*x*legendre(n-1,x) \n", " - (n-1)/n*legendre(n-2,x))\n", " D[(n, x)] = result\n", " return result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We test some simple cases." ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "legendre(0, x)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "legendre(1, x)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3/2*x^2 - 1/2" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "legendre(2, x)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5/2*x^3 - 3/2*x" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "legendre(3,x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, see if we can compute the Legendre polynomial of degree 50." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "12611418068195524166851562157/140737488355328*x^50 - 156050375086257748529223875175/140737488355328*x^48 + 226836112238787036521861509275/35184372088832*x^46 - 823773249709279237895181270525/35184372088832*x^44 + 4189728463575151392735706892025/70368744177664*x^42 - 7928255400303748020099876118755/70368744177664*x^40 + 5790298887862287879848224131675/35184372088832*x^38 - 6684039602901787158511168414725/35184372088832*x^36 + 24770264410753681822717859419275/140737488355328*x^34 - 18602568051449552212241926551825/140737488355328*x^32 + 1423900270604780539702468452115/17592186044416*x^30 - 712769410486857922635873160725/17592186044416*x^28 + 583174972216520118520259858775/35184372088832*x^26 - 194391657405506706173419952925/35184372088832*x^24 + 26248579962778792027330678575/17592186044416*x^22 - 5693353963757653481984400705/17592186044416*x^20 + 7838675747202566388239392275/140737488355328*x^18 - 1052956443654076082002306425/140737488355328*x^16 + 26998883170617335435956575/35184372088832*x^14 - 2052546673789621992207225/35184372088832*x^12 + 222078820442811559812585/70368744177664*x^10 - 8065816723104536070675/70368744177664*x^8 + 90048990529077755175/35184372088832*x^6 - 1067774591253886425/35184372088832*x^4 + 20146690401016725/140737488355328*x^2 - 15801325804719/140737488355328" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "legendre(50, x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Legendre polynomials are already defined in SageMath." ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5/2*x^3 - 3/2*x" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = polygen(QQ)\n", "legendre_P(3,x)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "legendre_P(50,x) - legendre(50, x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the test whether our Legendre agrees with the Legendre of SageMath is straightforward because of the canonical form in which the polynomials are stored." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the point $P = (1,1)$ on the curve defined by $xy - 2 x + 1 = 0$.\n", "\n", "Compute the slope of the tangent line to the curve at $P$ in two ways:\n", "\n", "1. with implicit differentiation,\n", "\n", "2. with a Taylor series." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 6" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With SageMath, the fastest way to solve this problem is via the Taylor series, truncated after the linear term. Let us do that first." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "x,y = var('x,y')\n", "f = x*y - 2*x + 1" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-x + y" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ts = taylor(f,(x,1),(y,1),1)\n", "ts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the tangent line is $y = x$, therefore, the slope is 1.\n", "\n", "Let us now do implicit differentiation." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x*Y(x) - 2*x + 1 == 0" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Y = function('Y')(x)\n", "F = f(y=Y) == 0\n", "F" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "diff(Y(x), x)" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dy = diff(Y,x)\n", "dy" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x*diff(Y(x), x) + Y(x) - 2 == 0" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dF = diff(F,x)\n", "dF" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[diff(Y(x), x) == -(Y(x) - 2)/x]" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol = solve(dF,dy)\n", "sol" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-(Y(x) - 2)/x" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope = sol[0].rhs()\n", "slope" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "the value of the slope is: 1\n" ] } ], "source": [ "print('the value of the slope is:', slope.subs({x:1}).subs({Y(x=1):1}))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 4: plotting and solving equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 7" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we want to plot the curve $x^4 + x^2 y^2 - y^2 = 0$\n", "for $x$ and $y$ both between $-1$ and $+1$.\n", "\n", "1. Sampling this curve as given in rectangular coordinates,\n", " how many samples do we need to take from the curve to\n", " obtain a nice plot?\n", " \n", "\n", "2. Convert the curve into polar coordinates and plot.\n", "\n", " How many samples of the curve are needed here?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 7" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x^4 + x^2*y^2 - y^2" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, y = var('x, y')\n", "f = x^4 + x^2*y^2 - y^2\n", "f" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAATsAAAEhCAYAAAAAkQc6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwiklEQVR4nO3deVxU9foH8M/IMojiZCJbIliuqSliAua+ILhk2vVqGqmZaWa5VF791S1tESs1rVyuXipz370VZmKKuGACF8rC7SYKKqiYDmgJIt/fH0+HCUFlmTPfOWee9+s1r3MYzzDPyOHDWb6LQQghwBhjOldDdgGMMWYLHHaMMYfAYccYcwgcdowxh8BhxxhzCBx2jDGHwGHHGHMIHHaMMYfgLLsAaysuLsb58+fh4eEBg8EguxzGmBUIIZCfnw8/Pz/UqFG1YzTdhd358+fh7+8vuwzGmAqysrLQoEGDKr1W1bBLSEjAhx9+iJSUFGRnZ2Pr1q144okn7vqavXv3YurUqfjll1/g5+eHadOmYfz48RV+Tw8PDwD0n1KnTp3qlM8YsxN5eXnw9/cv+f2uClXD7vr162jTpg1Gjx6NJ5988p7bZ2RkoG/fvhg7dixWrVqFAwcOYMKECahfv36FXg+g5NS1Tp06HHaM6Ux1Lk2pGnaRkZGIjIys8PZLly5Fw4YNsWDBAgBAixYtkJycjLlz51Y47BhjrDx2dTc2MTER4eHhpZ7r06cPkpOTcfPmzXJfU1BQgLy8vFIPGa5dA+LjgdRUKW/PmOoyMoBvvwV++012JVVjV2GXk5MDb2/vUs95e3ujqKgIubm55b4mOjoaJpOp5CHr5sT77wPduwMLF0p5e8ZUt3Ej0LcvMHas7Eqqxq7CDih7Tq4Mt3enc/UZM2bAbDaXPLKyslSvsTxhYbRMTJTy9oypTtm3lX1da+yq6YmPjw9ycnJKPXfx4kU4OzujXr165b7GaDTCaDTaory7Cgmh5YkTwOXLwB3KZUyThNB+2NnVkV1YWBji4uJKPbdz5060b98eLi4ukqqqmHr1gGbNaJ2P7pjeZGQAFy4ALi5AcLDsaqpG1bC7du0a0tLSkJaWBoCalqSlpSEzMxMAnYI+88wzJduPHz8eZ86cwdSpU3H06FF89tlniImJwauvvqpmmVbz2GO0PHBAbh2MWZuyT7drB7i5ya2lqlQNu+TkZAQFBSEoKAgAMHXqVAQFBeHNN98EAGRnZ5cEHwA0atQI27dvR3x8PNq2bYt33nkHH3/8sWaanXDYMb1S9mllH9cig94m3MnLy4PJZILZbLZ5o+Ljx4HmzQGjETCbacmYHrRqBfzyC7BlCzBokO3f3xq/13Z1zU7rmjYF6tcHCgqA5GTZ1TBmHb/9RkEHAJ06ya2lOjjsrMhgADp3pvV9++TWwpi17N9Py+bN6Y+5VnHYWVmXLrRMSJBbB2PWouzLyh9yreKwszIl7PbvB4qK5NbCmDXEx9Oya1epZVQbh52VPfIIcN99QH4+95Nl2mc2W/ZjDjs7sWjRIjz88MN49NFHpdbh5GQ5utuzR2opjFXbvn1AcTHw0ENAFcfMtBu6CbsXX3wR6enpSEpKkl0KunenJYcd0zplH1b2aS3TTdjZE2XH2LcPKCyUWwtj1bF7Ny179JBbhzVw2KmgdWvA0xO4fh04fFh2NYxVTW4u8GdPTw47Vr4aNSw7x65dcmthrKqUo7pWrYDbhpnUJA47lfTqRUsOO6ZVyr6r7Mtax2Gnkt69aXnoEN2+Z0xLhAB27qR1ZV/WOg47lQQGAk2aALdu8V1Zpj0nTwJnzgCurtpvX6fgsFNRnz60/O47uXUwVlnKPtupE1CrltxarIXDTkVK2O3YQacFjGnFjh20VPZhPeCwU1H37nQacPo0cOyY7GoYq5g//rBceqnEtM92TzdhZy/dxf6qVi3L9Y5vv5VbC2MVtXcvBd4DD1CzE73QTdjZU3exv+rXj5axsXLrYKyilH21Xz8ao1EvdBN29koJu4QEboLC7J8QwDff0Lqy7+oFh53KGjemKRaLiviuLLN/v/xC15iNRqBnT9nVWBeHnQ0MGEDLr76SWwdj96Lsoz176qfJiYLDzgYGDqRlbCxw86bcWhi7m//8h5bKPqsnHHY2EBZGE5VcvcpzUzD7df68ZZQe5WxETzjsbMDJybLzbN0qtxbG7mTbNlqGhgK+vlJLUQWHnY0MHkzLbdtomGvG7M2WLbSUMQm2LXDY2UivXoCHB3DuHPDDD7KrYay0y5cts4g9+aTUUlTDYWcjRiPQvz+tb9oktxbGbrdtG43Q06YNTa6jRxx2NjRkCC03buSBAZh92biRlso+qke6CTt77Bt7u4gIoHZtICuLBvVkzB7k5lpGJf773+XWoibdhJ299o39q5o1Le2X1q2TWwtjii1b6BQ2KIgGnNUr3YSdVgwbRssNG2gHY0y2tWtpqeybesVhZ2Ph4cD99wM5OTxcO5Pv7Fka0gkAhg6VW4vaOOxszNXVchF49Wq5tTC2di3dLOvcGQgIkF2NujjsJBgxgpabNwO//y63FubYVq6kpbJP6hmHnQSPPUazj+XnW7roMGZraWnAkSN0tqHnu7AKDjsJatQAnnmG1leskFsLc1xffknLxx8H6taVW4stcNhJooRdXBxdJGbMlgoLgVWraH3kSLm12AqHnSQPPQR06UIXh/nojtlabCxw6RLg40ON3R0Bh51Ezz5Ly88+45FQmG39+9+0fOYZwNlZbi22wmEn0ZAhQJ06wKlTwO7dsqthjiIryzIJ9pgxcmuxJd2EnRb6xt7O3d1yy3/ZMrm1MMehnEl07Qo0bSq7GtsxCKGv8Tfy8vJgMplgNptRp04d2eXc048/Am3b0qnE2bOAt7fsipieFRVRs6dz54A1a4CnnpJdUcVY4/daN0d2WtWmDQ2DXVRkuY7CmFq++YaCrn59y+jZjoLDzg5MmEDLf/2LQo8xtSxaRMtnn6UBZR0Jh50dGDIE8PSkC8c8tyxTy9GjNG5djRrA+PGyq7E9Djs74OYGjB1L6x9/LLcWpl+ffELLAQPoup2j4bCzExMm0JSLe/dSn0XGrOnKFUvj9UmT5NYiC4ednWjQwDL000cfya2F6c+yZTTCziOPAN26ya5GDpuE3eLFi9GoUSO4ubkhODgY+/btu+O28fHxMBgMZR7Hjh2zRalSTZlCy7Vr6Y4ZY9ZQWGi5PDJlCmAwyK1HFtXDbv369Zg8eTJef/11pKamonPnzoiMjERmZuZdX3f8+HFkZ2eXPJroeXD8P3XoQIMo3rzJ1+6Y9axdC5w/D/j6aqddnRpUD7v58+djzJgxeO6559CiRQssWLAA/v7+WLJkyV1f5+XlBR8fn5KHk5OT2qXahddeo+XSpYDZLLcWpn3FxcAHH9D6pEmO19zkr1QNu8LCQqSkpCA8PLzU8+Hh4Th48OBdXxsUFARfX1/07NkTe+4yWUNBQQHy8vJKPbSsXz/g4YeBvDxg8WLZ1TCt+/prID2d+mA7YnOTv1I17HJzc3Hr1i1439YHytvbGzk5OeW+xtfXF8uWLcPmzZuxZcsWNGvWDD179kRCQkK520dHR8NkMpU8/P39rf45bKlGDWD6dFr/6CMetp1VnRDA7Nm0PmECYDLJrUc2VfvGnj9/Hg888AAOHjyIsLCwkuffe+89rFy5ssI3HQYMGACDwYCvymlxW1BQgIKCgpKv8/Ly4O/vr5m+seW5eRNo1gzIyAAWLHDcpgKseuLiaDa7mjVpX9Jyv2u77xvr6ekJJyenMkdxFy9eLHO0dzehoaE4efJkuf9mNBpRp06dUg+tc3GxHN29/z5w44bcepj2CAHMmkXrzz+v7aCzFlXDztXVFcHBwYiLiyv1fFxcHDp27Fjh75OamgpfX19rl2fXRo0CGjYEsrOpzyxjlbFrF3DgAN2QmDZNdjX2QfUxSqdOnYqoqCi0b98eYWFhWLZsGTIzMzH+z6ulM2bMwLlz5/Dln7N/LFiwAIGBgWjZsiUKCwuxatUqbN68GZs3b1a7VLvi6gq8/jowbhwQHU3dydzdZVfFtEAI4M03aX38eMDPT2499kL1sBs6dCguX76Mt99+G9nZ2WjVqhW2b9+OgD9n5M3Ozi7V5q6wsBCvvvoqzp07h5o1a6Jly5aIjY1F37591S7V7owaBcyZQ9dbPv2U/0KziomNBQ4domt1yuUQxoN32r0vv6TZn+rWpeHb77tPdkXMnhUXA0FBwE8/0R/H99+XXZF12P0NClZ9I0ZQu7srV/Sz4zL1rF5NQWcyAf/4h+xq7AuHnZ1zcqJrdgA1Q+E5Ztmd3LgBvPEGrU+fDtx/v9x67A2HnQYMGEB9Zm/coJsWjJVn4UIgM5NG0OG2mWVx2GmAwQDMm0frX34JJCfLrYfZnwsXgPfeo/X33qObE6w0DjuNePRR4OmnaX3SJGpewJjijTeA/HwgONiyn7DSdBN2Wpw3trKio6mt3cGDNA0eYwAd6cfE0PrChdS/mpXFTU80ZvZsum7n6wscO0ajWTDHVVwMdOwI/PAD3blftUp2RergpicO6JVXgMaNqRuZ0kqeOa5//5uCzsPDMm4dKx+HncYYjZa5Pz/5BPjvf+XWw+S5eNHSQ+Ltt7lb2L1w2GlQeDgwdCidwjz/PE+s7aimTKHG5m3bAhMnyq7G/nHYadSCBdR1LCWF1pljiY2lm1Q1agDLlwPOqvdy1z4OO43y8bG0vfvnP4ETJ+TWw2zHbKbRcABg8mSgfXup5WgGh52GjR4N9O5NPStGjwZu3ZJdEbOFqVNpqs3GjYF33pFdjXZw2GmYwUCnMB4e1PZOOdJj+vXNN8Bnn9HP/rPPeIzDyuCw07iAAMs1uzfeAH78UWo5TEUXLwJjxtD61KnUX5pVHIedDoweDQwcSBP1DB/OM5LpkRAUdBcvAq1aAe++K7si7dFN2DlCd7E7UU5nfXxojtCpU2VXxKztk0/oFNZopDHr3NxkV6Q93F1MR3btojZ4QgDr1wN//7vsipg1pKRQl7DCQgo9R2xTx93FWCm9egEzZtD6mDHA8eNy62HVd+UKMGQIBd0TTwAvvii7Iu3isNOZWbOALl2Aa9eAJ5+kJdOm4mLgmWdowqVGjSx3YVnVcNjpjLMzsG4dXb/75Rfg2Wd57Dutevtty3W6TZto0iVWdRx2OuTrS78cLi7Axo2WOSyYdmzZQkfpAE2S3q6d3Hr0gMNOpx57jOaaBWj8u61b5dbDKi41FYiKovVJk2gqTVZ9HHY69vzzljt3I0YASUly62H3dvYs0L8/tZUMDwfmzpVdkX5w2OncRx8BERHAH3/QL9GpU7IrYndiNgN9+wLnz9NcwevX82gm1sRhp3POzsCGDUCbNtT6PiKClsy+3LgBDBoEHDlCN5diY2kIL2Y9HHYOwMMD2L6d+tGePEmBZzbLroopioqom9+ePZafVWCg7Kr0h8POQfj5ATt3AvXr0wXwfv2A69dlV8Vu3aLmQVu3UhOT//wHCAqSXZU+6SbsHLlvbEU1bUqBZzIBBw4AAwbwoAEyKcPqr1wJODnR5Ybu3WVXpV/cN9YB/fADDfqZnw9060YNV2vVkl2VY7l1i4Lus89oaPU1a2heEVY+7hvLqiQkBNixg64PxccDffrwNTxbunmT2s4pQbdqFQedLXDYOaiOHYG4OMspbffuwIULsqvSvz/+AP72NxqmydkZWLsWeOop2VU5Bg47BxYSQkd2Xl500+Kxx4D//U92Vfr122/UUPirr2g8uq1beRguW+Kwc3Bt2wL791NTh19/BcLCgMRE2VXpz6lTdDS9fz8dTX/3HTXyZrbDYcfQpAkFXLt2QG4undKuXi27Kv1ISKCj6OPHgQYNKPC6dJFdlePhsGMAqNX+3r00l0VBAfD008Brr1GDV1Y1QgBLlgA9e9IfkeBguhPeqpXsyhwThx0rUbs2DS2kjHY8dy5dY+IbF5V3/TowahQwYQL9wRg6lI7w/PxkV+a4OOxYKTVqALNnUwPXWrWoC1ObNjS/BauYn38GOnQAvvyS/j8/+IDuuvIcr3Jx2LFyDRkCJCcDLVvSkV14OJ3WFhTIrsx+FRcDH38MtG9Ps7z5+ADff0//bzycunwcduyOmjcHDh8Gxo2j609z59J1p+Rk2ZXZn4wM6pUyaRL9QYiMpAnLu3WTXRlT6CbsuG+sOtzdgaVLqYO6lxfNaxESQkcrPJAAXY/76COgdWtg926gZk0aITo2lv6/mP3gvrGswnJzgZdeogl9AKBhQ2D+fGDwYMc8TTtwgEaCTkujr7t0AWJigMaNpZalS9w3ltmUpyddaP/mGxobLzOTuj717An897+yq7Od06epi1enThR0desCy5bRzRwOOvvFYccqrV8/ugD/xhs0BtuePXQt76mngBMnZFennpwcuibXrBkd3RoMwHPPUWPhsWPpziuzX/zjYVXi7g688w79og8fTs+tWwe0aEENkn/+WW591pSVBUyeDDz4IN1tLSyko9mUFGD5choQldk/DjtWLQEB1LUsNZUGAy0upq9bt6Y7kt99R89p0eHDFNwPPggsXEgjloSG0mgxu3bxiMJaw2HHrKJtWxrNIyWFruMZDDRmXkQEjZA8Zw5w7pzsKu/t6lW6+9y+Pd11Xr2a7rh260bBffAg0KuX7CpZVfDdWKaKX38FPvkE+PxzIC+PnjMYgB49aFijJ56wn6YZeXk0yc2GDbRUGk67ulI3r0mT6Jokk8cqv9fCBhYtWiQCAwOF0WgU7dq1EwkJCXfdPj4+XrRr104YjUbRqFEjsWTJkgq/l9lsFgCE2WyubtnMCq5dEyImRohOnYSgpsn0MBiECAsTYtYsIQ4cEKKgwHY13bolRFqaEHPnCtG7txCurqVra9VKiHnzhLh0yXY1sbuzxu+16mG3bt064eLiIpYvXy7S09PFpEmTRK1atcSZM2fK3f7UqVPC3d1dTJo0SaSnp4vly5cLFxcXsWnTpgq9H4ed/Tp1SojoaCGCg0uHCyBEzZpCdO0qxGuvCbFmjRBHjghx40b13/PmTSFOnBBiyxYh/vlPISIihLjvvrLv37SpEDNmUAgWF1f/fZl1WeP3WvXT2JCQELRr1w5Lliwpea5FixZ44oknEB0dXWb7f/zjH/jqq69w9OjRkufGjx+PH3/8EYkVGFWST2O14exZ4NtvabazPXuAy5fLblOjBuDvTyOFODvTc0rj5WHDgDFjABcX6t2xcKHlRogQtH7xIrWJu3mz7PeuXRvo3Jm6ePXtS81JmP2yxu+1s5VrKqWwsBApKSmYPn16qefDw8Nx8ODBcl+TmJiI8PDwUs/16dMHMTExuHnzJlxcXEr9W0FBAQr+0js9T7lAxOxagwbUNm3sWAqmY8co+JYtA5S/c8XFwJkz9LhdQgKN+NuiBTBoUMXe08eH+vkOGEAjuTiruvcze6Pqjzs3Nxe3bt2Ct7d3qee9vb2Rk5NT7mtycnLK3b6oqAi5ubnw9fUt9W/R0dGYNWuWdQtnNnH5MvUn3b0b2LePGipX9DyjWzfqnlWvHjVmXrv23q/JyQFmzaKhlzp3pu/RuzcFL9M/m/xtM9zWcVIIUea5e21f3vMAMGPGDEydOrXk67y8PPj7+1enXKaiX38FNm+mU89Dh8q2wQsIoPZrrVvTqCtNmtBz9evfuf/tmjX0UAhBTUjOnKH3O3aMGjmnplIPj4wMenz5JW3fujWN0Dx4MDWhccR+vo5A1bDz9PSEk5NTmaO4ixcvljl6U/j4+JS7vbOzM+rVq1dme6PRCKPRaL2imdVdvkxhtHIlkJRU+t9atqR2a1270oQ0d9gtKsVgoP6qdetSeP2V2Uwhm5BAY80lJQFHjtDj3XepTeDTTwPPPEMhy3TESjdL7qhDhw7ihRdeKPVcixYtxPTp08vdftq0aaJFixalnhs/frwIDQ2t0Pvx3Vj7UFxMTUqGDy/dtMPJSYhevYRYvFiIzEzZVQqRmyvEypVCDBokhJtb6aYxERFCfP01NVVhcmmq6UlMTIxIT08XkydPFrVq1RKnT58WQggxffp0ERUVVbK90vRkypQpIj09XcTExHDTEw0pKhJiwwYhHn20dNOOoCAhFi4U4sIF2RXeWV6eECtWCNG9e+naGzemcP79d9kVOi5NhJ0Q1Kg4ICBAuLq6inbt2om9e/eW/NvIkSNF165dS20fHx8vgoKChKurqwgMDORGxRpQVERHSE2bWkLCaBTi2WeFSE6WXV3lnTwpxCuvlG6T5+0txPz5HHoyaKKdna1xOzvbEoL6xM6YYWkyUrcuDfI5caL2RwS5do26vM2dS+P3AdTu7623gGef5eYrtsKDdzKp0tJoQu0nnqCgq1uXZiY7c4aaeGg96ABqfPzSS8DJkzScU0AAcP48tddr04baBjJt4LBjlWY201FbcDBNrO3mRkd2GRm09PCQXaH1ubpaBupcsIDa96WnA336AE8+SWPeMfvGYccqZetW6rWwaBG1kRs6lAJg9mzq0aB3RiONgnLyJC2dnGhi8Ycfpol2tDp2nyPgsGMVcvky9VQYPBjIzqbGvt9/T6MTN2wouzrbq1uXjvBSU6l94LVrdLrbrRtw6pTs6lh5dBN2PJWienbtol4G69bRkcyMGcBPP9HYdI6udWvq6rZoEVCrFq23aQN88UXFu74x2+C7seyOioqAf/6TRhkGaGSQlSsB/ntSvowMYNQo6p0B0JHwv/6lz2uYtsZ3Y5lqsrPpyE0JunHjaLpEDro7a9SIBjWYPZuOgNeupeHd9TT5kJZx2LEyDh2iO6379tFRyYYNNC+Du7vsyuyfcpqfkECjqZw4QZP0bNokuzLGYcdKWbmSOuVnZ9MdxuRkYMgQ2VVpT8eOdCTcsydw/Tr9H86cydfxZOKwYwDol/DNN2m0j8JCaih86BCNAsKqpn59mmFtyhT6etYsGlHlL2PNMhvisGO4eRMYPZomvQaA6dNpzDm+sF59zs7A/PnAv/9N62vW0PSSZrPsyhwPh52D++MPGtZ8xQq63rR8ORAdTfM/MOsZM4bm3PDwAOLj6VLBhQuyq3IsvEs7sPx8IDISiI2lLl/btlGXKKaOXr3oxoW3N/DjjzQ0PHczsx0OOwdlNlO/zr17gTp1gLg4oH9/2VXpX9u2wP79NKDAyZM0j8bp07Krcgwcdg4oL4+uGyUmUren778HOnWSXZXjaNyYmvU0bkxB161b+TOoMevisHMw168D/frRnVYl6Nq3l12V4/H3p6PqJk0o6Hr0AM6dk12Vvukm7Lhv7L0VFFBH/v37aYSSXbtoJi8mh58fTRD+4IM0eEDv3kBuruyq9Iv7xjqI4mLqq7lhA3VYj4sDwsJkV8UAOpXt3Bk4e5a64+3eTYOGMgvuG8sqRAhq2LphA+DiQmPScdDZj8BA+uNTrx5N7ThkCLV9ZNbFYecAFiwAPv6Y1r/8kk6XmH1p3pyaANWsSb0uXnyRu5ZZG4edzn31FfDKK7T+4YfAsGFy62F3FhICrF9PDbqXL6dJfpj1cNjp2JEjwPDhdIQwbpwl9Jj9GjAA+OgjWv/HP4Cvv5Zbj55w2OnU5cvAwIHU1KRHD+CTTwCDQXZVrCJeegkYP57+SI0YYZmiklUPh50O3bpFR3QZGTSgpHJjgmmDwUDXWLt0oS59gwZRQ3BWPRx2OjRrFs1n6u5O/V3r1ZNdEassFxdg40YaAPT4cZqQm29YVA+Hnc589x3w7ru0vmwZ8MgjcuthVeflRSMcu7jQkFvKHXVWNRx2OnL+PBAVRUcA48fT9R6mbSEhwLx5tP7aazRyNKsaDjudKC6mUYYvXaKp/JQ7ekz7Jk6kbn43b1IvmPx82RVpk27CztH7xs6fT5363d1pflc3N9kVMWsxGGikY39/4H//AyZNkl2RNnHfWB346SfqU1lYSNfpxo6VXRFTQ0ICDQclBLBlC92ldRTcN5ahsJCu0xUWAo8/ziMN61mXLtTQGACefx64eFFuPVrDYadx775LR3aennRUxw2H9W3mTLrDnpsLTJgguxpt4bDTsLQ0mhwHABYvprkNmL4ZjTQ5krMzNUfZuFF2RdrBYadRRUV0ylpUBDz5JE9k7UjatgVmzKD1iROB336TWo5mcNhp1McfAykpwH33AZ9+KrsaZmuvvw60aEHX7aZNk12NNnDYaVBmJvDmm7T+4YeAj4/cepjtGY00DBQAxMTQBD7s7jjsNGjyZBrNpFMn6jPJHNNjj1maGU2YwKMb3wuHncZ89x0Nq+7kRDclavBP0KFFR9NADz//zJcz7oV/VTSksNDSev7ll4HWreXWw+SrVw+YM4fWZ84ELlyQWo5d47DTkE8/peF+vLyAt96SXQ2zF88+S3P/5uXRjQtWPt2End77xubmAm+/TeuzZ9O8r4wBdClj4UJa/+wzan/JyuK+sRrx0kt0ZNe2LQ3z4+QkuyJmb556igaB6NGDJkDXU28a7hvrIE6eBJYupfV58zjoWPmiowFXV5pke8cO2dXYHw47DXjjDeopERlJf7UZK09gIN24AmjAgOJiqeXYHQ47O5eSQhPmGAyWu26M3cmMGXQ998gRYM0a2dXYFw47O/fGG7QcMYLnk2D3dv/9lmGg3nqLGxr/FYedHTtwgK69ODtTGyrGKuLll6l50qlTwBdfyK7GfnDY2TGlLd3o0cBDD8mthWlHrVqWUVHefZcaozOVw+7KlSuIioqCyWSCyWRCVFQUrl69etfXjBo1CgaDodQjNDRUzTLt0v79NKeEiws3FGWVN24c4OtLg0asWCG7GvugatgNHz4caWlp2LFjB3bs2IG0tDRERUXd83URERHIzs4ueWzfvl3NMu3SO+/QctQoICBAailMg2rWtAz9NHs2X7sDAGe1vvHRo0exY8cOHDp0CCEhIQCA5cuXIywsDMePH0ezZs3u+Fqj0QgfBx63KCkJ2LmT2tMppyOMVdbzz1Pbu9OngbVraapNR6bakV1iYiJMJlNJ0AFAaGgoTCYTDh48eNfXxsfHw8vLC02bNsXYsWNx8S4zixQUFCAvL6/UQ+uUodZHjAAaNZJbC9Mud3dg6lRanzOH292pFnY5OTnw8vIq87yXlxdycnLu+LrIyEisXr0au3fvxrx585CUlIQePXqgoKCg3O2jo6NLrgmaTCb4+/tb7TPIcPw4sG0brStNCBirqhdeoHZ3R48CX38tuxq5Kh12M2fOLHMD4fZHcnIyAMBQTuc8IUS5zyuGDh2Kfv36oVWrVhgwYAC+/fZbnDhxArGxseVuP2PGDJjN5pJHVlZWZT+SXZk7l+YFffxx4OGHZVfDtK5OHQo8APjgA7m1yFbpa3YTJ07EsGHD7rpNYGAgfvrpJ1woZ3CtS5cuwbsS02D5+voiICAAJ0+eLPffjUYjjEZjhb+fPbtwAVi5ktZ5XgFmLS+/DMyfDxw8CCQmAmFhsiuSo9Jh5+npCU9Pz3tuFxYWBrPZjMOHD6NDhw4AgB9++AFmsxkdO3as8PtdvnwZWVlZ8PX1rWypmrN4MVBQAISE0JDbjFmDry/w9NM0/NO8ecCmTbIrkkO1a3YtWrRAREQExo4di0OHDuHQoUMYO3Ys+vfvX+pObPPmzbF161YAwLVr1/Dqq68iMTERp0+fRnx8PAYMGABPT08MGjRIrVLtwo0bwJIltP7KK3JrYfozZQott26lu7OOSNV2dqtXr0br1q0RHh6O8PBwPPLII1ipnKf96fjx4zCbzQAAJycnHDlyBAMHDkTTpk0xcuRING3aFImJifDw8FCzVOnWrgUuXQL8/QGd5zqToFUroFcvuiO7aJHsauTgwTvtgBBAcDCQmkpNBPguLFPD11/Tja+6dYGzZ6lpilbw4J06kZhIQefmBjz3nOxqmF717UvtNq9ccczhnzjs7MDixbQcNoxmi2JMDU5ONL8sQPucvs7p7o3DTrLcXGDjRlpXdkTG1DJ6NGA00pnE4cOyq7EtDjvJvviChuAJDgZ0OjEasyP16gF//zut/+tfcmuxNd2EnRanUhQCWL6c1p9/Xm4tzHGMG0fL9euBPxtCOAS+GyvRvn1Aly402GJ2NqDz1jXMTggBtGxJ/WWXLrWEnz3ju7EaFxNDy6FDOeiY7RgMlrv+yj7oCDjsJMnPt9yYGDNGbi3M8Tz9NM1tkpQE/Pyz7Gpsg8NOkk2bgN9/B5o2ddyO2UweLy+gf39ad5RJeTjsJFHmBRg1ik4rGLO1UaNouXo1TcKudxx2Epw5A+zdSyH39NOyq2GOKjKSmqLk5AC7dsmuRn0cdhKsXk3Lbt2o4z9jMri6Uq8dAFi1Sm4ttsBhZ2NCWMKOj+qYbMo+uG0bcP261FJUx2FnY0eOAOnp9Fd18GDZ1TBHFxICPPggBZ3e56jgsLOxtWtp2a8fcN99UkthDAYD8NRTtK7sm3rFYWdDQlAXHcByrYQx2ZR9cccO4OpVqaWoSjdhp4W+sSkpQEYGDZrYr5/sahgjLVvSTHaFhcBXX8muRj26CbsXX3wR6enpSEpKkl3KHSk9Jvr1o/6wjNkDg8EyEoqyj+qRbsLO3glhmdVpyBC5tTB2u7/9jZY7d+p3JBQOOxv58Ufg1Ckaer1vX9nVMFbaww8DzZvTqewd5qPXPA47G9myhZYREXwKy+yPwWBpCqXsq3rDYWcj27bRkqdJZPZK2Td37AD++ENuLWrgsLOBU6eoMbGTk2WkCcbsTXAw0KABNTD+/nvZ1Vgfh50NKLfzu3QB7r9fbi2M3YnBAAwcSOv/+Y/cWtTAYWcDStg9/rjcOhi7F2Uf/eYboLhYbi3WxmGnsqtXaa4JABgwQGopjN1T165A7do07FNKiuxqrIvDTmU7d9LAiM2bAw89JLsaxu7OaATCw2ldb01QOOxUpuwwfGOCaYWyr3LY2Sl77BtbXEy38QFuSMy0IzKSlsnJwIULcmuxJt2EnT32jf3vf4GLF+kayGOPya6GsYrx8QGCgmj9u+/k1mJNugk7e6TsKL160WCdjGmFcnSnnJnoAYedipSw69NHbh2MVZayz8bF6acJCoedSvLygMREWuewY1oTFgZ4eAC5uUBqquxqrIPDTiXx8dTk5KGHgEaNZFfDWOW4uADdu9N6XJzcWqyFw04lyjycvXvLrYOxqurVi5Z6mVOWw04lyg6i7DCMaY2y7+7fr49RUDjsVJCdDRw9Sh2rlVMBxrSmeXPAzw8oKAAOHpRdTfVx2Klg925aBgXxKCdMuwwGoEcPWlf2aS3jsFPBnj20VHYUxrRK2YeVfVrLOOxUEB9Py27dZFbBWPUp+3BSEnDtmtRSqk03YWcvfWPPngV+/RWoUQPo1ElqKYxVW2Ag0LAhNaPS+nU73YSdvfSN3buXlkFBgMkktRTGqs1gsBzdKfu2Vukm7OxFQgItu3aVWwdj1tKlCy2VfVurOOysTBmVWNlBGNM6ZV8+fBi4cUNuLdXBYWdFubnUvg7gIZ2YfjRuDHh70wTadjSCWqVx2FnRgQO0bNEC8PSUWwtj1mIwWG627d8vt5bq4LCzIiXs+KiO6Y0Sdso+rkUcdlak3JrnsGN607EjLRMTtTu+naph995776Fjx45wd3fHfffdV6HXCCEwc+ZM+Pn5oWbNmujWrRt++eUXNcu0ioICGrMfsOwYjOlFUBDg5gb89htw4oTsaqpG1bArLCzEkCFD8MILL1T4NR988AHmz5+PTz/9FElJSfDx8UHv3r2Rn5+vYqXVl5ZGgefpCTRpIrsaxqzLxQVQ2usrg9JqjaphN2vWLEyZMgWtW7eu0PZCCCxYsACvv/46Bg8ejFatWmHFihX4/fffsWbNGjVLrTZlBwgNpQu6jOlNWBgttRp2zrIL+KuMjAzk5OQgXJmlF4DRaETXrl1x8OBBjBs3rsxrCgoKUFBQUPJ1Xl6eTWq9Xbt2wIQJlr9+jOlNeDhw6RIQESG7kqqxq7DLyckBAHh7e5d63tvbG2fOnCn3NdHR0Zg1a5bqtd1Lly7ckJjpW8+e9NCqSp/Gzpw5EwaD4a6PZOVKfRUZbjsPFEKUeU4xY8YMmM3mkkdWVla13psxpk+VPrKbOHEihg0bdtdtAgMDq1SMj48PADrC8/X1LXn+4sWLZY72FEajEUajsUrvxxhzHJUOO09PT3iq1D2gUaNG8PHxQVxcHIL+nJK8sLAQe/fuxfvvv1+h7yGEACDv2h1jzPqU32fl97sqVL1ml5mZid9++w2ZmZm4desW0tLSAACNGzdG7dq1AQDNmzdHdHQ0Bg0aBIPBgMmTJ2P27Nlo0qQJmjRpgtmzZ8Pd3R3Dhw+v0HsqTVT8/f1V+UyMMXny8/NhquLYaaqG3ZtvvokVK1aUfK0cre3Zswfd/hwk6/jx4zCbzSXbTJs2DX/88QcmTJiAK1euICQkBDt37oSHh0eF3tPPzw9ZWVnw8PAouc6Xl5cHf39/ZGVloU6dOhX6Po8++milxsZTc/vKbMuf1Xrfv7Lbq/m9Hf2zCiGQn58PPz+/Cr/n7VQNuy+++AJffPHFXbe5/bDUYDBg5syZmDlzZpXes0aNGmjQoEG5/1anTp0K7yhOTk4V3lbt7Sv7vQH+rNb6/mrWzp/17m7/rFU9olNw39g7ePHFF+1m+8p+78riz2qd7dX+f6wsR/qsFWEQ1bnipxF5eXkwmUwwm82V/uuiNfxZ9Yk/a/U5xJGd0WjEW2+95RBNVPiz6hN/1upziCM7xhhziCM7xhjjsGOMOQQOO8aYQ+CwY4w5BA47xphD0G3YOdL8F1euXEFUVBRMJhNMJhOioqJw9erVu75m1KhRZYbmCg0NtU3BlbB48WI0atQIbm5uCA4Oxj5lFvI72Lt3L4KDg+Hm5oYHH3wQS5cutVGl1VeZzxofH1/u8GrHjh2zYcWVl5CQgAEDBsDPzw8GgwHbtm2752us9TPVbdg50vwXw4cPR1paGnbs2IEdO3YgLS0NUVFR93xdREQEsrOzSx7bt2+3QbUVt379ekyePBmvv/46UlNT0blzZ0RGRiIzM7Pc7TMyMtC3b1907twZqamp+L//+z+8/PLL2Lx5s40rr7zKflbF8ePHS/0Mm9j5BCjXr19HmzZt8Omnn1Zoe6v+TIXOff7558JkMt1zu+LiYuHj4yPmzJlT8tyNGzeEyWQSS5cuVbHC6klPTxcAxKFDh0qeS0xMFADEsWPH7vi6kSNHioEDB9qgwqrr0KGDGD9+fKnnmjdvLqZPn17u9tOmTRPNmzcv9dy4ceNEaGioajVaS2U/6549ewQAceXKFRtUpw4AYuvWrXfdxpo/U90e2VXWvea/sFeJiYkwmUwICQkpeS40NBQmk+medcfHx8PLywtNmzbF2LFjcfHiRbXLrbDCwkKkpKSU+nkAQHh4+B0/V2JiYpnt+/Tpg+TkZNy8eVO1WqurKp9VERQUBF9fX/Ts2RN79uxRs0wprPkz5bD7093mv1D+zR7l5OTAy8urzPNeXl53rTsyMhKrV6/G7t27MW/ePCQlJaFHjx6lJi+SKTc3F7du3arUzyMnJ6fc7YuKipCbm6tardVVlc/q6+uLZcuWYfPmzdiyZQuaNWuGnj17IiEhwRYl24w1f6Z2NeHOvcycOfOek+skJSWhffv2VX6Pysx/oaaKflagbM3AveseOnRoyXqrVq3Qvn17BAQEIDY2FoMHD65i1dZX2Z9HeduX97w9qsxnbdasGZo1a1bydVhYGLKysjB37lx00dnMT9b6mWoq7Oxt/gs1VfSz/vTTT7hw4UKZf7t06VKl6vb19UVAQABOnjxZ6VrV4OnpCScnpzJHNnf7efj4+JS7vbOzM+rVq6dardVVlc9antDQUKxatcra5UllzZ+ppsLO3ue/sKaKftawsDCYzWYcPnwYHTp0AAD88MMPMJvN6NixY4Xf7/Lly8jKyioV9DK5uroiODgYcXFxGDRoUMnzcXFxGDhwYLmvCQsLw9dff13quZ07d6J9+/ZwcXFRtd7qqMpnLU9qaqrd/Pysxao/00rf0tCIM2fOiNTUVDFr1ixRu3ZtkZqaKlJTU0V+fn7JNs2aNRNbtmwp+XrOnDnCZDKJLVu2iCNHjoinnnpK+Pr6iry8PBkfocIiIiLEI488IhITE0ViYqJo3bq16N+/f6lt/vpZ8/PzxSuvvCIOHjwoMjIyxJ49e0RYWJh44IEH7Oqzrlu3Tri4uIiYmBiRnp4uJk+eLGrVqiVOnz4thBBi+vTpIioqqmT7U6dOCXd3dzFlyhSRnp4uYmJihIuLi9i0aZOsj1Bhlf2sH330kdi6das4ceKE+Pnnn8X06dMFALF582ZZH6FC8vPzS34XAYj58+eL1NRUcebMGSGEuj9T3YbdyJEjBYAyjz179pRsA0B8/vnnJV8XFxeLt956S/j4+Aij0Si6dOkijhw5YvviK+ny5ctixIgRwsPDQ3h4eIgRI0aUaZLw18/6+++/i/DwcFG/fn3h4uIiGjZsKEaOHCkyMzNtX/w9LFq0SAQEBAhXV1fRrl07sXfv3pJ/GzlypOjatWup7ePj40VQUJBwdXUVgYGBYsmSJTauuOoq81nff/998dBDDwk3NzdRt25d0alTJxEbGyuh6spRmszc/hg5cqQQQt2fKY9nxxhzCNz0hDHmEDjsGGMOgcOOMeYQOOwYYw6Bw44x5hA47BhjDoHDjjHmEDjsGGMOgcOOMeYQOOwYYw6Bw44x5hD+H/JZGfhdUcC7AAAAAElFTkSuQmCC", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pltf = implicit_plot(f, (x,-1,+1), (y,-1,+1), plot_points=5000)\n", "pltf.show(figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even with 5000 points (which already takes some time), we do not get a nice plot around (0,0) which is a singular point on this curve in rectangular coordinates." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us now convert to polar coordinates." ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "r^4*cos(t)^4 + r^4*cos(t)^2*sin(t)^2 - r^2*sin(t)^2" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r, t = var('r,t')\n", "pf = f.subs(x=r*cos(t),y=r*sin(t))\n", "pf" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{r: -sin(t)/(sqrt(cos(t)^2 + sin(t)^2)*cos(t))},\n", " {r: sin(t)/(sqrt(cos(t)^2 + sin(t)^2)*cos(t))},\n", " {r: 0}]" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s = solve(pf, r, solution_dict=True)\n", "s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before plotting, let us simplify the expressions." ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-sin(t)/cos(t)" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rt0 = s[0][r].full_simplify()\n", "rt0" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "sin(t)/cos(t)" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rt1 = s[1][r].full_simplify()\n", "rt1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As the $\\cos(t)$ appears in the denominator, we have to avoid plotting around multiples of $\\pi/2$, otherwise we would divide by zero." ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAHWCAYAAADKGqhaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABT+klEQVR4nO3dd3RUVdcG8GeSkIQShh6pAZEWirQQioKChCpNpYeigiK8SJEmKmChqCioCSBVkCYdKQEUpCYgJSKGLlKEEEBSaAGS8/2xv1ADpMzMuXfm+a01Cxkzc/clM7Nnn3vOPhallAIRERHZhJvuAIiIiJwJEysREZENMbESERHZEBMrERGRDTGxEhER2RATKxERkQ0xsRIREdkQEysREZENMbESmYBSCvHx8WA/FyLjY2IlMoGEhARYrVYkJCToDoWInoCJlYiIyIaYWImIiGyIiZWIiMiGmFiJiIhsiImViIjIhphYiYiIbMhmiTU2FkhMtNWzEREROUZiIhAXZ7vns0liPX8eeOopYMUKWzwbERGR46xYAfj6AjExtnk+myRWX1+gYkVgwQJbPBsREZHjLFgAPPssUKCAbZ7PZkPBHToAa9bYtpwmIiKyp7g4yV3t2tnuOW2WWNu2BW7eBJYts9UzElFISAj8/f0REBCgOxQip7RihVxjbdvWds9pUTbs6l2vHpA1KxAWZqtnJCIAiI+Ph9VqRVxcHHLmzKk7HCKn0bQpcOUKsGWL7Z7TpsttOnQAfvkFuHDBls9KRERkexcvAhs2AO3b2/Z5bZpYX3lF/ly0yJbPSkREZHtLlwLJycCrr9r2eW2aWPPnBxo25OxgIiIyvoULgQYNbDcbOIXNOy916ABs3QqcPm3rZyYiIrKNc+eATZtsOxs4hc0Ta6tWgJeXfBMgIiIyosWLAQ8PoHVr2z+3zRNrzpxAs2YcDiYiIuNasABo1AjIk8f2z22XJvwdOgB79gBHj9rj2YmIiDLu5Elgxw7bzwZOYZfE2qwZkCMHMH++PZ6diIgo4+bNk54LLVrY5/ntklizZpVrrfPnA7ZrP0FERJQ5SgGzZwNt2gA+PvY5ht32Y+3QATh0CNi/315HICIiSp/duyU3deliv2PYLbE2bAjkzQv8+KO9jkBERJQ+s2cDBQvK+lV7sVtizZIF6NhREuvt2/Y6ChERUdrcvCmzgTt3Btzd7XccuyVWAOjaFYiOll6MREREOoWFSX9gew4DA3ZOrFWrAhUqALNm2fMoRM6L28YR2c7s2UCVKpKX7Mmm28al5ssvgeHDpXLNndueRyJyXtw2jihz/vtPrq2OHQv072/fY9m1YgWATp2ApCS2OCQiIn1++klyUYcO9j+W3RNrwYLSNorDwUREpMvs2ZKLnnrK/seye2IFgG7dgJ07Ze0QERGRIx09CoSH23/SUgqHJNaXXwZy5QJ++MERRyMiIrrrxx9lgxh7tTB8kEMSq7e3jGvPmSNj3ERERI6QnCzDwG3bSrtdR3BIYgVkOPjff4Fff3XUEYmIyNVt3Qr88w8QHOy4YzossQYEAGXLcjiYiIgcZ/p04JlngOefd9wxHZZYLRapWpcuBeLiHHVUIiJyVbGxwKJFwJtvSg5yFIclVkD6M968KeuJiIiI7GnePODWLWmv60gOTayFCwONGwNTpzryqES2FxoaihIlSsDb2xvVqlXD1q1bH/vziYmJGD58OPz8/ODl5YWSJUtixowZDoqWyDVNmwY0b+6Ytav38nDs4YCePWUT9H37pGcjkdksXLgQ/fr1Q2hoKOrUqYMpU6agSZMmiIqKQrFixVJ9TNu2bXH+/HlMnz4dzzzzDGJiYnCb2z4R2c3evZJnPv7Y8ce2e6/gB92+DRQrJsk1NNSRRyayjcDAQFStWhWTJk26c1+5cuXQqlUrjBkz5qGfDwsLQ/v27fH3338jT548GTomewUTpc877wArVgAnTwIeDi4hHToUDMgJvvGGLNi9etXRRyfKnJs3b2LPnj0ICgq67/6goCDs2LEj1cesXLkS1atXx+eff47ChQujdOnSeO+993D9+vVHHicxMRHx8fH33Ygoba5dA+bOBbp3d3xSBTQkVkAS65UrbMxP5nPx4kUkJSXB19f3vvt9fX0RHR2d6mP+/vtvbNu2DQcOHMCyZcswYcIELF68GL17937kccaMGQOr1XrnVrRoUZueB5EzW7wYiI8HXn9dz/G1JNbixaUZ8vff6zg6UeZZHpi7r5R66L4UycnJsFgsmDt3LmrUqIGmTZviq6++wqxZsx5ZtQ4bNgxxcXF3bqdPn7b5ORA5q2nTgAYNgKef1nN8LYkVkElMO3cCf/yhKwKi9MuXLx/c3d0fqk5jYmIeqmJTFCxYEIULF4bVar1zX7ly5aCUwpkzZ1J9jJeXF3LmzHnfjYie7K+/pNvSm2/qi0FbYk2ZAs2lN2Qmnp6eqFatGjZs2HDf/Rs2bEDt2rVTfUydOnVw9uxZXLly5c59R44cgZubG4oUKWLXeIlczaRJgK8v0KaNvhi0JdYsWWT8e84cudBMZBYDBgzAtGnTMGPGDBw8eBD9+/fHqVOn8PbbbwOQYdwu9+xP1bFjR+TNmxfdu3dHVFQUtmzZgkGDBuH1119HVkd1BSdyAVeuSMP9N98EPD31xaEtsQJy8gkJ7MRE5tKuXTtMmDABH3/8MSpXrowtW7ZgzZo18PPzAwCcO3cOp06duvPzOXLkwIYNGxAbG4vq1aujU6dOePnll/HNN9/oOgUipzR3rqw26dlTbxwOX8f6oMaNpZ9jRITOKIiMjetYiR5PKaByZaBECWD5cr2xaK1YAaBXL5nEtHu37kiIiMisduwA9u+XnKKb9sTavLl0YvruO92REBGRWU2aBJQsCTRsqDsSAyRWd3dpPbVgAXDhgu5oiIjIbGJiZHu4Xr0AN+1ZzQCJFbi7V9706bojISIis5kxQxJqt266IxGGSKx58wIdOkhTfm74QUREaZWUBEyeDLRvL7nECAyRWAGgTx/g9Gng5591R0JERGYRFiY72Bhh0lIK7ctt7lWnDuDtDfz6q+5IiIyFy22IUteoEfDff8CuXXJJ0QgMU7ECUrVu3Ci9HomIiB7nr7+A9euBfv2Mk1QBgyXWV16R/sEhIbojISIio5s4EShYEHjtNd2R3M9QidXTE3jrLen1GBurOxoi/UJCQuDv74+AgADdoRAZysWL0mu+d2+9fYFTY6hrrABw7hzg5weMHg28957uaIiMgddYie43ejTwyScy6TVfPt3R3M9QFSsgZX2nTlLi37qlOxoiIjKamzelW19wsPGSKmDAxAoAAwYAZ85w1xsiInrYokUyuvnuu7ojSZ3hhoJTNG4sbar27DHWbC8iHTgUTCSUAmrUAPLkAdat0x1N6gxZsQLAwIHAvn3Apk26IyEiIqPYsUN2Q+vXT3ckj2bYijVlb70iRYDVq3VHQ6QXK1Yi8eqrwIEDQFSUMRrup8agYcnw78CBwJo18g9IRESu7dgxYNkyqVaNmlQBAydWQJoqFyoEfPWV7kiIiEi3r76SWcBdu+qO5PEMnVg9PYG+fWURcHS07miIiEiXmBhg5kzJCVmz6o7m8QydWAHpxOTpyTaHRESu7NtvAXd3Y+1i8yiGT6y5cgFvvCF7tV67pjsaIiJytCtXpLjq2VOW2Rid4RMrIBeqY2NlGICIiFzL9OlAQoKxl9jcy7DLbR7UsaOsXzp6FMiSRXc0RI7F5Tbkqm7dAp55BqhXTzZoMQNTVKwA8P77skv83Lm6IyEiIkdZuBA4dQoYNEh3JGlnmooVAFq1Ag4elHWt7u66oyGyv5CQEISEhCApKQlHjhxhxUouJSkJqFABKFkSWLVKdzRpZ6rE+vvv0iNy4UKgbVvd0RA5DoeCyRUtWiSf9RERQGCg7mjSzlSJFQCCgmQ90759bM5ProOJlVxNcjJQpQrg6wusX687mvQxzTXWFMOHA3/8wf7BRETO7Oefgf37gQ8/1B1J+pmuYlUKeP55GXvfsYNVK7kGVqzkSpQCAgKAHDmA337THU36ma5itVikao2I4JZyRETOKCxM9uL+4APdkWSM6SpWQL7NVK8uXZl+/VV3NET2x4qVXIVSQJ068qdZRyVNV7EC8g/9wQfAxo3Ali26oyEiIlvZuBEID5drq2ZMqoBJK1ZAvs1Uqwb4+MgYvFl/AURpwYqVXMULL0hv4N9/N+/nuikrVkD+wT/+WCpWDgcTEZnf1q3A5s0yImnWpAqYuGIFpGqtVUt+AWYdiydKC1as5OyUAl58Ebh8WfoUuJm27DNxxQrcrVojIoC1a3VHQ0REGfXLL1KtfvqpuZMqYPKKFZBvOfXqAVevArt3s2ol58SKlZyZUkDNmvL5HR5u/s9xk38vkF/AJ58Ae/cCy5frjoaIiNLr55+BXbukWjV7UgWcoGJN8dJLwPnz0u7Q7MMIRA9ixUrOKqUncJ48stTGGRKr06SgTz4BDhyQ3RCInEVISAj8/f0REBCgOxQiu/jpJ+kJ7CzVKuBEFSsANG0KnDghCZb7tZIzYcVKzuj2baB8edlvdc0a3dHYjtNUrIDMED50CJg3T3ckRET0JHPmAEeOSLXqTJyqYgWA1q3lOuvBg4CXl+5oiGyDFSs5m8REoEwZ6fu+eLHuaGzLqSpWAPjsM+DkSWDyZN2REBHRo0yZApw6JSONzsbpKlYA6NkTWLoUOHZMdsAhMjtWrORMYmPlumrr1sC0abqjsT2nq1gBYNQo4Pp1YOxY3ZEQEdGDRo8GbtxwzmoVcNLEWrAg8N57wIQJMtRARETG8M8/wMSJwKBBQKFCuqOxD6ccCgaAhATgmWeAxo2BH37QHQ1R5nAomJxFx47Apk3A0aNAjhy6o7EPp6xYAdmnddQomc4dGak7GiIi2rULmD9fGvo4a1IFnLhiBYBbt4CKFYFixYD163VHQ5RxrFjJ7FI2TLl8WYodZ27i47QVKwBkyQKMGwds2MDESrYVGhqKEiVKwNvbG9WqVcPWrVvT9Ljt27fDw8MDlStXtm+ARAazYoVsZP7ll86dVAEnr1gB+ZZUty4QHy874Dj7L5Tsb+HChQgODkZoaCjq1KmDKVOmYNq0aYiKikKxYsUe+bi4uDhUrVoVzzzzDM6fP4/IdFyjYMVKZnbrlrQuLFECWLdOdzT25/SJFQB27pS9/qZOBd58U3c0ZHaBgYGoWrUqJk2adOe+cuXKoVWrVhgzZswjH9e+fXuUKlUK7u7uWL58ORMruYzvvgP69pUh4EqVdEdjf049FJwiMBDo3Bl4/30gLk53NGRmN2/exJ49exAUFHTf/UFBQdixY8cjHzdz5kwcP34cI0aMsHeIRIYSFweMHAl07+4aSRVwkcQKSLOIa9ecd0EyOcbFixeRlJQEX1/f++739fVFdHR0qo85evQohg4dirlz58LDwyNNx0lMTER8fPx9NyIzGj1aGvZ88onuSBzHZRJr4cJSsX7zDXD4sO5oyOwsD2wcqZR66D4ASEpKQseOHTFq1CiULl06zc8/ZswYWK3WO7eiRYtmOmYiRztyBPj6a2DwYOdtBpEal7jGmuLGDcDfHyhb1rn2/iPHuXnzJrJly4ZFixahdevWd+5/9913ERkZic2bN9/387GxscidOzfc75k1l5ycDKUU3N3dsX79etSvX/+h4yQmJiIxMfHO3+Pj41G0aFFeYyXTUApo1kx2GouKArJm1R2R46RtXMpJeHsD48cDbdoAq1fLL50oPTw9PVGtWjVs2LDhvsS6YcMGtGzZ8qGfz5kzJ/7888/77gsNDcXGjRuxePFilChRItXjeHl5wYv7HpKJrVoFrF0LLFvmWkkVcLHECgCtWgH16wP9+wMNGwKenrojIrMZMGAAgoODUb16ddSqVQvff/89Tp06hbfffhsAMGzYMPz777+YPXs23NzcUKFChfseX6BAAXh7ez90P5GzuHED6NcPCAoCUvm+6fRcLrFaLNIAunJl4KuvgKFDdUdEZtOuXTtcunQJH3/8Mc6dO4cKFSpgzZo18PPzAwCcO3cOp7j7A7mw8eNlA5TVq+Uz19W41DXWew0cCEyaJOP///95SGRYXMdKZnHqlMxj6d0b+OIL3dHo4bKJNSFBfvkBAcDy5bqjIXo8JlYyi9deA7Ztk9UXrvpSdZnlNg/y8ZH9WlesAH7+WXc0RETmt2YNsHixDAW7alIFXLhiBWQ6eJMm8s3qr7+AbNl0R0SUOlasZHTXrgEVKgAlS8qmJ654bTWFy1asgPziv/sOOHcO+Owz3dEQEZnXZ58BZ88CoaGunVQBF0+sAPDMMzIz+IsvgEOHdEdDRGQ+UVHyGfr++0CpUrqj0c+lh4JT3LghQxhFiwIbN/LbFhkPh4LJqFI2MD9/Hti/H2BfE1asAKQj06RJwG+/AdOn646GiMg8Zs2SDcwnTWJSTcGK9R6vvw4sWSLDGoUL646G6C5WrGRE0dHSf71ZM2DOHN3RGAcr1nuMHy8zg3v1kuENIiJ6tP/9D/DwkB1s6C4m1nvkzi0z2n7+GVi4UHc0REBISAj8/f0REBCgOxSi+yxdKmtWv/sOyJdPdzTGwqHgVLz2mlxvPXiQLxgyBg4Fk5H8958MAdesKbvXcMLn/VixpuK774CkJODdd3VHQkRkPAMGyGoKrllNHRNrKnx9pd3hvHmypyAREYl164AffpA5KYUK6Y7GmDgU/AhKAU2bAn/+Ke0OrVbdEZEr41AwGUFCgqz5L1UK2LCB1eqjsGJ9BIsFmDwZiIsDhgzRHQ0RkX5DhgAXLwJTpzKpPg4T62P4+QFjxwJTpgC//qo7GiIifdatkyYQn38OlCihOxpj41DwEyQnAw0aAMePy7Awh4RJBw4Fk07//SdDwBUrAmFhrFafhBXrE7i5Scuu2FjOEiYi1/TOO8D168CMGUyqacHEmgZ+fsDEiTITbtky3dEQETnO/PnSMCc0lK1e04pDwWmkFNCqFbBjB3DggCzJIXIUDgWTDmfOyPBvo0bAggW6ozEPVqxpZLEA338vf3bvLtdeiYicVXKybEySLZtUq5R2TKzp4Osr11vXrgW++UZ3NERE9jNpkqxVnTEDyJNHdzTmwsSaTk2bAv37A4MHA3v36o6GiMj2Dh4EBg2SSUuNGumOxnx4jTUDEhOBWrWAK1ckuebIoTsicna8xkqOcuMGEBgI3LwJ7N4NZM+uOyLzYcWaAV5eciH/7FnZj5DIXrhtHDna4MHA4cMyE5hJNWNYsWbCDz8A3boBc+cCHTvqjoacGStWcoSVK4GWLWWHr969dUdjXkysmaAU0LmzbIy+bx9QsqTuiMhZMbGSvf37L1CpEvDcc8Dy5WwEkRlMrJkUHw9UqQLkzQts2wZ4euqOiJwREyvZ0+3bwEsvAceOAX/8IZ9nlHG8xppJOXNKZ5J9+4CPPtIdDRFR+o0YIYXBvHlMqrbAxGoDNWoAn30GjBsHrF+vOxoiorRbuxYYPVo+w+rW1R2Nc+BQsI0kJ8sa1z17ZAlO0aK6IyJnwqFgsodTp+RSVq1aMnHJjaWWTfCf0Ubc3IAff5T2X6++KmtdiYiM6uZNoF07WYf/ww9MqrbEf0obypcPWLJELv5zizkiMrLBg2WE7aefeF3V1phYbax6dSAkBJgyBZg5U3c0REQPmzdPtsIcP166LJFt8RqrnfToAcyZI9vMVa2qOxoyO15jJVv54w+5pvrqqzIEzPWqtsfEaic3bshC60uXpN8mh1ooM5hYyRb++09G1axWYPt2mRNCtsehYDvx9pbrrfHxQKdOQFKS7oiIyJUlJUnr1bg4YOlSJlV7YmK1Iz8/aR6xfj0wfLjuaIjIlb3/vuyvOn8+UKKE7micGxOrnQUFAV98Ic0jZs/WHQ0RuaJZs4DPPwe+/FI+k8i+eI3VAZQC3nhDdsHZtAmoXVt3RGQ2vMZKGbV9O1C/PhAcDEydyslKjsDE6iCJidLk+sgRYNcuGSYmepKQkBCEhIQgKSkJR44cYWKldPnnH2m5Wq6cDANzkxDHYGJ1oAsXgICAuzPycuTQHRGZBStWSq+EBBkdu3pVvszny6c7ItfBa6wOlD+/7N36998yLJOcrDsiInJGSUmyGuHkSfnMYVJ1LCZWB6tYUbqerFgBDBumOxoickbvvw+sXg0sWACUL687GtfDxKrByy/L7LzPPwcmTdIdDRE5k8mT784AbtpUdzSuyUN3AK6qf38ZpunTByhcGGjRQndERGR2K1YAvXsDffsC/frpjsZ1cfKSRklJwGuvAWFhwObNMrGJKDWcvERPEh4ONGggVerChYC7u+6IXBcTq2bXr8ub4fhxeWM8/bTuiMiImFjpcY4ckRnA/v7S6c3bW3dEro3XWDXLmhVYuVKW4DRpAly8qDsiIjKT6GigcWOgQAFg+XImVSNgYjWAfPmAtWuB2FgZxklI0B0REZnBlStA8+bSgCYsDMiTR3dEBDCxGkbJkvLGOHQIaNNG3ihkXKGhoShRogS8vb1RrVo1bN269ZE/u3TpUjRs2BD58+dHzpw5UatWLaxbt86B0ZIzunEDaN1ahoHXrAGKFdMdEaVgYjWQKlVkWHjLFmkgwa3mjGnhwoXo168fhg8fjn379uH5559HkyZNcOrUqVR/fsuWLWjYsCHWrFmDPXv24MUXX8TLL7+Mffv2OThycha3bgHt2gHbtkkDiGef1R0R3YuTlwxo+XLg1VcluU6fDrjx64+hBAYGomrVqph0zyLkcuXKoVWrVhgzZkyanqN8+fJo164dPvroozT9PCcvUYqkJKBzZ9nveeVKub5KxsKPbANq1Qr44Qe59e0ru+OQMdy8eRN79uxB0AN7bwUFBWHHjh1peo7k5GQkJCQgDy+IUTolJwNvvQUsWiRLaphUjYkNIgyqUyfg2jWgZ08ge3Zg7Fhu92QEFy9eRFJSEnx9fe+739fXF9HR0Wl6jvHjx+Pq1ato27btI38mMTERifdcaI+Pj89YwOQ0lJLGMjNmyN7OrVvrjogehYnVwHr0kJ0p+veXnXA+/FB3RJTC8sC3HKXUQ/elZv78+Rg5ciRWrFiBAgUKPPLnxowZg1GjRmU6TnIeH34IfPONtCzs3Fl3NPQ4HAo2uH79gE8/BT76CBg/Xnc0lC9fPri7uz9UncbExDxUxT5o4cKFeOONN/DTTz/hpZdeeuzPDhs2DHFxcXdup0+fznTsZF5jxwKffSafAW+9pTsaehJWrCYwfLgMC7/3ngwHDxigOyLX5enpiWrVqmHDhg1ofc9Y3IYNG9CyZctHPm7+/Pl4/fXXMX/+fDRr1uyJx/Hy8oKXl5dNYiZzmzBBdsIaOZLvfbNgYjWJTz+VaywDB8rf+QbTZ8CAAQgODkb16tVRq1YtfP/99zh16hTefvttAFJt/vvvv5g9ezYASapdunTBxIkTUbNmzTvVbtasWWG1WrWdBxnfF18AgwcDQ4fKqBWZhCLTSE5WatgwpQClvvxSdzSuLSQkRPn5+SlPT09VtWpVtXnz5jv/r2vXrqpevXp3/l6vXj0F4KFb165d03y8uLg4BUDFxcXZ8CzIyD79VN7rH34o730yD65jNRmlgA8+AEaPlusuQ4bojogcgetYXYdSwKhRcvv4Y05aNCMOBZuMxSLDwu7uMjx09aq8AbkUh8j87v3iPGaMvMfJfJhYTchikW+yOXJIxXrliswWZHIlMi+l5Hrql1/K+5nzKMyLidXEBg+W5Nq7tyTXSZO4uTGRGd2+DfTqBUybJmtV//c/3RFRZjCxmtw770hnptdfl+Q6axbg6ak7KiJKqxs3gI4dpe/vDz8AXbrojogyi4nVCXTtKsm1UyfgwgVg6VLAx0d3VET0JPHxQMuWQESEbL7RvLnuiMgW2HnJSbz6KrBuHbBrF/DCC8D587ojIqLHOX9e3quRkcCGDUyqzoSJ1Ym88AKwdStw7hxQuzZw9KjuiIgoNSdOAM89B0RHA5s3y3+T82BidTKVKgHh4UCWLJJcf/9dd0REdK89e4A6dWQW8Pbt8p4l58LE6oT8/OQN+8wzUsWuXas7IiICZIJS3bpAkSLyHi1RQndEZA9MrE4qb17g11+B+vWBl18GvvtOd0SUESEhIfD390dAQIDuUCgTlJJm+q1aAU2aAL/9BjxhMyQyMbY0dHK3b8uuOBMnytKciRMBD84FNx22NDSv27eBd98FQkNl7fmYMYAbSxqnxo9YJ+fhId+Uy5WTRhLHjgELFwK5cumOjMj5JSQA7doB69cD338P9OihOyJyBH5vchFvvXV3OU6tWsDx47ojInJuR48CNWvKtdS1a5lUXQkTqwtp0EAWoiclATVqyNo5IrK9tWuBgAAZBt65E2jYUHdE5EhMrC6mTBlJrgEBQOPGsvUcr7IT2YZS8p5q1kzWpu7aBZQtqzsqcjQmVheUJw+wejUwbJjcXnlFWqsRUcZdvQq0by/vqeHDZWmN1ao7KtKBs4Jd3IoVQHAwUKgQsGyZTHIi4+GsYGM7fFjaip44AcyeDbRpozsi0okVq4tr2VK6M7m7y3XXxYt1R0RkLvPnA9Wr372eyqRKTKx057prkybAa6/Jetfr13VHRWRsN27IHqodO979glq+vO6oyAiYWAmAbDO3cKEsYp8xAwgMBA4e1B0VkTEdOya9uGfOlPWpc+YAOXLojoqMgomV7rBY5Bv4rl0yrFW9uiRZXoUnEkrJZuRVqkjzh4gIWZ9qseiOjIyEiZUeUqmSDGt16AC88YZsoM5Zw+TqYmPlPdGtm0xU2rsXqFxZc1BkSEyslKrs2YFp02RixqpVkmx/+013VER6bN0KPPssEBYGLFggQ8A+PrqjIqNiYqXHat8e+OMPoHhx4MUXgf79ObGJXMe1a/Kar1cPKFZM3gvt2umOioyOiZWeqEQJYONG4KuvgEmTgKpVuYG6o3DbOH22bZMqdfJk4MsvZcTGz093VGQGbBBB6XLwINClC7BvHzB0KPDBB4C3t+6onB8bRDjOtWvSOWniRGmiP3OmLEkjSitWrJQu5coBO3YAI0YAX3zBa6/kXLZvlwlJkyfL63vrViZVSj8mVkq3LFmADz8EIiMBX1+59vrmm8B//+mOjChj4uKAPn2A558H8uaV1/bAgdKRjCi9mFgpw8qVAzZvBqZMkVaI5crJjEleXCCzUEpes2XLyvrU8ePl2iqrVMoMJlbKFDc3oGdPufZar56s82vWDPj7b92RET3e0aNAo0bymq1TR17D/fuzSqXMY2IlmyhYEPjpJ9kq68ABqV6HDZPuNERGcvUq8NFHQMWKklxXrZIRlyJFdEdGzoKJlWzq5ZeBQ4ckqU6cCJQuDcyaBSQn646MXJ1SwNy5Msw7bhwwYADw118ywkJkS0ysZHPZsgEjR0qCfeEFoHt3aeq/Y4fuyMhV7dolTfM7d767wcTo0fJaJbI1Jlaym2LFpCXi1q1SsdapI1tsnTypOzJyFceP302m164BmzYBS5YATz+tOzJyZkysZHfPPSedmqZPlw5OpUoB//sfcO6c7sjIWZ09Kzs1lS0rr7kpU6Rp/gsv6I6MXAETKzmEmxvw+uuyj+WoUcCPPwIlSwKDBwOXLumOjpzFpUvymipZUibTjR4tr7mePTnblxyHLQ1Ji9hYWTP49deSdAcMkBu79aWOLQ0f78oVeS19+SWQlCSvpYEDAatVd2TkiphYSasLF4CxY4GQENmqbuBA4J13gFy5dEdmLEysqfvvP+C772QG+pUr8toZNgwoUEB3ZOTKmFjJEP79F/jsM2DGDMDTE3jrLVmsX6iQ7siMgYn1fufOSYU6aRJw+7a01Bw0SCbMEenGxEqGEh0t1UdoKHDjhuykM2iQrId1ZUys4tgx2b5wxgzAywvo3Rvo148VKhkLJy+RoTz1FDBmDHDqFPDJJ9IVp2xZ4JVXZBcdV/sayP1YZanWunXSyKF0aemSNGKEvEZGj2ZSJeNhxUqGduMGMGeOVCmHDgHly0uV0rkz4OOjOzrHccWKNT5eGuN/9x1w5AhQpYos02rfHsiaVXd0RI/GipUMzdsb6NEDiIoCfvlFKpY+fYDCheVD9tAh3RGSrR06JL/bwoXlOnvlytJkZM8e6eLFpEpGx4qVTOfUKVnwP3WqzCp+7jmgWzegbVvnrWKdvWJNSJB1pzNnymbj+fPLBLa335YES2QmTKxkWomJwNKl0uR/wwapZF55RZLsCy/I+lhn4YyJNSlJKtFZs4BFi4Dr14GGDaUqbdVKRiuIzIiJlZzC6dNyLXbWLNkKrGhRqWDbtgUCAgCLRXeEmeMsiTU5WSrSRYtkEtK5c9K3t3t3mQHO5TLkDJhYyakoBYSHS8vEJUuAmBjAzw947TVJstWrmzPJmjmxJifL7+SnnySZnj0rw7uvvSa3mjWda3SBiImVnFZSErBli3ygL1ki12OLFpVlG82bA/Xrm2cijNkS67VrsjxqzRpg+XJpAFKo0N1kWqsWkyk5LyZWcgm3b0uSXbkS+Pln4O+/JanWry+Jtn59mXFs1GrWDIn16FFJpGvXSlJNTARKlABeflmSae3aTKbkGphYyeUoBRw+LM0nVq0Ctm2T6rZQIZn0VL8+8OKLkhSMkmiNlliVki8nW7bIbfNm4MQJaUdZrx7QpAnQtKmxv6wQ2QsTK7m8hARJrps2yW3vXrkuWKwYULeubJJdowbw7LPSRk8H3Yk1MRE4cADYtetuMj17VpJm5crA88/LjN4XX5TNFIhcGRMr0QNiY2UZyKZNMoM1MhK4eRPIkkWSSEAAUK2adIHy93fM2llHJtZr16RJw969wO7dctu/H7h1C/DwkAlgdevKrU4d7kRE9CAmVqInSEyUxLJr193b4cN3+xYXLSpJNiXRlisnG23nz2+7YVBbJ1alZDLX4cPAwYOSSFP+PHlS/r+bm5xP9epyq1ZNqnazTPgi0oWJlSgDUqq6v/6SW1SU/HnixN2EmyOHLPUpVAgoWFD+vPeWP79s7J4zpwwxPy4JpyWxKiW9lePi5Hb5sszG/fdf4MyZu3+m/HdiojzOzU2+CJQtK18KUv6sWJHDukQZkabEqpRCQkKCI+IhMrWrV4Hjx6XqO3FCklh0tDRCSPnz1q2HH+fmJhN/smSR4VY3NwWl1D27+cTj8mU/5Mt3Gm5uOZGcfDeBJyXJUPWNG3Jt+EGenrJutGBB+TMlsRcuDDzzjDRo0HXtmMhsfHx8YHnCUFSaEmvKt2UiIiJXlpbLMXapWAMCAvD777+n+edt8dj4+HgULVoUp0+fzvA1KB1xZ/axmT1vXXFn5vFGOOdr12SWbFTU3WuThw5JVZrC2xsoXlyuwT71FODrC/zww+f48svB8PWVfUTz5pXJT+7u9x8nMTERiSljtQDOnTuHGjVqICoqCoUf05U+MRG4ckWGgi9dAs6fl+5Tn3wyBS1bvoXTp6WSPnXq/sq5YEEZAk65+fsDFSoA2bLZ7t8svXT+njP7eH6OpY8R3tNplZaK1SMtT2SxWNJ1su7u7hl+UWTmsQCQM2dOLcfWec5Axs9bZ9y6ftfpPe7NmzJ5afdu4Pz50XjuuZz46y8ZdnVzk+HU8uWBN96Q/y5ZUm5PPfVwQ4TVq5eiS5dP0x1zCh8fnyfGnj//w/d9//1MzJgx6M7fk5NlucyJE7Ie9ehRuUb8yy+yc1BysiR7f3+ZtBQXF4yoqJwZmrxklt+zLR/Pz7GMMePnWGrSlFjTq3fv3loem1m64nbFc7bF4+113IsX767V3L5dkurNm3Lt86mn6qJmTdkvtFo1STzp2YXFKOfs5gYUKSK355+//2dv3JBKfM8e+TKxZw9w8eIw1KolybZqVdmqL+VWoED6ju0oOl+frviedsVzfhSnmRWsewG9Lq543rY+55S+tuvWScUWFSX3lyghiaNGDVluonOpyZkzZ+4MlRUpUsThx09pEPH77/JlY+tWmaAFAJUqAY0aAUFB8u9lq+3e+Np2jXMGnO+87VKx6uDl5YURI0bAy8WmN7rieWf2nJWSYc+wMEmmW7dK4vDzk+5B778vzQ+KFrVx4JmQcq66fs9eXlKhV6smm48DslXfli3A+vWyZd8XX8gXj3r17ibacuUyvpaXr23X4Wzn7TQVK9HjpGxdtmSJbI5+8qQkgRdekCTQqBFQpoxx+9oa/Ru9UsCff8oXlfXr735ZKVpUNi1v3VqGnD2c5qs80aMxsZLTun1bmsMvWQIsWyYzdgsWlA/5li2lKrXVsKW9GT2xPujaNalmV6+WbePOnJEZ0C1bAh06SE/hB2dCEzkLJlZyKkrJdcDZs4EFC2TJiZ8f0KYN8Mor5t0H1GyJ9V5KySSopUtlo/Njx+QLTvv2QOfOQJUqxh0pIMoIJlZyCqdPAz/+KAn10CHpLNSpE9C2rVwXNPsHt5kT671SvvjMnStffGJiZN1sp05yK1FCd4REmcfESqaVmChV0PTpwMaNMqzbpg3QpQvQoIFzDTU6S2K91+3bwK+/SpJdulTaQTZsCPTsCbRoIa0YiczIhINid12+fBnBwcGwWq2wWq0IDg5GbGzsEx938OBBtGjRAlarFT4+PqhZsyZOnTpl/4BtIKPnnOKtt96CxWLBhAkT7BajrT14zq1bD0C/fjdQtCjQsaN0EJo+Xa6h/vgj8OKLt/D++0NQsWJFZM+eHYUKFUKXLl1w9uxZ3adC9/DwAI4fD8XWrSVw61Ye+Pl9hOjoOLz2mkx6GjpU+i6nWLp0KRo2bIj8+fMjZ86cqFWrFtatW6fvBDIoNDQUJUqUgLe3N6pVq4atW7em6XHbt2+Hh4cHKleubN8A7SC955yYmIjhw4fDz88PXl5eKFmyJGbMmOGgaG1AmVjjxo1VhQoV1I4dO9SOHTtUhQoVVPPmzR/7mGPHjqk8efKoQYMGqb1796rjx4+rVatWqfPnzzso6szJyDmnWLZsmXr22WdVoUKF1Ndff23fQG2ocePGqnz5Cuqrr6JUrVr/KSBJeXgkqHffVergwYd/PjY2Vr300ktq4cKF6tChQyo8PFwFBgaqatWqOTx2W4mLi1MAVFxcnO5QbGbBggUqS5YsaurUqSoqKkq9++67Knv27GrDhrOqb1+lcuVSClCqeXOlfvlFqb5931Xjxo1Tu3btUkeOHFHDhg1TWbJkUXv37tV9Kmn2qHM+efLkYx8XGxurnn76aRUUFKSeffZZxwRrIxk55xYtWqjAwEC1YcMGdeLECbVz5061fft2B0adOaZNrFFRUQqAioiIuHNfeHi4AqAOHTr0yMe1a9dOde7c2REh2lxGz1kppc6cOaMKFy6sDhw4oPz8/EyTWPftO6iAN1SJElcVoFSVKkoNG3ZMAdmeeM732rVrlwLwxA8wo3LGxFqjRg319ttv33df2bJl1dChQ5VSSl27ptSMGUpVqiQJtkIFpaZNk/tT+Pv7q1GjRjky7Ex50jk/Srt27dQHH3ygRowYYbrEmt5zXrt2rbJarerSpUuOCM8uTDsUHB4eDqvVisDAwDv31axZE1arFTt27Ej1McnJyVi9ejVKly6NRo0aoUCBAggMDMTy5csdFHXmZOScATnv4OBgDBo0COXLl3dEqJn233/AqFFA3bp+AL5HpUrZsHmztNcbPbokrNYsjz3nB8XFxcFisSBXrlx2i5nS7ubNm9izZw+CgoLuuz8oKOjO7zVrVqB7dyAyUq6hlygB9OgBFCsmr42LF5ORkJCAPHnyaDiD9EvLOadm5syZOH78OEaMGGHvEG0uI+e8cuVKVK9eHZ9//jkKFy6M0qVL47333sP169cdEbJNmDaxRkdHo0AqTUoLFCiA6Hu3GLlHTEwMrly5grFjx6Jx48ZYv349WrdujTZt2mDz5s32DjnTMnLOADBu3Dh4eHigb9++9gzPJi5eBIYPl11hxo0D/P3/gp9fIyxfLutOU2b3Pumc73Xjxg0MHToUHTt2NN3En5CQEPj7+yMgIEB3KDZ18eJFJCUlwdfX9777fX19H/q9Wiyy7nXlSuDwYaBdO2DsWKBw4duIiRmK559v78jQMyw955zi6NGjGDp0KObOnQsPE3bXyMg5//3339i2bRsOHDiAZcuWYcKECVi8eLHWXsTpZbjEOnLkSFgslsfedu/eDQCpbt2jlHrklj7J/78LdMuWLdG/f39UrlwZQ4cORfPmzTF58mT7ndQT2POc9+zZg4kTJ2LWrFlP3OrIkR4+5wKwWMYhf/4rGD36ChISxmH58ki0aLEeXl4PTyx73Dnf69atW2jfvj2Sk5MRGhpqj1Oxq969eyMqKipT258Z2YO/wyf9XkuVAr77Dvj666VITv4aHh49UKNGPvzvf7KJvBmk9ZyTkpLQsWNHjBo1CqVLl3ZUeHaRnt9zcnIyLBYL5s6dixo1aqBp06b46quvMGvWLNNUrYb7CtSnTx+0b//4b6DFixfH/v37cf78+Yf+34ULFx76dpQiX7588PDwgL+//333lytXDtu2bct40Jlkz3PeunUrYmJiUKxYsTv3JSUlYeDAgZgwYQL++eefTMWeUSnnfOGCO2bMyIsFC3LB3V2hc+fL6Nr1MnLnbonixYvjzJm96T7nFLdu3ULbtm1x4sQJbNy40XTVqjPLly8f3N3dH6paYmJinvh7XbhwIQYM6I7lyxfh+eez4JtvgPHjgWnTgHfeAYYMefKOOzqk95wTEhKwe/du7Nu3D3369AEgSUcpBQ8PD6xfvx7169d3SOwZlZHfc8GCBVG4cGFYrdY795UrVw5KKZw5cwalSpWya8w2oe/ybuakTOTZuXPnnfsiIiKeOJGnVq1aD01eatWqlerQoYPdYrWVjJzzxYsX1Z9//nnfrVChQmrIkCHpmvxjaxcuKNW/v1Le3krlzKnUhx8qldpchYz+nm/evKlatWqlypcvr2JiYuxxCg7lrJOXevXqdd995cqVe+xEnnnz5ilvb2+1bNmy++6/fFmpjz5SysdHqezZlRo6VKmLF+0QdCal55yTkpIeeu/26tVLlSlTRv3555/qypUrjgo7U9L7e54yZYrKmjWrSkhIuHPf8uXLlZubm7p278w1AzNtYlVKlmFUqlRJhYeHq/DwcFWxYsWHlp6UKVNGLV269M7fly5dqrJkyaK+//57dfToUfXtt98qd3d3tXXrVkeHnyEZOecH6ZwVfPWqUp99JsnUx0epkSPlQ/Fx0nvOt27dUi1atFBFihRRkZGR6ty5c3duiYmJdjoz+3LGxJqyDGP69OkqKipK9evXT2XPnl39888/Simlhg4dqoKDg+/8/Lx585SHh4cKCQm573caGxt752cuXlRq2DBJrj4+So0apZSR8k96z/lBZpwVnN5zTkhIUEWKFFGvvvqq+uuvv9TmzZtVqVKl1JtvvqnrFNLN1In10qVLqlOnTsrHx0f5+PioTp06qcsPfEoDUDNnzrzvvunTp6tnnnlGeXt7q2effVYtX77ccUFnUkbP+V46Euvt27JUolAhpbJkUapvX6XSWkim95xPnDihAKR627Rpk03Py1GcMbEqpVRISIjy8/NTnp6eqmrVqmrz5s13/l/Xrl1VvXr17vy9Xr16qf5Ou3bt+tDznj+v1IABSnl6KlW4sFI//KBUUpIDTigN0nPODzJjYlUq/ed88OBB9dJLL6msWbOqIkWKqAEDBpimWlVKKbY0JLtbvx4YMED2QG3XDvjsM6BkSd1RmYsztjR0hOPHpYPT4sVA1arAV1/JfrFE9mS4WcHkPI4dk23CGjUC8uQBdu2SxutMquQoJUsCixbJ/rDu7rL/bps28tokshcmVrK5a9eA998HypcH9u0DFi6UfVGdbCkmmchzzwEREdLwf/dueW1+8IG8VolsjYmVbCosDKhQQYbchg2TLdzatjX/tm1kfm5usmnD4cPy2vzyS8DfH1ixQrazI7IVJlayieho2bi6SRPg6aeBP/8ERo4EsmXTHRnR/bJmldfmgQOSWFu1km3qTp/WHRk5CyZWypTkZGDyZNms+tdfgTlzgA0bpEMOkZE98wywerXsBbt3rwwPT5kir2mizGBipQz780+5dtWrF/DqqzLs27kzh33JPCwWoHXruzPW334baNCAk5soc5hYKd1u35YlM9WqAZcvy8SkadOAvHl1R0aUMblyAVOnAr/8Apw8CVSsKNdgb9/WHRmZERMrpcuRI1KlfvQRMHCgbOlVt67uqIhso0EDGYnp1QsYPBioXVv+TpQeTKyUJsnJwLffApUrA5cuybrAMWMALy/dkRHZVvbsMqt9xw7g6lWgenWpXnntldKKiZWe6NQpoGFDoG9f4PXXpUqtXVt3VK7BWfdjNYOaNWVSU9++Ur2+9BJnDlPasKUhPdaPPwK9ewM5cwIzZ8qHCzkeWxrqtXEj0LUrcOWKzBxu21Z3RGRkrFgpVdeuAW+8AQQHAy+/LNeZmFTJVdWvD+zfDwQFyezhLl2A+HjdUZFRMbHSQw4eBGrUAObPB2bMkLWpuXLpjopIr9y5pdf17NnA8uVy7TUyUndUZERMrHSf2bPlAyM5Gfj9d6B7d65LJUphscgozp49MsmpZk1pkMILanQvJlYCIEO/r78u15Fee02SavnyuqMiMqZSpYDwcLlc0qsX0KEDh4bpLiZWQlSUDP0uWCATlGbNkm/jRPRo3t5ASIjs3rRmjbyHDh3SHRUZAROri/vpJ9nOTSmpUrt10x0Rkbm0bStb0bm5SXJduVJ3RKQbE6uLSkqSrbPatZPNyHft4tAvUUaVLg3s3Ckz51u2lN1z2FDCdTGxuqDYWFlC8/nnwBdfyObPHPolyhwfH2DxYuDTT4GPP5bt6OLidEdFOrBBhIs5eFC+UV+8KNdUg4J0R0RpwQYR5rJ6NdCpE/DUU7I0p2xZ3RGRI7FidSHr18vyAE9PuZ7KpEpkH82ayXvM3V2uu65frzsiciQmVhcxdSrQtKnsTBMeDpQsqTsiIudWqhQQESHvuaZNZWtFcg1MrE4uORkYOhTo2VM2cV6xQq4FEZH9+fjILOGePYEePWTCICc1OT8P3QGQ/Vy/Lj1NlyyRbbD69WMXJSJH8/CQ9a4lSwKDBgEnTshacW9v3ZGRvbBidVIxMdI4fPVqYOlSoH9/JlUz4rZxzsFiAQYOBBYtklGjBg1kAiE5J84KdkLHjsnEpGvXgFWrpPcvmRtnBTuPnTuBFi1kmHjtWrkWS86FFauTiYwE6tSRmb87dzKpEhlNYKBMasqSRSY27dunOyKyNSZWJ7J1K1CvHlCsmPy3n5/uiIgoNSVKyHu0WDHghReALVt0R0S2xMTqJFatkuHf6tWBjRuB/Pl1R0REj5Mvn7xXq1cHGjUCfv5Zd0RkK0ysTmDOHGmf1rSpTFbichoic/DxkfdskyZA69ayHzKZHxOryX3zjSyp6dZNdqrhFH4ic/H2lvdut26yH/LEibojosziOlYT++ILYPBgWRs3bhyX0xCZlYeHdEfLnVvWm9+8Ke9rMicmVpMaN046Kn34ITBqFJMqkdlZLLLjlJeXfGFOSpL3OJkPE6sJjR4NDB8OjBgh+z4SkXOwWIBPPpHm/cOGSXIdPlx3VJReTKwm88knwEcfSZX60Ue6oyEiW7NY5P3t7g588IEkV77XzYWTl0wkJZmmJFdyvMuXLyM4OBhWqxVWqxXBwcGIjY195M/funULQ4YMQcWKFZE9e3YUKlQIXbp0wdmzZx0XNJnSRx/JpukpI1PskWceTKwmoNTdN9fo0fItlvTo2LEjIiMjERYWhrCwMERGRiI4OPiRP3/t2jXs3bsXH374Ifbu3YulS5fiyJEjaNGihQOjJrMaPhwYM+bul2omV5NQZHijRikFKDV2rO5IXFtUVJQCoCIiIu7cFx4ergCoQ4cOpfl5du3apQCokydPpvkxcXFxCoCKi4tLV8zkHL74Qj4DPvlEdySUFrzGanATJ0q1+tlnwJAhuqNxbeHh4bBarQgMDLxzX82aNWG1WrFjxw6UKVMmTc8TFxcHi8WCXLlyPfJnEhMTkZiYeOfv8fHxGY6bzO+994DERBmtypFDluSQcTGxGtisWfIGGjRIZgiSXtHR0ShQoMBD9xcoUADR0dFpeo4bN25g6NCh6Nix42N3qRkzZgxGjRqV4VjJ+bz/PnDlimwBmT27bJxOxsRrrAa1dCnwxhvy5mHzB/saOXIkLBbLY2+7d+8GAFhS+UUopVK9/0G3bt1C+/btkZycjNDQ0Mf+7LBhwxAXF3fndvr06YydHDkNi0XmWPTpA7z1FjBvnu6I6FFYsRrQL78AHToAr70GTJrEpGpvffr0Qfv27R/7M8WLF8f+/ftx/vz5h/7fhQsX4Ovr+9jH37p1C23btsWJEyewcePGJ+6p6uXlBS8vrycHTy7FYpHLQ1euSCvTbNmkTzgZCxOrwYSHyxulQQNpyO3urjsi55cvXz7ky5fviT9Xq1YtxMXFYdeuXahRowYAYOfOnYiLi0Pt2rUf+biUpHr06FFs2rQJefPmtVns5Hrc3IBp04CrV4H27YH164G6dXVHRfeyKMUJ3Eaxf7/sp1qxIhAWJt9GyViaNGmCs2fPYsqUKQCAnj17ws/PDz/fs+dX2bJlMWbMGLRu3Rq3b9/GK6+8gr1792LVqlX3VbZ58uSBp6dnmo4bHx8Pq9WKuLi4J1a75BoSE2VHqz17ZG/XihV1R0QpeI3VII4fl/1Un35a9mVkUjWmuXPnomLFiggKCkJQUBAqVaqEOXPm3Pczhw8fRlxcHADgzJkzWLlyJc6cOYPKlSujYMGCd247duzQcQrkJLy8gGXLgOLFZdu5U6d0R0QpWLEawMWLQO3acv1k2zZuUk4PY8VKj3LunHx+ZM0qnx958uiOiFixanb9OtCiBRAXB6xdy6RKROlTsCCwbh1w4QLw8svymUJ6MbFqlJQEdO4MREYCq1bJMDARUXqVLi2fIZGRMqHp9m3dEbk2JlaN3nsPWL4cWLgQCAjQHQ0RmVlgILB4MbB6NTBwoO5oXBsTqyYTJsjt229l+IaIKLOaNJHPlG++ASZP1h2N6+I6Vg1WrgQGDAAGDwbeeUd3NETkTHr1Ag4elA5NpUrJmnhyLM4KdrD9+2UGX+PGwE8/yWJvoifhrGBKj9u3gebNgZ07gYgIII37Q5CNMLE6UEyMXEvNm1cWdGfPrjsiMgsmVkqvuDigVi1JshERXIbjSKyXHCQxEWjTRv5csYJJlYjsy2qVZjP//Sd9x2/d0h2R62BidQClgLffBnbvlqRatKjuiMgsQkJC4O/vjwBOG6cMKFlSdsraulWuuXJ80jE4FOwAX34pe6r++CPQqZPuaMiMOBRMmTFjhmxDOWEC8O67uqNxfpwVbGcbNgBDhgBDhzKpEpEer78uM4UHDAD8/YGGDXVH5NxYsdrRyZNAtWpA9eqyaJtbwFFGsWKlzEpKkpnCv/8O7N0LFCumOyLnxWusdnLjBvDKK4CPDzB3LpMqEenl7i6Xo3LkAF59VSZSkn0wsdpJnz7AX38BS5bI8hoiIt3y5pW2h3/8AfTrpzsa58XEagdTpwLTpwOTJgFVq+qOhojorurVge++k5aHs2frjsY58Rqrje3aBTz/vEwWmDRJdzTkLHiNlWxJKZklPH++NI949lndETkXJlYbunBBJisVKgRs3gx4eemOiJwFEyvZ2vXr0l41IUHW2OfKpTsi58GhYBtJTpa9VW/ckGsYTKpEZGRZs8ockEuXgC5d5DOMbIOJ1UY+/1zWrM6dCxQpojsaIqIne/ppYM4caX04bpzuaJwHE6sNbN8OfPABMGwYF14Tkbk0by6fXx98AGzZojsa58BrrJn0339A5cqy2Pq33wAP9rIiO+A1VrKn27dl39a//wYiI7lEMLNYsWaCUkD37sDVqzK7jkmViMzIw0MuY127JrOFWW5lDhNrJnzzDbByJTBrFnesISJzK1JEmvWvWAGEhuqOxtyYWDMoMhIYPFi6l7z8su5oiIgyr2VL6Ro3cKB0Z6KM4TXWDLh2TbqXeHnJ4mourSF7CQkJQUhICJKSknDkyBFeYyW7u3EDqFlTegnv2QNky6Y7IvNhYs2APn2kZeGePbIFE5G9cfISOdKhQ9KO9fXXpf0hpQ+HgtNpzRogJAT44gsmVSJyTmXLymdcSAgQFqY7GvNhxZoOMTFAxYrStnD1asBi0R0RuQpWrORoSgFNmsi11j//BPLl0x2RebBiTSOlgDfflLZfM2YwqRKRc7NY5LPu5k3grbe4BCc9mFjTaOZMafs1fTrw1FO6oyEisr9ChYDvvweWLuUWc+nBoeA0OHUKqFABaNNG1qwSORqHgkmnbt0kue7fDxQvrjsa42NifQKlgKAg4OBB4MABbq1EejCxkk7x8VJclCkDrF/PS2FPwqHgJ5g8GfjlF7nWwKRKRK4oZ05g6lT5LJw2TXc0xseK9TH+/huoVEn2WZ08WXc05MpYsZIRvPkm8NNPMnpXrJjuaIyLifURkpOBF1+U66v79wM+ProjIlfGxEpGEBcHlC8vw8Jr13JI+FE4FPwIU6bI3oQzZjCpEhEBgNUqs4TXrZOVEpQ6VqypOH1avpW1by8vIiLdWLGSkXTrBixfLkPCRYrojsZ4mFgfoJTsVrNvH/DXX5ywRMbAxEpGcvmyFB9VqgCrVnFI+EEcCn7AggXSrjA0lEmViCg1uXPL5bI1a4Aff9QdjfGwYr3HxYtAuXJA/frAwoW6oyHitnFkbB06yBKcQ4eAvHl1R2McTKz3CA6WavXgQcDXV3c0RHdxKJiMKDpaipHWrWWiJwkOBf+/X3+VIY3x45lUiYjS4qmngHHjZIbwb7/pjsY4WLECSEyURhBPPSUvDl6IJ6NhxUpGlZwM1K0rl9L++APw8tIdkX6sWAF8/rl0WQoNZVIlIkoPNzeZyHT8ODB2rO5ojMHlE+uxY8BnnwHvvSfTx4mIKH3KlwcGDwZGj5aJTK7OpYeClQKaNJEXQlQUkC2b7oiIUsehYDK669eBihWlh/Cvv7r26J9LV6yLF0trrm+/ZVIlIsqMrFnls3TTJvlsdWUuW7HGxwNlywKBgcCyZbqjIXo8VqxkFi1aAJGRMhLoqgWLy1asI0bITg0TJ+qOhIjIeXz9NXD+vGtPZHLJxBoVJUMWH33EPQWJiGypZEmZDJqy2sIVuVxiVQro1w8oUUL+JEqPy5cvIzg4GFarFVarFcHBwYiNjU3z49966y1YLBZMmDDBbjES6fb++0D+/MCAAboj0cPlEuvKlcCGDcBXX3EhM6Vfx44dERkZibCwMISFhSEyMhLBwcFpeuzy5cuxc+dOFCpUyM5REumVPTvw5ZfAihUyQdTVuNTkpRs3ZL1VqVLA2rWuPR2c0u/gwYPw9/dHREQEAgMDAQARERGoVasWDh06hDJlyjzysf/++y8CAwOxbt06NGvWDP369UO/dAyZcPISmY1SwIsvSj/h/fsBT0/dETmOS1WsX38NnDolfzKpUnqFh4fDarXeSaoAULNmTVitVuzYseORj0tOTkZwcDAGDRqE8mnsQpKYmIj4+Pj7bkRmYrEA33wDHD0qf7oSl0msZ89Kh6U+fWQ3BqL0io6ORoECBR66v0CBAoiOjn7k48aNGwcPDw/07ds3zccaM2bMneu4VqsVRYsWzVDMRDpVqgS88w4wahRw7pzuaBzHZRLr0KGypmrECN2RkNGMHDkSFovlsbfdu3cDACypDHUopVK9HwD27NmDiRMnYtasWY/8mdQMGzYMcXFxd26nT5/O2MkRafbxx4C3t3wGuwoP3QE4QkQEMGcO8P33QK5cuqMho+nTpw/at2//2J8pXrw49u/fj/Pnzz/0/y5cuADfR+w1uHXrVsTExKDYPeu6kpKSMHDgQEyYMAH//PNPqo/z8vKCF2fXkRPInVt6CPfsKdXrPVdSnJbTT15SCqhdW/pY7tkDuLvrjojMKmXy0s6dO1GjRg0AwM6dO1GzZs1HTl66dOkSzj0wBtaoUSMEBweje/fuj53wdC9OXiIzS0oCqlYFcuYEtmxx/jkuTl+xLlokFesvvzCpUuaUK1cOjRs3Ro8ePTBlyhQAQM+ePdG8efP7EmTZsmUxZswYtG7dGnnz5kXevHnve54sWbLgqaeeSnNSJTI7d3dpGNG4sSzBadVKd0T25dTXWBMTZVy/WTOgQQPd0ZAzmDt3LipWrIigoCAEBQWhUqVKmDNnzn0/c/jwYcTFxWmKkMiYGjUCGjYEhgwBbt3SHY19OfVQ8Pjx8kv880/OBCZz41AwOYM//gCqVAFCQoBevXRHYz9Om1hjY4GnnwbatgUmT9YdDVHmMLGSs+jWTRr0HDsG+PjojsY+nHYoeNw4GQrm8hoiIuP45BPZtvPzz3VHYj9OmVj//Ve2g+vfHyhYUHc0RESUomhR2QBl/Hj5rHZGTplYR42SZhCDBumOhIiIHjR0qDTq/+gj3ZHYh9Ml1kOHgOnTgeHDAatVdzRERPQgq1WS6syZsj+2s3G6yUtt2gB79wKHD3NbOHIenLxEzubmTaB0aaB6dWDxYt3R2JZTVawREcCyZXJxnEmViMi4PD1lcumSJVIMOROnqlgbNABiYoDISHZZIufCipWc0e3bQIUKQMmSwOrVuqOxHaepWDdvBjZulIlLTKpERMbn4SGf2WvWAI/Z0th0nKZifeEFIC5OGu27Oc3XBXJ1ISEhCAkJQVJSEo4cOcKKlZxOcrJ0Y8qTR4ojZ2jQ7xSJddMmoH59ae7cooXuaIhsj0PB5MxWrgRatgQ2bABeekl3NJln+sSqFFC3rmwL9/vvzvFth+hBTKzkzJQCataUz+/wcPN/jpt+0PSXX4Bt22Sc3uy/DCIiV2SxAJ99BuzcCaxapTuazDN1xaoUUKeObKIbEcHESs6LFSs5O6WAF18ELl8G9u0z91wZE4cOrFsnwwasVomIzM1iAT79FNi/X+bLmJlpK9aUMXl3d2D7diZWcm6sWMlV1K8v237u2WPez3XTVqxr1gC7dgEff2zef3wiIrrfhx/KULCZG0aYsmJVCggIALJmBbZsYWIl58eKlVxFykqPxESZzGTGz3dTVqyrVskwAatVIiLnYrHIzje//y7zaMzIdBVryrVVLy+pVolcAStWciVKAbVrS5I14xwa01WsmzbJtdX339cdCRER2UNK1RoeLm0OzcZ0FetLLwH//WfuGWNE6cWKlVyNUkCNGkC2bLLJipmYqmLdtQv49VepVplUiYicl8UCDB8ul/zMtvONqSrW1q2BqCi5cWs4ciWsWMkVJScD5csDpUubq2mEaSrWv/4Cli8Hhg5lUiXXERISAn9/fwQEBOgOhcjh3NyAIUNk95u//tIdTdqZpmLt0gX47Tfg2DHA01N3NESOxYqVXNXNm0DJktKR6YcfdEeTNqaoWE+cAObNA957j0mViMiVeHoCAwdKDjh5Unc0aWOKxDp+PJA7N/Dmm7ojISIiR3vzTSBnTuCrr3RHkjaGT6wXLwIzZgB9+8q0ayIici05cgD/+x8wdSpw4YLuaJ7M8Il10iT5s1cvvXEQEZE+ffrInyk5wcgMnVhv3AC+/Rbo3h3Il093NEREpEu+fEC3bkBIiOQGIzN0Yp0zR4aC+/fXHQkREenWrx8QEyMTmYzMsMttkpMBf39ZHLxkie5oiPTichsi0aIF8PffwJ9/GrcDn2Er1lWrgMOHZYkNERERAAwYIM0iNmzQHcmjGbZifeEFWRhsth6RRPbAipVIKAVUqwYUKACEhemOJnWGrFgjI2U3A15bJSKie1ksUrWuWwccOKA7mtQZMrF+8w1QtKg03SciIrpX27ZAoULAhAm6I0md4RLrhQsy46t3b8DDQ3c0RERkNJ6e0jDixx+B8+d1R/MwwyXWKVNkR4MePXRHQkRERtWzp+x0ZsSGEYZKrDdvAqGhQHAwkCeP7miI9OO2cUSpy5NHGkaEhgKJibqjuZ+hZgXPmwd06iQXpMuX1x0NkXFwVjDRww4elH4Hc+YAnTvrjuYuQyXWwEDZwcDI65OIdGBiJUrdSy8BV68C4eG6I7nLMEPBERHArl3Au+/qjoSIiMyiTx/JH7t3647kLsMk1m+/lV3imzbVHQkREZlF8+ZAsWLSnN8oDJFYY2KARYtkiY2bISIiIiIz8PAA3n4bmD8fuHRJdzTCEGlsxgyZNt21q+5IiIjIbN58U1odTp+uOxKhPbEmJQGTJwPt23OJDRERpV/+/JJDQkMlp+imPbGGhQEnTwLvvKM7EiIiMqs+fSSXrF6tOxIDLLdp1kxaUhlpRheR0XC5DdGT1agB5M4tDfp10lqxnjgBrF0L9OqlMwoiInIGvXoB69dLbtFJa2KdMkUaQnTooDMKorS7fPkygoODYbVaYbVaERwcjNjY2Cc+7uDBg2jRogWsVit8fHxQs2ZNnDp1yv4BE7mQtm0lp+iexKQtsSYmymzgbt2AbNl0RUGUPh07dkRkZCTCwsIQFhaGyMhIBAcHP/Yxx48fx3PPPYeyZcvit99+wx9//IEPP/wQ3t7eDoqayDVkzy5tcWfMAG7f1heHtmusCxfKLK6oKKBcOR0REKXPwYMH4e/vj4iICAQGBgIAIiIiUKtWLRw6dAhlypRJ9XHt27dHlixZMGfOnAwfm9dYidJm3z6galVg+XKgZUs9MWirWKdNA557jkmVzCM8PBxWq/VOUgWAmjVrwmq1YseOHak+Jjk5GatXr0bp0qXRqFEjFChQAIGBgVi+fLmDoiZyLVWqANWqAVOn6otBS2I9cQL45RdZ1EtkFtHR0ShQoMBD9xcoUADR0dGpPiYmJgZXrlzB2LFj0bhxY6xfvx6tW7dGmzZtsHnz5kceKzExEfHx8ffdiChtevSQibFnzug5vpbEOmOGXGB+9VUdRye638iRI2GxWB572/3/68EsFstDj1dKpXo/IBUrALRs2RL9+/dH5cqVMXToUDRv3hyTJ09+ZExjxoy5M0HKarWiaNGiNjhTItfQoQPg7S25RgcPRx/w9m1g5kygY0e50EykW58+fdC+ffvH/kzx4sWxf/9+nD9//qH/d+HCBfj6+qb6uHz58sHDwwP+/v733V+uXDls27btkccbNmwYBgwYcOfv8fHxTK5EaZQzp8zhmT4dGD5cWuY6ksMTa1gY8O+/HAYm48iXLx/y5cv3xJ+rVasW4uLisGvXLtSoUQMAsHPnTsTFxaF27dqpPsbT0xMBAQE4fPjwffcfOXIEfn5+jzyWl5cXvLy80nEWRHSvnj2lYt2wAWjc2LHHdvhQ8LRpQOXKMmuLyEzKlSuHxo0bo0ePHoiIiEBERAR69OiB5s2b3zcjuGzZsli2bNmdvw8aNAgLFy7E1KlTcezYMXz33Xf4+eef8Q77eBLZTY0aQMWKeiYxOTSxnjsHrFol1eojLkkRGdrcuXNRsWJFBAUFISgoCJUqVXpoGc3hw4cRFxd35++tW7fG5MmT8fnnn6NixYqYNm0alixZgueee87R4RO5DItFJjGtXCltcx16bEeuYx07Fhg1Cjh7Vvo5ElHacB0rUfpdvgwUKgSMHAkMGeK44zqsYlVKhoFfe41JlYiI7C93bll9Mm2a5CBHcVhi3bIFOH4ceOMNRx2RiIhc3RtvAMeOAdu3O+6YDkuss2YBTz8N1K3rqCMSEZGrq1sX8PMDfvjBccd0SGK9ehVYvBjo0oWTloiIyHHc3CT3/PQTcP26g47piIMsXQpcuSInR0RE5EhdugDx8dKY3xEcMiu4YUPg5k3gMa1RiegxOCuYKHOeew7IkUOaFNmb3SvW06eBX38Funa195GIiIhS17WrdGE6e9b+x7J7Yv3xR2mGzIb7RESkS9u2gKen5CR7s2tiVUpmYrVpI02RiYiIdLBagVatJCfZ+wKoXRPrrl3A4cMcBibKqJCQEPj7+yMgIEB3KESm17UrEBUF7Nlj3+PYdfLS//4nM4JPnXL8tj1EzoSTl4gy7/ZtoFgx4JVXgG+/td9x7Fax3r4NLFwoG84yqRIRkW4eHkCnTsCCBcCtW/Y7jt0S66+/AhcuyIbmRERERtChA3DxouQoe7FbYp03DyhTBqhSxV5HICIiSp8qVYDSpaVqtRe7JNbr14Fly6RaZQtDIiIyCotFqtZly4AbN+xzDLsk1tWrgYQECZ6IiMhI2reXFodr1tjn+e2SWOfNAwICgFKl7PHsREREGVe2LFC5MjB/vn2e3+aJNTZWKlZOWiIiIqPq0AFYtUoqV1uzeWJdulSmMbdrZ+tnJiIiso327eUa68qVtn9umyfWefOAF18ECha09TMTERHZRrFiQJ069hkOtmliPXcO2LiRw8BERGR8HToA69cDly7Z9nltmlgXLgSyZJGm+0REREb26qtAcjKweLFtn9emiXX+fKBJEyB3bls+KxERke35+gINGti+WYTNEuuJE7KbDdeuEhGRWXToAGzeDPz7r+2e02aJddEiIGtWoFkzWz0jEXHbOCL7at1aLmH+9JPtntNm28ZVrw6UKCEJlohsi9vGEdlPq1bA2bMy6moLNqlYY2KAAweAtm1t8WxERESO06EDsH+/5DJbsFnFGhcHeHkB3t62eDYiuhcrViL7SUyUZhFWq22ez8M2T2O7gIiIiBzJy0tutmK3/ViJiIhcERMrERGRDTGxEhER2RATKxERkQ0xsRIREdkQEysREZENMbESERHZkM0aRBCR/SilkJCQAB8fH1gsFt3hENFjMLESERHZEIeCiYiIbIiJlYiIyIaYWImIiGyIiZWIiMiGmFiJiIhsiImViIjIhphYiYiIbOj/AFBzCaU7XgnyAAAAAElFTkSuQmCC", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p0 = polar_plot(rt0, (t, -pi/4, pi/4))\n", "p1 = polar_plot(rt1, (t, -pi/4, pi/4))\n", "(p0+p1).show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the plot looks good with the default number of plotting points." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us try 100 points total for the plots." ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAHWCAYAAADKGqhaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7Y0lEQVR4nO3debjN5f7/8ddmyzyURpIi54jiIAqZihQlVOxolqQUpTTqnNJcSpJEMoQkYygi8zxGplIZMpPMw57W74/3b3+pZNg+a92fz2c9H9fl2nWd6+z1Dnu91vv+3Pf7TohEIhEBAABPZHFdAAAAYUKwAgDgIYIVAAAPEawAAHiIYAUAwEMEKwAAHiJYAQDwEMEKAICHCFYAADxEsAIA4CGCFQAADxGsAAB4iGAFAMBDBCsAAB4iWAEA8BDBCgCAhzwL1l27pEOHvPpuAADExqFDlmFe8SRYt26VzjtP+uorL74bAACxM2qUdP75lmVe8CRYzztPKlNGGjzYi+8GAEDsDB5sGXbeed58P8+WgpOSpK+/lnbv9uo7AgAQXbt3W3YlJXn3PT0L1iZNpMOHraUGACAIRo6UkpMtw7ySEIlEIl59s2rVpHz5pLFjvfqOAABEz403Svv3S9Omefc9PT1uk5Qkffut9PvvXn5XAAC8t2OHNGGCt8vAksfBetttUnq6NHy4l98VAADvDRsmRSKWXV7ydClYkurUsXD97jsvvysAAN6qVUvKls1WWr3k+eSlpCRp8mRp82avvzMAAN7YtEmaOtX7ZWApCsHaqJGUmCgNHer1dwYAwBtffmlZ1bix99/b86VgSbr5ZmnnTmnmTK+/MwAAp69yZencc6NzRDQqQ/iTkqRZs6R166Lx3QEAyLw1a6Q5c6KzDCxFKVgbNJBy5pQ+/zwa3x0AgMwbPNgy6uabo/P9o7IULEnNmklLlkjLlkkJCdF4BQAATk0kIpUuLZUrJw0cGJ3XiNp9rHffLa1YIS1eHK1XAADg1CxaJK1cKd11V/ReI2rBWru23RTQv3+0XgEAgFPTv79dEVe7dvReI2rBmphoy8Gffy6lpETrVQAAODkpKZZJzZpZRkVL1IJVsuXgbdu8n2oBAMCpGj9e2r7dsimaorZ5SbKHxGXK2INiLkEHALjUtKk9X126NLqvE9WONSHBPhmMHMkF6AAAd3btsmEQ0e5WpSgHq2Rr2cnJjDgEALgzdKg9Y23WLPqvFdWl4Ax16th/0JQp0X4lAAD+rkYNKUcOe84abVHvWCVrvadOldaujcWrAQBwxNq10rRpsVkGlmIUrI0aSblzc6YVABB7n31mGdSwYWxeLybBmieP1KSJ1KePXYIOAEAspKdb9jRpYuEaCzEJVklq0cLa8cmTY/WKAIB4N2WK3WbTokXsXjMmm5ckO9N62WVS+fLSoEGxeEUAQLxr3lxauNDOr8bqQpiYdawJCdL990vDh0t//BGrVwUAxKs//pCGDbNuNZa3rMUsWCXbkZWaSscKAIi+gQOltLTY7QbOELOl4AwNG0rr19vVPQAAREu5ctIll9hKaSzFtGOVbDl48WLuaQUARM+iRdL338d201KGmAdrvXp2F96nn8b6lQEA8aJ3b6lQIalu3di/dsyDNTFRuuceW/s+dCjWrw4ACLuDBy1j7r03uveu/pOYB6tky8F//MFgfgCA94YNsxvV7r/fzevHfPNShjp1pP37pVmzXLw6ACCsqlWTsmWTJk1y8/pOOlZJevhhafZsNjEBALyzbJk0Y4bUurW7GpwF6803S4ULSx995KoCAEDY9OhhG2RjNXD/WJwFa2Ki9OCD9oB51y5XVQAAwmLfPrtF7YEHbCnYFWfBKtl/fHIy18kBAE7foEG2d6dlS7d1ONu8lKFJE+mHH6QVK2I7yxEAEB6RiF3yUqSI9NVXbmtx2rFKtolp1Sq72gcAgMyYN88mLbnctJTBeccaiUilS9uvL790WQkAIKjuvVeaNk36+Wcpi+OW0XnHmpBgXeuIEdKmTa6rAQAEzc6d0hdfSK1auQ9VyQfBKkl33SXlyCH16uW6EgBA0PTta9fDuZq09FfOl4IzPPSQNHq0tHat223SAIDgSEuTSpSQqlSRBgxwXY3xRccq2QPnTZvc7+YCAATHmDHSmjVS27auKznCNx2rJF1zjS0JT5zouhIAQBBce63dlOanufO+6Vgl28T03Xd2/AYAgONZulSaPNlf3arks2C99Vbp3HOlbt1cVwIA8LuuXW3mfOPGriv5M18Fa/bs1rX26WPbpwEAOJbt222z0iOP+G/Dq6+CVbJgTU+3GwoAADiWnj1tDsKDD7qu5O98tXkpQ6tWtjt47VrrYgEAyJCSIl18sVS/vgWs3/iuY5Wkxx+XtmyxmwoAADja0KF2PPOxx1xXcmy+7FglqUED6ddf7eYbbr0BAGS4+mopTx7/Hs30ZccqSe3bS8uXS+PHu64EAOAXc+ZIc+f674jN0XzbsUYiUsWK0plnShMmuK4GAOAHjRtb07VypT8G7h+LT8uy5d/27a3VX7LEdTUAANd+/FEaOVJ66in/hqrk445Vsp1fxYtLNWtK/fu7rgYA4NKDD9plLWvW2Phbv/Jx5tuh33btpM8/lzZudF0NAMCVLVukfv3s2aqfQ1XyebBK0gMPSLlySR984LoSAIArXbtKZ5xhV4z6ne+DNV8+qWVLm8S0d6/ragAAsbZ3r9S9uw0PKlDAdTUn5vtglewQ8L590qefuq4EABBrvXpJ+/fbo8Eg8PXmpaM1by7NnCmtXu2/gcsAgOhITrZNrNddJ/Xt67qakxOIjlWSnn5aWreOMYcAEE8GD5Y2bJCefNJ1JScvMB2rJN1yi12CvmKFlDWr62oAANGUni6VKWMD98eMcV3NyQtMxypJzz8v/fSTDWAGAITbiBE2Zem551xXcmoC1bFKUt26dqvBkiX+nrwBAMi8SEQqX14qWNC/w/b/SeCi6YUXpGXLbPoGACCcxoyRvv9e6tjRdSWnLnAdqyTVqCEdOCDNm8eVcgAQNpGIVKmSlDOnNG2a62pOXeA6Vsm61gULpG+/dV0JAMBr48fbe/yLL7quJHMC2bFGInbRbbZs0vTpdK0AEBaRiHTNNbYjeNasYL6/B7JjTUiwrnXmzGAuEwAAjm3yZAvUjh2DGapSQDtWyT7VlCsnnX128HaMAQCOrWZNG2E7f35wgzWQHatkv+EdO0rffUfXCgBhMH26NHVqsLtVKcAdq2Rr8FdeKeXJY38YQf6DAIB4FolItWpJu3ZJixcH+/08sB2rZAMiOnWyTzksBwNAcE2YYA3SK68EO1SlgHeskn3KqVLFutc5c4L/BwIA8SYSkSpWlLJnl2bMCP77eKA7Vsn+AF55xYZFBGlIMwDADB8uLVwovfZa8ENVCkHHKv15bX7RImYIA0BQpKVJl18uFS0qjRvnuhpvhCKCEhLsWeuSJfbJBwAQDAMG2HWgr77quhLvhKJjzXDDDdL69dIPP3BfKwD43eHD0r//bac7wnQdaCg61gydOkkrV0qff+66EgDAifTqJf32m713h0moOlZJathQWrrUAjZ7dtfVAACOZf9+qXhx6cYbpT59XFfjrVB1rJLtKlu3TvroI9eVAAD+Sdeu0s6d0n//67oS74WuY5WkBx+0TUw//ywVKOC6GgDA0f74QypWTLrrLgvYsAldxypJ//ufdPCg9OabrisBAPzV229LycnS88+7riQ6QhmshQpJTzwhdeliD8YBAP6wZYv0/vtSu3bSeee5riY6QrkULEl79kiXXirddJP06aeuqwEASPaobuhQ6ZdfpDPPdF1NdISyY5WkfPnsoXjfvnauFQDg1tKlUu/e9rgurKEqhbhjlaSUFKlUKalECenrr11XAwDxKxKR6tSxx3PLlknZsrmuKHpC27FK9gf3+uvSN9/YhegAADfGjrX34XfeCXeoSiHvWCX7lFS1qh1GXrhQSkx0XREAxJeUFOmKK6QLL7R7V8Nwg83xhLpjlewP8P33bW3/k09cVwMA8adHD+mnn6TOncMfqlIcdKwZ7rtPGj1aWr063A/NAcBP/vjDTmg0bmyzgeNB6DvWDK+9ZjcpvPSS60oAIH506mTDIMI2aP944iZYL7hAeuEF6cMPbUA/ACC6Vq+WunWTnn1WOv9819XETtwsBUvWsWYcv/nmm/hY6wcAVxo0kJYssYvMc+Z0XU3sxE3HKtk1cp07S+PHc64VAKJp9Gj79e678RWqUpx1rNKRQ8rr19sh5TPOcF0RAITLwYNS6dK2OjhuXPytDsZVxyrZH3CXLtKvv9oNCwAAb735prRhg/TBB/EXqlIcdqwZnn7azrcuW2ZbwQEAp++XX6xbbd9eevVV19W4EbfBun+/dPnlUvHi8TEJBACiLRKxG8WWLZNWrJBy53ZdkRtxtxScIXduqXt3m105cKDragAg+EaPto2hXbrEb6hKcdyxZkhKsnBdtUoqWNB1NQAQTAcO2HHGyy6zcI3nVcC47VgzdOliA6Kfesp1JQAQXG+8IW3eLHXtGt+hKhGsOv9828HWp480ZYrragAgeFatsvfRp56yIzbxLu6XgiUpPV2qXl3ascOmhGTP7roiAAiG9HSpRg1pyxa7RSzehkEcS9x3rJKUJYv08ce2TfyNN1xXAwDB0auXNGOG1LMnoZqBjvUozz9vt9svWSKVLOm6GgDwt02bbLPSbbdJvXu7rsY/CNajHDx45Jb7yZN5AA8Ax3PrrdatrlwpnXWW62r8g6Xgo+TMaTfdT50q9e3ruhoA8K+RI6Xhw20XMKH6Z3Ssx3DXXXYOa9Uq6ZxzXFcDAP6ye7edWS1XzoZCsLr3Z3Ssx9C5s31t395tHQDgRx06WLh2706oHgvBegznnms333z2mTRxoutqAMA/xo+3HcDvvCNddJHravyJpeB/EIlItWpJa9faLuH8+V1XBABu7dpll5eUKmUBS7d6bHSs/yAhwaYx7dwpPfaY62oAwL22baW9e+1oDaH6zwjW47jkEruot39/6csvXVcDAO589ZW9F77/vlSkiOtq/I2l4BOIRKQmTewGnB9+kAoXdl0RAMTWjh22BFyxogUs3erxEawn4fffbXBE6dL2XCELfT6AONK0qTRhgrR8uXTBBa6r8T8i4iQULGgDIyZOtKVhAIgXQ4bYr+7dCdWTRcd6Ctq2tWH9CxbYsggAhNnWrbZSV6uWhStLwCeHYD0FBw9KV14pZcsmzZ3L9XIAwisSkW6+WZo3z5aAmUJ38lgKPgU5c0oDBkgrVkgvvui6GgCIng8+kMaOtcdghOqpoWPNhDfflJ591m7AqVHDdTUA4K3vv5euukpq3Vrq0sV1NcFDsGZCWpp07bVHpjIVKOC6IgDwxv799sgrRw5pzhweeWUGS8GZkDWrHZTetUtq08Z1NQDgnXbtpPXrpc8/J1Qzi2DNpKJFpQ8/lAYOlAYPdl0NAJy+IUOkTz6xO1ZLlnRdTXCxFHwaIhEpKUn69ltp6VLGfAEIrrVrpf/8R6pb15oFjtZkHsF6mnbulMqUkYoXt7GHiYmuKwKAU5OaahsxN260jUvsGzk9LAWfprPOkgYNkmbOlF54wXU1AHDqXnrJzuYPGkSoeoFg9UD16tLrr9sxnK++cl0NAJy8MWOkV16RXn5ZqlLFdTXhwFKwRyIRqVEjaepUaeFCqVgx1xUBwPGtXm031tSsKQ0fzgUjXiFYPbRrl1Shgi2lzJxp58AAwI/27ZMqV5aSk21sYf78risKDz6feKhAAWnoUJur2bat62oA4NgiEemBB6Q1a6QRIwhVrxGsHitXzs639uxpQyQAwG/ee0/64gupTx+pVCnX1YQPS8FREIlI999vf3HnzrVL0gHAD6ZMkWrXlp54QnrrLdfVhBPBGiUHDtjzi337LFzPPtt1RQDi3YYNUvny9mF//HjO3UcLwRpFa9bYDRElS0oTJjB3E4A7hw/b0cDNm+3kAlfBRQ/PWKPokkukkSNtx12rVrZEDAAuPPaY3cY1fDihGm0Ea5RVqSJ9+qnUr58NkACAWOvd2zZUdu9uV8Ihulhhj4FmzaQff7TL0UuUkG691XVFAOLF/PnSI4/Yqtn997uuJj7wjDVGIhHpjjts5OG0aXxqBBB9GzdKlSrZzVtTp7LPI1YI1hg6eFCqVcsuEZ43T7rwQtcVAQir/ftts9L27fZ+c/75riuKHzxjjaGcOW0zU7Zs0s0321EcAPBaerp0113STz9Jo0cTqrFGsMbY+efbX/Sff5buvNN+AADAS889Zx/iP/9cKlvWdTXxh2B1oEwZafBgC9hnnnFdDYAw6drVTiC8+650002uq4lPBKsj9etLnTtLb79tW+EB4HQNHiy1ayc99ZR9hRtsXnIoEpFat7ZgnTDB7kQEgMyYOFGqV89OH/Tpw92qLhGsjqWk2A/DwoU2U7hECdcVAQiahQvtg3m1atKoUbZBEu4QrD6wa5cN7E9LswvSGTcG4GStXi1VrSoVKyZ9952UO7frikCw+sQvv9gPR+HC0uTJUr58risC4HebN9v7Rvbs0owZUsGCriuCxOYl3yhe3K5x+uUXqUEDGyYBAP9k927pxhul5GR77yBU/YNg9ZGyZaWxY21KStOm9vwVAP7q0CGpYUNp3Tpp3DjpootcV4SjEaw+U7WqXes0bpzUogUDJAD8WVqaDZeZM8fOwl9+ueuK8FcEqw/dcIP02WfSgAF2Fo2n4AAkey9o00YaMUL64gvpmmtcV4Rj4do4n2ra1HYLP/SQbWTq1ElKSHBdFQBXIhGpY0epRw/pk09sLwb8iWD1sVatpD17pA4d7LD3Sy8RrkA8ygjVV1+V3nrLHhPBvwhWn3vqKfuhevpp+3fCFYgvkYj0wgvSa6/ZCNQnn3RdEU6EYA2ADh3sK+EKxJdIxG6qeeMNmy3+xBOuK8LJIFgDgnAF4kskYrdfvfWW3VTz+OOuK8LJIlgDhHAF4kPG45+335a6dJHatnVdEU4FwRowhCsQbpGI7a3o3NnuVn30UdcV4VQRrAFEuALhFIlI7dtL770ndesmPfKI64qQGQRrQB0drvv2Se+8w/2LQJBFIvYc9f33pQ8/lB5+2HVFyCyCNcA6dLAroh59VNqyRerbVzrjDNdVAThVqanWnfbsKXXvLrVu7boinA6ujQuBoUOl5s2l6tWlYcO4cg4IkoMHpTvukMaMsYlK997ruiKcLoI1JKZMkW65xa6f+/pr6fzzXVcE4ER27rTRhIsXS19+KdWr57oieIFgDZGlS22Af44cdj9jiRKuKwLwT377zX5et2616yKvusp1RfAK211CpEwZafZse85atao0f77rigAcy4oVUpUq0v790syZhGrYEKwhU7So/aAWLy7VqmWdKwD/mDnTrns780xp1izp3/92XRG8RrCGUMGC0sSJUs2a0k032d2uANwbNUqqXdtWl6ZNkwoVcl0RooFgDancuaWRI6W777Zfb7/NhemAK5GInU1t3FiqX18aN04qUMB1VYgWzrGGWGKibd+/4AI787pxo41Jy5rVdWVA/Dh40M6l9utnM3/5GQw/dgXHie7dbZBEnTrS55/b8x0A0bVunXWpK1bY8Ie77nJdEWKBYI0jEydKTZtaqI4aJZUu7boiILwmTpSSkqS8eaXhw6Vy5VxXhFjhGWscqV3bjuDkymXb+4cPd10RED6RiN2hWreuVKGCtGABoRpvCNY4U6yYnXWtV0+69VbpxRel9HTXVQHhsHevdPvtdjnGM8/YFLSCBV1XhVhj81Icyp1b+uIL+xT9/PM2Tm3AACl/fteVAcH1449So0bShg22GtSokeuK4Aoda5xKSJCefdYGf0+fbkvDP/7ouiogmEaNkipVsmXgefMI1XhHsMa5evXsjSBLFqliRdsxDODkpKVJHTtKDRtK110nzZ0rlSzpuiq4RrBC//qXNGeOTWlq1ky67z67PB3AP9u5035mXn1Veu01rmzEERy3wf+JRKT+/e3C5UKFrHutUMF1VYD/LFli51N37bKfk+uvd10R/ISOFf8nIUG65x5p0SI7e1e5sk2JYdcwYDI+fFaubN3pggWEKv6OYMXf/OtfdiTnscekJ5+UbrxR2rLFdVWAW9u22RG1e+6RmjSxW2ouucR1VfAjloJxXOPH2xtJJGKzTm+4wXVFQOwNHWrzfiWpRw8LWOCf0LHiuOrWtedJ5ctb5/rEE9KhQ66rAmJj507b0Hf77VK1atLy5YQqToyOFSclPV16/32bKFOsmNSrl73RAGE1ZozUsqV9kOzWzQI2IcF1VQgCOlaclCxZpMcftylNZ50lVa8utWpluyKBMNm9W7r/funmm22lZvlyqXlzQhUnj2DFKSldWpoxwz7BDxoklSrFMH+Ex8SJ0hVX2DPVTz6xrrVQIddVIWgIVpyyLFnsrOuKFdKVV9ozp8aNpU2bXFcGZM6ePdLDD9t9xSVKSD/8ILVoQZeKzCFYkWlFitiM1CFDpFmzpMsusx2TnHtFUEQiNuChZEnb9d6tmzRhglS0qOvKEGQEK05LQoLtmFy50r62bi3VqCGtWuW6MuD4Vqyw+b7NmklVqtjf4UcesRUZ4HTwVwieOPNMeyY1ebK0datUtqwNJ9+zx3VlwJ/t22e728uWtSvexo2zZ6oXXeS6MoQFx23guYMHbTB5585Snjx252vr1lL27K4rQzyLRCxAn3hC2rHD/l4+9RR/L+E9OlZ4LmdO6ZVXpNWr7Tqt9u2lf//bZqympbmuDvFo7lw7d92kiR2hWblSeuEFQhXRQbAiai680AZJLF9ut+Tcc49UrpwdYWCdBLGwZo2UlCRdfbW0d69tTBo1Srr4YteVIcwIVkRdyZJ2V+WcOVLBgnbwvnp120kMRMOaNVKbNvZ3b9o06dNP7dam2rVdV4Z4QLAiZq66Spo0SfrmG+seqlaVbrnFOlrAC0uXSnfeaWdRBw+2DXSrV0v33Sdlzeq6OsQLghUxlZBgN+QsWiQNHGgH8cuUsTe+9etdV4egmj5dql/fdvpOny699560bp09R82d23V1iDcEK5zIksXOD65aZcP9v/7a7oFt3176/XfX1SEI0tOl0aNt5aN6dQvS/v2ln3+WHn2UQIU7HLeBL+zbJ737rvTOO9bVtmljh/WZ04q/SkmxZd4337THCFWqSM8+K9Wrx3AH+APBCl/Zvt3eMHv2tOu6kpLsVp1y5VxXBtf275d697bz0evX29LvM89I11zjujLgzwhW+NLu3fYm2rWrLfHVrGkH++vXpyuJN7//Ln34of1d2LVLuuMOqUMHu4UG8COCFb6WmiqNGGGbUWbPtt2ebdtK997LM7Sw++03+3Pv2dMGizzwgD2D5wwq/I5gRWDMnm1vtMOGSfnz207ili3trCLCY+VK6a23pAEDbCRmmza2Gencc11XBpwcghWBs3atXe/Vt68tE1arZh1sw4bSWWe5rQ2Zs2uX9MUXdnXb7Nm2aa19e/vglDev6+qAU0OwIrAOH7Zl4l697FadrFltss7ttxOyQZCaaiMG+/WTRo603b5169roy4YNmeOL4CJYEQqbN0vDh0tffmkj7LJmtbs2M0K2YEHXFUKyGdGLF9vl4gMGSFu2SKVLW5g2b87xKoQDwYrQ2bLlzyGbJYt07bUWso0aEbKxlpoqLVhgw++HDJF+/dX+DO64wwK1QgU7uwyEBcGKUNu69UjITp1qb+BHd7Jnn+26wvBJT7eZvZMm2a9p02w29FlnSY0b2+99rVpStmyuKwWig2BF3DhWyF57rV0EUL26LUlyRvbURSLSjz8eCdLJk6WdO6UcOWzc4LXX2q8KFQhTxAeCFXFp27Y/LxenpkpnnmlBUK2aTfO58krpjDNcV+pPa9ceCdJJk+wZd2Ki3WCUEaRXX23hCsQbghVx78ABae5cuxVl+nQ77rF/v4XCVVdZ0FarJlWuHL9HPzZvtk40I0jXrLGOv3z5I0F6zTV27hSIdwQr8BepqdL33x8J2hkzbIZxlizSf/5jQ9+vuMKWjkuXlgoUcFywh9LTpU2b7A7Tn36yZ6WTJ9vQBsn+ezOCtEYN6/IB/BnBCpxAJGIhkxG08+ZZ8KSl2f9euLB0+eVHgrZkSemii6QLLvDn5dqRiLRjh/03ZQRoxteff7YOXrIPEpdeagF67bW24ei889zWDgQBwQpkwuHDtmFn+XJp2bIjX3/91YJLslAtXFgqUsSC9uivBQvaWMb8+a3jzZPn9I6cpKdLe/bY5QVH/9q1y75u327hmRGgu3cf+f8WKWIzmP/1rz9/veQSnjEDmUGwAh46cED65RcbIL9+/d+/bthgE4b+KiHBJg1lz247ZxMSjgR0xv8eidjXrFntn1NTpeRkC/nk5H+uKVs2C/JLL/17gBYvLuXK5f3vAxDPEl0XAIRJrlzSZZfZJqccOY6EZHq6hd+BA3bs568SEo6EZ3q67bDNWEb+ayebmnokdBMTj4RtSor9b0fLm9c6z2LF7OvR/3zxxYQqEA10rEAm7N9vS7+//GI7ZI/+tX79keevCQk2pi8j1I7+VbiwLQXnzWudqhfTh9LSpD/+sOMwGfX8+uuRf1637s/d7Xnn/b2ujOAtUsSCG8CpIViBE0hNtd2x06dLc+bYrNuffjrSNRYseOzgvOQSqWhRfw2Tz9j1+9cPAxm/Nmz48zPiSy+186iVK9uv0qX9uSEL8BOCFfiLgwdt5+/R51r37rWArFDBzm6WK2dHby69VMqXz3XF3jl82DrujKBdutQ+TCxZYt1wnjxSpUpHgvbqq5m9DPwVwYq4t2fPkRCdPt0GxicnW2BmTGKqVs0mMcXrJKH9++33ZfZsC9rZs216lSSVKmXHcerWtWv74vX3CMhAsCIu7dwpffWVNHSo9O23tvHn/POPhGi1ajYEgmXPY4tErKOdPdvmLn/3nT3LzZNHql/fhu3feGP8TqpCfCNYETd27LALtYcOtSBIS7OO9LbbpHr1bFmX68syb+VKm788fLi0aJEtndepYyHboAFLxogfBCtCbetWacQIC9MpU6zTql7dwrRRIy7Wjpa1a+33ffhwaeZMm+JUo4b9njdsKF14oesKgeghWBE6mzbZG/rQofbMNCHBxvHddpu9qTOWL7a2brVLzocPt5WC1FS73KBxYwvaEiVcVwh4i2BFKPz225EwnTnTno3WqSPdeqvdt8qF5v6wa5c0Zoz9WY0bZzuwL7/cPvTcf7+dnQWCjmBFYKWlSWPHSt26SRMm2FzbunXtTfrmm7l5xe8OHJDGj7eQHTnS/v2mm6TWraXrr+fSeQQXwYrA2bFD6t1b+ugjmyR01VX2ZtyoUbjOlMaTvXulgQPtz3TpUpv+1KqVdbGsNiBoCFYExvz50ocfSoMH27/fcYf0yCN2vhThEInYOdmPPpKGDLF/v/12++BUpQq7thEMBCt87dAhe4P98EObhlS0qPTww3Qy8WDHDqlvX6lHD5vJfMUVUps20t13M4QC/kawwpd+/92enXbrZm+w119vb6r16jG0Id6kp0sTJ1oXO2qU7ep+/HHpoYdY+oc/EazwlU2bpHfekXr2tDfUFi2kRx+1u0OB1ault96S+vWTcue2D1uPPSadc47ryoAjCFb4wp490ttvS50728Qe3jBxPBs2SO++K338sT2HbdlSevJJjuvAHwhWOJWSIvXqJf3vf7YztF076Zln7J5S4ER+/1364AOpa1dp3z5bHn7hBencc11XhnhGsMKJSMTOLj7zjC3v3X231KkTHQcyZ+9eex7/5pt2vrl9e/vFJQBwgSPYiLnZs+32mMaN7TLwxYtt9yehiszKm1d69lnbPfzQQ9Ibb0jFi1sne/iw6+oQbwhWxMzq1TYVqUoVW7b79lsba1e2rOvKEBYFC9qz+tWrbfrW449LJUtKAwbYZjggFghWRN327bazt1QpO4var5+0cKHN8gWioUgRm871ww/Sf/4j3XWXVK6cXQIARBvBiqg5cEB67TVbkuvfX3rlFenHH+15KmdREQulStn1dbNm2ZnX2rWl5s2lLVtcV4YwY/MSPJeWZkHasaO0bZtNSnrhBSYlwa1IxP5ePvmk7UZ//XXpwQf5kAfv0bHCM5GIPTMtV85GDlatKq1cKXXpQqjCvYQE6Z57pFWrpCZN7ANflSq2eQ7wEsEKT2zcaPee3nijVKCADVL/4gtbBgb8pGBBm+w1Y4bdB3vllXZ+es8e15UhLAhWnJZIxDaJlCplt88MHSpNnWpXuQF+VrWqbaJ7800bUnLZZfb3l4djOF0EKzJt7Vobjv/AA9Ktt0orVthXrvZCUGTLZs9cV66UKla0K+rq15d+/dV1ZQgyghWnLD3dptxcfrnt8h03Tvr0U+nMM11XBmTORRfZJLCRI6Vly6TSpW1He0qK68oQRAQrTslPP0k1ati51LvvtjehunVdVwV445ZbbOXl0UelF1+Urr5aWr7cdVUIGoIVJyUtzSbalC0rbd4sTZkide/OfZgInzx57Gq6OXNsc1OFCnbrUlqa68oQFAQrTmjNGutSn37ajigsXWr/DoTZlVfa5qaHH5aeekq69lr7WQBOhGDFcQ0aZF3qhg2227dzZylXLtdVAbGRM6fd+zppkrRunZ3RHj3adVXwO4IVx5ScbM+ZmjeXGjSQliyxG2mAeFSzpvT99/a1QQPpueek1FTHRcG3GGmIv9m40Y4dLFggvf++XcPFERrAzri+/bZdUVezpvT551yqjr8jWPEnU6fauLds2eyw/NVXu64I8J8pU6SkJJszPGSIDZsAMrAUDEn2SbxzZ+m66+wM36JFhCrwT2rWtJ+RYsXsn99/n4lNOIJghfbulZo2tQk07dvbBeQsbwHHV6iQbWp67DGbNZyUZD9LAEvBcW7VKqlxY+m336S+fW0kIYBTM3So3ehUqJA0fLjNzkb8omONY8OG2XxUyQboE6pA5tx2m/0MJSZKlSrZpibEL4I1DqWmSh062JvBjTdK8+ZJJUu6rgoItn//W5o7V2rYUGrWzI6rJSe7rgousBQcZ7Zts2dB06bZ2LbHH+coDeClSETq0UNq29bGIQ4ZIhUp4roqxBLBGkfmzLEuNSXFftgZSwhEz7x59vN28KBNMKtTx3VFiBWWguPERx9J1atLRYvaMQFCFYiuSpXsZ61CBemGG+xnEPGBYA259HQ7QvPww1KrVtLkyVLhwq6rAuLD2WdLY8dKbdocGeafnu66KkRbousCED3JydK990qDB0sffGA/3ABiK2tWGyBRrJjtaVi3TurfX8qRw3VliBaesYbUnj12fGbaNGngQHvWA8CtESPsYoty5aRRo6yjRfgQrCG0ZYtUr57066/2w8vzVMA/5s6Vbr5ZKlBA+vpr6dJLXVcEr/GMNWR++kmqUkXaulWaPp1QBfzmqqtsh36WLFLlytLixa4rgtcI1hCZN89u2cieXZo9W7riCtcVATiWYsWkmTOlSy6RatWSZs1yXRG8RLCGxDff2A9oiRLSjBnSRRe5rgjA8RQsKE2cKJUtK11/vfTdd64rglcI1hDo18+e2dSubT+oBQu6rgjAyciXzz4UV6sm1a8vjR7tuiJ4gWANuG7d7EhNixY2VD9XLtcVATgVuXJJI0dasDZubMfjEGwEa4B16WKDvtu3t9mkiZxKBgIpe3bpiy+kO+6wAf69e7uuCKeDt+KA6tzZLiZ/+mnp9dcZpA8EXWKi3YmcO7f0wAPSvn02yB/BQ7AG0JtvSs88Iz3/vNSpE6EKhEWWLFL37lLevFK7dtLevfZzzs94sBCsAfPqq9ILL0j//a/94gcOCJeEBPvwnC+f1LGjda6sSgULwRogL70k/e9/0ssv2w8cgHBKSLAP0Lly2R6KnDntgzSCgWANgEjEfqg6dZJee0169lnXFQGIhSeekA4dsuXgvHnt3+F/BGsAvPGGhepbb9m1UwDix3PP2bPW9u2lPHmkBx90XRFOhGD1uU8+sR+sl14iVIF49dpr9qz1oYcsXJs1c10Rjodg9bGRI+1y8kce4ZkqEM8SEuxO1337pLvvtiM5t9ziuir8E66N86lp02x+6C23SIMG2WXJAOJbaqoNkfjqK2nMGKlOHdcV4VgIVh9askSqXl2qWFEaO9amsgCAJCUnSw0bSlOnSt9+azdawV8IVp/59Vf7QSlcWJo82XYCAsDRDh6UbrzR7nKdPFkqX951RTgaweojW7dK11xjz1NmzJDOPdd1RQD8au9eu9Hql1/s0VGpUq4rQgaC1Sf27JFq1pS2bDlyATIAHM/OnXYP8/bt9mG8WDHXFUEiWH3h8GGpXj1p4UL75FmmjOuKAATF1q12n6skzZ7Nfcx+wLVxjqWlSXfeKc2aZZccE6oATsV550njxkm7dtmmpkOHXFcEgtWhSERq00YaPtzuYsz41AkAp6JYMTuCs2CBdN99Unq664riG8Hq0Esv2QXlvXpJDRq4rgZAkF19tTRggH1IZ6CMWwSrIx99ZMH6+uvS/fe7rgZAGNx6q/T22zYC8ZNPXFcTv9i85MCoUVKjRlLbttK773LPIgDvZDxi+vhj6euvbYIbYotgjbFFi+xZar16tmSThTUDAB5LTbVxqNOn2/G9K65wXVF8IVhjaOtW6corpfPPt3FkuXK5rghAWO3bZ6NRt2+X5s+39x3EBv1SjKSkSE2a2NeRIwlVANGVJ48N6k9NlZo2tfcexAbBGiPt29tZ1aFDbQ4wAERboULSl1/ae88zz7iuJn4QrDHQr5/0wQdS1642CxgAYuWaa6TOnW2j5BdfuK4mPvCMNcoWLrTbapo3t+3v7AAGEGuRiE14GzVKmjtXKl3adUXhRrBG0fbtUoUKtmlg2jQpRw7XFQGIV/v3S5Ur22zyefOk/PldVxReLAVHSWqqbVY6fNhGFhKqAFzKndvei7Zule69l7GH0USwRslTT9k1Tl9+KV14oetqAEC69FLps8/sZMKbb7quJrwI1igYMEDq0kV67z07RwYAfnHzzdILL9ivCRNcVxNOPGP12OLFUpUqdm6sTx82KwHwn7Q0qX59uw1n4UKpaFHXFYULweqhHTtsstLZZ9sosZw5XVcEAMe2c6dtrixY0B5bsQ/EOywFeyQ9XbrrLtt5N3w4oQrA3846y96rli+3of3wDsHqka5dpXHjbGPARRe5rgYATqxcObvCsndvuxca3mAp2ANLlkiVKkmPPGLTTQAgSFq3lj791B5hVarkuprgI1hP04ED9lz1jDNsokn27K4rAoBTc/iwVKOGtGmTbWY65xzXFQUbS8GnqX17ae1aadAgQhVAMGXPbheEHDwotWhhIxCReQTraRg5UurRw86rlirluhoAyLwLL7RnraNHSx9/7LqaYGMpOJM2bpTKlLEBEMOHc14VQDi0bm03ci1aJJUs6bqaYCJYMyE9XapTR1q1Slq61M6BAUAYHDgglS8v5colzZlj+0dwalgKzoR33pEmT5b69ydUAYRLrly2Z2TZMqljR9fVBBPBeooWLJCef17q0EG67jrX1QCA98qXl159VXr7bWnSJNfVBA9Lwadg3z77C5c/vzRzJkskAMIr45HXjz/aI6+zznJdUXDQsZ6Ctm3tnNegQYQqgHDLksU2MR04ID34IEdwTgXBepK+/NImk3zwgVSihOtqACD6LrxQ6tlTGjbMbuvCyWEp+CSsXy+VLStdf700eDBHawDEl/vus3BdtoxZ6CeDYD2BSESqXVtavdpmAp95puuKACC2du+WSpeWLr9c+uYbmosTYSn4BD791HbFffIJoQogPuXPb0vC48dLffu6rsb/6FiPY9MmG1XYqBHPFwDgnnukUaPsDtfChV1X418E63E0bizNmiWtWMFWcwD44w9bEi5f3mYKsyR8bCwF/4Nhw6QRI6Ru3QhVAJDscViPHtLYsdKAAa6r8S861mPYudOWgCtXZsA+APxV8+a2iWn5cumCC1xX4z8E6zHcf78F6ooVUqFCrqsBAH/5/fcjzceIETQff8VS8F9MnGgbld55h1AFgGMpWFD66CPbyDR4sOtq/IeO9Sj799s5rUsukb77jk9hAHA8TZvae+XKldI557iuxj/oWI/SsaO0ZYud1yJUAeD4PvjAhvU//bTrSvyFYP3/5s6V3n9f6tRJuvRS19UAgP+de670+uv2+GzGDNfV+AdLwZKSk6UKFaTs2aU5c6TERNcVAUAwpKdLVarYo7RFi6Rs2VxX5B4dq6Q33pBWrZJ69yZUAeBUZMliG5lWrJC6dHFdjT/Efce6cqXdXNOhg/TKK66rAYBgatdO6tXL3lPj/QacuA7WSMSugluzxq5DypHDdUUAEEx79kiXXSZVqmRnW+NZXC8Fjxxp51a7dCFUAeB05Msnvfeeva+OGeO6GrfitmM9eNAmh5QqZXMvAQCnJxKRbrhB+uknG3eYK5frityI24717beljRvtExYA4PQlJEgffiht3hzfe1biMljXrbOzV088If3rX66rAYDwuPRS6bnnbCzsypWuq3EjLpeCb79dmjlT+vFHKW9e19UAQLgcPixdcYV04YXxOR427jrWSZOkoUNtKZhQBQDvZc8ude0qTZ5sN4XFm7jqWFNSpHLlpAIFpOnT4+9TFADE0k032VHGlSulnDldVxM7cdWxZkwH+eADQhUAou2996RNm+x5azyJm4512zbbqJSUJPXo4boaAIgPHTpI3brZnpYiRVxXExtxE6wtW0rDhtn5qrPPdl0NAMSHPXusqalTR/rsM9fVxEZcLAUvWGAD9jt1IlQBIJby5ZNeflkaMEBauNB1NbER+o41EpGqVrUrjRYu5PYaAIi11FS77OScc2yncNj3uIS+Yx0xQpo9W3r3XUIVAFxITLQNTFOnSqNHu64m+kLdsaakSJdfLl1yiTRunOtqACB+Zdwmtn69HcEJ84Xooe5YP/1UWr1aevNN15UAQHxLSLCudfVqqWdP19VEV2g71v37bWZl7drxsxMNAPzu/vttOfjnn6X8+V1XEx2h7VjffVfaudN2AgMA/KFTJ+nAAbsIJaxCGazbtklvvSW1aSNdfLHragAAGQoXlp58UurSRfrtN9fVREcol4IffdSWf3/5RSpY0HU1AICj7d0rFSsmNWwo9erluhrvha5j/eUXG1n4zDOEKgD4Ud68dmdrnz42DS9sQtexJiVJM2bYzrN4uk0BAILk0CEbdVi5svTFF66r8VaoOtb58+0P6OWXCVUA8LMcOaT//lcaMkRatMh1Nd4KVcd63XXS1q3SkiVS1qyuqwEAHE9q6pEhPt9847oa74SmY508WZo0SXr1VUIVAIIgMdGO34wbJ02b5roa74SiY41EpBo17GzU/PnhH/AMAGGRni5VrGhLwzNmhOP9OxQd66RJ0vTp0v/+F44/FACIF1mySK+9Js2aJY0d67oabwS+Y41EpGrVpORkae5cghUAgiZj1fHgQWnevOC/jwe+Y504UZo5k24VAIIqIcF2CC9YEI6byALdsUYiUpUq9nX2bIIVAIIqEpGuuUZKSwv++3mgO9bx46U5c+hWASDoEhKkF1+0R3oTJriu5vQEtmONRKSrr7YH37NmEawAEHQZ7+uJicHeIRzYjvWbb+wh98svB/c3HwBwREbXOmuWzSYIqkB2rJGIVKmSlD27HbMhWAEgHCIRO9eaO7c0darrajInkB3r2LG2e+yllwhVAAiTjK512rTgBmvgOtaMNfgzzrDfeIIVAMIlEpHKlbOrP7/7znU1py5wHeuUKfZs9fnnCVUACKOMrnXSJNvEFDSB61ivv17atk1avJhgBYCwSk+XypaVChWyo5VBEqiOdeFCO9/0zDOEKgCEWZYsUseO0rff2ryCIAlUx9qkiV2Iu2qVnXMCAIRXWpp0xRV2X2uQBvQHpmNdvVoaOlR66ilCFQDiQdastkL59dfSsmWuqzl5gelYW7aURo+W1q61e/sAAOGXnCwVLy5dd53Ut6/rak5OIDrWjRulfv2kxx8nVAEgnpxxhtSunTRwoLRhg+tqTk4ggrVLFylXLql1a9eVAABirWVLm8T0/vuuKzk5vg/WP/6QevSQHn5YypfPdTUAgFjLl88aq48/lnbvdl3Nifk+WD/8UEpJkdq2dV0JAMCVxx6TDh+2cPU7X29eOnBAKlpUuv12qXt319UAAFx64AHbIbxmjV3C4le+7lj795d27pSefNJ1JQAA1558Utq8WRo0yHUlx+fbjjU9XSpVSipdWho2zHU1AAA/uOUWm2uwbJlNZ/Ijn5ZlsyF//NGO2AAAIEkdOkgrV9qSsF/5tmO9/npbBp4/n7nAAIAjqlSx861Tpriu5Nh82bEuX27D9h9/nFAFAPxZu3Z2CfrSpa4rOTZfdqwtW9rA5bVr7VMJAAAZUlKkYsWkunWlTz5xXc3f+a5j3b5d+uwzqU0bQhUA8HfZstnQoIEDpR07XFfzd74L1p49bfn3wQddVwIA8KuWLaVIxJ8dq6+WgpOTpYsvlm66yQIWAIB/0qKFXYT+66/WxfqFrzrWIUPs8G+7dq4rAQD43WOP2Y03I0e6ruTPfNOxRiJSxYpSwYJ2hhUAgBOpWVNKS5OmT3ddyRG+6VjnzpUWLrRPIAAAnIzHHpNmzJAWLXJdyRG+6Vjvvtt+c1avlrJmdV0NACAIUlOl4sWlWrWkvn1dV2N80bHu2GHPVx96iFAFAJy8xEQ7nvn559K2ba6rMb4I1j597Bnrffe5rgQAEDQtWlhT5pfTJM6XgtPTpRIlbPbjZ5+5rAQAEFQtWkgTJ9rRG9crn8471owzSA8/7LoSAEBQtWolrV/vj1MlzjvWBg3sN2PxYgbuAwAyJxKRypeXLrpIGjXKbS1OO9Z162zY/sMPE6oAgMxLSLANsGPGSL/95rYWp8Has6eUJ4/UrJnLKgAAYdCsmZQrl9S7t9s6nAVrcrINT777bgtXAABOR968Fq6ffGLnW11xFqzDh9uZo9atXVUAAAibVq2kjRulr792V4OzzUu1atlRm6lTXbw6ACCsKlWSzj7bXbg66VhXr5amTOHOVQCA91q1ksaNk9audfP6ToK1d2+pQAGpcWMXrw4ACLOkJHve6uoS9JgHa0qKDUq+6y4pZ85YvzoAIOxy55buvNOauJSU2L9+zIN17Fhp61bpgQdi/coAgHjRqpW0ZYs0enTsXzvmm5fq15e2b5fmzYvlqwIA4k3FitIFF0hffRXb141px7phgz1Qbtkylq8KAIhH991nO4O3bInt68Y0WPv3l3LkkJo2jeWrAgDi0R132E03AwfG9nVjthQciUglS9r5Iq6HAwDEQtOm0vLl0g8/xG4mfcw61rlzpZ9+ku65J1avCACId/fdZ8G6cGHsXjNmwdqvn3ThhTZxCQCAWKhTRypUSOrTJ3avGZNgPXRIGjzYzq66vtkdABA/sma1y14GDbIsioWYBOtXX0m7drEMDACIvfvuswyK1bGbmGxeql9f+v13ac6caL8SAAB/V7WqlC+f9M030X+tqHesW7ZI48dL994b7VcCAODY7r1X+vZbu1Iu2qIerAMHSomJnF0FALjTpImUPXtsjntGdSk4EpHKlrXzq0OGROtVAAA4sTvvlObPl1atiu6Z1qh2rN9/b4dy2bQEAHDt3nttnsL8+dF9nagGa79+0nnnSXXrRvNVAAA4sVq1bCh/tEccRi1Yk5Ot+ObN7RkrAAAuZc1q84MHD5ZSU6P3OlEL1m+/lXbssIO5AAD4QfPm0rZt0sSJ0XuNqAXr4MFSqVJSmTLRegUAAE5NuXK2oTaay8FRCdYDB6SRI63ljtVtAgAAnEhCgnWtI0ZI+/dH5zWiEqxjxljBSUnR+O4AAGRes2aWUaNGRef7R+Uca+PG0m+/RX9LMwAAmVG1qlSggDR2rPff2/OOdfdu6euvbRkYAAA/at7cxu1u3+799/Y8WEeOtKM2TZp4/Z0BAPDG7bfb1+HDvf/eni8F33CDdPCgNHWql98VAABvXX+9lJIiTZ7s7ff1tGPdvt3OBrFpCQDgd0lJ1gRu3uzt9/U0WIcOta+33ebldwUAwHuNGtlkwIzs8oqnS8E1akg5c0rjxnn1HQEAiJ6bbpJ27ZJmzPDue3rWsW7YIE2fzm5gAEBwJCVJM2faEVGveBasQ4ZIZ5whNWzo1XcEACC6GjSwC9C9vDPcs6XgihWlIkWis3UZAIBoadzYVl3nzfPm+3nSsW7daheaswwMAAiapCRp6VK79cYLnnWsu3dLOXJYSw0AQFAcPiwdOiTlz+/N94vKrGAAAOJV1O5jBQAgHhGsAAB4iGAFAMBDBCsAAB4iWAEA8BDBCgCAhwhWAAA8RLACAOAhghUAAA8RrAAAeIhgBQDAQwQrAAAeIlgBAPAQwQoAgIcIVgAAPESwAgDgof8HopWQL705yHYAAAAASUVORK5CYII=", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p0 = polar_plot(rt0, (t, -pi/4, pi/4), plot_points=50)\n", "p1 = polar_plot(rt1, (t, -pi/4, pi/4), plot_points=50)\n", "(p0+p1).show(axes=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Conclusion: in rectangular coordinates, even with 5000 plot points we do not get a nice enough plot around the origin. In polar coordinates, 100 plots suffice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 8" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve $x^2 a^2 - 2x^2 a - 3 x^2 - x a^2 + 4 x a - 3 x + a^2\n", " + 2 a - 15$ for $x$ for all values of the parameter $a$.\n", "\n", "Be as complete as possible in your description of the solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 8" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "x, a = var('x,a')\n", "p = x^2*a^2 - 2*x^2*a - 3*x^2 - x*a^2 + 4*x*a - 3*x + a^2 + 2*a - 15" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[x == 1/2*(a - sqrt(-3*a^2 - 26*a - 19) - 1)/(a + 1), x == 1/2*(a + sqrt(-3*a^2 - 26*a - 19) - 1)/(a + 1)]" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sx = solve(p,x)\n", "sx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that this solution will not work for $a = -1$." ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "a^2 - 2*a - 3" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "px2 = p.coefficient(x,2)\n", "px2" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[a == 3, a == -1]" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sa2 = solve(px2,a)\n", "sa2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another problem is the value 3 for $a$." ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p3 = p.subs(a = 3)\n", "p3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For $a = 3$, the entire polynomial vanishes and x can be anything." ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-8*x - 16" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p4 = p.subs(a = -1)\n", "p4" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[x == -2]" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(p4,x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For $a = -1$, there is one solution for $x$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 9" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Find the point with real coordinates on the curve\n", "$xy - 2 x + 3 = 0$ closest to the origin." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 9" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Find the point with real coordinates on the curve $x y - 2 x + 3$\n", " closest to the origin." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We apply the technique of Lagrange multipliers." ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y) |--> x^2 + y^2" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, y, Lambda = var('x,y,Lambda')\n", "target(x,y) = x^2 + y^2\n", "target" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The target expresses the square of the distance of the point $(x, y)$ to the origin $(0,0)$." ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y) |--> x*y - 2*x + 3" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "constraint(x,y) = x*y - 2*x + 3\n", "constraint" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y) |--> (2*x, 2*y)" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dt = diff(target)\n", "dt" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x, y) |--> (y - 2, x)" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dc = diff(constraint)\n", "dc" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-Lambda*(y - 2) + 2*x, -Lambda*x + 2*y, x*y - 2*x + 3]" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sys = [dt(x,y)[0] - Lambda*dc(x,y)[0], \\\n", " dt(x,y)[1] - Lambda*dc(x,y)[1], constraint(x,y)]\n", "sys" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{x: 0.493535367624353 + 1.8118700339534315*I,\n", " y: 1.5801426285062603 + 1.541382927845156*I,\n", " Lambda: 2.0261972405615523 - 1.1922959243930134*I},\n", " {x: 0.493535367624353 - 1.8118700339534315*I,\n", " y: 1.5801426285062603 - 1.541382927845156*I,\n", " Lambda: 2.0261972405615523 + 1.1922959243930134*I},\n", " {x: -2.1655756642279553, y: 3.385312783318223, Lambda: -3.126478616924477},\n", " {x: 1.1785049215216814, y: -0.5455981941309255, Lambda: -0.9259158751696065}]" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sols = solve(sys,[x,y,Lambda],solution_dict=True)\n", "sols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two real solutions and we have to evaluate those two solutions in the target\n", "to see which one is closest to the origin." ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{x: 0.493535367624353 + 1.8118700339534315*I, y: 1.5801426285062603 + 1.541382927845156*I, Lambda: 2.0261972405615523 - 1.1922959243930134*I} has value \n", " -2.9183064646721366 + 6.659653628869352*I\n", "{x: 0.493535367624353 - 1.8118700339534315*I, y: 1.5801426285062603 - 1.541382927845156*I, Lambda: 2.0261972405615523 + 1.1922959243930134*I} has value \n", " -2.9183064646721366 - 6.659653628869352*I\n", "{x: -2.1655756642279553, y: 3.385312783318223, Lambda: -3.126478616924477} has value \n", " 16.150060598394123\n", "{x: 1.1785049215216814, y: -0.5455981941309255, Lambda: -0.9259158751696065} has value \n", " 1.6865512394897515\n" ] } ], "source": [ "for sol in sols: \n", " print(sol, 'has value \\n', target(sol[x],sol[y]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the last value is the one closest to the origin." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How many real solutions does the system defined by $x^2 - 2 y^2 - 1 = 0$\n", "and $x y - 2 x - 3 = 0$ have?\n", "\n", "Let us first try to see if we can work over ``QQbar``." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 10" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [], "source": [ "R. = PolynomialRing(QQbar, order='lex')\n", "p = x^2 - 2*y^2 - 1\n", "q = x*y - 2*x - 3" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Ideal (x^2 + (-2)*y^2 - 1, x*y + (-2)*x - 3) of Multivariate Polynomial Ring in x, y over Algebraic Field" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "J = Ideal([p,q])\n", "J" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{y: -0.4808124295021981?, x: -1.209281267831274?},\n", " {y: 2.747679206408034?, x: 4.012415985744022?},\n", " {y: 0.8665666115470820? - 1.068362531678099?*I, x: -1.401567358956374? + 1.321102825438973?*I},\n", " {y: 0.8665666115470820? + 1.068362531678099?*I, x: -1.401567358956374? - 1.321102825438973?*I}]" ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" } ], "source": [ "J.variety()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see there are two real and two complex conjugated solutions.\n", "Let us verify with a lexicographic Groebner basis." ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Ideal (x^2 - 2*y^2 - 1, x*y - 2*x - 3) of Multivariate Polynomial Ring in x, y over Rational Field" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "RQ. = PolynomialRing(QQ, order='lex')\n", "JQ = Ideal([RQ(p),RQ(q)])\n", "JQ" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[x - 2/3*y^3 + 4/3*y^2 - 1/3*y + 2/3, y^4 - 4*y^3 + 9/2*y^2 - 2*y - 5/2]" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gb = JQ.groebner_basis()\n", "gb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The triangular structure of lexicographic Groebner basis shows there are 4 solutions." ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "y^4 - 4*y^3 + 9/2*y^2 - 2*y - 5/2 has type \n" ] } ], "source": [ "py = gb[1]\n", "print(py, 'has type', type(py))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To find the real roots of py, we have to cast the polynomial in a real ring.\n", "We cannot do this directly from Singular, but convert to a Sage expression." ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z^4 - 4*z^3 + 9/2*z^2 - 2*z - 5/2 has type \n" ] } ], "source": [ "z = var('z')\n", "pz = py.subs(y = z)\n", "print(pz, 'has type', type(pz))" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(-0.480812429502198, 1), (2.74767920640803, 1)]" ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rpz = RR['z'](pz)\n", "rpz.roots()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now see there are two real solutions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 11" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider $y'' + 6 y' + 13 y = 0$, with $y(\\pi/2) = -2$ and $y'(\\pi/2) = 8$.\n", "\n", "1. Find an exact solution to this initial value problem\n", " and use this to create a function $f$ which returns\n", " a numerical 10-digit floating-point approximation of the solution.\n", " \n", "\n", "2. Solve this initial value problem numerically.\n", " \n", " Compare the solution with the value for $y(2)$ and\n", " also with $f(2)$ obtained in the first part." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 11" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "13*y(x) + 6*diff(y(x), x) + diff(y(x), x, x)" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x,y = var('x,y')\n", "y = function('y')(x)\n", "ode = diff(y,x,2) + 6*diff(y,x) + 13*y\n", "ode" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2*cos(2*x)*e^(3/2*pi) - e^(3/2*pi)*sin(2*x))*e^(-3*x)" ] }, "execution_count": 85, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol = desolve(de=ode, dvar=y, ics=[pi/2, -2, 8])\n", "sol" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [], "source": [ "def funsol(z):\n", " return sol(x=z).n(digits=10)" ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-2.000000000" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" } ], "source": [ "funsol(z=pi/2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The system from $y\" + 6 y' + 13 y$ is\n", "\n", "$y' = z$,\n", "\n", "$z' = -6 z - 13 y$,\n", "\n", "in the ``rhs()``, ``w[0] = y`` and ``w[1] = y' = z``.\n", "\n", "So, for a numerical solver, we define the right hand side" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import odeint\n", "from numpy import linspace" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [], "source": [ "def rhs(w,t):\n", " return [w[1], -6*w[1] - 13*w[0]]" ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[8, -22]" ] }, "execution_count": 90, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rhs([-2,8], pi/2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because we are interested in the value at 2,\n", "the end of the span for x ends at 2 and begins at pi/2.\n", "Because we work numerically with scipy,\n", "we have to use the pi from scipy." ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [], "source": [ "from numpy import pi as scipypi\n", "xspan = linspace(scipypi/2,2,100)" ] }, { "cell_type": "code", "execution_count": 92, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "the numerical solution at 2: [-0.15189484 1.65169878]\n" ] } ], "source": [ "numsol = odeint(rhs,[-2,8],xspan)\n", "print('the numerical solution at 2:', numsol[-1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we compare with f(2):" ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.1518948078" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "funsol(2)" ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.15189483999619088" ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "numsolyat2 = numsol[len(numsol)-1][0]\n", "numsolyat2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us compute the difference." ] }, { "cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.216155164e-8" ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs(funsol(2)-numsolyat2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A 5-by-5 variable Toeplitz matrix has the following form:\n", "\n", " [t0 t1 t2 t3 t4]\n", " [t8 t0 t1 t2 t3]\n", " [t7 t8 t0 t1 t2]\n", " [t6 t7 t8 t0 t1]\n", " [t5 t6 t7 t8 t0]\n", "\n", "for the symbols in the list ``[t0, t1, t2, t3, t4, t5, t6, t7, t8]``.\n", "\n", "For general dimension $n$, the $(i,j)$-th element\n", "of the Toepliz matrix $T$ is\n", "\n", " $T_{(i,j)} = \\left\\{\n", " \\begin{array}{lcl}\n", " j - i & {\\rm if} & j \\geq i \\\\\n", " j - i + 2 n - 1 & {\\rm if} & j < i.\n", " \\end{array}\n", " \\right.$\n", "\n", "Define a variable Toeplitz matrix, for any dimension $n$.\n", "\n", "Test your construction for $n = 5$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 12" ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [], "source": [ "n = 5" ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[t0, t1, t2, t3, t4, t5, t6, t7, t8]" ] }, "execution_count": 98, "metadata": {}, "output_type": "execute_result" } ], "source": [ "T = [var('t' + '%d' % k) for k in range(0,2*n-1)]\n", "T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The indices in a Toeplitz matrix:" ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [], "source": [ "idx = lambda i,j: (j-i if j >= i else 2*n-1+j-i)" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0 1 2 3 4]\n", "[8 0 1 2 3]\n", "[7 8 0 1 2]\n", "[6 7 8 0 1]\n", "[5 6 7 8 0]" ] }, "execution_count": 100, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Matrix(5,5,idx)" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [], "source": [ "toeplitz = lambda i,j: (T[j-i] if j >= i else T[2*n-1+j-i])" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[t0 t1 t2 t3 t4]\n", "[t8 t0 t1 t2 t3]\n", "[t7 t8 t0 t1 t2]\n", "[t6 t7 t8 t0 t1]\n", "[t5 t6 t7 t8 t0]" ] }, "execution_count": 102, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Matrix(5,5,toeplitz)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### question 13" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maximize $x_1 + x_2$ \n", "\n", "subject to \n", "\n", "$-x_1 + 2 x_2 \\leq 8$,\n", "$4 x_1 - 3 x_2 \\leq 8$, $2 x_1 + x_2 \\leq 14$,\n", "$x_1 \\geq 0$, and $x_2 \\geq 0$.\n", "\n", "Define this problem and then solve it.\n", "\n", "What are the values of $x_1$ and $x_2$ at the optimal solution?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### answer to question 13" ] }, { "cell_type": "code", "execution_count": 103, "metadata": {}, "outputs": [], "source": [ "reset()" ] }, { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Maximization:\n", " x_0 + x_1 \n", "\n", "Constraints:\n", " - x_0 + 2.0 x_1 <= 8.0\n", " 4.0 x_0 - 3.0 x_1 <= 8.0\n", " 2.0 x_0 + x_1 <= 14.0\n", "Variables:\n", " x_0 is a continuous variable (min=0.0, max=+oo)\n", " x_1 is a continuous variable (min=0.0, max=+oo)\n" ] } ], "source": [ "p = MixedIntegerLinearProgram(maximization=True)\n", "x = p.new_variable(nonnegative=True)\n", "p.set_objective(x[1]+x[2])\n", "p.add_constraint(-x[1]+2*x[2], max=8)\n", "p.add_constraint(4*x[1]-3*x[2], max=8)\n", "p.add_constraint(2*x[1] + x[2], max=14)\n", "p.show()" ] }, { "cell_type": "code", "execution_count": 105, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10.0" ] }, "execution_count": 105, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p.solve()" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{1: 4.0, 2: 6.0}" ] }, "execution_count": 106, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p.get_values(x)" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 10.3", "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.10.14" } }, "nbformat": 4, "nbformat_minor": 4 }