{ "cells": [ { "cell_type": "markdown", "id": "733da50c", "metadata": {}, "source": [ "---\n", "format:\n", " html: default\n", "---\n", "\n", "# A3: Newton, secant & interpolation {.unnumbered}\n", "\n", "Assignment 3: Solving nonlinear equations in 1d and polynomial interpolation
\n", "Published: Wed 7 Oct\n", "\n", "
\n", "Note:\n", "\n", "
\n", "\n", "
\n", " \n", " Download (.ipynb)\n", " \n", "
" ] }, { "cell_type": "markdown", "id": "d5ef604d", "metadata": {}, "source": [ "The assignment is worth **100 pts**. The symbols ✍️, πŸ’» indicate theory ($\\LaTeX$), practice (code), respectively." ] }, { "cell_type": "code", "execution_count": null, "id": "81b62b99", "metadata": {}, "outputs": [], "source": [ "Name = \"YOUR NAME HERE\"" ] }, { "cell_type": "markdown", "id": "f30b81c4", "metadata": {}, "source": [ "The following cell contains code that you should use throughout the assignment (you do **not** need to implement Newton's method or the secant method yourself): ```Newton(f, f_prime, x1)``` and ```secant(f, x1, x2)``` are the functions from Lecture 5 (each returns the vector of iterates $x_1, x_2, \\dots$) and ```order_table(x, ΞΎ; Ξ±)``` (from A2) prints the errors $e_n := |x_n - \\xi|$ together with $\\frac{e_{n}}{e_{n-1}^\\alpha}$ (approximating $\\mu$ if the order is $\\alpha$) and $\\frac{\\log e_{n}}{\\log e_{n-1}}$ (approximating $\\alpha$). " ] }, { "cell_type": "code", "execution_count": 1, "id": "60461a73", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "order_table (generic function with 1 method)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using Plots, LaTeXStrings, Printf\n", "\n", "# Newton's method x_{n+1} = x_n - f(x_n)/f'(x_n) (from Lecture 5)\n", "function Newton(f, f_prime, x1; N=100, tol=1e-10)\n", " x = [x1]\n", " for n ∈ 2:N\n", " push!(x, x[n-1] - f(x[n-1])/f_prime(x[n-1]))\n", " if abs(f(x[end])) < tol\n", " return x\n", " end\n", " end\n", " @warn \"Newton did not converge within $N iterations.\"\n", " return x\n", "end\n", "\n", "# secant method x_{n+1} = x_n - f(x_n)/f[x_n, x_{n-1}] (from Lecture 5)\n", "function secant(f, x1, x2; N=100, tol=1e-10)\n", " x = [x1, x2]\n", " for n ∈ 3:N\n", " push!(x, x[n-1] - (x[n-1] - x[n-2])/(f(x[n-1]) - f(x[n-2])) * f(x[n-1]))\n", " if abs(f(x[end])) < tol\n", " return x\n", " end\n", " end\n", " @warn \"secant did not converge within $N iterations.\"\n", " return x\n", "end\n", "\n", "# prints the errors e_n = |x_n - ΞΎ| together with \n", "# e_n / e_{n-1}^Ξ± (\\to ΞΌ if the order of convergence is Ξ±) \n", "# log(e_n) / log(e_{n-1}) (\\to order of convergence)\n", "function order_table(x, ΞΎ; Ξ±=1)\n", " e = @. abs(x - ΞΎ)\n", " @printf(\"%4s %22s %12s %16s %20s\\n\", \"n\", \"x_n\", \"e_n\", \"e_n/e_{n-1}^Ξ±\", \"log e_n/log e_{n-1}\")\n", " @printf(\"%4d %22.16f %12.3e\\n\", 1, x[1], e[1])\n", " for n ∈ 2:length(x)\n", " @printf(\"%4d %22.16f %12.3e %16.6f %20.6f\\n\", n, x[n], e[n], e[n]/e[n-1]^Ξ±, log(e[n])/log(e[n-1]))\n", " end\n", "end" ] }, { "cell_type": "markdown", "id": "9d367d9b", "metadata": {}, "source": [ "The following theorem from lectures may be useful:\n", "\n", "
Theorem (Newton). \n", "\n", "Suppose $f: [a,b] \\to \\mathbb R$ is twice continuously differentiable with $f(\\xi) = 0$ and $f'(\\xi) \\not= 0$ and that $\\left|\\frac{f''(x)}{f'(y)}\\right| \\leq A$ for all $x,y \\in [a,b]$. Then, $x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_n)}$ converges at least quadratically to $\\xi$ for all $x_1 \\in [a,b]$ with $|x_1 - \\xi| \\leq A^{-1}$, and the asymptotic error constant is $\\mu = \\frac12 \\left| \\frac{f''(\\xi)}{f'(\\xi)} \\right|$.\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "13066b80", "metadata": {}, "source": [ "
\n", "\n", "## A: Newton's method in practice [30 pts] {.unnumbered}" ] }, { "cell_type": "markdown", "id": "f202ffa8", "metadata": {}, "source": [ "***Exercise 1. (Newton's method in practice)*** Consider the following functions:\n", "\n", "\\begin{align}\n", " f_1(x) &= x^3 - 3x - 5, \\nonumber\\\\\n", " f_2(x) &= x - \\log(1+x), \\nonumber\\\\\n", " f_3(x) &= \\arctan x, \\nonumber\\\\\n", " f_4(x) &= \\mathrm{sign}(x) \\sqrt{|x|}. \\nonumber\n", "\\end{align}\n", "\n", "Here's a plot of these functions:" ] }, { "cell_type": "code", "execution_count": null, "id": "b2208f7e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzddXgUV9sG8Hs3G3d3I4JDcAtSrEigRYoUtxQrRdpC6UcJNaDlxQlaCkFKcW+A4MFdAiTE3d2z8v2xk2wSAkR2d1ae39Wr18zZzc7TJnDnzBzhiEQiEEIIIeqKy3YBhBBCCJsoCAkhhKg1CkJCCCFqjYKQEEKIWqMgJIQQotYoCAkhhKg1XsM/IjEx8dixY0+fPtXV1R06dGj//v0b/pmEEEKIfEihR7hp06bHjx937drVzc1t7Nixf/31V8M/kxBCCJEPTsMn1AuFQi6XCdS1a9eeOnXq+vXrDS6MEEIIkQcp9AgrUhBAamqqpaVlwz+TEEIIkQ8pPCOs8OLFC39//2vXrr3vDV999dXbt29NTEzEp1paWlu2bNHW1n7f+8vKyjQ1NSXnhYXc0FAA0NMTNm4sparJR1T/LhA20HdBESj+dyE+npOWxgHg4CCytFTN5TPr+l3Q0dHh8T6SdFILwoiIiEGDBm3evLlt27bve096enrTpk179uwpPjU0NDQ1Nf3AZ+bl5enp6UnOY2O5PXoAQJMmwpAQqZRNPqr6d4Gwgb4LikDxvwurV3P27OEA2LFDNG2aagZhXb8Lle9Zvo90gjAmJqZv377Lli2bOHHiB95maGjYsWPHUaNG1fJjuVxulf+GitTMzq7NfxuRiurfBcIG+i4oAsX/LmRmMgeWlhwul8NqLbIii++CFD4uPj6+T58+ixYt8vX1bfinfUj5PVXk5Mj2QoQQooQyMpgDc3NW61A2UgjCFStWxMTErFu3zs3Nzc3NreLOp/Tp6kL8QLGoCCUlsroKIYQoJwrC+pHCrdHff//9hx9+kHzixx5LNoiREdLSACAnB1ZWMrwQIYQom/R05oCCsE6kEFqWlpbymzJhYsIEYXY2BSEhhFQQiZCdDQAcDszM2K5GqSj0g98aVIyXycpitQ5CCFEs2dng8wHAyAiKPctD4ShbEFaMl6EgJISQSui+aL0pWxBWmkHBah2EEKJYaKRMvSltEFKPkBBCKqkIQgsLVutQQsoWhHRrlBBCakI9wnpTtiCkW6OEEFITCsJ6U7YgpB4hIYTUhIKw3pQtCCtmx1AQEkJIJTRqtN6UNggrFpclhBBCPcIGoCAkhBBVQKNG642CkBBCVAH1COuNgpAQQlQBBWG9KVsQGhhASwsACgpQXMx2NYQQoigoCOtN2YIQtLgMIYRUV9E10NWFnh7b1SgbJQzCit92Kn7/IYQQ9UbdwYagICSEEKVHkwgbgoKQEEKUHvUIG4KCkBBClB5NImwICkJCCFF61CNsCApCQghRehSEDaHMQVjxdJgQQtQbBWFDKGEQVtwCpyAkhBAANGq0YSgICSFE6VGPsCEoCAkhROnRqNGGoCAkhBClpxY9wtJsxJ3QfvMLMu5J94N50v04eTAzA48HPh85OSgtZdbgJoQQNaayQSjiI+0Oki8iOQgZDyASaAGI3o4RqeBqS+siShiEHA7MzJCaCpEIGRmwtWW7IEIIYVNZGfLyAIDHg7Ex29VIRd5bJF9C0kWkXEFZXvVXNXTB0ZDi1ZQwCAFYWiI1FQDS0igICSFqLiMDIhEAmJmBw2G7mnorzUbKZSRdQvJF5EfV8AaOBszalZr10Go6ExxphpdyBqGVFUJCADBxSAghakyJ74uK+Mi4j6SLSLqIjPsQCWp4j0Ej2PSDbT9Y94aWaUlenpaBoXSrUM4gtLRkDigICSFqT/kmEeZHIjmI+ae0pp1lefqw6AKbvrDpC7N2si5HOYPQyoo5SEtjtQ5CCGGfcvQI+QVIv8OEX+ajGt7A0YCpFxN+Vj3Ald9ASOUMwooeIQUhIUTtKe4kQpEAWU+Z8Eu9AWFpDe/RsYFtP9gPgU0faJnJvURAWYOwokeYksJqHYQQwj6F6xEWJSH5EhLOIvkySjNreANPDxZd5Xbn86OUPAjpGSEhRO0pRBB+/M4nF6Ztyu98dpfiLMCGU84gtLZmDqhHSAhRe6wFoUiIzEdIuoDki0i/C2FZDe/Rc4Rtf9j2h3UfaCtIj7U65QxC6hESQkg5eY8aLUlHylUkByHhLIoSa3hDlTufbQFFn9uonEFY0SNMTma1DkIIYZ88eoQiATLuIfE/JAYi6zFEwupvEN/5tO0Pm/6w7CrPMZ8Np5xBaGQEXV0UFaGoCPn5MDBguyBCCGGNDEeNFiUi6QISA5F8qeYJfzrWsBsA2wGw6QttRRu0WlvKGYQAbGwQFQUAKSkUhIQQdSblHqGIj7RbSPwPSYHIeg6Iqr+Bw4NlF9gOgO0AmLWR653PvDze2bPo3VtyX1AalD8Ik5Ph5sZ2NYQQwg6hEFlZAMDhwNS0AR9UkoaUa0g4g4QzKM2u4Q061rDqAXsf2A+BVkOuVF/Xr2PSJN2YGLi4IDwcGlJbd1uZg1CMHhMSQtRYdjYEAgAwNoamZh2/WCRE5gMknEPieWQ+rqHzx9WCVXfYfgrbATBpKZWC66O4GP/3f1i3DkIhAKSkoKyMgpCCkBBCgPrdFy3NRNIFJJxHUiBKatrhXN8ZdoNgOwA2vcFj+9nTkyeYMIHZaAEQmZlx9uyBjo4Ur6C0QUgDRwkhpE5BmBOChLNIDkLq9Rrm/ImX+hTf+VSQOQ98PlauxC+/oKy82kGDCtavN/DwkO51lDYIK7YhTKxpFgshhKiHj0wi5Bci/TYSziDuBArjaniDtiWse8HeB/ZDoWUiuzrr7PVrTJqEBw+YUwMDrFkDX19Rfr7UL6X8QZiUxGodhBDCppp7hPkRSDiLhLM1L3XN0YB5R9gPhu1AeQ/7rA2RCDt3YuFCFBQwLV26YO9eSLsjWEFpg9DOjjmgHiEhRI1VBKGVJR+pt5BwDglnkPumhrdqm8P2U9gNhu2nCrvaGWJiMGUKrl5lTnV04OeHb7+V4tCYd1EQEkKIEivMzhzb9cKQNmeGtQ1E0Ltz3jkwawPbgbAfDPOO4MgwTqQgIABz5yIvjzlt1Qp798LLS9aXVdogtLYGjwc+H+npKC2FljIt50MIIQ2VH4mEM0g4+0Pj6xpN3hn5oqELy26w94HjCOg5sFFfHaWkwNcXp08zpzweFi3Czz/L5+92pQ1CLhdWVkhMhEiE5GQ4ObFdECGEyJhIgPQ7SDiL+NPIfS1u06j8gE/HGrb9YT8EdgPZn/ZQe0eOYNYsyU3eZs2wdy/at5fb9ZU2CAHY2TH3RRMTKQgJISqrNBOJ5xF/GkkXUJZb7UWhiHs/ouOZx0OGzhzcaUBrVgqsv+xsfP019u9nTjkczJiBdeugpyfPKpQ8CMXoMSEhRPXkRyL+FBLOIPUmRPzqr2oawqY/7H0+GTvoxn0rAKNWsFBjgwQGYto0yV/gLi74+2/06iX/QigICSFEkeSEIPYIEs7WvM+7gSts+sHeB7b9xZu8h8Uyr0h/6wnZyc3Fd99hxw5Jy4QJ2LIFhoaslKPMQUhz6gkhqkFQjJQrTP+v6J250RwuzDvBYSjsh8K4WeVXRCJkZjLH8t6evt6uXMGUKYgtD3B7e+zahQEDWKxImYPQ3p45SEhgtQ5CCKmX0kwkX0bCGcSfevfhHzR0YOkNex84fQFdu5q+Hrm5KC0FAEND6a6+KRtFRVixAn/+yaydDeCLL7B1K+sZrsxB6FA+JjiupnWDCCFEMeVHIO4kEk4j7RZEguqvalvCYQjsh8KmH3gfGTNSsb6aEtwXvXcPkyYhNJQ5tbTE1q0YMYLVmhjKHISOjsxBfDyrdRBCSC2I17xOOIO0WzW8atCIWfDauhc4tf2b+SMLjSqIsjL89ht+/ZXZLwrAwIHYtUsyzoNtKhGEcXEQicBRsOXyCCFEJEDabcSfQPxJ5EdVf5WjAYsuzMM/o8b1+Hgl6BG+fImJE/HkCXNqZIQ//4SvL6s1VafMQWhoCCMj5OaiuBiZmYr9GxEhRJ0IS5B8GXEnkHAaxanVX9XQhW1/OHwGex9oWzbkOgodhAIB1qzB8uUoKWFa+vbF7t2SPozCUOYgBODggFevACAujoKQEMIyfiFSLiP2COJPoyyn+qtaJszMB4dh0JTOPAHFDcLISEyejJs3mVNdXSxfju++A5fLalk1U/IgdHRkgjA+Xg4LsxJCSA1KMpB4DrFHkHQJwpLqr2pbwG4gnL6A7afgSnnlzPpsTy9r726i1Lkz9uxB4/rc+5UPJQ/CioGjNF6GECJnhQmIP47YY0gLrmHwp6EHHIfB4XOYdwJHVt0ghesRxsRg6lRcucKcampi6VIsWybTTZQaTlWCkGZQEELkoyAGcccQewwZdyESVn/VrC0chsHhc5i0kEMtihWEu3djwQLklk+I9PJCQABatmS1plpR8iCsWGu7YpECQgiRhfwIrfADSD2LjIeAqMpLHA1YejP9P31neRZVcWuU5SBMScFXX+HUKeZUvpsoNZySB6Fz+c8cBSEhRBbyo5BwGrFHkHZb+938s+gMpy/gNAq6tu/5etlSiB7hkSOYPVtSipsb9uyBtzd7BdWZqgRhTAyrdRBCVMsHVr6uyD/nMdCxZqM4CZaDMCsLS5ZI1s4Wb6K0di309dmopv6UPAgdHcHlQihEQgL4fPCU/D+HEMIucf7FHK7Y9lZCQ4dv0YvXaCzsh0LLhI3iqqu84raZmdwvf/w4Zs1CavksSfY2UWo4JU8ObW3Y2CAxEXw+EhIkHURCCKk9cf5F/4O8sOovaejCpg+cvoDD50XFHEOW9gmqUVYW+HwAMDaW78O4jAx8/TX++Yc55XAwbRrWrmVrE6WGU/IgBODkxGzDFBtLQUgIqYPs54g5jNjDyHtb/SVNQ9j5wGkEbAdKVr4uzpNzgR/Gzn3R8+fh6yvZ88fWFtu2YehQOVYgfcofhM7OuHsXAGJi0L0729UQQhRezivEHkbMv8h9U/0lLRPYD4XjCNj2h4aib2sk7yDMzsbixVV20/3iC2zbxsZtWSlTiSAUi45mswxCiIKTjP98Z/MHnh6se8N1IhyGird9VwpynTtx9ix8fZFUvmmwnR127MDgwbK/sDwofxC6uDAHUe+s7E4IIXnhiD2M2MPIelb9JZ4B7IfAeRRsByh+/+9dcuoR5uTg+++rdwQVYDddKVL+IHR1ZQ4oCAkhFUoyEHcMUQFIu119/ruGDmz6wukLOA4Hz4Cl+qRAHkF44QKmT5esYWltja1bMWyYzK7HDuUPwkaNmIPISFbrIIQogNIsxB1H9EGkXq++/qeGLuwGwmkU7H3AU7KJbjWS7Yrbubn47jvs3AlR+a8RKtcRrKD8QejszEwljI9HaamyrOhDCJEmQRHiTyPmHyQGVt//gasNuwFwGgX7IdLa/EhByHB7+gsXMGOGZA1na2ts24bPP5f2ZRSF8gehtjbs7BAfD4EAcXFwc2O7IEKIvIgESL+DqH2I+QdlVec2cLiw6AKnL+AyDtqKsCK19Mnk1qj4iWC1jqC/P9uLmcqW8gchgEaNmFvYkZEUhISohcxHiApAzL8oTqn+knEzuE6E60S21v+UG+kHYWAgfH0lHUErK/j7Y8QIKX264lKVILxxA6DHhISoOvESMFH7kR9R/SWjpnAeBZcvYejJRmUskGYQZmVh/nwEBEhaRo/Gpk2wtGzwRysBVQlCMQpCQlRSfiSiDyD6nxqWANV3gvMYOH8J09ZsVMYmqQXh6dOYOVMyR9DaGv7+GD68YR+qTFQiCCtuh0a880siIUR5leUg/hSi9iH5cvUpEFqmsPeB60TY9AE4LNXHJqEQ2dkAwOE0YGmXGheLUfUngu9SiSB0d2cOwsNZrYMQIg3CMiQFIioACWchKK7ykqYhHD6H81jY9gNHJf76qq/MTAgEAGBqWt99d86exVdfMWs1A7Cxgb+/6s0RrA2V+Eny8GAOwsMhEoGjjr8eEqIKMu4jaj9iDqEkrUo7VxO2n8JlPByGQkOXpeIUS4PmTlTbRxCqPEewNlQiCE1NYWaGzEwUFCApCXZ2bBdECKmLwnjEHUPk3zWsgiYeAtpoMutb4Cqa+j8gPHIEc+YgrfxXDVtbbN2Kzz6TYm1KRyWCEIC7O+7fB4DwcApCQpRDWS7iTyL2CBL/q74KjJ49HEeg0RSYerFUnKKrTxCmpmLOHBw9KmlR+O0jCssK43Pjk/OTE/MSUwpSUvJTUnJTxrQe069RPyleRVWC0MNDEoQ9erBdDSHk/UR8JF1E1H7En4SgqMpLmoZwHAHXibDqCQ6XpfqUQ52D8OBBzJsnWZbNyQk7d6J/f1nUVifF/OKk/KTEvMSKf5LykhLyEsT/zi3JffdLDr0+lPF9hg5Pauukq0oQVoyXCXtng2lCiILIC0P0QUTuQUFMlXbxKjCuE+HypVKvgi1PdQjCpCTMnInTp5lTDge+vvjjDxgZya68ysqEZSn5KeKOXeV/i9Muoyjj4x9RlYWehaaGphQrVJUg9CyfQhsaymodhJB3lOUh9jAidiP9dvWXTFvDdSKcx6r8KjBSV9sVtwMCMH8+srKYU1dX7NyJPn2kW4xAJEgtSK2xS5ecn5ySnyKqNvul1nR4OrYGtnaGdraGtrYGtlb6VmaaZsNaDtPgaEixflUJwsaNmQMKQkIUhQipNxD5N2KPgl9Q5RVdW7iMg+tEmLRkqTal9/FdeRMTMWcOTp5kTjkczJiBNWtgWM+Vx0UiUXJBclxOXFxuXHxufEx2jPggNic2JT9FUO0pb61pcjWtDawdjRytDawdjBxsDGzE/7Y3tLc1tDXXrZ7zeXl5hgZSXjxdVYLQ0xMcDkQiRESAz6/vtBpCiDQUJSJqHyJ2Ia/q1F6OBqw/gbsvHD4HV5q3ttTQh26NikTYtg2LFyOvfCFyNzf89Rd69qzNJ+eV5sVkx0RnR8flxsXlxMXnxsfkxMTlxCXkJZQKSutRqgZHw0rfys7QTvyPraGtOOTsDe1tDGys9a05bM95U5XAMDSEnR0SElBaiuhoySNDQojcCEuQdBFR+xB3AiJ+lZeMmqLRJDSaAh0rlopTNe8NwrdvMX06s/wyAC4X8+bht9+gp1ftE4r5xYl5iZFZkRX/JOYlJuUnRWVF1eNOpqmOqa2hrZ2hXcWdzIpjJ2MnHlehs0ahi6sbT08kJABAWBgFISFylROCqH2I2F19IrymMZxHw3UCLL1Zqkxl1TChns/Hli348UcUlN+IdnfHrl05nb2is9/GxMVEZ0dX/BOTE5NZlFnXi5rrmjsYOTgZOzkZOzkYOTgaOzoZOzkaOdoZ2mlpKPFesCoUhI0b4+pVAHjzBoMGsV0NIWqgNAtR+xC5u/pEeA4X1n3gNgUOw6AhtTHupLLqPcLnz0XTpiWEPYywRIQnIiy4EV2aRNjrRtwdlnUt6/0fUwMdno6zsbOziXOVzDNydDJ20tOs3q1UDSoUhE2aMAdv3rBaByFqIC0Y4TsQe6T6WqAGjdBoMlwnQd+JpcrUQmFJaTYvCh4RHLPIX++9jXx8KTz1TdQAUbFPxVuEwCskfeAzoMPTcTVxdTZxdjFxcTZ2djFxER/bGqjdCF4VCsJmzZiDV69YrYMQ1VWaicgAhO+ovh0STw+OI9BoKqx7qudeELJTKiiNyIp4k/4mLCMsIjMiIisiIjMiPjdeNEcAQARseAgAeP9sQl2erqupa+WoEyefjYGNnP4bFJ4qBmFICKt1EKKKUm8gfAfijlXvApq1h/sMOI+BppxmZ6u2hLyE0PTQsIywsIwwcfhFZ0fXfmaClY65m4Wnm5mbm6lbxb+t9WmZ1o9QoSC0t4exMXJykJ2N5GTY0C87hDRYaTZiDyNsM7JfVGnXNITzWLj7wqwdS5UpvfzSfHHghWaEVoRfXmnex78S4IJrm63RJLPMLQtumXAr1Hb74iu3BT8b6hjLumyVpEJBCKBJE9y7BwCvXlEQEtIgmY8QvgPR+8EvrNJu1g7uvrQWWl1lF2e/SnsVkhYSkhryKu3V6/TX8bnxtflCLofrbOzsae7pae7Z2KKxm5aN297TTpv2a/PLmHd8+ikObIezswyrV3WqFYTNmkmCsHdvtqshRAmVZiLib0TsRG7VRZo0jeAyDu4zYNqGpcqUSUXsvUp7JU6+hLyE2nyhqY6pp7lnE4sm4tjzNPf0MPOQrC595gxmzWLmiQGZMDvZc93UwIky+q9QH6oVhC1aMAcvX7JaByFKKOsJ3m5D9IHqy6FRF/Bjckty32a+FQeeOPxqMyddk6vpaOzYyLRRM8tmzS2bNzJtJP6n5nenpuKbb3DoUEXDUYyci82TO9PzPymgICREvQlLEXcMoZurr4itaQyXcXD3hWlrlipTUCKRKCo76mny02cpz54lP3uW8iw6O/qjX6Wtod3Usqk485pZNmth1cLFxKW2663s3YtFiySri9rY7OmwZcqZ4ajf9vTkHSoahCEhEInA9vp1hCi0wgSEb0f4ThQnV2k3aweP2XAeA55qzp6uq2J+8cvUl0+Tnz6Ie/A68/WzlGc1bpJXmTj2mlo0bWHVoqll0xZWLRqZNqrPhgmRkZg5E5cuMaccDiZNwtq1V+ebihvqvD09qYlqBaGdHczNkZGB7GzEx8PRke2CCFFImY8QugExhyAskzRyteDwGdx9YdOXvcoUQmpB6rOUZ0+Snoj7fKEZoXwh/wPv19LQamrRVBx4TS2atrRuWc/Yq0woxK5dWLQI+flMi6srtm0Tb6Vbn+3pyfupVhACaNEC168DwIsXFISEVCEoRuxhvFlbfUU0XVu4ToTnXOg5sFQZyxLyEh4mPnyQ8OBR0qNnyc+S8j+4IgtgoWfR2rq1l41Xa5vWra1bN7VsqindzTSePsX06Xj0iDnl8fDNN/j554qFsykIpUt1g/D5c1pxlBBGfgTCdyJiF0qq7gZu1g6N58F5rLptipRemP4w8eGDxAfi/Ptw8nE5XHcz99bWrZuaNu3o3LG1dWsHI5n9xlBUhBUr8L//gV/eB/Xywq5daFdlviYFoXSpXBC2asUcPHv2wfcRogZEQqRcQfgOxB1H5dVJNHTg9AWafguTVu//YpWSX5r/NPnpo6RHjxIfPUp69Drt9QdGdWppaLmbubezbdfOrl0723ZeNl4GWgYQbwlb311ta+XmTfj6SlZL1tXF999j6VJoVd/Y4eO78pK6ULkgbF0+wo2CkKizshzE/Is366svCmroDrfpcJ8BLTOWKpMToUj4Ku1VcGzwnfg7DxIehGaECkXC973ZUMuwnV279nbt29u197Lx8jDz4HK48qwW2dlYvBg7d0JUHs/du2PHDsleApWUlSE3FwB4PJiYyLFI1aVyQdiqFTQ0IBAgLAyFhe/uRUmIist5hdD1iD5QZUUYDhd2g+E5F7b9VHhR7IKygvsJ92/F3rodd/t23O2ckpz3vVOHp+Nl49XBrkN7u/Yd7Ds0Nm8s7+SrrOo0eZiYYPVqzJjxvnHvGRlMXJqa0tB46VC5INTVhbs7QkMhECAkBB06sF0QIfKSHIQ3a5EYiMo3/bTN0WgqPGbBwJW9ymQoKT/pYeLDW7G3gmODHyY+LBGU1Pg2Hpfnae5Zcbezg30HbQ1tOZdag6QkfP01jh2TtPj4YNs22Nt/4IvSyjc/trSUZW3qROWCEEDr1ggNBYCnTykIieoTliL+FF6vQcb9Ku2mbeAxEy7jVWw6oEAkeJn6Mjg2+E7cneDY4JicmPe909bAtptTt26O3To7dPay8ZIsVKYIhEJs24YffmDucgJwcIC/P4YM+eiXUhBKnSoGYZs2OHwYAJ48YbsUQmSpJAPh2xG2GUWVBj1yuHD4DE0WwtKbvcqkrLCs8E78neDY4Ntxt+/G333ffHYuh9vcsrm3k3dXx67dnLq5mihqJzgkBF99hVu3mFMuF7Nm4fffYVSrrawoCKVOCkGYl5d3+fLlhw8f5uTkbNq0qeEf2FAV44wfP2a1DkJkJj8SoRsQ8VeVdUHFY0GbL4VRDSMslI5AJHia/DQoMigoMig4NriYX1zj2/Q19b1svLydvMWdPzNdxR4EVFyM337DH3+gtJRpadYMO3agW7fafwbNnZA6KQTh06dPN2zYYGNjc/r0aYUIwrZtmYNnz1BWBk31miBFVFxaMEI3Vp8OoWMNj5nw/Brayr30pFAkfJ7y/HLU5StRV27E3Mgvza/xbfaG9hXdvtbWrWu7YifrbtzAV19JZkdoamLhQqxYAe26Pa2kHqHUSeEHqHv37levXn306NHp06cb/mlSYG4OJyfExqK4GG/eoGVLtgsipMFEfMT8g9drkfmwSrupF5oshPNocKtPNVMiYRlh4vC7Fn0tvTC9xvc0t2zey6VXV8eu3k7eTsZOcq6wodLSsGgR9u+XzI7w9sb27WjWrB4fRj1CqVOS36Tqqm1bxMYCwKNHFIREuQlLEXNI//kvKAiv1MqBTR80ngd7HyWdDpGcn3wz9mZQZNCF8AvvG/Bia2Dr7eTdt1HfgR4DHY2Uc8VEkQh79+LbbyVz4E1NsXo1pk+v99QH6hFKnVyDMCws7P79+3v37hWfampqHjhwQPv9twXy82u+MfJRWi1bap88CaDszp3iESPq9yFErN7fBdJAHEGBZvRurchNnOJkyRw3Dd0yhzGljeYIDTwBIE+ZvjslgpLbCbcDIwIvRV8Kzwqv8T1WelY9HHv0dOrZ06mni7FLRXteXp6cqny/uv5Z4EZG6ixYoHH1akULf8CA4vXrRXZ2aMAfq5QUPUADgL5+UV7eh1YDV0l1/S7o6OhofuwBmVyD0MHBwdvbu3///uJTQ0NDi4/17eu5oJE3M15O8+lTTZkuiaQeZLusFHlXaSZCNyF0I0ozJY3aFmj8NTxmaWpbKtdz78S8xPNvz59/e/5S5KUaH/sZaxv3dOnZx7VPb9fezS2bcxR4lnht/yyUlWHtWvj5obh8jI+rK/z9eQMGNHx348zyHwonJ131/KMp9QYR8jQAACAASURBVL+R5BqEenp6zZs379tX9pu8tG8PDgciEZ49Q0lJXZ9FE8Ka4lS89ceb9SirtCqKjnWp81StVkugWavh9YpAKBI+SX4SFBl0JvTM7bjb7y7sqcvTbWvbVnzns6dLTynv3sCu4GB89RVevWJOeTzMno3ffoNBw0MQoFujMqCizwjNzNCoESIiUFKCly+rLdxOiCIqiMabdQjfCUGRpNHAFY2/gftXJYVlWppK8Mt/TknOxYiL58LO/Rf+X2pB6rtv8DT3HOwxeLDnYG8nb4VY20W6srOxfDk2b4awfF3Ttm2xY4cU/woSiSRPG2l7emmRQhAmJSU1b95cIBAUFRWZmZk5Ozs/UYSZ7B06ICICAO7fpyAkCi0vHK//RMRuiCo97zFpgabfweVLcMR/SMve88UKITIrMigy6EzYmYsRF0sFpdVe1eBodHboPKTxkCGeQ5pZ1mecpHI4cAALFyK1PP6NjfH775g5E1xprmKak4OyMgAwNISOIi2Vo9SkEITW1tYR4sgBAHCl+l2vv06dcOgQANy9i1mz2K6GkJpkPcWbtYg+WGVSoGkbNP8BTiMVfDioQCS4GXPzxJsT58LORWRFvPsGGwObwR6DB3kM6ufWz1BLCbqz9RcZidmzceGCpMXHB/7+stgbnO6LyoIUgpDL5Zqamjb8c6Ssc2fm4O5dVusgpCZpwXi1GgnnqiyQbdkNzRbD/uOrTbKoVFB6JerK8dfHT745mVaYVu1VLofbzradj6fPYM/BbW3aKvKwF+koLcWff+K331BUfkPbyQlbtsDHR0YXpEmEsqCizwgBtGkDbW2UlODtW2Rk0N10ohhESDiHkN+RfqdSIwcOQ9F8Kcw7slbXxxTxiwLDA4+/Pn427Gx2cXa1V420jfq79R/sMXigx0BrfWtWKmTB1auYPVuyUgyPh3nzsGKFtAbF1Ih6hLKgukGorY02bXD3LkQi3L2LwYPZLoioveQgPP2hytIwHC7sBqHlCpi1ff+XsalEUBIYHng45PDp0NPvznywNbAd1nTY8KbDezj3UKlhnx+VkoJvv8WBA5KVYtq3x/btkvUdZYaCUBZUNwgBdOnC3BcNDqYgJGxKPI8Xfsh4IGnhaqPRZDT7HgaN2CvrvcqEZZciLv0b8u+pN6fe3d7W1cR1eNPhw5sO7+zQmc39bFkhFGL7dixdiuzybrGxMX79FbNmQUNDDtenW6OyoNJB6O2NdesAIDiY7VKIukq6gOfLkXFP0sLTg/tMNF0EXTv2yqqZUCS8GXtz//P9x18fzyzKrPZqE4smI5uNHN50eBubNqyUxzru8+f49tsqww58fLB1Kxwc5FYDBaEsqHQQdu/OTKt/8ADFxTTWmMhVchCeL0f6bUmLhi48ZqLZYugo3FO01+mv9z/fv//5/tic2GovuZm6jW4xelTzUa2tW7NSm0IoKMAvv+ivWQNB+fhed3ds3oxPP5VzIXRrVBZUOggtLeHpidBQlJTgwQN07852QUQ9ZNzDsx+RfFnSwtVCo8louVzReoFZxVlHQo4EPAt4d/EXRyPHYU2HfdHsi26O3VR/8OeHnTmDOXMQF8ec6uhg8WIsWcLK79YUhLKg0kEIwNsboaEAEBxMQUhkLvs5nv0fEs5IWrjacJ+BZkugZ89eWdUV84tPh57e93zfhfALZcIqU/Ut9CzGtBjzZcsvO9t3Vvf8AxARgblzERgoafnkE/j7owlrWx/TrVFZUPUg7N4df/0FADdv4ocf2K6GqK68cLxYjphDEJWvrcXVhNs0NF8KPQXaP+hh4sNdj3cdenmo2hAYHZ6Oj6fPhFYTBrgP0NJQ4q0NpebdVbNtbYv9/HR8fVkti3qEMqHqQVi+DQVu3YJAIJ9hXUS9FMbj5S+I/BsVXSsOF85fopUfDNxYrUwipyTnwPMDux7vepJcZflDDjjdnLpNaDVhVPNRJjombJWncK5dw+zZeP2aOeVyMX06/vyzjMNhfaAB9QhlQdWD0M0NdnZITERuLl68gJcX2wURFcLPx+s1ePVHlWWybfqizZ8wVZSftODY4F2Pdx15daSwrLByu6e55/hW48e3Gu9q4spWbYooKQkLFzKrM4p16oStW9GmDQCwvSdiSQlTgqYmjI3ZrUWlqHoQAujeHf/+CwBXr1IQEukQluLtNrz8BSXpkkabPmj9G8w7sVeWRHZx9uGQw1sebHme8rxyuw5PZ4jnEN92vn1c+9AjwCr4fPj746efkFN+09jUFL//Dl9f6a6a3RAV90UtLOq9vz2pgRoEYZ8+TBBevowFC9iuhig7EWKP4tlS5FXaYN2sPdqshnVv9qqSCI4N3vJgy4nXJ0oEJZXb29i0md52+rhW44y1qSvxjuBgzJ2LZ8+YUw4H48djzRpYWbFaVnX0gFBG1CAIe5f/9XTjBsrKoKlOC0ER6Uq/iyffIa3S+gx6jmjxf3CbDrYXWCksKzzw4sCW+1uepTyr3G6kbTS2xdgZ7Wa0s6XNyGqSkoLvv8e+fZLF0po2hb8/evVis6r3oAeEMqIGQejmBhcXREcjLw8PHqBrV7YLIkooNxTPlyH2qGSzCC1TNFuMxt9Ag+XxExFZETsf7dz1eFdGUUbl9na27Xzb+X7Z8ksDLRmuAa3EhELs34+FCyUb3erp4bvv8MMP0FbQHYOpRygjahCEAHr3xu7dAHD5MgUhqZvSLDxfjrdbJbvmaujAcy6aL4UWm7uPiUSiS5GX1t1ddzHiorBizgagr6k/vtX42R1mt7JuxWJ5iu7BA8yZgweVVn/18cHmzXB2Zq+mj6MeoYyoRxD26SMJwmXL2K6GKAmRAOE78HwZSsp7DOJ5Ea1/hT6bf12WCEr2P9+//u76l6kvK7d7mnvOaj9rstdkmgjxIRkZ+PlnbN4MYflvDx4e2LgRAwawWlatUI9QRtQmCMWLjt65g4IC6OuzXRBReClX8Wg+sisNubTpA68/2N0vKa0wbfeT3ZvubUrIS6ho5HK4vV17+7bzHd50uAaHZsq+nxLeC62GeoQyoh5BaG2N5s3x8iVKS3HzplL86kdYUxiPZ0sRtU/SoueI1r/CdSJ7NSEsM2zbtW37n+8v5hdXNBpqGU5tM3Vep3mNTBVxLyfFcvs25szB06eSlpEjsXYtHBVo3Z+Poh6hjKhHEALo3x8vXwLAf/9REJKa8Qvx+g+8Wg1Bedjw9NH0WzRbwuKImEdJj1beXHnizYnKDwKdjJ2+7vj1jHYzaC7Ex2VmYsWKKvdC3d2xcSMGDmS1rPqgHqGMqE0QDhqEtWsB4Px5bNjAdjVE0YgQfRBPF6Ow4pYjB40movVK6NqyVVNwbPDqW6vPhZ2rvC9EW9u233T6ZmzLseq1I3z98PnYsgV+fpJNdPX18eOPWLQIWkq5nir1CGVEbYKwe3cYGCA/H+HhePsWHh5sF0QURt5bPJiD5EuSFrP2aL8BFuwMMBaJRGffnl15c+Wd+DsVjRxwBnsOXuK9pJtjN1aqUj6XL+ObbxASImkZMQJr18LJib2aGoqCUEbUJgi1tNCnD06dAoD//qMgJAAgKMKr1QhZBWH5Iiy6tmjpx9YEeaFIeO7tuRXXVjxKelTRyOVwB3kM+r7j993daB+x2omOxqJFOH5c0tK4MTZuRP/+7NUkBUIhMjMBgMOBuTnb1agWtQlCAAMHSoJw3jy2qyFsSzyHh18jP4o55Wqi8Xy0/Ak8FqafC0SCgy8O/nz95/BMycpt2hrak70mf9ftOzdTtzy2l3tWDkVFWL0af/yBovJl0PX18e23SjQu9AOysiAQAICxMS2QJWVqFoRi166hsBB6eqxWQ9hTlIinS6qMC7X0Rgd/mLSUfy1CkfBwyOEV11e8SX9T0WioZfhV+68WdF5gZ6hYO9ortCNH8O23iI1lTjkcTJqElSthY8NqWVJD90VlR52C0MkJLVrg5UsUFyMoCEOHsl0QkTsRH6Eb8Xw5+PlMi7YFvFbDbQog78X8RSLR8TfH/a75VZ4Xb65rPq/TvLkd55rpmsm5HiX24gXmzcO1a5KWjh2xcSM6KcROINJCQ0ZlR52CEMBnnzGTKE6epCBUO9kvcHcqMh+Wn3PgOh5t/wdtFn7BDooMWhK0pPKzQEMtw9kdZi/xXkLrwtRBVhb8/ODvD375Anjm5li2DHPnqt4u3NQjlB31C8LffgOA06fB54OnZv/5aktYghe/4PUfkk3kTVqiw1ZYsjAC87/w/5ZfXf4gUbLKpaGW4fzO8xd2WUgRWAfiZWK++w6pqUyLpiZmzcLPP6vqlrUUhLKjZknQvj0cHREXh4wMBAcr5k4rRMrS7+DedOS8Yk41dNHSD00WQO7z8B4kPvj+0vfXoq9VtOhr6k9rO21p96XW+tZyLka53b+Pr7/G/fuSlt69sXEjmjdnryaZo1ujsqNmQcjh4LPPsHkzAJw6RUGo4gRFeLECr9dAJGBaLLuh0y4YNZFzIW8z3y69vPTYq2MVU+N1ebqzOsxa3G2xlb5ibf2q6BITsWQJ9u+XbB/o6Ihff8VENhfAkw/qEcoOy7uJsuDzz5mDEyckf5aI6km9ifNeeLWaSUGeHrxWoe8NOadgakHq3PNzm29pfvTVUXEKamlofd3x64hvIv7X/3+UgnVQWooNG9CkiWQTXT09LF+OsDB1SEFQj1CW1KxHCKBHD5iaIisLMTF4/BjtaNtulcPPx+NFCN8p2UTXdgA6boe+XJcUKSwr3HR/08qbK3NKcsQtHHBGNhv5e5/f3c3c5VmJKjh6FIsXIzKSOeVwMGYM/vgDDg6sliVX1COUHfULQk1NfPYZ9uwBgH//pSBUNRn3cHs88sqnpWuZou06NJokzxKEIuHeZ3t/vPxjUn5SReMnLp+s7re6g10HeVaiCh49woIFuHlT0uLlhY0b0V3t1tmhHqHsqN+tUQCjRzMHhw/T3VHVIeLjxQpc9JakoONw+LyScwreirvVaVenqaemVqRgS6uW5748d2XSFUrBuklMxJQp6NhRkoLm5ti6FQ8fqmEKgnqEsqR+PUIAffvCygqpqYiJwd276NKF7YJIgxVE4/YEpAUzp5pGaPMn3H3lWUJcbtziS4sPvTxUMSLG3tD+l96/TGo9icvGyqVKrKgIGzfit99QsbCcpiamTMFvv6lzb6iiR0hBKHVqGYQ8Hj7/HDt2AMC//1IQKr2oADyYI1ksxqILuu6DgZvcrl/EL9p4b+OvN37NL2Vq0OXpzus078cePxpqGcqtDFUgEuHoUXz/PaKjJY19+2L9etWeGvFRBQUoLAQAbW0YsLAaropT119UK+6OHjki2a6TKJ2SdNwYhjuTmBTk8NByOfrdlGcKHg453HhT4yVBSypScFTzUa/nvl7VdxWlYN3cuYMuXTBqlCQFW7fG5cu4dEnNUxDUHZQxtewRAujZEzY2SE5GYiIuX0a/fmwXROou+RLuTpFspWvUFF33w6yt3K4fnhk+9/zcCxEXKlra2rZdP2B9dyd1fILVILGxWLwY//4reWZvbY1ff8WUKaq3Ulr90ANCmVLXHqGGBsaMYY4DAlgthdSdiI+ni3Hl0/IU5MBzLgY+klsKlgnLVt9a3XJry4oUNNc1Xz9g/f0Z9ykF66agAH5+aNIEhw4xKailhXnzEBqK6dMpBSvQkFGZUtceIYBJk7B+PQAcPw5/fxjSXSwlUZyM4DFIvc6c6tqi827YDpDb9YMig+acnxOWESY+1eBozOk455dPfjHSNpJbDapA/Diw8sZJAHx8sGEDGjVirywFRT1CmVLXHiEALy+0bg0AhYU4epTtakjtpN7Af20lKWjvg0HP5ZaCyfnJE09M7L+vf0UKtrVte3va7Q0DNlAK1s3162jXDqNGSVKwQwfcvIkzZygFa0RBKFNqHISAZGUmujuqBEQI3YArfVGUBAAcDbRcjh6noC2PW0VCkTDgWUBz/+b7nu8Tz44w0TER3wvtaN9RDgWojtBQjBqFXr3w5AnTYmeH7dtx9y68vVmtTKFlZDAH5uas1qGi1PjWKIBx47B4Mfh8XL+OiAi4yW+oIambslzcm4bY8o67tiW6HYCNnIY4PUl+MuP0jIq9AzngjGs1bk3/NbRlRN0kJ8PPD3/9Jdk7UE8P332H776Dvj6rlSmBis2mrGh5WhlQ7x6htTUGDQIAkQg7d7JdDXmPrGcIbCdJQUtvDHoqnxQsEZT835X/67SzU0UKNrFocnnS5X3D9lEK1kF+Pvz84OGB7duZFORyMXEiQkPh50cpWBsVQUi3RmVBvYMQwIwZzMHu3SgpYbUUUpOofbjUVbJqmrsv+lyBrp0crnwv4V677e1+u/lbmbAMgC5P99fevz6b+ewTl0/kcHUVIRQiIACenlixAvnlKx707YuHD7F3r1otmd1AFc8IqUcoC+p9axTAoEFwdkZMDNLScPo0vviC7YJIOWEZHs5F+A7mVNMInXfDcYQcrlzEL/rp6k/r7qwTlG9k2NO5566hu2jXiLoJCsLChXjxQtLSrBn8/OhPWT3QYBmZUvseIZeLqVOZ4+3bWS2FVFKSjiv9JClo0hIDHsgnBW/H3W67ve2a22vEKainqbeq76ork65QCtbBvXvo2RP9+klS0MEBf/+NFy8oBeuHbo3KlNoHIYBp08DjAcCVK3j9mu1qCJATggudJHMkXMah/10Yesr6skX8oiVBS3r83eNN+htxSw/nHs9mPlvcbTGtml1bEREYPRpduuDGDabF2Bi//46wMEyeDC79b6yP0lJm+XEeD6ambFejiujnErC3x5AhACASYdMmtqtRe0mBuNgV+ZEAwOHCaxW67gdPT9aXvRp9tYV/i9W3Vos7gsbaxjuH7Lw26Rp1BGsrPR3z56NZM8nuZuI1YsLD8cMP0NVluz4llprK/B+1sACHw3Y1qoiCEAAwbx5zsHcvMjNZLUW9he/EtSEoywUATUN0P4Fmi2V9zWJ+8ZKgJX0D+kZmMRugD3Qf+GL2i+ltp3Pob53aKC3Fhg3w8MCGDSgtZRp9fBASgg0baE2whqORMrJGQQgA6NVLssrM33+zXY16EuHpD7jvCxEfAPSd0e8WHIbK+qpPk5+239F+9a3VQpEQgJmuWcCwgPPjzjsaOcr60qpAIMDff8PdHfPnIzubaezeHXfv4swZuFNnWjpopIysURCW+/pr5mDzZggErJaifoSluDUOr1Yxp2bt8eldmLSU7TVFwj9u/dFpV6eQtBBxi4+nT8jskAmtJsj0uipCJMLx42jZElOnIi6OaWzSBCdO4MYNdOrEanGqhkbKyBoFYbkvv2R+yqKjceQI29WoE34+rg9BzD/Mqf0Q9L0GHRuZXjM2J7ZPQJ/FQYtLBaUAdHm66wesPz3mtI2BbK+rIoKC0KkTRoyQDC6ztsbWrXjxAp9/zmplqolujcoaBWE5XV3Mncsc//knq6Wok5J0XO6DpIvMqecc9DgBnmyXGjny6ojXNq9r0dfEpx3sOjyd+fSbTt/QE8GPu38fffuiXz88eMC0GBvjl18QHo6ZM5nR10Ta6NaorFEQVjJ3LgwMAODxY1y6xHY1aqAwDpd6IOM+c9rqF7TfDI4Mt6DLLMocdWTUqCOjsoqzAPC4vOU9l9+edtvTXOZzM5TemzcYNQqdO+PyZaZFSwu+vggLw//9H/MHh8gGBaGsURBWYmYmmVy/atUH30oaLD8Sl3og9zUAcDTQcTta/J9ML3g/4X77He2PvGLue7uauF6ddNWvlx+PS/2YD4qJwZQpaNECR44wo/h5PPj6IjIS27fT3To5oGeEskZBWNXChdDUBIArVxAczHY1qisvDEE9URANAFxteP8Ld1/ZXU0oEv5+8/duu7tFZUeJW3zb+b6Y/cLbifb9+aDUVHzzDRo3xp49zAgyDgejR+PVK2zfDnt7tutTF/SMUNYoCKtydsaE8kGDv/zCaimqK/cNgj5BYTwA8PTQ87RM105Lzk/+dP+nP175kS/kAzDXNT815tR2n+36mrTpwfvl52P1anh4YONGyWL0ffviwQMcOgQPD1aLUzt0a1TWKAjfsXQp88z/4kXcucN2NSonJwSXeqAoEQB4Buh5Drb9ZXe1K1FX2m5vGxQZJD7tZN/poe/DoY1lPj1RiRUU4M8/4eKCJUuQm8s0du2K69dx6RLatWO1ODVFt0ZljYLwHW5uGDeOOV62jNVSVE5eOK70Q0kaAGga4ZNAWPeS0aX4Qr7fNb9++/ol5ScB4IAzr9O8m1Nvupi4yOiKSq+wEGvXolEjfP+9ZEP0Vq1w+jRu3UKPHqwWp75ooVE5oCCsyY8/Mp3Cy5cRFMR2NaqiIBZX+qIoCQA0jdD7Eiy7yehScblxn+z9ZMX1FeL1Yqz0rf4b/9+GARs0uZoyuqJyKyrC+vVwc8OiRZLeh5sb9u/HkyfMSryEJbTQqBxQENbEw0MyfPTHH5kfQ9IQxam4+ikKYgDxc8EzMO8oo0udf3vea5tXcCwz1mmA+4AXs1586vapjC6n3IqLsWED3NywYAGSk5lGZ2ds347XrzFuHO0XwToaKSMH9FP+HsuWMevl37+P48fZrkbJlWbj6qfIfQMAXC14H4WVTO6zCUXC5deWD/lnSGZRJgBNrubqvqvPf3neSp/+CnlHcTE2bYKbG+bPR1IS0+joiK1bERYGX19m+DRhG42UkQMKwvdwcMCcOczxkiWSNfVJXfHzcaUfsp4CAFcT3kdgN1AW18koyhh8cPDP138W3w51Mna6MeXG992+p/ViqispwZYtcHfHvHlITGQaHRywZQvevsXMmdDSYrU+UgWNlJEDCsL3W7oU5uYAEB6OzZvZrkY5iQS4NRaZDwGAo4EuATLaUOJJ8pMOOzoEhgeKT3u59Lo/435nh86yuJYSKytDQACaNcPcuUhIYBqtrLBqFcLCMHs2tLVZrY/UgG6NygEF4fuZmkpGjf76K9LTWa1GOT1eiISzzHGHrXAeI4uLBDwL6PYXM1meA87ibouDJgZZ61vL4lrKqrQU27fD3R2TJiGS2XYRdnbYsAExMVi8mDbOVVh0a1QOKAg/aPZseHoCQFYWvv+e7WqUTfgOhG5kjpv/APcZUr9CiaDE94zvpJOTivhFAIy0jY6OOrqq7yoNWS5YqmTKyrBzJzw9MXMmYmOZRhsbrF+P8HDMmwcdHVbrIx9BQSgHFIQfpKmJdeuY4z17cPMmq9UolcT/8KD8IavTSLT6VepXiMuN6/F3j52Pd4pPW1u3fvzV4+FNh0v9QsqqrAx//QVPT/j6IiaGabS2xtq1iIzEN99QL1Ap0DNCOaAg/JhBgzBiBACIRJg1C2VlbBekDHJCcGsss9e8WXt03guOlH/SLkZc9NrmdT+B2bliYuuJt6fddjN1k+5VlFVZGf7+G02aYPp0REczjVZWWLMGkZFYsIAiUInQM0I5oCCshU2bYGwMACEhWLuW7WoUXkkarg5EWQ4AGLii1znw9KR7hXV31w06MEg8R0JbQ9t/sP/ez/fqaUr5KkqppIQZDjN1quRZoIUFli/H27dYtAh69H9JydCtUTmgIKwFW1ssX84c//yz5Fds8i7xMNHCOADQMkHPs9CR5u+xJYKSyScnL7ywUCASAHA0crwx5cas9rOkeAlllZuL1avh7IxJkxAezjRaWGDVKkRFwc8PRkas1kfqiW6NygEFYe18/TW8vACgsFCykT151/NlSL4MABwuuv0D42ZS/Oyk/KRee3rtfbZXfNrVsev9Gfc72stqhRqlkZ6OZcvg7IwlS5CSwjRaWWHlSkRFYfFi2jVXedFCo/JBQVg7PB62bWOWmzp3Drt2sV2QQko8h5DyDY1b/ATbAVL87CfJTzrv6nw3/q74dHyr8ZcnXrYxsJHiJZRPfDwWLICLC379FdnZTKOTEzZuRFQUliyhCFR2tNCofFAQ1lqnTpg3jzlesEDyAIaIFSXhzhRABAB2A9FCmht3HHp5qNtf3WJzYgHwuLxVfVftG7ZPh6fG4/7DwzFjBtzcsH49CgqYxsaNsXs3wsPx9df0LFA10EgZ+aAgrItVq9CyJQDk52PCBGbPbgIAItydzOyvpOeALvukNUxUIBIsCVoy9thY8UxBM12zwPGBi7stlsqHK6UXLzBxIpo2xa5dkpX/WrXC3r0ICcGUKbRGqCqhkTLyQUFYF9ra2LOHWYnx9m388QfbBSmMN+uQdBEAOFx0CYC2uVQ+Nbs42+egz+pbq8WnLaxaPJjxoI9rH6l8uPK5dQtDhqB1a+zbBz6fafT2xrlzePoUEydCg5YRUDU0UkY+KAjrqG1byQhSPz88esRqNYohNxTPfmSOm34P60+k8qlR2VFd/+pasXzoYI/BwVODG5k2ksqHKxM+H4cP6/XuDW9vnD0r2RRswABcv46bNzFoED0+UlV0a1Q+KAjrbvFidO0KAKWlGDlSspe3ehIJcW86BMUAYNYerX6Wyqfejb/beVfn1+mvAXDA+annT2fGnjHWNpbKhyuN3FysWwcPD4werfHwIdPI5WLkSDx6hP/+o13jVR7dGpUPHtsFKCENDQQEoF075OQgOhpjxiAwUH3vSoVtRlowAHC10Hk3pLEF/JFXRyadYJYP1eXp7h2294tmXzT8Y5VJbCw2bsTOncjNlTRqa2P8eHz7LZo0Ya8yIlcUhPJBPcJ6cXPDvn3MbIqgIPzwA9sFsaQgRnJTtPlSmLRs+EduuLdhzNEx4hQ01zW/OOGieqXgkyeYOBHu7vjf/yQpaGlZumABIiKwaxeloFqhZ4TyQUFYX0OGSDZpWrMG//7LajUsebwI/HwAMGmJ5g39bYAv5M8+N3t+4HzxzrruZu53pt/xdvJueJlKQCjEmTPo1w9t22LfPsmStu7uWL8e0dElK1bA3p7VEgkL6BmhfNCt0QZYvhwvXuD4cYhEmDIFHh5o25bt6GRt/gAAIABJREFUmuQo5SrijgEAOOi4A9wGbWueV5o35uiY82/Pi0+7OXY7OeakhZ5Fg6tUeMXFOHwYK1fizZsq7d264ZtvMHw4c9ddvL4IUTPUI5QPCsIG4HCwezdevkRYGIqKMHo07t2DmRnbZcmFSIBH85lj1wmwaNBe8Al5CT4HfZ4mPxWfTvaavN1nu5ZGg5JVCSQmYtMm7NiBzExJo5YWRo/GwoXMkn5EvVGPUD7o1mjDGBvjxAkYGgJAeDh8fCTLfKi28O3Ifg4APAN4rWzIJz1Nftp5V2dxCnLAWd5z+e6hu1U8BcXT/lxdsWqVJAXNzLBkCaKiEBBAKUgAlJRIFho1MWG7GpVGQdhgzZphzx5m4MydOxgxQrLeh6ri5+PFCua4+VLo2tX7k06+Odltd7f43HgA2hraAcMC/Hr5cVR1VlxZGQ4fxiefoE0b7Nsn+Tlxd8fmzYiNxcqVsKv//0yiYtLSmFmjlpY0U1S2KAilYfhwbNzIHF+4oPqrr71Zj+JUANB3QdOF9f6YDfc2jDg8orCsEICZrtmFCRfGtxovrRoVS3w8fvoJzs4YPRrXrknau3fH8eMIDcWcOdDXZ608opBo7oTc0DNCKZkzBxkZzKIzhw/D1BRbt6rmb3FlOXhTvjtxy+XgatfjMwQiwfzA+Zvvbxafupm6nRt3rrF5Y2nVqChEIgQFwd8fZ85U+d1IUxPDh2PRInTowF5xRNHRSBm5oSCUnp9+QkYG0zXcvh1mZvj9d7ZrkoGwzSjNAgBDT7hOqMcH5Jfmjz029mzYWfFpF4cup8aestRTrT/rmZnYswfbtuHt2yrt9vbw9cX06XQLlHwUjZSRGwpCqVq/HllZ2LcPAFauhFCIlStVql8oKMKbDcxx86Xg1Hk9ncS8xEEHBj1LeSY+Hddy3F+f/aWtUZ9upYJ68ABbt+LQIRQVSRo5HPTpg1mzMHQoePSHjtQK3RqVG/ozKVUcDv7+GwUFOH4cAFavRlIS/vpLdf7ui9rH7LWk7wLXcXX96tCM0E/3fRqTEwOAA86ynsv8eqrK0JiSEpw+jQ0bcOtWlXZjY4wejXnz0Lw5S5URZUVBKDeq8he04tDQwIEDGDMGp04BQEAAioqwfz+zeZNyE+HNeuaw8Txw6vbDcy/hns9Bn/TCdABaGlo7h+yc2Hqi1EtkwZMn2LsXAQHIyqrS3r49Zs3CmDG0Ry6pH3pGKDcUhDKgo4OjRzFtGgICAODIEeTk4PhxpR8WmHIdua8BQNMYbtPq9KVBkUHD/x2eV5oHQF9T/8ioIwPdB8qiRvlJTsbBg9izBy9eVGnX1cXo0Zg9mwbCkAaiZ4RyQ0EoGzwe9uyBmRnWrweAixfRpw9OnICtLduVNUD4dubAdQI0jWr/dfuf7596amqZsAyAma7Z2S/PdnHoIosC5UF8CzQgAIGBkt1xxTw8MHMmJk9Wl9WFiIzRrVG5oSCUGQ4H69bBzAw//QQA9+6hfXscO4bODVqNjDWlWYg7wRy7+9b+6zbc27DwwkLxOtouJi6B4wOVdZrEvXvYuxeHDlW/BaqnhxEjMGkSevdWqYFRhG3p6cyBhRqsucsuCkIZW7YM5ub45hvw+UhMRK9e2LQJM2awXVbdxR2DsAQAzDvWcrslkUi05PKSP279IT5tYdUicHygvaGybaEQH499+xAQUH1RbA4H3btj8mSMHMmssUeIVFU8I6Rbo7JGQSh7s2ejRQt88QVSU1FSAl9fPHyITZuUbPhM9EHmwKVWg0VLBaVTTk05+IL5ql4uvU6OOalMW8wXF+PMmZpvgTo64ssvMX063N1ZKo6ovuJi5OQAgKYmTE3ZrkbVURDKRY8euHcPw4bh6VMA2LEDz55h3z54eLBdWe0UpyD1BgBwNOD08W1yC8oKRh4eGRgeKD79vMnn/4z4R4enI9MapUMkQnAw9u7FkSNVdocHYGSEkSMxeTK8vekWKJG1yt1B+nGTNQpCeXFxwa1bmDEDBw8CwL17aNMGa9bgq6+U4Mc88T+IBABg6Q3dj4z3ySjKGHJwyJ34O+LTOR3mbBy4kctR+FVtQ0Jw9Cj27UNERJV2Lhd9+mDiRAwfThMhiNykpDAHdF9UDigI5UhPDwcOoH17LFmC0lIUFGDWLJw5g127FH00aSKzXy7sB3/4jVHZUQP2DwjLCBOfLu62eFXfVTItraEePcKxY8zK19U0aYKJEzFhAhwc2KiMqDV6QChPFIRyt2AB+vXDhAnMbdLz59GsGfz8MHcusxe5ohHxkXyJObYb9IE3vkx9OWD/gIS8BAAaHA3/wf6+7eowvlR+hEI8eYIzZ3DwYPW1QAEYG2PoUEyciD59lKCzTlQUBaE8URCyoUUL3LmDH3/E+vUQCpGdjfnzcfAgtm5F27ZsF1edRvZjlGYDgL4zjN+7Tti16GufH/o8pyQHgLaG9v7h+0c2Gym3ImuFz8e1azh+HCdPIimp+quGhhg0CCNHYsgQaKvQ2qdEOVEQyhMFIUt0dPC//8HHB76+CA8HgPv30bEjZs7ETz8p1M++RsZt5si69/vec+LNiS+PfVnMLwZgqmN6euxpbydv+ZT3ccXFuHQJx4/j9GnJdvAVzMwwdCiGD0e/ftBRhuE8RD3QsjLyREHIqk8+wYsXWLUKq1ejuBgCAbZsQUAAFi3CokUwMGC7PgDQyLrHHFl2q/ENWx5smfffPPGUeVsD28Dxga2sW8mtvPcqKMD58zh+HOfOIS+v+qs2Nvj8cwwfjk8+UZ0l0YkKoR6hPNFfAWzT0YGfH8aPx9y5uHABAPLy4OcHf3/88ANmzGB9hVJuDrNlEiyqr4kjEolWXF+x4voK8WlTi6aB4wOdjJ3kWV4VfD4ePEBQEC5fxp07KC2t/gZnZwwbhhEj0LUruAo/kJWoMRo1Kk8UhIrB3R2BgTh3DkuX4vlzAEhNxYIF+O03zJ2LuXNhbs5OYSUZ3KJ4ANDQhVGVpdEEIsHsc7N3PNohPu1k3+nsl2ct9NhYDCokhAm/69erT/4Ta9wYw4dj+HC0by/34gipD+oRyhMFoSIZPBgDB+LAAfz0E6KjASA9HX5++PNPTJqEadNYGEqT/Zw5MGlRed+lEkHJuGPjjr0+Jj4d4jnk0MhDeppynGYXF4fLl5n8S06u4Q0cDlq3Zvp/tBcgUTYUhPJEQahguFxMmIDRo7F7N9asYSZ3FxTA3x/+/mjdGlOmYNw4+a3Cm1s+u67SeNGs4qwhB4fcimN2oJ3aZup2n+08rux/ljIzce0aE35hYTW/x8UFffow/9BfIUQ5iUS09YRc/X97dx4XVdX/AfwzC8M67DvigqAoipCIK4pLoplLJO6k6ONeWT8te6zUNEtTn8xS08oNTNFwSc2lzJVUENkFU0QRREBWYQYYmPn9MeOAiooyc+8M9/t+9co5x+Ger4B+uPeeew4FoU4SiTBrFqZPx759+OYbxMer+hMT8cEH+PhjDB+OsDAMGaL1Rw/Lb6leiFXrauaW5w6JGJKUpzpT/DTg0+X9l2txl/nycly6pDr5i49HbW0D77G1Rf/+qvCj9T+J/isuhkwGAObmMDZmuxoOoCDUYQIBxo3DuHE4cwZbtyIqChIJAFRXIyoKUVFwcsLw4QgKwqBBMH+JDQJfgiRL9cK0NYCM4ozB4YNvFd8CwOfx1w1Z957/exoesbYW164hJgaXLiEmBqmpDYefiQkCAlTh5+NDM19Ic0LXRRlGQagPAgMRGIgffkBkJLZtw0XVMp7IzcWWLdiyBQYG6NULQ4YgKAg+PppcD0V6T/XC2Dm1IDUoPEi5cIyQL9wyfEuYT5gGhqiuRmoqkpKQlISrV3HlCsrLG36nUIhu3VTh17MnPfZOmit1EDo4sFoHZ2gmCDdu3Lh3715TU9MPP/xw0KBBGjkmeZK5OaZPx/TpSE/Htm0ID69bHkUmw9mzOHsW//0vHB0RGAgfH/j6wsenqT9SVqruVCSV5g04GFIoLQRgYmCyL2TfGx7PW27tmWQy3LmDGzeQkoLERCQlIT1ddRmoQUIhvLzQrx8GDkRgoLZOfAnRJXRGyDANBOGOHTtWr169Y8eOnJyct99++/Lly56enk0/LHkmT0+sWoWvvsLFizh+HMePIz4ecrnqd+/fx5492LNH1XRygo9P3X/u7i93FVFWqvw1+MDUwkopAEsjy8PjDzdq4RiZDJmZuHkTN27g5k3Vf7dvP7m939NcXeHvjx494O+Prl1Zf5KSEIbRQ4QM00AQfv/990uXLu3bty+AU6dObd68+dtvv236YckLCATo0wd9+uDLL5GfjxMncPw4Tp7EgwePvS03F7m5OHZM1TQ2hrMznJzq/nN2hqMjnJ3h4ABLSxgYPPbhtZXKX4uqpQAczRyPTzrexaELJBJUVUEiQWEhcnORn4+CAty7h4IC5OcjN1f14oWZB4DHg5sbunSBtzd8fNCtG5ydm/7pIUR/0fpqDGtqEMrl8qSkpO7duyubPXr0iIiIaHJV5CXZ2yM0FKGhkMuRkIArV5CQoLrw+MT9NqkUGRlPbrn3BIFAdQXS1LTmqxKhAACqFWhTJji5Q+r+aXdUVb16qTweXF3h7g4PD/j4wNsbnTtDLH71AxLS7KgvjdKzE8xoahAWFxfLZDIrKytl09raOk99Vv8UmxJn/hmL3ReOKJs8Hs/d3V05895lmLWpq2ruQ/aRIkl2FQC5XM6vdx2Px0ebSQ5CEz4AKJCxI6+2Uo6nGIgFbSaqfo6SldZm7smHooFizNyMnAeryi6/VXnvZHGDNdv6i61fU635+eDyw6L4hudxNFj/E5iq3wwIRJtAtAFGwa611E6SIEhO5icl5VW7F7gENli/25mvze8lKF/fDFxU7tAJwH0zOFyHcuLNT+VwL6xtdX/uoxTkJYf8UmvQwEVLkeRBx4NzAYDPl3r4pgcug6ERRCKFoSEMDRWGhsq5PI/Vn1jQjD7/WqlfLpcLhHz9rR96/vlX1m//1qPlf7Vc/0Aev+sEHgArmdnDhyaaql/fP//K+svLy1+qfpvO5q2DHBsct+7gCkVDf8pGq66uNjIyysjIaNOmDYDw8PANGzZcunSpwTeHTZzauW2XgIAAZdPU1LRlS9W6lKZORjyBaq5jVYlMVl4DQCKRmNTbE5wn4Jk61e0PIC2oqq1q4BMhNBYY2YhUDQUqcisV8gb+jCJzA5G56ucAeY1Ccr+ywZqNbEVCI9WzejWS2sqip5avfHb9T9CJ+rNLkJ+PoqLH/p+XZ3o7npeXi5IS1NZWmTlUiW1W9caBDohpCQByoHw5LCQy0wePdu8zMpI6eNSaWcHSEtbWsLWFpSVsbWFtLbSzMGpjCwcH2NnBQESff43UL5FITMWm+ls/9Pzzr6xfbiYTP7p6odX6p09XPTz8Y4QocBB9/h+r/+HDh2KxuPH1i8yFBmYvOOVrahACcHZ23rlzp3Ky6NKlS9PS0iIjIxt855QpU/r37z958uRGHln5B25ieeRlVVdWvLM/NPLGAWMeJMrH0/ki+KfByAjGxqr/E2bR3wVdwNhXwdMT168DQGoqOnZkYEB9oo2vggYeQx4/fvymTZsUCkVZWdmuXbvGjx/f9GMStlTIKkZGjY68cQCAWP3dYWABNzc4O8PKilKQEG2jWaMM00AQLlq06N69e25ubu7u7oGBgSNGjGj6MQkriiuLB4cPPn7zuLL5rs87qt8woKf3CGFIdTVKSwFAIIC1NdvVcIMGHp+wsbG5ePHi3bt3jY2NbRlbDJpoWm55blB4UHJ+srK5sPfCz31G4eROADCkv46EMCQ/H8obVnZ2tHQgQzS2xJqrq6umDkWYd6v41uDwwRnFGQAEPMGGYRtmdp2JuwdUv230gjlXhBBNoWVlmEdrjRJcuXfljV1vFEgKABgKDCOCI0Z3HA0AFZmqd5i2Yq86QriFgpB5FIRc9+etP4Mjg8urywGYicwOjD0wyO3RarHlt1UvKAgJYQqtuM08CkJO+zX517BDYdW11QBsTWyPTjjq7+Jf99tl6aoX4nZsVEcIF9EZIfPoVix3rb+8PvRAqDIFW1m0Oh92/rEUBFCaonph0YHx6gjhKFpfjXkUhFykUCiWnlk67/g8uUIOwMvO6/zU8562j+8ZUpkHaS4AhcAEZm1ZqZMQDqIzQubRpVHOqZHXzDoy65f4X5TNHi16HJlwxMbY5sn3FcYof5VbeAt49AMTIQyhIGQeBSG3VMgqxuwb88eNP5TNke1H7h6921jY0GIxBdHKX2utugkYq48QzqMgZB4FIYcUSYuG7x7+z91/lM0wn7Atw7cI+c/4Hsg/p/y11ro3M+URQkCzRtlAQcgVt0tuD4kYcr3wurK5sPfClYNWPvPd1SUojAUAHr/WpicjBRJCANqVlw0UhJyQWpA6JGJIdlk2AAFP8MMbP8zym/W8D7j/FxQ1AGDtpzCwYqRGQghKS1FZCQBmZqi3Bx3RLgrC5u/snbMjd48srSoFYCgwDA8OD+kY8oKPyX60uJrzUC1XRwipQ/tOsIKCsJk7mH5wQtQEaY0UgKWR5aFxh/q26vuCj5FXIeeI6rXr21oukBBSh2bKsIKmxTdnW+O3huwLUaagk5nT6cmnX5yCAHL/hKwMAMTusOys5RoJIXUoCFlBQdhsrYpeNe33aTXyGgCetp4X/3PRx9GnUR+ZGa56QaeDhDCLgpAVdGm0GapV1M49Ondz3GZls5tzt6MTj9qZNG69puoi5Pyuet16gnYKJIQ0jIKQFRSEzU1VbdWk/ZN+u/absvm62+tRY6PEInFjPz4zArWVAGDjD0tv7dRICGkYPTvBCgrCZqWksmTE7hHns84rm6Heob+M/MWAb9DYj1fI8e/3qtdtp2qhQELI89CsUVbQPcLmI7ssu8/WPuoU/KTPJztG7XiJFASQfQgPbwKAyBqtJ2mhRkLI89ClUVbQGWEzkZiXOGzXsJyHOQD4PP7/gv43r/u8lz5K+lrVC4/ZEJpqtEBCyItRELKCgrA5OJFxImRvyMPqhwBEAtGOUTvGdRr30kcpuqJaaJtviHZzNV0jIeTFKAhZQUGo97YlbJt5eKZMLgMgFon3jdkX1DboVQ6U8qXqRevxMHbSXIGEkEapqUFxMQDw+bB5amM0oj0UhHpMoVB8cfaLL85+oWy6iF2OTjzaxaHLqxyr8DKyfwcAHh+e/6e5GgkhjVVQALkcAGxsIKR/mxlEn2x9VV1bPe33aRFJEcqmt4P30QlHW5i3eMXDJSwCFADQciytJkMIK2gDJrZQEOql4sri4MjgM7fPKJuD2w7eF7LP3ND8FQ93/0/k/Q0APCE6L9FMiYSQl6R+iNCucatfEE2hINQ/t4pvvfnrm2kP0pTN2X6zv3/jewHvVbeRV9QifqHqddtpMG+viRoJIS9N/RChoyOrdXAPBaGeuZB1ITgyuEBSAIAH3spBKz/u/XGTjnjjRxTHA4DQFJ0+10SNhJBXoQ5CujTKMApCfRKeFD799+lVtVUAjIRGO0btGOM1pklHrCpE8qNroZ0+g4lLk2skhLwiCkK2UBDqB+UE0WVnlymgAGBrYrt/7P6AlgFNPW78R6gqBABxO5osSgi7KAjZQkGoB6pqq6YdmrYreZey2cm+0+Hxh1tbtm7qcfP+xq3tqtd+34MvauoBCSFNQEHIFgpCXfdA8uCtyLcuZF1QNoPaBkWGRFoYWjT1uDXluDRN9ciE69twGtzUAxJCmkYbQbh9+/Y89XH1mUAgCAsLs9HOQgMUhDotMS9xxO4RWaVZyub73d//X9D/Xn2CaH3xH6HiNgAY2qDbDxo4ICGkabQRhDNmzPjggw/4fL3fX2H//v0eHh4jR47UxsEpCHVXVFrU5AOTK2QVAIR84fqh62f7zdbMoe//hRuqbXvRdT2MaLI2ISxTKFTPEfJ4Gn6OcMWKFQYGL7MLjU66fv26QqHQ0sEpCHWRQqFYenbp8rPLlVNjLAwt9obsHdxWQ1cvqwpwcbLqomiLkbQNPSG6oKgIMhkAWFrC0JDtajiGglDnlFeXTz44eX/afmXT3dr90LhDHe06aujwClyaCuk9ADBygP/mF72fEMIEminDIgpC3XK37O6oPaOu5l5VNoPaBu0evdvKyEpjA6R/h5wjAAAeeu6AEf2dI0QnUBCyiIJQh5zPOj967+j8CtXKuzO6ztjwxgYhX3Nfo4JoJDxahqbDfDi90m5NhBAtoCBkEQWhrtgSt+XdP95VbitoKDD88c0fp/hM0eQA0lxcCIFcBgA2/uiyQpMHJ4Q0DQUhiygI2VdZUzn3j7lb47cqm85i5wNjD/i7+GtyDHk1LoyBNBcADO0Q8Bs9Pk+ITuFmEB49erSoqKhdu3bdu3dnsQy9f7hE390uud1nax91CnZ36R47PVbDKQgFLv8HBRcAgCdEnz0wcdXo8QkhTaUOQnt7VutgkEQi2bRpk0wmk0gk7FZCZ4RsOpFxYmLUxEJpobI5xWfKpmGbjIRGGh4meRkyw1WvfVbCYYCGj08IaTIOnhFeuXLF09Nz6tSpbBdCZ4QsUSgUq6JXDds1TJmCQr5w5aCV20Zu03wK3tmD5C9Ur9tOQ4f5Gj4+IUQTuLY9/ZUrV7Zv3y6RSGJiYtiuhc4I2VBWVRZ2KEz9pKCL2GXfmH09W/TU/Ej553FxiurZeYf+6LZR80MQQjSBa2eEfn5+Mpls/vz5Xl5ebNdCQci49AfpwZHB6v3l+7bqGzk60tFMC4uclWfgfDDkVQBg0REB+2mCDCE6S31GqO17hPfvg8lbclZWsHrGg9Dp6ent27dnrpRnoyBkVERSxMwjMyUyCQAeeAt6Lfhq4FeafFJQTZKNU4NQ9QAAjJ0QeAwiS82PQgjRhNJSVFYCgFgMExMtDrRsGZYsefHbNOvjj7Fq1ZOdVVVVIpFIKNSJDKJ7hAyRyCTTfp8WeiBUmYJikXjfmH3fvP6NVlKwMg9/D1JtLiE0Rb/DMG2p+VEIIRrC2HXRgwe1e/zGD5qSktKxo6ZWjmwqnUjjZi+1IHXMvjHXCq4pmx1sO+wfu9/T1lMrg1UV4u/XUXYdAPiGCIiCdVetDEQI0RDGgnDJEixapDr7ZIZIhKVLG+hPSEjo2fPJiRESicREq2fEz0BBqHU7E3fOOTpHuZsSgFDv0I3DNpqJzLQymKwMZ4aiJBkAeAL0Cqd11AjRfYwF4ciR0M6Ofi+hoKAgOzs7ISHhs88+U3cWFRUdOXKksLDwww8/ZL4kujSqReXV5ZP2T5p8ULWnoLHQeN2QdTvf2qmtFKypwJlhKIwFAB4fPXeiZYhWBiKEaBSnpoyeOHFi7969bm5uDvX+tNbW1r6+vtrbcfD56IxQW1ILUsfuG5takKpsdrDtEBkS2dm+s7bGq6nA2RGq5WPAQ7eNtNEgIfqCU0E4adKkwsJCGxsbtgupQ0GoFTsTd84+Ols5LwZAqHfopjc3mRqYamu86hKceQMPLqqaXdfBfaa2xiKEaBqnghCATqUgKAg1rqSyZPrh6b9d+03ZFIvEm4dvHt9pvBaHrCrE6SAUxQEAeHhtDdq/r8XhCCGaxrUgbJBcLpfL5awMTUGoSf/c/Sf0QOit4lvKpq+jb2RIpIe1hxaHrLyPv4NQkgQA4KHbD/CYo8XhCCFawNkgVM6a6dmzp1AovH37tqur661bt9zc3Bgug4JQM6pqq5acXrL6n9VyheonmlDv0M3DNxsLjbU4ankG/g5CeQYA8ATw34K27C9fSwh5WdwMQrlcPm/evICAgNzc3OnTp49kbz4rBaEGJOUlTdo/KTk/Wdm0NbH9afhPozxHaXfU4gScHorK+wDAE6LndrSeqN0RCSHawbUVt5VSU1OtrKzmzp1rYGDAbiX0+ESTKBSK7y5/5/+TvzoFB7cdnDgrUespmH8WfwWqUpBviD6RlIKE6KnyclRUAICxMcy082iVDsrOzj548KCFhUVOTg7btdAZYRNklWZNOTjl9O3Tyqax0PjrQV+/7/8+j8fT8sD7cHEyaqUAILJEv8Ow66PdEQkhWqO+LuqohbX3dVZtbW1iYuKIESPYLgSgIHxl+67tm3VkVpG0SNns5twtPDi8vY32V1JPXYHEz1U7K5m4IPA4LDtpfVBCiNZw8wZhq1at8vPzg4ODzXTgLJiC8KUVSgvnHJ2zN3WvsinkC//b57+f9/vcgK/ly9xyGWJm4tY2VdPCC4FHYdpKu4MSQrSM0SD86y8sXcr0YqOffophw57oVigUUqlUF1IQFIQvKyotau7RuXkVqu9cd2v38LfCe7ToofWBqwpxYQzy/lY1HQeiz2+0sxIhzQCjQfjJJ4iL0/4wj1uw4OkgvHnzZtu2bZmu5Bloskxj5Vfkv3PgndF7R6tTMNQ7NH5mPBMpWJqKE93rUtAtjPYXJKTZYDQIAwO1P8ZT+vd/ui8hIaFrV13ZGIfOCBtld8ru94+9/0DyQNlsYd5i85ub3/B4g4mxc47gn4mQlQEAjw/v5fBaxMS4hBBGMBqEa9bgP/9RzVJlhpUVnnpAXiaTxcfHjxs3Tt2jUCgOHz4sk8kqKysnTmR6DjwF4QvkVeTNOTpnf9p+dU9Ix5Af3/zR2tha+4MrkPoVkhZD+ZC+0Ay9wtFCyw9mEEKYpQ5Ce3tGxvPUzk6ojbZ9+/bs7OzKykpvb291Z2xsbE5OzuzZs19//fXRo0cbGhoyWRIF4TMpFIptCdsWnFxQXFms7Gll0eqnET+97vY6E8NXF+Gfd3DvqKpp1gZ9D8FSa5tXEEJYwrVZoyNHjrx06ZK/v3/9Tj8/vw4dOtTW1ioUCpFIxHBJFIQNS3uQNuvIrHN3zimbPPBmd5u9ctBKsUjMxPBFcbgQgvJMVdMhEH32wdCWiaEJIcziWhBaWVkNHTr0iU4+ny8Wi3/55Zf169dr/VHsp1AQPqmypvIP6P61AAARU0lEQVSr81+til5VXVut7Glr1faXkb/0a9WPoQoydyJmlup5efDQ/j34roG2n80ghLCEm+urPe3y5cuDBg2SSqUymYzhRdcoCB9z7s65mUdmpj9IVzaFfOGcbnNWDFihrT3ln1BbiSvvIeNnVdPAHD22wvVtJoYmhLBBKkVZGQAYGsLCgu1qGCeTydLT0728vPLz8zdu3GhoaCiXy3/++ecXf6RGURCqFEmL/nvqvz/F/aRQLtoCdHXqunn45q5OTE3wfXgD50c/2lAJsPJBwG8w05XnbAgh2qA+HbS3B+NXBNk3d+5cT0/P6OjoWbNm7dixg60yKAghV8gjkiIWnFxQIClQ9lgYWnzR/4t3/d8V8AQMFXEnEjEzVM9IAGg7DX4/QGDE0OiEEJZw7QZhfZmZmQ8fPgwODmZ9w3quB2H03ej3/ngv/n68umes19h1Q9Y5mjG1/K2sDHHv49ajH4UExui2AW5hDI1OCGEVl4MwJibG19e3devWbBfC4ZVlcstz3znwTsDWAHUKulu7/zHxjz2j9zCXgg8u4phvXQqK3TH4IqUgIdzB2SC8cuVKeHi4VCqNiYlhuxZOnhHWyGs2xG5YfHpxWZXqUqSJgclHvT76pM8nRkKmrkYqapC2FkmfQy5T9bQMgf8WWjiNEE7hbBD6+fnJZLLZs2fbM7SOwPNwLghPZpz84PgHaQ/S1D1jvMasGbzG1dyVuSIq7uBiKPLPq5oGFui2gXbWJYSD6k+WYWK4ivxdybvUz4YxQMgXju803lns/PRvFRYW6kIKglNBmFqQ+tHJj47dPKbu8bLzWj90/YA2AxitI2sfYmagukTVtO2JXhEwe3ItPkIIF9y/r3rBzBnhOwfeOZFxgomR6jmYfvB82PknOnNycpycnBiu5Fk4EYSF0sJlZ5dtjN1YI69R9pgamC7otWBRwCKRgMG1fCrzEDsXd6NUTb4I3svQ4SPwuHunlhCOUwehzoQCQxITE319fdmuQqWZB6G0RvrtxW9XXlj5sPqhskfAE0z1nbp8wHIHU2Yvyd/ehbh5qCpUNc090SsC1rqyCwkhhBUMnxHufGvntoRtxdJiJgYDAFgYWUzuMvnp/vj4+AEDBjzRU1xczOfzAxnfK6rZBqFCofgt7beFfy7MLMlUdw5yG7Rm8JouDl0YLUV6DzGzkHP4UZsHj5nwXQuhCaNlEEJ0jzoIHRmZq25var+w90ImRnq248ePp6Sk5OTk9O3bV90pl8u3b9/+3XffBQQEnD//5HVUbWueQfjXrb8W/rXwau5VdY+nrefq11e/2e5NpkvJ2ofY2XUngqat0f0nOA5iugxCiO6pv76aJWcmjLdu3bqsrGzw4MF8ft1dIT6f/+233169enXgwIHMl9TcgvDcnXOLTi2Kvhut7nE0c1zWf9lU36nMLROjJM1F7Cxk//6ozYP7dPiugQEj+1cQQnRebq7qhaMjh9ZX8/T09GxoT0SZTMbj8SoqKuRyef2MZEDzmaZxNffq0F1D+23vp05BUwPTz/p+duO9G9Nfm85sCipwcwuOeNaloJkbBv4F/82UgoQQNfVDhMxcF9Vld+7ciYqK8vX1vXnz5r///svw6M3hjPB64fUV51bsSt4lV+7kDhjwDcJ8w5YGLnUyY3wmVnkGLs9A3t+qJk+A9u/D+0u6I0gIeQLDNwh1mYmJSUlJyfXr183Nzd3d3RkeXb+D8EbRjS/Pfflr8q/q5yIEPEFol9ClgUtbWbRiuhp5Fa59g9SvH20lCJh7osdW2PZkuhJCiD5geMqoDoqIiCgrK2vfvv3AgQOnT59+//79bdu2MXxdFPobhMqzwF+Tf61V1Cp7eOAFdwhePmB5B9sOLBSUdwaxc1D2aMEanhDt5qDL13QiSAh5Fo5fGi0rK4uMjAwNDTU3NwdgYGDg6srgCl/16F8Qpj1IW3FuxZ6UPeoIBBDUNujLAV/6OfuxUJA0F1fn487uuh6bbuj2I6xfY6EYQoj+4PgZYWxsbKdOncaMGcN2IXoVhNcKrq28sLL+WSCA3q69lw9Y3r91fxYKUshxOwJxH6K6SNUjskTnpWj3LhieoUoI0UNcvkd44cKFH3/80dzc/OzZs/369WO3GP0IwrjcuK/Pf30g/YB6OgyAoe5DF/db3KNFD3ZqKo5HzCwU1ttApGUI/L6HESd/tCOEvDwuB2GfPn02bNjw8ccft2/fnu1adD4I/878e+WFlX/e+rN+5zCPYYv7LfZ38WenpuoiJH6Gm5uhTmVzT3TbCAc2zkoJIXqLnSAsS0ONhLnhBEaw6Ag08JjkjRs3PDw8mKvk2XQ3CM/dOffhsQ+v5tWtDsMDb3j74Z/3/Zyde4EAFLW4uRlJi+tWihGawOtTdFgAPoOLdxNCmgX1HkzM3SOMX4C0tUwN9ojHLHTb9ERfRUWFsbEx8xNEG6SjQahQKN6KfKtIqrr3JuQLx3iN+aTPJ53tO7NWU94ZxM1DSVJdj/Mw+H0PszaslUQI0VslJZBKAUAshqkpU6PmnWFqpBcMmpyc3Lkze/+eP05Hg5DH4zmaORZJi4yFxmG+YfN7znezYm/HvvJMxH9Ut30SALO2eG0tWoxkrSRCiJ5j57qozyokL0FtJXMj8g3RadHT3QkJCb179366v6CgwM7OTvtlPUZHgxDAmSlnTqafHOw52M6E6U9KnZpypK5E+tq67xuhGTp9Cs8PwTdkrSpCiP5jJwgdB8KRhVWt68vKyrp161ZSUtKKFSue+K3bt2+vWbPmhx9+YLgk3Q1COxO7ER4jxCYsLc6pfDQi4RNIHy2LCx5ajobvGpi2ZKckQkgzwtkpo1euXImNjfXz87OysqrfL5fLCwsLZTIZ8yXpbhCyKfcE4j9CSXJdj40//NbDpjt7NRFCmhXOBmFwcPCoUaOeniYTHR3dq1cvVkqiIHxcSRLiP0LuyboeExd4r4DbOw1O/yWEkFejXl+Ng8vKPJ2Ct2/f/vfff6uqqnJycrKyslq2ZPTCGwXhI9J7SFqMW9uhXrZGaIaOH8NzPq0XSgjROI4vNPqE1q1bT5s2rayszMHBgeEUBAUhAFQX4fp6pP8PsoeqHh4frSfCZxWMGd/FiRDCDZy9NPocsbGxAQEB2dnZLVq0YHJcbgdhTTnS1yFtDWSldZ3OQ+G7GhZe7JVFCGn+KAifNnAgOzNauRqEtZW48SOufY3K/LpO667wWQnHQeyVRQjhCo5vPaFTuBeECjnuRiHhE5Tfqus094T3MrQcTTNiCCEMkMtRUAAAPB7s7dmuhvM4FYQKZP2GxM/w8N+6PtOW8PoUbaeCx6lPBSGETQUFqKkBABsbiGihYrZx5l//+38hYSGK6pbwhqEdOsxH+3kQGLFXFiGEi+gGoU7hQBA++AcJi5B/tq5HZI3278FzPgxYWraGEMJt9OyETmnWQViSjJTlyNpX1yM0Rbt30fETiCzZK4sQwnWMnREaGhq2b9+ex9P72Q95eXkzZszQ0sGbaRAWxiL1S2QfBhSqHr4hPGbBaxGM6MY0IYRljE0ZzcjIKC8v1+4YjDAwMNDew4XNLggLLiDlS+SeqOvhCeE2GZ2XwMSVvbIIIaQOY+ur2dvb29O01BdpRkGY9zdSlj+2AySPj5Zj4P0FxO1Yq4oQQp5Ck2V0SrMIwoILSFqMvNN1PTw+XN+G9xcw78BeWYQQ0jAKQp2i10GoQM4RpHyJwpi6Pr4BWo2D16cwb89eYYQQ8jwUhDpFP4NQIce9o0j+AkVxdZ18EVqNRafFELuzVxkhhLwYBaFOeXJTKN2hUChWr179VK8cWftwtBPOjqhLQb4h3GdgRAZ67qQU1LhvvvmG7RIIVq9erVAoXvw+ok1r166tra198ftepLoaxcUAIBTCxqbpx+OWdevWaXwXe17T/3bl5eXt2rUrLi6utLT0yJEjz3nnlClT+vfvP3ny5MYctqqqysLCorKyUtWWVyNzJ1JXojyj7k0GYnjMhud8eihCe3g8DXyTkCYyMjIqLS01NDRkuxBOs7CwyMrKsrCwaOJx7t6Fcsc9Z2fk5GigME6xt7dPSUnR7FRYDVwavXv3bkpKipOT08GDB5t+tAbUlCPjF6StheRuXaeBGB5z0PFjiKy1MighhGgHXRfVNRoIQj8/v61bt8bFxW3atKnpR6vPTqxA4me4sRHVxXW9hrbw/ADt3oVBU38uI4QQ5lEQ6hrdnSwj+Hfdv6urkbqirsvYCZ7/B49ZEJqxVxchhDQJBaGuadTtn+Li4tjY2Kf7e/fubWpqqnwdFxfXt2/fioqK5xynZcuW9+7dU696x+PxnJ2dn7UI3hfDC6uk5WZmqsxLyhbtTzCrrtH7FfP0TlZWVkvlDQ3CnqysLFdX12awYqReu3v3bosWLZr+VZBK+0okQwAYGZ03NT2midI4JDs729nZmc9v7EzPCRMmrFix4vnvadQZYXZ29nffffd0v6enpzoIG+Py5cvp6enqP4CxsbGtre1z3l9VVaWeHdAdmN74kYjm1P8qELbQV0EXaOGr0B1YoNEDNn8v+1VwcnJ64Xs0NiGwMWeEhBBCiK7RwHOEcrm8uLj44cOHCoWiuLi4tLS06cckhBBCmKGBM8Lc3FwvLy91s1WrVvHx8U08JiGEEMIMelaaEEIIp+nuEmuEEEIIA3T3OcL6MjIy4uLiSkpKxo8fLxaL2S6HK2pqarZv337t2jUvL6/JkycLhfrx3dLM3Llz58qVK4WFhcHBwc+fZU20RKFQxMTEnDp1qqioqHPnzuPHjxeJRGwXxUVRUVGJiYnl5eUdOnSYOHGiiYmJpo6sB5dGs7Ozvb29fXx8Tp8+nZmZ2bp1a7Yr4oqwsLDr16+HhoaGh4d7enpu3bqV7Yo4p6KiwsXF5bXXXjt37lxsbKyvry/bFXFRZmbmgAEDQkJCXFxcfv31V2Nj41OnTgkEArbr4pyRI0d2797dzMzswIEDZWVlly9f1tRP53oQhHK5nM/nS6VSExMTCkLGZGVltWvXLisry97ePj8/v2XLljdu3HB1dWW7Lm5R/vXk8XgmJibR0dEUhKyoqanh8XjK5CstLbWzs4uLi+vcuTPbdXFXZWWlubn51atXO3XqpJED6sE9wsavIEA0KDo62svLS7nEu729vZeXV3R0NNtFcQ6Px6PVZFgnFArV538ymay2tla94hVhxaVLl8zMzDT4cznd9SENu3//vp2dnbrp4OBw7949FushRBfMmzcvJCSkTZs2bBfCUWPHjj179mxFRUVUVFTT98NS04mTrYKCAmFDoqKi2C6Nu4RCYf09SGUyGU0QIBz32WefpaSkaHybHdJ4P//8c1xc3Ndffz1hwoQcze3lqBNBaGdnV9OQt99+m+3SuMvFxaX+91lOTo6zszOL9RDCruXLlx86dOjUqVNWVlZs18JdYrHYxcXl3XffdXNzO3HihKYOqxNBSHTQgAEDsrKyUlJSAKSkpGRlZQ0cOJDtoghhx9q1ayMiIk6ePElPsLClsrJSLpcrXxcWFmZmZrZq1UpTB9ePe4SBgYFlZWUARowYIRKJoqOjaSV+bbO0tFy8ePGQIUOCgoJOnDixZMkSDV6RJ403fPjw3Nzcqqoq5YNTx44dq3/vljDg2rVrCxYscHNzGz58uLLn22+/DQgIYLcqrrl48eK0adP8/Px4PN7p06fffPPNAQMGaOrgevD4BICEhIT696t8fX1pKikzEhMTlQ/Ue3t7s10LRyUnJ1dXV6ub3t7eBgYGLNbDQRKJJC0trX6Pu7s7/VzIvPT09GvXrgHo2LGjp6enBo+sH0FICCGEaAmdVxFCCOE0CkJCCCGcRkFICCGE0ygICSGEcBoFISGEEE6jICSEEMJpFISEEEI4jYKQEEIIp1EQEkII4TQKQkIIIZxGQUgIIYTT/h/NZtzYo4RjawAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f1(x) = x^3 - 3x - 5\n", "f2(x) = x - log(1 + x)\n", "f3(x) = atan(x)\n", "f4(x) = sign(x) * sqrt(abs(x))\n", "\n", "plot( f1, -1, 3; ylims=(-1,2), label=L\"f_1\", lw=3, color=:blue )\n", "plot!( f2, -0.99, 3; label=L\"f_2\", lw=3, color=:red )\n", "plot!( f3, -1, 3; label=L\"f_3\", lw=3, color=:green )\n", "plot!( f4, -1, 3; label=L\"f_4\", lw=3, color=:orange )\n", "hline!( [0], linestyle=:dash, primary=false )" ] }, { "cell_type": "markdown", "id": "f2f2953c", "metadata": {}, "source": [ "
\n", "\n", "(i) ✍️ Show that $f_1$ has a root $\\xi \\in [2,3]$ and that $f_2(0) = f_3(0) = f_4(0) = 0$.\n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "facb43a4", "metadata": {}, "source": [ "
\n", "\n", "(ii) $f_1$: \n", "\n", "πŸ’» Run Newton's method for $f_1$ with $x_1 = 2$ using ```Newton``` and ```order_table``` from the code cell at the start of the assignment. If the iteration converges, estimate the order of convergence $\\alpha$ and the asymptotic error constant $\\mu$ numerically.\n", "\n", "✍️ Find $A$ such that $\\left|\\frac{f_1''(x)}{f_1'(y)}\\right| \\leq A$ for all $x, y \\in [2,3]$ and compare your numerical estimate of $\\mu$ with the theoretical value $\\frac12 \\left| \\frac{f_1''(\\xi)}{f_1'(\\xi)} \\right|$." ] }, { "cell_type": "code", "execution_count": null, "id": "f87d311b", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "ccd6a4b8", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "a4e8e78e", "metadata": {}, "source": [ "
\n", "\n", "(iii) $f_2$: \n", "\n", "πŸ’» Run Newton's method for $f_2$ with $x_1 = 0.5$ using ```Newton``` and ```order_table``` from the code cell at the start of the assignment. If the iteration converges, estimate the order of convergence $\\alpha$ and the asymptotic error constant $\\mu$ numerically.\n", "\n", "✍️ Why does the convergence behaviour for $f_2$ not contradict **Theorem (Newton)** above?" ] }, { "cell_type": "code", "execution_count": null, "id": "64f05362", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "eabdf36e", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "44827b80", "metadata": {}, "source": [ "
\n", "\n", "(iv) $f_3$: \n", "\n", "πŸ’» Run Newton's method for $f_3$ with $x_1 = 0.5$ using ```Newton``` and ```order_table``` from the code cell at the start of the assignment. If the iteration converges, estimate the order of convergence $\\alpha$ and the asymptotic error constant $\\mu$ numerically.\n", "\n", "✍️ For $f_3$, the convergence (from $x_1 = 0.5$) is faster than quadratic. Explain why. What happens if you instead start at $x_1 = 2$? Explain.\n", "\n", "*Hint:* see the discussion of Kepler's equation (Example 5.5)." ] }, { "cell_type": "code", "execution_count": null, "id": "6b2c133f", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "9f9480ba", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "ce5fce15", "metadata": {}, "source": [ "
\n", "\n", "(v) $f_4$: \n", "\n", "πŸ’» Run Newton's method for $f_4$ with $x_1 = 0.5$ using ```Newton``` and ```order_table``` from the code cell at the start of the assignment. If the iteration converges, estimate the order of convergence $\\alpha$ and the asymptotic error constant $\\mu$ numerically.\n", "\n", "✍️ Compute $x_{n+1}$ explicitly in terms of $x_n$. Hence show that Newton's method does not converge for any $x_1 \\not= 0$. Why does this not contradict **Theorem (Newton)** above?" ] }, { "cell_type": "code", "execution_count": null, "id": "93ecab4d", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "2d06567f", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "a19e1665", "metadata": {}, "source": [ "
\n", "\n", "## B: Efficiency of Newton and secant [20 pts] {.unnumbered}" ] }, { "cell_type": "markdown", "id": "6f185e05", "metadata": {}, "source": [ "The *efficiency index* of a method with order of convergence $\\alpha$ requiring $\\theta$ function evaluations per iteration is $E := \\alpha^{1/\\theta}$. (Here, evaluating $f$ or $f'$ each counts as one function evaluation.)\n", "\n", "***Exercise 2. (Newton vs. secant)*** Let $f_1(x) = x^3 - 3x - 5$ and $\\xi \\in [2,3]$ be as in Exercise 1." ] }, { "cell_type": "markdown", "id": "843b59f6", "metadata": {}, "source": [ "(i) ✍️ What is the efficiency index of Newton's method and of the secant method?" ] }, { "cell_type": "markdown", "id": "254cff31", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "e61607df", "metadata": {}, "source": [ "
\n", "\n", "(ii) πŸ’» Apply the secant method (```secant``` from the code cell at the start of the assignment) to $f_1$ with $x_1 = 2, x_2 = 3$. Estimate the order of convergence and compare with the theoretical values $\\alpha = \\frac{1 + \\sqrt5}{2}$ and $\\mu = \\left| \\frac12 \\frac{f_1''(\\xi)}{f_1'(\\xi)} \\right|^{\\alpha/(1+\\alpha)}$ from Lecture 5." ] }, { "cell_type": "code", "execution_count": null, "id": "ebd500fe", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "8019d0c7", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "50437a24", "metadata": {}, "source": [ "
\n", "\n", "(iii) πŸ’» For both Newton's method (with $x_1 = 2$) and the secant method (with $x_1 = 2, x_2 = 3$), count the total number of function evaluations needed to obtain $|x_n - \\xi| < 10^{-12}$. Which method is cheaper? Is this consistent with the efficiency indices of the two methods?\n", "\n", "*Hint:* use ```Newton``` and ```secant``` from the code cell at the start of the assignment. Both stop once $|f(x_n)| <$ ```tol``` (default $10^{-10}$), so you may need to pass a smaller tolerance, e.g. ```Newton(f, f_prime, x1; tol=1e-14)```." ] }, { "cell_type": "code", "execution_count": null, "id": "d1bcbcd0", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "3fb00905", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "68582bfb", "metadata": {}, "source": [ "
\n", "\n", "## C: Barycentric formula & divided differences [30 pts] {.unnumbered}" ] }, { "cell_type": "markdown", "id": "d530c319", "metadata": {}, "source": [ "Recall from Lecture 8 that, for distinct nodes $X = \\{x_0, \\dots, x_n\\}$, the divided differences satisfy the recursion\n", "\n", "\\begin{align}\n", " f[x_0, \\dots, x_n] = \\frac{f[x_1, \\dots, x_n] - f[x_0, \\dots, x_{n-1}]}{x_n - x_0}, \\nonumber\n", "\\end{align}\n", "\n", "the polynomial interpolant can be written in *Newton form* \n", "\n", "\\begin{align}\n", " I_Xf(x) &= \\sum_{k=0}^n f[x_0, \\dots, x_k] \\, \\pi_k(x), \\nonumber\\\\\n", " \\pi_k(x) &:= (x - x_0) \\cdots (x - x_{k-1}), \\nonumber\n", "\\end{align}\n", "\n", "and in *Barycentric form* \n", "\n", "\\begin{align}\n", " I_Xf(x) &= \\left. \\sum_{j=0}^n \\frac{\\lambda_j f(x_j)}{x - x_j} \\middle/ \\sum_{j=0}^n \\frac{\\lambda_j}{x - x_j} \\right., \\nonumber\\\\\n", " %\n", " \\lambda_j &:= \\frac{1}{\\prod_{k\\not=j} (x_j - x_k)}. \\nonumber\n", "\\end{align}" ] }, { "cell_type": "markdown", "id": "90c94da9", "metadata": {}, "source": [ "
\n", "\n", "***Exercise 3. (Divided differences)*** Let $f(x) = \\frac1x$.\n", "\n", "(i) ✍️ Compute the divided difference table of $f$ on $X = \\{ x_0, x_1, x_2 \\} = \\{1, 2, 4\\}$ and write down $I_Xf$ in Newton form. Now add the node $x_3 = \\frac12$: update your table and write down the new interpolant. Which terms did you have to compute again?" ] }, { "cell_type": "markdown", "id": "cb60c1bf", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "14494d35", "metadata": {}, "source": [ "
\n", "\n", "(ii) ✍️ Show that, for any distinct nonzero nodes $x_0, \\dots, x_n$,\n", "\n", "\\begin{align}\n", " f[x_0, \\dots, x_n] = \\frac{(-1)^n}{x_0 x_1 \\cdots x_n}. \\nonumber\n", "\\end{align}\n", "\n", "Check that this agrees with your table in (i)." ] }, { "cell_type": "markdown", "id": "fe628b33", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "14da7fa8", "metadata": {}, "source": [ "
\n", "\n", "(iii) ✍️ Let $X = \\{x_0, \\dots, x_n\\}$ be distinct nodes and $x \\not\\in X$. By adding $x$ as an extra interpolation node, use the Newton form to show that \n", "\n", "\\begin{align}\n", " f(x) - I_Xf(x) = f[x_0, \\dots, x_n, x] \\, \\ell_X(x) \\nonumber\n", "\\end{align}\n", "\n", "where $\\ell_X(x) = (x - x_0) \\cdots (x - x_n)$ is the node polynomial. Hence, using (ii), write down an exact formula for the error $\\frac1x - I_Xf(x)$ when $f(x) = \\frac1x$ (and $x, x_0, \\dots, x_n \\not= 0$)." ] }, { "cell_type": "markdown", "id": "667c2952", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "17786875", "metadata": {}, "source": [ "
\n", "\n", "***Exercise 4. (Newton form vs. Barycentric formula)*** \n", "\n", "(i) ✍️ Compute the Barycentric weights $\\lambda_0, \\lambda_1, \\lambda_2$ for $X = \\{1, 2, 4\\}$ and write down the Barycentric formula for the interpolant of $f(x) = \\frac1x$ from Exercise 3 (i)." ] }, { "cell_type": "markdown", "id": "61b6d220", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "42170255", "metadata": {}, "source": [ "
\n", "\n", "(ii) πŸ’» The code below implements ```divided_differences(x, y)```, which, given the nodes ```x``` $= (x_0, \\dots, x_n)$ and the data values ```y``` $= \\big( f(x_0), \\dots, f(x_n) \\big)$, returns the vector of Newton coefficients $\\big( f[x_0], f[x_0,x_1], \\dots, f[x_0,\\dots,x_n] \\big)$ by computing the divided difference table column by column (using $O(n^2)$ operations), and ```newton_eval(c, x, t)```, which evaluates the Newton form at $t$ using $O(n)$ operations by nested multiplication: $c_0 + (t - x_0)\\Big( c_1 + (t - x_1)\\big( c_2 + \\cdots \\big) \\Big)$.\n", "\n", "Run it, and use it to check your answer to Exercise 3 (i)." ] }, { "cell_type": "code", "execution_count": null, "id": "48c6169f", "metadata": {}, "outputs": [], "source": [ "# Newton coefficients (f[x_0], f[x_0,x_1], ..., f[x_0,...,x_n]): the divided difference\n", "# table is computed column by column, overwriting c in place. Here y[j] = f(x[j]) are the data values\n", "function divided_differences(x, y)\n", " c = float.(copy(y))\n", " n = length(x)\n", " for k ∈ 1:n-1 # column k of the table\n", " for j ∈ n:-1:k+1 # bottom to top, so that c[j-1] still holds column k-1\n", " c[j] = (c[j] - c[j-1]) / (x[j] - x[j-k])\n", " end\n", " end\n", " return c\n", "end\n", "\n", "# evaluates the Newton form c_0 + (t - x_0)(c_1 + (t - x_1)(c_2 + ...)) by nested multiplication\n", "function newton_eval(c, x, t)\n", " p = c[end]\n", " for k ∈ length(c)-1:-1:1\n", " p = c[k] + (t - x[k]) * p\n", " end\n", " return p\n", "end\n", "\n", "# test: g(x) = x^3 on 4 nodes has divided differences g[x_0,...,x_3] = 1 and is reproduced exactly\n", "g(x) = x^3\n", "X = [1, 2, 4, 1/2]\n", "c = divided_differences(X, g.(X))\n", "println(\"Newton coefficients of x^3 on X: \", c)\n", "println(\"p(3) = \", newton_eval(c, X, 3), \", g(3) = \", g(3))" ] }, { "cell_type": "markdown", "id": "d21b38d3", "metadata": {}, "source": [ "
\n", "\n", "(iii) πŸ’» Let $X_n := \\{ \\cos\\frac{j\\pi}{n} \\}_{j=0}^n$ be the Chebyshev nodes of the second kind. The code below computes (an approximation of) $\\max_{x\\in[-1,1]} |f(x) - p(x)|$ for $f(x) = \\frac{1}{1 + 25x^2}$ and $f(x) = e^x$ and $n = 10, 20, 40, 80, 160$, where $p = I_{X_n}f$ is evaluated using *(a)* the Newton form from (ii), and *(b)* the Barycentric formula (```ChebyshevInterpolant``` from Lecture 8). Run it." ] }, { "cell_type": "code", "execution_count": null, "id": "bde8d5d7", "metadata": {}, "outputs": [], "source": [ "# Barycentric formula for Chebyshev nodes of the second kind (from Lecture 8)\n", "function ChebyshevInterpolant( f, n )\n", " x = @. cos( Ο€*(0:n)/n )\n", " y = @. f( x )\n", " Ξ» = @. (-1.0)^(0:n)\n", " Ξ»[1] /= 2\n", " Ξ»[end] /= 2\n", "\n", " return function ( t )\n", " num = 0.0\n", " den = 0.0\n", " for j ∈ 1:(n+1)\n", " if t == x[j]\n", " return y[j]\n", " end\n", " num += Ξ»[j] * y[j] / (t - x[j])\n", " den += Ξ»[j] / (t - x[j])\n", " end\n", " return num / den\n", " end\n", "end\n", "\n", "Runge(x) = 1/(1 + 25x^2)\n", "egrid = range(-1, 1, length=2001)\n", "\n", "for (name, f) ∈ [(\"1/(1 + 25x^2)\", Runge), (\"exp(x)\", exp)]\n", " println(\"f(x) = \", name)\n", " @printf(\"%5s %16s %16s\\n\", \"n\", \"error (Newton)\", \"error (Bary)\")\n", " for n ∈ [10, 20, 40, 80, 160]\n", " x = @. cos( Ο€*(0:n)/n )\n", " coeffs = divided_differences(x, f.(x))\n", " p = ChebyshevInterpolant(f, n)\n", " err_newton = maximum( abs(newton_eval(coeffs, x, t) - f(t)) for t ∈ egrid )\n", " err_bary = maximum( abs(p(t) - f(t)) for t ∈ egrid )\n", " @printf(\"%5d %16.2e %16.2e\\n\", n, err_newton, err_bary)\n", " end\n", " println()\n", "end" ] }, { "cell_type": "markdown", "id": "e3b69812", "metadata": {}, "source": [ "
\n", "\n", "(iv) ✍️ Describe what you observe in (iii). Explain why the two methods give different results even though, in exact arithmetic, they compute the same polynomial. " ] }, { "cell_type": "markdown", "id": "f92471b4", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "7c1e5a90", "metadata": {}, "source": [ "
\n", "\n", "(v) πŸ’» ✍️ The code below reorders the nodes $X_n$ using the *Leja ordering*: $x_0$ is a node of largest absolute value, and each subsequent $x_k$ is chosen among the remaining nodes to maximise $\\prod_{i x_1 > \\cdots$? What about with the Leja ordering?" ] }, { "cell_type": "code", "execution_count": null, "id": "b3d84f21", "metadata": {}, "outputs": [], "source": [ "# Leja ordering: start at a node of largest |x|, then repeatedly choose the remaining node that\n", "# maximises the product of distances to the nodes already chosen (sums of logs avoid overflow)\n", "function leja_order(x)\n", " x = float.(collect(x))\n", " i = argmax(abs.(x))\n", " order = [i]\n", " s = zeros(length(x)) # s[j] = sum of log|x[j] - x[i]| over chosen nodes i\n", " for k ∈ 2:length(x)\n", " s .+= log.(abs.(x .- x[i]))\n", " s[order] .= -Inf\n", " i = argmax(s)\n", " push!(order, i)\n", " end\n", " return x[order]\n", "end\n", "\n", "for (name, f) ∈ [(\"1/(1 + 25x^2)\", Runge), (\"exp(x)\", exp)]\n", " println(\"f(x) = \", name)\n", " @printf(\"%5s %16s %16s %16s\\n\", \"n\", \"error (Newton)\", \"error (Leja)\", \"error (Bary)\")\n", " for n ∈ [10, 20, 40, 80, 160]\n", " x = @. cos( Ο€*(0:n)/n )\n", " coeffs = divided_differences(x, f.(x))\n", " xL = leja_order(x)\n", " coeffsL = divided_differences(xL, f.(xL))\n", " p = ChebyshevInterpolant(f, n)\n", " err_newton = maximum( abs(newton_eval(coeffs, x, t) - f(t)) for t ∈ egrid )\n", " err_leja = maximum( abs(newton_eval(coeffsL, xL, t) - f(t)) for t ∈ egrid )\n", " err_bary = maximum( abs(p(t) - f(t)) for t ∈ egrid )\n", " @printf(\"%5d %16.2e %16.2e %16.2e\\n\", n, err_newton, err_leja, err_bary)\n", " end\n", " println()\n", "end" ] }, { "cell_type": "markdown", "id": "e5a2c7d6", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "30ba0bb8", "metadata": {}, "source": [ "
\n", "\n", "\n", "## D: Reflection [20 pts] {.unnumbered}\n", "\n", "**Goal:** Reflect on your learning strategies, track your workload, and help us improve future course material. \n", "**Time:** take 5-10 mins (this shouldn't take too long)\n", "\n", "Estimate the approximate hours you spent: \n", "* **Assignment:** ____ \n", "* **Reading notes/materials:** ____ \n", "* **Exercises:** ____ \n", "\n", "(Reminder: you should be spending around 8 hours per week on this course outside of the 4 hours of lectures)\n", "\n", "
Any other comments: \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "Take a moment to reflect on your learning over the past two weeks:\n", "\n", "What was the easiest topic or task and why? Where did you struggle, and how did you overcome those challenges? Which resources (e.g. lecture notes, practice exercises, peer discussions, ...) helped you the most? How did you manage your time and workload (e.g., all at once, spread out over several days, or working with peers)? Looking back, what would you do differently if you were starting these two weeks over today?\n", "\n", "
Comments: \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(Optional) Any changes that would improve the lectures, exercises, or assignments for the future? You can also suggest things anonymously \n", " \n", " here\n", " (opens in a new tab)\n", " \n", "\n", "
Comments: \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.12.6", "language": "julia", "name": "julia-1.12" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.12.6" } }, "nbformat": 4, "nbformat_minor": 5 }