{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "2d69121f",
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: verify download of index files...\n",
"└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/MatrixDepot.jl:139\n",
"┌ Info: reading database\n",
"└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:23\n",
"┌ Info: adding metadata...\n",
"└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:67\n",
"┌ Info: adding svd data...\n",
"└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:69\n",
"┌ Info: writing database\n",
"└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:74\n",
"┌ Info: used remote sites are sparse.tamu.edu with MAT index and math.nist.gov with HTML index\n",
"└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/MatrixDepot.jl:141\n"
]
}
],
"source": [
"using FundamentalsNumericalComputation\n",
"FNC.init_format()"
]
},
{
"cell_type": "markdown",
"id": "468aabaf",
"metadata": {},
"source": [
"(section-localapprox-integration)=\n",
"# Numerical integration\n",
"\n",
"In calculus you learn that the elegant way to evaluate a definite integral is to apply the Fundamental Theorem of Calculus and find an antiderivative. The connection is so profound and pervasive that it's easy to overlook that a definite integral is a numerical quantity existing independently of antidifferentiation. However, most conceivable integrands have no antiderivative in terms of familiar functions.\n",
"\n",
"(demo-int-antideriv)=\n",
"```{proof:demo}\n",
"```\n",
"\n",
"```{raw} html\n",
"
\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"The antiderivative of $e^x$ is, of course, itself. That makes evaluation of $\\int_0^1 e^x\\,dx$ by the Fundamental Theorem trivial."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "846ccf24",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.718281828459045"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"exact = exp(1)-1"
]
},
{
"cell_type": "markdown",
"id": "ac214035",
"metadata": {},
"source": [
"```{index} ! Julia; quadgk\n",
"```\n",
"\n",
"The Julia package `QuadGK` has an all-purpose numerical integrator that estimates the value without finding the antiderivative first. As you can see here, it's often just as accurate."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "396c3a49",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Q = 1.718281828459045\n"
]
}
],
"source": [
"Q,errest = quadgk(x->exp(x),0,1)\n",
"@show Q;"
]
},
{
"cell_type": "markdown",
"id": "d1607fa0",
"metadata": {},
"source": [
"The numerical approach is also far more robust. For example, $e^{\\,\\sin x}$ has no useful antiderivative. But numerically, it's no more difficult."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "0b81d925",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Q = 1.6318696084180515\n"
]
}
],
"source": [
"Q,errest = quadgk(x->exp(sin(x)),0,1)\n",
"@show Q;"
]
},
{
"cell_type": "markdown",
"id": "6706865d",
"metadata": {},
"source": [
"When you look at the graphs of these functions, what's remarkable is that one of these areas is basic calculus while the other is almost impenetrable analytically. From a numerical standpoint, they are practically the same problem."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "8c03d245",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAFACAIAAADF7xm+AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1QU9/0//vfuIiAX5SpoJcCi4gVM5HIE1GgRUxI1GC8xp37SqEm9JN+YpK01tYnW5AQlpw3xEiNaNalNvEBroib1QpWfgmiE9YrxBuyFy4LAAju7y+6y8/79MbqhLPcwLMw+HyfHMzO8d/bFZngyvGfm/RZRSgkAAAiL2N4FAABA70O4AwAIEMIdAECAEO4AAAKEcAcAECCEOwCAACHcAQAECOEOACBACHcAAAFCuAMACFC/Dne1Wl1eXt5pM4vF0gfFCAmlFMNOdBcOsx5gWdbeJQwwLMv21s9mvw73L774YseOHZ020+v1fVCMkJhMJrPZbO8qBhgcZj2AD627jEZjb/1GdOqVvQAAQI/VNJELavb/U1OTWbw9oXf2iXAHALADJUMvqGleFT2vprc1j/pinETi9DgikfTC/hHuAAB9gaXkdv3jQK+kKl0bfeu+rr32dgh3AAC+GC2koIbmqmluFZtXRTXGTtq7SohTL10J5SvcjUbjnTt3zGbz2LFjPTw8bBuo1WrrxRYnJ6cnnniCp0oAAPoSYyaXqh+leV4VNTTbpwxewj07O3vx4sVBQUEuLi4lJSUHDhxITk5u1WblypX5+fmenp6EkICAgIsXL/JRCQBAH1AbyJWHbF4VzVXTHx5Scz+4BZSXcA8ODpbJZMHBwYSQ7du3r1q1Si6X2zb7+OOPly5dykcBAAC8slBSpKF5VTS/iuZW0VJtv3twhJdwHz16tHU5JiampqaGUioSiVo1Yxjm/v37wcHBzs7OfJQBANCLGkwkv5rmV7H51fRSNdX272dFeL+g+tlnny1YsMA22QkhaWlp27dvLy8v37hx49q1a9t8uVwuP3PmDLfs5eUVHR1t24ZlWTwI1y0sy4pEInxo3YLDrAcE8KHdbaCXqsnFanqpmtyupyz/J+jcj2fHbcTizq+68hvun3zyyQ8//JCfn2/7pb179/r5+RFCCgoKZsyYERcXN23aNNtmly9frqqq4pYjIiLGjRtn28ZgMEh65b5Qh2E0GkUiUXOznS70DEw4zHrAYDB0JYb6Fb1FJKsVXa4VX34o+qFWVGvsJGd7F8uyen1Tp0eam5tbpx8sj+G+a9euHTt25OTk+Pr62n6VS3ZCSExMzMyZM3Nzc9sM98WLF2/evLnjN6KUtnlDDrRn0KBBIpEIvWHdgsOsZ/r/h0YJuVtPLz+kl6rppWp6S0Ob7ffHhlgs9vDw6JXTCL7Cff/+/Zs3bz537lyn9zhSShUKxezZs3mqBACgFY2RXKqmlx+yl6vp5Yed334+EPES7qdPn37ttdeWLVuWmZnJbVmzZs3gwYPT09PPnj17/PhxhmHefPPNxMREV1fXzMzM2trahQsX8lEJAAAhxMySW5pH5+aXq+m9BuEPi8pLuA8ZMoS7QKrRaLgt3CCWkyZN4v5Gc3FxGTt27OnTp1mWjYyM3Llzp4+PDx+VAIDDetBIf3hIrzykPzykV2v/92EiESFCT3dewj0uLi4uLs52+4wZM2bMmEEIGTRo0Lp16/h4awBwWA+byA8P6Q/V7JUa+kM1re2gs0XoyU4wtgwADFxaM5HV0IKaR6fn/fBJIjtCuAPAgKFvJtdqaUENLXhIC2ro3Ya+uPF8gEK4A0D/ZWLJ9RZpfrvenvcpDiwIdwDoR0wsuVlHC2uorIYW1tAbddSENO8RhDsA2FOThdx4nOayWnqzrl8MqSgACHcA6FP6ZnLdmuY1tAg9LfxAuAMAvxpMj3paCmvo7Xp6Ez0tfQLhDgC9TKWjV2vo1VpSUDXoZkOzgsEdLXbAS7gXFRUdOHDg1q1b7u7u8+fPX7x4sW2bxsbGjz76SCaTjR8//r333vP39+ejEgDgm4WSew30Wi29yv1X0/LpIbFDPC/UL/ES7l999ZVYLF69enVtbe2aNWsaGhpWrFjRqs0rr7xCCHn//ff379//wgsv5Obm8lEJAPS6BhO5Xkdv1NHrtfRaLS2qt9s0odABXsI9NTXVulxeXn706NFW4V5aWvr999+r1Wpvb++4uLiAgIArV67ExsbyUQwA/ByUkFItvVZLr9fSG3Xkeh0eBB0YeO9zLy4uDgoKarXx6tWr4eHh3t7ehBBnZ+eYmJjCwkKEO0B/0GgmN+vozTp6o47eqKM362hjq/nkHGDULQHgN9yzs7OzsrKuX7/eantVVRWX7BxfX1+1Wt3mHr788suTJ09yy0899dS2bdts2zAM00v1OgpuJiZM1tEtQj3MLJQUa0W36kVFDeKiBnFRvUih62zuISQ7b1iW1Wp1XZmJqdM2PIb75cuXlyxZ8q9//Ss4OLjVlzw8PJqamqyrer3e09OzzZ0kJye/8cYb3LKnp2d7zdrbDm1ydnZGuPeAMA6zagPhTslvaeiNOlqkoU0We9cEj4nFYk9Pz349E1NhYWFKSsr+/ftnzpxp+9WgoCC5XM6yLDcNYGlp6csvv9zmfgICAtqcFBsAuqLBRIo09JaG3tLQIg29WUcfNnX+KhAAXsL9+vXrc+bMycjIeO6551pu//bbbwMDAydPnjx16lRnZ+djx47NmzcvPz9fqVQ+++yzfFQC4FCaLOTHenpLQ2/VcWlOcI+5w+Lrbpmampply5Zxq+PHj+fudNy5c2d8fPzkyZOdnJw+//zzZcuWbdmy5f79+5999ln/n0UXoL8xWsiP9ZRL8x/ryS0NLWmkllZhjoufjoqXcN+3b9+uXbusq9b+o2PHjnH9MISQOXPmyOXyBw8ehIaGDhkyhI8yAITEzBKVjhZp6G0NKdLQwhp6t8Emym0h2R0VL+Hu7u7u7u5uu93FxaVVsyeffJKPAgAGOl0zuVNPf6yntzX0TgMp0tASLQbYgm7A2DIA9qcxkh/r6e16eqeeFmnonQai0FKcc8PPgXAH6GsqHb1bT+42PErz2xqqNti7JhAchDsAj4wWcr+R3qmndxvIj/XcAmXMnb8Q4GdCuAP0modN5G49vdNAuX/v1JNSbReueQLwAOEO0BNGC7nXQO810HuN5G49vdtA7zXQOmPnLwToGwh3gE6wlCgZeqNKUqZi7zZwOU4UDGVxSg79GMId4H9U6On9BnK/kd5voPcbyb0G+qCRGi2EkEGEYBAWGDB4DHeLxdLY2Ojp6enk1Ma7MAxjNj+6riQWi4cOHcpfJQBtethEHjTS+w30XgO930juN9D7jbjaCQLBV7hPnz69sLBQp9OdPXv2l7/8pW2DJUuW5Ofnc8PsBQQEXLx4kadKAAghVYZHOf6gkT5oJA8a6YNG2mCyd1kAvOEr3N99991JkyZFRUV10Objjz9eunQpTwWAY6KElOtocSMp1tIHDT/luBbn4+BgOgp3jUbTckqNbunKKI96vV6pVI4YMaLNfhuAjhktpFRLS7SkuJEWa2lxIy1uJKVajE4OQEgH4b5u3br6+vqMjIzMzEypVMrHoOoffvjhX//615qamg8++ODtt99us41CocjOzuaWvb29Mba7Y6ozkhItLWmkxVpS3PhooUyH+1UA2tVuuMfExEilUkLIokWLMjIyej1Vd+/eHRAQQAi5dOlSUlJSbGzslClTbJtdunTJOgNfREREeHi4bRudTicSdTYxGLTQb6fZM7FEqRPJGZFcJ5IzRK4TyxmRnCENZvz/BYfAsizD6LsyzZ51hN32tBvudXV1W7ZsIYQkJSX1oMROcclOCImLi0tMTDx//nyb4b548eLNmzd3vCtKKYaD75ZBgwbZN9wpIZV6WqolpdpH/5ZoaSlOxsHhicViDw8PfqfZ8/HxKSwsJIRkZmb+/LfpAKVUpVLNnj2b13cBe3nYREq1VK6lpVoiZ2iplsq1RM5wd44DAF/aDXcvL6+MjIyVK1d6eXnJZLLu7nf//v3V1dUMwxw8ePCHH3549dVX/fz8tm7dmpOTc/ToUYZh3nnnncTExMGDBx8+fLi6unrhwoU/7xsBO3vYRORaqmConCEKLZUzVK4lpVqqa7Z3ZQAOqd1wnzVrlkajIYT4+Pj0oGemsbFRo9G8/vrrhBCNRsOyLCEkMjKS6wpwcXEJCQk5duyYyWSaMGHCp59+6uvr2/NvAvoKJUStJwqGKhiqYB6nuZbIGapHiAP0JyLatSkBCgsLrddUWy7zasuWLQ0NDZ32uWu1Wu5hKOiiTi+oGi1EpaNKhigZqmCIgqHcgkqH7hQAHgV7kOJFYn773FtKS0uTSqVHjhxJS0sjhNTV1fVZvgOvHjYRFUNVOqpguIVHaV6pxzRAAANbl8JdKpUuWrRIKpUWFxeHhYXNmjUrMzMT4T5Q6JuJgqEqHSnTUSVDlQxRaEVlOpFKbzagLwVAoLoU7gUFBYsWLYqOjj5z5kxYWBghxMvLi+fCoHuaLKRMR8t1RKmjZTpSpqMqhqh0VMXQ2jYGGcdt4wAC16Vwf/HFF6Ojo999913rFplMNmvWLN6qgrYZmkmZjlboiUpHy3SkXEeVOlKmo2U6Wo1JOAGghS6Fe3R09O7du5OSkurr67nHVrnnm4AP9SZSrqMqHanQUxVDyvW0QkeVOlKha/McHACgDV0dsSs6Olqj0RQXF5eUlMTExPR4QDEghJhZUmWgKh2p1D86Aa/QkzIdVRuIEvcUAkBv6N5wjGFhYVyfO3Sq2kDUBlqme/yvnpbpSIWeVuipWk9wLwoA8Apj7fZcTROpMtByPVHrabmeVOppeYsoN7H2rg8AHBgv4a7T6bKzswsKCqqrq3ft2tXmkI0Mw6Smpl69enX8+PHr16/vh0+ospRUNz0K7ioDLdeRagMt0z0K9CoDHucBgP6Ll3C/e/duenr6E088ceDAgc8//7zNcF+6dKnJZFq7du0XX3zxwgsvnD9/no9KOsaYSYWecv0nlXpS3UQrdKTKQNUGwq024+wbAAYmXsI9KioqJyfn7t27Bw4caLOBXC4/fvx4ZWWlj4/P1KlThw0bxscjryaWVBtohZ5UGWiVgVToyUMDLdeT6kfxjUuXACBY9ulzl8lkY8aM8fHxIYQ4OzvHxsYWFBT0ONwL68RKNVtlIJV6Wt1EKnS0ykCqDLhxEAAcl33Cvaqqikt2jq+vr3W6pVYOHz58+fJlbnnChAmpqam2bTZclVyoQf83AAx4LMtqtbquzMTUaRv7hLu7u3tTU5N1Va/XtzeVUlxc3PLly7llb2/vNkd/1LI4RQcAIRCLxZ6enn03KmSvCwoKUigUlFLuWqtcLv+///u/NlsGBwd3Opr8KA96ra73iwQAGLg6mWK1d3333XcFBQWEkGnTpkkkkhMnThBCLl++LJfLn3322b6sBABA2Hg5c29sbAwJCbFYLIQQPz8/Pz+/e/fuEUK2bdsWHx8fExPj5OS0c+fOZcuWhYeH37lzZ/v27ZhtAwCgF/ES7p6ensXFxdZV633uR48etfYlpaSkJCYm3rt3TyqVYqQaAIDexUu4i0SiNvPazc2t5aqnpydm/AAA4EOf9rkDAEDfQLgDAAgQwh0AQIAQ7gAAAoRwBwAQIIQ7AIAAIdwBAASIx7Fl6urqSktLw8LCvLy8bL/KMIzZbOaWxWLx0KFD+asEAMDR8HXm/tVXX40ePfrtt98eNWpUVlaWbYMlS5YEBQVxM24nJCTwVAYAgGPiJdz1ev2bb77573//+8KFC19//fUbb7xhMplsm+3YsaOurq6urq6oqIiPMgAAHBYv4X769Gk/P7/p06cTQp555hlXV9dz587ZNjMajSqVimUxUSkAQC/jJdyVSmVoaKh1NTQ0VKlU2jZ777334uPjfXx8duzY0d6uFApF9mOFhYV8VAsAIDy8XFDV6/UuLi7W1cGDB+t0ulZtdu7c+Ytf/IIQkpeX98wzz0RHR8fHx9vu6tKlS9YZ+CIiIsLDw23bWCwiQnph4hIAAPtiWZZh9F2ZZk8s7uTUnJdwDwgIqKv7aW6k2trawMDAVm24ZCeETJkyZebMmTk5OW2G++LFizdv3tzx20kkTR03AAAYEMRisYeHR69Ms8dLt0xMTMzVq1cNBgMhRKvV3rx5s4OhfSml5eXlGNIdAKAX8RLukZGRCQkJy5Yt++9//7t06dKkpKTRo0cTQrZv375w4UJCCMMwr7/+elZW1vHjx3/zm99UVFQsWLCAj0oAABwTXw8xZWZmpqampqenT5w48U9/+hO3cdy4cdysTM7OzgEBAQcPHmxubo6IiJDJZP7+/jxVAgDggPgKdy8vr48//rjVxqSkpKSkJEKIs7Pzxo0beXprAADA2DIAAAKEcAcAECCEOwCAACHcAQAECOEOACBACHcAAAFCuAMACBDCHQBAgPgK97t3765YsWLevHm7d++mlNo20Ol0GzdunDdv3h//+EeNRsNTGQAAjomXcG9oaHj66adHjBjx29/+Nj09fevWrbZtli9fXlBQsHr16rKysvnz5/NRBgCAw+Jl+IEDBw6MHTv2L3/5CyFEIpGsWrVqzZo1LUcfVigU33zzTXl5uZ+f34wZM4YNGyaTyaKiovgoBgDAAfFy5l5QUDB16lRueerUqQqForq6umUDmUw2ZswYPz8/QoiLi0tsbOyVK1d69l537949derUzywYAKA/qK2t/fLLL3tlV7ycuVdVVU2cOJFb9vDwcHV1raysbDlfh1qt9vHxsa76+vpap1tq5fDhw5cvX+aWIyIiUlNTWzXQ6/WDqu/98pfJvfkNCF1xcbFEIgkJCbF3IQPJxYsXo6OjW04xBh3TaDRyuXzSpEn2LmQgKb1+zyg2MgzTcTO7zcTk5ubW1PRodiSLxWIymTw8PFo2cHd3NxqN1lWDweDu7t7mruLi4pYvX84te3t7t9oP916+l/ecPbC216p3AO+++w93d/f3/9/79i5kIBm56uXDr18aOXKkvQsZMM6evf3B3g/Ovp9j70IGkt8c+bvLjBm2QdcDvIR7UFCQXC7nlpVKpUgkGjFiRMsGI0eOVCgUlFJueHeFQhEUFNTmroKDg7lRggEAoOtEbd6n+DPl5eUtWLDg1q1bfn5+GzZsuHbt2rFjxwghJ0+eHDZsWFRUlNlsDg4O3rNnz+zZswsKChITE8vKyoYMGdJqP4cOHUpNTQ0ICOjgvYxG44MHDyZMmNDr34WAlZeXSyQS24ltoQM3b94cO3bsoEGD7F3IgKHVaisrK8eMGWPvQgYSuVzu6enp6+vbcbNdu3aFhYV13IaXM/cpU6bMnz9/4sSJoaGhKpXq5MmT3Pb09PT4+PioqKhBgwbt2LHjlVdeiYiIuHXrVnp6um2yE0JeeumlYcOGsSzLR5EAAAPUsGHDOm3Dy5k7Ry6X19TUREZGWq9BNTY2Ojk5ubm5cav19fX37t2TSqXcbTMAANBbeAx3AACwF4wtAwAgQAh3AAAB4uWCau+yWCw3btxwdXUdN25ce21KSkpqa2snTpzY8hkTo9F448YNPz+/0NDQPqm0H9HpdEVFRcOHD2/vHtOHDx/K5fLhw4e3vHG7vr7e2k3n7Ozc3sMHQvXgwYP6+vonn3yyzVtiWg5v5+LiYr10RCktKiqyWCyRkZGdPlciMCzL3rhxw9nZefz48bZfNRgM1uddON7e3oQQhmHMZjO3RSwWDx06tA9K7T+ampqampq8vLzaa1BTU1NSUhIWFtbynpmeHGa0f6uuro6MjJw0adKoUaOSk5Obmpps26xYsWLEiBEJCQlBQUG3b9/mNhYVFY0cOTIhIWHEiBGrVq3q26rt7OLFiwEBAdOmTfP399+wYYNtg5SUFB8fn8mTJ/v6+r744otms5nb7ubm9sQTT0ilUqlU+oc//KFvq7YnlmVffvnloKCguLi40NDQ4uJi2zYikSg4OJj7cN577z1uo06ne/rpp8eOHTtx4sSYmBiNRtO3hdtTTU3NU0899dRTT40ZMyYpKclgMLRqkJqaKn3Mz89v6NChLMtSSufOnevv789tT0hIsEft9iGTySIjIyUSibu7e3ttDhw44OPjM336dB8fn6+//prbqNPppk+f3t3DrL+H+9q1axcsWEApNRqNkyZN2rdvX6sGFy5cCAwMrKmpoZSuX79+3rx53Pbnn39+/fr1lNKamprAwMDc3Ny+LdyeJk+evHXrVkqpQqHw9PS8f/9+qwanTp0ymUyUUo1GExwc/OWXX3Lb3dzcVCpVH1fbH5w8eTI4OLihoYFS+tZbby1ZssS2jUgk4g6zlrZu3ZqQkGA2m1mWffbZZz/44IO+KLd/WL9+fUpKCsuyJpMpNjY2IyOjg8avvPKK9Rxr7ty5X3zxRZ/U2L+Ul5efP38+JyenvXDX6/Xe3t45OTmU0rNnz/r5+XG/Mrdt22Y9zJ577rlNmzZ15e36e7gHBwd///333PLf/va35OTkVg3WrFmzcuVKbvnBgwdOTk56vZ5hGIlE8uDBA277ypUr33rrrT6r2b7kcrlEIuE6WCilKSkpmzdv7qD9r371q7S0NG7Zzc3t2rVrFRUVvFfZzyxfvvx3v/sdt3z9+nVXV9fm5uZWbUQi0Z07d9RqdcuNU6ZMsYbakSNHIiIi+qDafmLUqFHffvstt7xt27aZM2e211Kr1Xp4eFy5coVbnTt37s6dOxUKhfVPRodSUFDQXrgfO3YsLCzMuhoaGvrdd99RSqdOnbpr1y5uY2Zm5oQJE7ryRv26i5Bl2fLycmuPeWhoqFKpbNVGqVRaG4SEhHAvqaioYFk2ODi4gxcKlVKp9PX1tfZjdvy93759Oz8/f86cOdyqWCxOSUl58sknQ0NDc3Jy+qDafqLlURQaGtrU1NRqHFNCiEQiSU5OjoiIGDVqVF5eXpsvdJzDjFJaVlbWxe/966+/DgkJiYmJsW7ZtGnTtGnTfHx82pzswWEplUqpVGpdDQkJUSgUxOYw4zZ2ql9fUDWZTM3Nzc7Oztyqq6urTqdr1cZgMFgvokokEicnJ66Nk5OTk5NTBy8UKr1eb/3ECCGurq51dXVttnz48OGCBQs2bNhgvRp2//79wMBASml6evqLL74ol8utlw2FreVR5OrqSgixPWBUKhX34WzevHnRokWlpaUuLi6tXqjX6+njEZOEzWKxGI3Gjn82rfbt2/faa69ZV3fv3s0NfZGfn5+UlBQbG5uQkMB3wQNCqx/ewYMH6/V6YnN8cpceOz3M+vWZu6ur69ChQ63ZVFtbO3z48FZtAgICrA20Wq3JZBo+fHhAQIDZbG5sbLS+0HHGUQkMDOSut3Cr7X3vdXV1s2bNWrhw4e9///uWryWEiESit99+m7vfpm9qtruWR1FtbS0hxPZIs344a9eurampuX//vu0LAwMDHSHZCSFOTk6+vr4d/2xyioqKrl69+utf/9q6xXpAxsfHJyYmnj9/nu9qB4rAwMCWp2Lc9UJic5gFBAR05TDr1+FOCImNjc3NzeWWc3NzY2NjWzWIiYlp2SAkJMTf33/YsGHBwcHWv53bfKFQjR492snJ6erVq9xqXl6e7ffe0NCQnJycmJj44YcftrmThw8fNjU1tRxzX9haHWbjx4/v4DZQtVptNpu5Dyc2NvbChQvWFzrOYUa68LPJ2bt377x58/z9/W2/xPXtOM5h1qmYmJhr165xg7lrtdobN25wfVld/Khb6/l1gT5x4sQJPz+/zMzMjIyMIUOGcHc6GgwGf3//oqIiSqlGo/H399+0adOJEyfCw8M//fRT7oWffPJJeHj4iRMnNm3a5O/v71D3qK1bty42NvY///nPO++8I5VKuRtjzp49a71WM3369DFjxmQ8xt1KlJOT8/7773/zzTf//Oc/o6KiZs+ezd245gjUarW3t3daWtrx48dDQ0P37NnDbZ89e/bOnTsppadPn964ceO33377j3/8Y+LEifPnz+caFBYWDhkyZN++fQcPHrTe5+AgTp065ePjc/jw4T179gwZMuTGjRuUUrPZ7O/vf+3aNa6N0Wj09/c/efKk9VWNjY0rVqw4dOjQ0aNHX3rppREjRtjegyRUDMNkZGT8+c9/dnFxycjIOHjwILfdephRSpOTkxcsWHDq1KkXXnhhzpw53EaZTNbyMDt37lxX3k7CzXTab40ZMyY8PPzgwYMVFRVbt261zupSVVWVmJjo6enp6uo6d+7cU6dOXbhwYenSpatXr+b+YImLi3N1dT106BDLsnv27HGoORZmzJih1+uzsrI8PDz+/ve/c0+OmEwmi8Uyc+ZMQkhZWVlYWJjmMT8/v/DwcEJIXl7euXPn5HJ5SkpKWlqa9aKF4Hl4eDz33HPff/99Xl7e6tWrly1bxm2vra2NiIgICQmhlObm5p49e1alUi1cuPCjjz6SSCSEkOHDh8fHx2dlZf34448ffvihQ809EBYWNn78+EOHDqlUqvT0dOvppFqtnjlzJjfOq0qlEolEr776qvW5G7FYfPPmzfPnz8tksnHjxu3du7fNk3pBMhgMx48fZ1k2Li5Oo9GYTKYpU6aQFocZIeT555//8ccfv/vuu4iIiPT0dK6rvWeHGQYOAwAQoP7e5w4AAD2AcAcAECCEOwCAACHcAQAECOEOACBACHcAAAFCuAMACBDCHQBAgBDuAAAChHAHABAghDsAgAA5yshQAD2j0WiOHDlCCKmvr4+KiiopKZHJZBkZGfauC6ATCHeAjuzevXvdunWEEG9v7yNHjkil0t27d9u7KIDOYVRIgI4UFhZGR0cTQkQi/LDAQII+d4COcMl+5swZhxqrHQQA4Q7QuezsbGu4FxYW2rcYgK5AuAO0KzMzc9GiRYSQrKwsqVRKCDlz5oyXl5e96wLoHC6oArRLKpVKpdK0tLTTp09z11G9vLzCwsLsXRdA53CNCABAgNAtAwAgQAh3AAABQrgDAAgQwh0AQIAQ7gAAAoRwBwAQIIQ7AIAAIdwBAAQI4Q4AIEAIdwAAAUK4AwAIEMIdAECAEO4AAAKEcAcAECCEOzx+YA8AAB8xSURBVACAACHcAQAECOEOACBACHcAAAHq1+GuVqvLy8s7bWaxWPqgGCGhlGJ6xe7CYdYDLMvau4QBhmXZ3vrZ7Nfh/sUXX+zYsaPTZnq9vg+KERKTyWQ2m+1dxQCDw6wH8KF1l9Fo7K3fiP063AEAoGcQ7gAAAoRwBwAQIIQ7AIAAOfG0X6PReOfOHbPZPHbsWA8PD9sGarXaerHFycnpiSee4KkSAAAHxEu4Z2dnL168OCgoyMXFpaSk5MCBA8nJya3arFy5Mj8/39PTkxASEBBw8eJFPioBAHBMvIR7cHCwTCYLDg4mhGzfvn3VqlVyudy22ccff7x06VI+CgAAcHC89LmPHj2aS3ZCSExMTE1NTZu35TMMc//+fZPJxEcNAACOjK8+d6vPPvtswYIFIpHI9ktpaWnbt28vLy/fuHHj2rVr23y5XC4/c+YMt+zl5RUdHW3bhmVZPAjXLSzLikQifGjdgsOsB/ChdRf3ibUZmC2JxZ2fl/Mb7p988skPP/yQn59v+6W9e/f6+fkRQgoKCmbMmBEXFzdt2jTbZpcvX66qquKWIyIixo0bZ9vGYDBIJJJeLVzgjEajSCRqbm62dyEDCQ6zHjAYDF2JIbAyGAzNzc2dHmlubm6dfrA8hvuuXbt27NiRk5Pj6+tr+1Uu2QkhMTExM2fOzM3NbTPcFy9evHnz5o7fiFLa5g050J5BgwaJRCJnZ2d7FzKQ4DDrGXxo3SKRSJydnXvlNIKvX6r79+/fvHlzdnZ2p/c4UkoVCkWbvwAAAKBneDlzP3369GuvvbZs2bLMzExuy5o1awYPHpyenn727Nnjx48zDPPmm28mJia6urpmZmbW1tYuXLiQj0oAABwTL+E+ZMgQ7gKpRqPhtnB3y0yaNIn7G83FxWXs2LGnT59mWTYyMnLnzp0+Pj58VAIA4JhE/Xlc7y1btjQ0NHTa567VarmHoaCLuAuq6HPvFhxmPcAwDPrcu6W60eDr3jt97rzfCgkAAC3pm4mCoUqGqHRUxVAFQ1Q6WqYjKh19wt2paH7vvAvCHQCg97GUqA1UwRAlQ1U6omKonCEqhqp0tKap3VcVa3utAIQ7AEDPMeZHp+FKHVUxVMkQBUNVOlKmo+buP7/1C7deKwzhDgDQuYdNXIhTBUMUDJVriZKhCobWGe1dWTsQ7gAAj7CUVOqpnCFyLVXqiEJLubPyUoYaBtoD3Qh3AHA4Fkoq9FSuJXKGyrXcmfij/nGTUMbCQbgDgGBZKCnXPToTL9USOUMVWipnetghPrAg3AFgwOO6U7j4LtUSuZZyp+QqBwjx9vAS7kVFRQcOHLh165a7u/v8+fMXL15s26axsfGjjz6SyWTjx49/7733/P39+agEAASmykBKtVSupaWPzsepnCEKrXC6U3oLL+H+1VdficXi1atX19bWrlmzpqGhYcWKFa3avPLKK4SQ999/f//+/S+88EJubi4flQDAAKUxktLH2f0ozbUD8sKmvfAS7qmpqdbl8vLyo0ePtgr30tLS77//Xq1We3t7x8XFBQQEXLlyJTY2lo9iAKA/a7KQUi0t0T46DS/VPsr0ekzR9vPw3udeXFwcFBTUauPVq1fDw8O9vb0JIc7OzjExMYWFhQh3AAGzUKJiaCnzKLutIV6pt3dlAsVvuGdnZ2dlZV2/fr3V9qqqKi7ZOb6+vmq1us09fPnllydPnuSWn3rqqW3bttm2YRiml+p1FBg4rAdwmHVdjVGk0InkDLmnsVSamxQ6kZwRlelFDntts+tYltVqdV2ZianTNjyG++XLl5csWfKvf/3LOlm2lYeHR1PTT8Mr6PX69sbbS05OfuONN7hlT0/P9pphuL5ucXZ2Rrj3AA6zVgzNpERLS7Xcv49Oxku0lDFbm+AY6x6xWOzp6dmvR4UsLCxMSUnZv3//zJkzbb8aFBQkl8tZluWmASwtLX355Zfb3E9AQECbk2IDQJ9hKSnT/RTiJVyIN1K1wd6VQft4Cffr16/PmTMnIyPjueeea7n922+/DQwMnDx58tSpU52dnY8dOzZv3rz8/HylUvnss8/yUQkAdAt3j0rL8/ESLW40HJD4ulumpqZm2bJl3Or48eO5Ox137twZHx8/efJkJyenzz//fNmyZVu2bLl///5nn32GEf0B+pKJJQrrOXiLKNf012GwoLs6CffCwsKCgoL6+novL6+YmJgu9pDs27dv165d1lVr/9GxY8e4fhhCyJw5c+Ry+YMHD0JDQ4cMGdKj4gGgc2oDKWls0Z2ipaVaUqajbP+dhA16QbvhnpmZmZ2dLZVKo6KiCCElJSVHjhzZvXt3UlLSokWLOt6pu7u7u7u77XYXF5dWzZ588skelQ0AremaH3WFPz4NJyVaWqLFUz8Oqu1wz8zMlEqlGRkZbX4pIyNj5cqVPBcGAG2z3jDOnY9be1SqcHkTWmg73JOSklreh95Sp6ftANBbapp+6kixLigZxx0MC7qu7XBvmewajaa9oAeAXtFkIT8leONPUa41d/5agDZ1ckF13bp19fX1GRkZXEcNbjkH+Dksj28YL/3fp37wCD70uk7CPSYmRiqVEkIWLVqUkZGBcAfoomrD4wRvMZoKelSgz3QS7nV1dVu2bCGEJCUl9Uk9AANMo7nVcIaPolyHe1TArjoJdx8fn8LCQkJIZmZmn9QD0E81WUh5vXWun0fT/ZRqaS2e+oF+qZNw9/Ly4m589PLykslk3dq1xWJpbGz09PR0cmrjXRiGMZsfXS0Si8VDhw7t1s4BeGJiiYp5FOLWiTdLtbRS70oIzsZhwOgk3GfNmqXRaAghPj4+3eqZmT59emFhoU6nO3v27C9/+UvbBkuWLMnPz+eG2QsICLh48WJ3ygb4uZpZUq6nci0p1VLrxJulWlKhpxY8ugkDX+djy3D3QXb3Uuq77747adIk7unW9nz88cdLly7t1m4Buou7QUXxaL7NR+fjcoaU6Wgzrm2CcLX7hGoHDyt15QnVrozyqNfrlUrliBEj2uy3AegW7kycC3E59+/jEMcNKuCA2k7VRYsWrVu3zsfHZ+HChWFhYYQQjUZTUlKSnZ1dV1eXlpbWK+/94Ycf/vWvf62pqfnggw/efvvtNtsoFIrs7Gxu2dvbG/digpklqsdn4gqGlmqJgqEKnIkD/K92T5nT0tKKi4uzsrKs2RoVFbVixQou63++3bt3BwQEEEIuXbqUlJQUGxs7ZcoU22aXLl2yzsAXERERHh5u20an04lEol6pykEMiGn2DM1EpRcpdaIyvUipIyq9WKEjKp240kDQJw5CxbIsw+i7Ms2edYTd9nTUHxIWFrZu3bp169Z1u8Au4JKdEBIXF5eYmHj+/Pk2w33x4sWbN2/ueFeUUgwH3y2DBg3qP+FebyJKhioYKn98Ds6tYhgscEBisdjDw6Ovp9njaZAZSqlKpZo9e3av7xn6D5aSSj2VM0TJUCVDlDqqZKhcS5QMbcTwKQA86CTcNRrNkSNH6uvrCSHZ2dlnzpzp4n73799fXV3NMMzBgwd/+OGHV1991c/Pb+vWrTk5OUePHmUY5p133klMTBw8ePDhw4erq6sXLlz4c78V6Ad0zUTBUBVDlAzlesYVDFXiqiZAn+sk3N99992oqChueBkvL6+u77exsVGj0bz++uuEEI1Gw7IsISQyMpLrCnBxcQkJCTl27JjJZJowYcKnn37q6+vb828C+paFErWeyhmiYqhKR1Q6KtcSlY4qGVqHxzUB+gcRpR1dnGp5T2Tfj/27ZcuWhoaGTvvctVot9zAUdFEXL6g+bCJlOqp6fC9KmY7rTiEVetyXAsCLYA9SvEjcR33u3GC/hJAjR4701k2Q0H80mh+dgJfpKPewT5mOqhii1GF6NoABrEvdMtxySUkJ//UALxgzUepomY6U66iSIfJGUYWelOmbVTpMBwEgTJ2E+5EjR6zPDRUXF/NfD/Rcg+lR50m5nqoYUq6nXJSX6Wi9qVVb7rEA3C4OIFhth7u1e93Ly8ua6VlZWTzd8w5dpzaQSj0t15EyHa3QU+WjECcqHWVwDg4Aj7Ud7jExMQUFBd7e3s8884x1MMiCggKEex8wWkjFT/FNVDpaoSflOlqmI5V6asKVTADogrbD3Xq2fvr0aet4A+iW6S2UkKrHJ+Dlelqpp2U6otZTlY6oDbQaT2YCwM/WSZ+7NdkLCwsxaFe31BlJpZ6W64laT8v1pGWCVxnwRA8A8KuTcF+5cuWKFSvefffdpKSk7OxsdMu0pDGSSgOt1JMK/aN/K3Q/LTdZ7F0fADiwTsJ94cKFUqm0pKRk3bp1XZ9GVafTZWdnFxQUVFdX79q1q80hGxmGSU1NvXr16vjx49evX98Pn1DlOk+qDbRM9/jfJlquI2oDrdATNeIbAPqxTsK9pKREJpNxA790/T73u3fvpqenP/HEEwcOHPj888/bDPelS5eaTKa1a9d+8cUXL7zwwvnz57tb+s+nayYVOlplIFUGWmkg1QZaoSdVBqrWk0oDqTLgOUwAGKg6CXdrb0xGRkbXdxoVFZWTk3P37t0DBw602UAulx8/fryystLHx2fq1KnDhg3jo0/fxP6U11UGUqlvtUp1eAITAASq8wuq3DVVbl69lhH8c+JYJpONGTPGx8eHEOLs7BwbG1tQUNDjvRXWiZVqlsvr6qafTsZrMYgVADiqboznnpaWJpVKrSPM1NXV9Tjfq6qquGTn+Pr6WqdbauXw4cOXL1/mlidMmJCammrbZsNVyYUa9H8DwIDHsqxWq+vKTEydtulGuEul0kWLFkml0uLi4rCwsFmzZmVmZvYs3N3d3Zuamqyrer2+vamU4uLili9fzi17e3u3OfqjlsUpOgAIgVgs9vT07JVRITuZha+lgoICQkh0dLT1ymq3RnhvKSgoSKFQWEcblsvlQUFBbbYMDg5Oeqy9XySjPDBGCgDA/+hGuL/44ovR0dGZmZncxEyEEJlM1q03++6777jfENOmTZNIJCdOnCCEXL58WS6XP/vss93aFQAAdKAb3TLR0dG7d+9OSkqqr6/nRnjfsmVLmy0bGxtDQkIsFgshxM/Pz8/P7969e4SQbdu2xcfHx8TEODk57dy5c9myZeHh4Xfu3Nm+fTtm2wAA6EXdCHdCSHR0tEajKS4uLikpiYmJaW9iJk9Pz5YD0Vjvcz969Ki1LyklJSUxMfHevXtSqbSPJ3gCABC87oU7x3p/ZHtEIlGbee3m5tZy1dPTE+PVAADwoRt97gAAMFAg3AEABAjhDgAgQAh3AAABQrgDAAgQwh0AQIAQ7gAAAtST+9y7qK6urrS0NCwsrM0haBiGMZvN3LJYLB46dCh/lQAAOBq+zty/+uqr0aNHv/3226NGjcrKyrJtsGTJkqCgIO55qISEBJ7KAABwTLyEu16vf/PNN//9739fuHDh66+/fuONN0wmk22zHTt21NXV1dXVFRUV8VEGAIDD4iXcT58+7efnN336dELIM8884+rqeu7cOdtmRqNRpVKxLCYqBQDoZbyEu1KpDA0Nta6GhoYqlUrbZu+99158fLyPj8+OHTva25VCoch+rLCwkI9qAQCEh5cLqnq93sXFxbo6ePBgnU7Xqs3OnTt/8YtfEELy8vKeeeaZ6Ojo+Ph4211dunTJOgNfREREeHi4bRuLRURIL0xcAgBgXyzLMoy+K9PsicWdnJrzEu4BAQF1dXXW1dra2sDAwFZtuGQnhEyZMmXmzJk5OTlthvvixYs3b97c8dtJJE0dNwAAGBDEYrGHh0dfT7PXdTExMVevXjUYDIQQrVZ78+bNDob2pZSWl5djSHcAgF7ES7hHRkYmJCQsW7bsv//979KlS5OSkkaPHk0I2b59+8KFCwkhDMO8/vrrWVlZx48f/81vflNRUbFgwQI+KgEAcEx8PcSUmZmZmpqanp4+ceLEP/3pT9zGcePGcbMyOTs7BwQEHDx4sLm5OSIiQiaT+fv781QJAIAD4ivcvby8Pv7441Ybk5KSkpKSCCHOzs4bN27k6a0BAABjywAACBDCHQBAgBDuAAAChHAHABAghDsAgAAh3AEABAjhDgAgQAh3AAAB4ivc7969u2LFinnz5u3evZtSattAp9Nt3Lhx3rx5f/zjHzUaDU9lAAA4Jl7CvaGh4emnnx4xYsRvf/vb9PT0rVu32rZZvnx5QUHB6tWry8rK5s+fz0cZAAAOi5fhBw4cODB27Ni//OUvhBCJRLJq1ao1a9a0HH1YoVB888035eXlfn5+M2bMGDZsmEwmi4qK4qMYAAAHxMuZe0FBwdSpU7nlqVOnKhSK6urqlg1kMtmYMWP8/PwIIS4uLrGxsVeuXOnZe929e/fUqVM/s2AAgP6gtrb2yy+/7JVd8XLmXlVVNXHiRG7Zw8PD1dW1srKy5XwdarXax8fHuurr62udbqmVw4cPX758mVuOiIhITU1t1UCv1w+qvvfLXyb35jcgdMXFxRKJJCQkxN6FDCQXL16Mjo5uOcUYdEyj0cjl8kmTJtm7kIGk9Po9o9jIMEzHzew2E5Obm1tT06PZkSwWi8lk8vDwaNnA3d3daDRaVw0Gg7u7e5u7iouLW758Obfs7e3daj/ce/le3nP2wNpeq94BvPvuP9zd3d//f+/bu5CBZOSqlw+/fmnkyJH2LmTAOHv29gd7Pzj7fo69CxlIfnPk7y4zZtgGXQ/wEu5BQUFyuZxbViqVIpFoxIgRLRuMHDlSoVBQSrnh3RUKRVBQUJu7Cg4O5kYJBgCArhO1eZ/iz5SXl7dgwYJbt275+flt2LDh2rVrx44dI4ScPHly2LBhUVFRZrM5ODh4z549s2fPLigoSExMLCsrGzJkSKv9HDp0KDU1NSAgoIP3MhqNDx48mDBhQq9/FwJWXl4ukUhsJ7aFDty8eXPs2LGDBg2ydyEDhlarraysHDNmjL0LGUjkcrmnp6evr2/HzXbt2hUWFtZxG17O3KdMmTJ//vyJEyeGhoaqVKqTJ09y29PT0+Pj46OiogYNGrRjx45XXnklIiLi1q1b6enptslOCHnppZeGDRvGsiwfRQIADFDDhg3rtA0vZ+4cuVxeU1MTGRlpvQbV2Njo5OTk5ubGrdbX19+7d08qlXK3zQAAQG/hMdwBAMBeMLYMAIAAIdwBAASIlwuqvctisdy4ccPV1XXcuHHttSkpKamtrZ04cWLLZ0yMRuONGzf8/PxCQ0P7pNJ+RKfTFRUVDR8+vL17TB8+fCiXy4cPH97yxu36+nprN52zs3N7Dx8I1YMHD+rr65988sk2b4lpObydi4uL9dIRpbSoqMhisURGRnb6XInAsCx748YNZ2fn8ePH237VYDBYn3fheHt7E0IYhjGbzdwWsVg8dOjQPii1/2hqampqavLy8mqvQU1NTUlJSVhYWMt7ZnpymNH+rbq6OjIyctKkSaNGjUpOTm5qarJts2LFihEjRiQkJAQFBd2+fZvbWFRUNHLkyISEhBEjRqxatapvq7azixcvBgQETJs2zd/ff8OGDbYNUlJSfHx8Jk+e7Ovr++KLL5rNZm67m5vbE088IZVKpVLpH/7wh76t2p5Yln355ZeDgoLi4uJCQ0OLi4tt24hEouDgYO7Dee+997iNOp3u6aefHjt27MSJE2NiYjQaTd8Wbk81NTVPPfXUU089NWbMmKSkJIPB0KpBamqq9DE/P7+hQ4eyLEspnTt3rr+/P7c9ISHBHrXbh0wmi4yMlEgk7u7u7bU5cOCAj4/P9OnTfXx8vv76a26jTqebPn16dw+z/h7ua9euXbBgAaXUaDROmjRp3759rRpcuHAhMDCwpqaGUrp+/fp58+Zx259//vn169dTSmtqagIDA3Nzc/u2cHuaPHny1q1bKaUKhcLT0/P+/futGpw6dcpkMlFKNRpNcHDwl19+yW13c3NTqVR9XG1/cPLkyeDg4IaGBkrpW2+9tWTJEts2IpGIO8xa2rp1a0JCgtlsZln22Wef/eCDD/qi3P5h/fr1KSkpLMuaTKbY2NiMjIwOGr/yyivWc6y5c+d+8cUXfVJj/1JeXn7+/PmcnJz2wl2v13t7e+fk5FBKz5496+fnx/3K3LZtm/Uwe+655zZt2tSVt+vv4R4cHPz9999zy3/729+Sk5NbNVizZs3KlSu55QcPHjg5Oen1eoZhJBLJgwcPuO0rV6586623+qxm+5LL5RKJhOtgoZSmpKRs3ry5g/a/+tWv0tLSuGU3N7dr165VVFTwXmU/s3z58t/97nfc8vXr111dXZubm1u1EYlEd+7cUavVLTdOmTLFGmpHjhyJiIjog2r7iVGjRn377bfc8rZt22bOnNleS61W6+HhceXKFW517ty5O3fuVCgU1j8ZHUpBQUF74X7s2LGwsDDramho6HfffUcpnTp16q5du7iNmZmZEyZM6Mob9esuQpZly8vLrT3moaGhSqWyVRulUmltEBISwr2koqKCZdng4OAOXihUSqXS19fX2o/Z8fd++/bt/Pz8OXPmcKtisTglJeXJJ58MDQ3Nycnpg2r7iZZHUWhoaFNTU6txTAkhEokkOTk5IiJi1KhReXl5bb7QcQ4zSmlZWVkXv/evv/46JCQkJibGumXTpk3Tpk3z8fFpc7IHh6VUKqVSqXU1JCREoVAQm8OM29ipfn1B1WQyNTc3Ozs7c6uurq46na5VG4PBYL2IKpFInJycuDZOTk5OTk4dvFCo9Hq99RMjhLi6utbV1bXZ8uHDhwsWLNiwYYP1atj9+/cDAwMppenp6S+++KJcLrdeNhS2lkeRq6srIcT2gFGpVNyHs3nz5kWLFpWWlrq4uLR6oV6vp49HTBI2i8ViNBo7/tm02rdv32uvvWZd3b17Nzf0RX5+flJSUmxsbEJCAt8FDwitfngHDx6s1+uJzfHJXXrs9DDr12furq6uQ4cOtWZTbW3t8OHDW7UJCAiwNtBqtSaTafjw4QEBAWazubGx0fpCxxlHJTAwkLvewq22973X1dXNmjVr4cKFv//971u+lhAiEonefvtt7n6bvqnZ7loeRbW1tYQQ2yPN+uGsXbu2pqbm/v37ti8MDAx0hGQnhDg5Ofn6+nb8s8kpKiq6evXqr3/9a+sW6wEZHx+fmJh4/vx5vqsdKAIDA1ueinHXC4nNYRYQENCVw6xfhzshJDY2Njc3l1vOzc2NjY1t1SAmJqZlg5CQEH9//2HDhgUHB1v/dm7zhUI1evRoJyenq1evcqt5eXm233tDQ0NycnJiYuKHH37Y5k4ePnzY1NTUcsx9YWt1mI0fP76D20DVarXZbOY+nNjY2AsXLlhf6DiHGenCzyZn79698+bN8/f3t/0S17fjOIdZp2JiYq5du8YN5q7Vam/cuMH1ZXXxo26t59cF+sSJEyf8/PwyMzMzMjKGDBnC3eloMBj8/f2LiooopRqNxt/ff9OmTSdOnAgPD//000+5F37yySfh4eEnTpzYtGmTv7+/Q92jtm7dutjY2P/85z/vvPOOVCrlbow5e/as9VrN9OnTx4wZk/EYdytRTk7O+++//8033/zzn/+MioqaPXs2d+OaI1Cr1d7e3mlpacePHw8NDd2zZw+3ffbs2Tt37qSUnj59euPGjd9+++0//vGPiRMnzp8/n2tQWFg4ZMiQffv2HTx40Hqfg4M4deqUj4/P4cOH9+zZM2TIkBs3blBKzWazv7//tWvXuDZGo9Hf3//kyZPWVzU2Nq5YseLQoUNHjx596aWXRowYYXsPklAxDJORkfHnP//ZxcUlIyPj4MGD3HbrYUYpTU5OXrBgwalTp1544YU5c+ZwG2UyWcvD7Ny5c115Owk302m/NWbMmPDw8IMHD1ZUVGzdutU6q0tVVVViYqKnp6erq+vcuXNPnTp14cKFpUuXrl69mvuDJS4uztXV9dChQyzL7tmzx6HmWJgxY4Zer8/KyvLw8Pj73//OPTliMpksFsvMmTMJIWVlZWFhYZrH/Pz8wsPDCSF5eXnnzp2Ty+UpKSlpaWnWixaC5+Hh8dxzz33//fd5eXmrV69etmwZt722tjYiIiIkJIRSmpube/bsWZVKtXDhwo8++kgikRBChg8fHh8fn5WV9eOPP3744YcONfdAWFjY+PHjDx06pFKp0tPTraeTarV65syZ3DivKpVKJBK9+uqr1uduxGLxzZs3z58/L5PJxo0bt3fv3jZP6gXJYDAcP36cZdm4uDiNRmMymaZMmUJaHGaEkOeff/7HH3/87rvvIiIi0tPTua72nh1mGDgMAECA+nufOwAA9ADCHQBAgBDuAAAChHAHABAghDsAgAAh3AEABAjhDgAgQAh3AAABQrgDAAgQwh0AQIAQ7gAAAuQoI0MB9IxGozly5AghpL6+PioqqqSkRCaTZWRk2LsugE4g3AE6snv37nXr1hFCvL29jxw5IpVKd+/ebe+iADqHUSEBOlJYWBgdHU0IEYnwwwIDCfrcATrCJfuZM2ccaqx2EACEO0DnsrOzreFeWFho32IAugLhDtCuzMzMRYsWEUKysrKkUikh5MyZM15eXvauC6BzuKAK0C6pVCqVStPS0k6fPs1dR/Xy8goLC7N3XQCdwzUiAAABQrcMAIAAIdwBAAQI4Q4AIEAIdwAAAUK4AwAIEMIdAECAEO4AAAL0/wOn3oqvAbwujgAAAABJRU5ErkJggg==",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {
"filenames": {
"image/png": "/Users/driscoll/repos/fnc-julia/_build/jupyter_execute/localapprox/integration_8_0.png",
"image/svg+xml": "/Users/driscoll/repos/fnc-julia/_build/jupyter_execute/localapprox/integration_8_0.svg"
}
},
"output_type": "execute_result"
}
],
"source": [
"plot([exp,x->exp(sin(x))],0,1,fill=0,layout=(2,1),\n",
" xlabel=L\"x\",ylabel=[L\"e^x\" L\"e^{\\sin(x)}\"],ylim=[0,2.7])"
]
},
{
"cell_type": "markdown",
"id": "4b8079a6",
"metadata": {},
"source": [
"```{raw} html\n",
"
\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"```{index} ! numerical integration\n",
"```\n",
"\n",
"```{index} see: quadrature; numerical integration\n",
"```\n",
"\n",
"Numerical integration, which also goes by the older name *quadrature*, is performed by combining values of the integrand sampled at nodes. In this section we will assume equally spaced nodes using the definitions\n",
"\n",
"```{math}\n",
" :label: nc-nodes\n",
" t_i = a +i h, \\quad h=\\frac{b-a}{n}, \\qquad i=0,\\ldots,n.\n",
"```\n",
"\n",
"::::{proof:definition} Numerical integration formula\n",
"A **numerical integration** formula is a list of **weights** $w_0,\\ldots,w_n$ chosen so that for all $f$ in some class of functions,\n",
"\n",
"```{math}\n",
" :label: quadrature\n",
" \\begin{split}\n",
" \\int_a^b f(x)\\, dx \\approx h \\sum_{i=0}^n w_if(t_i) = h \\bigl[ w_0f(t_0)+w_1f(t_1)+\\cdots w_nf(t_n) \\bigr],\n",
" \\end{split}\n",
"```\n",
"\n",
"with the $t_i$ defined in {eq}`nc-nodes`. The weights are independent of $f$ and $h$.\n",
"::::\n",
"\n",
"Numerical integration formulas can be applied to sequences of data values even if no function is explicitly known to generate them. For our presentation and implementations, however, we assume that $f$ is known and can be evaluated anywhere.\n",
"\n",
"```{index} interpolation; by piecewise polynomials\n",
"```\n",
"\n",
"A straightforward way to derive integration formulas is to mimic the approach taken for finite differences: find an interpolant and operate exactly on it. \n",
"\n",
"## Trapezoid formula\n",
"\n",
"One of the most important integration formulas results from integration of the piecewise linear interpolant (see {numref}`section-localapprox-pwlin`). Using the cardinal basis form of the interpolant in {eq}`plbasis`, we have\n",
"\n",
"```{math}\n",
"\\int_a^b f(x) \\, dx \\approx \\int_a^b \\sum_{i=0}^n f(t_i) H_i(x)\\, dx = \\sum_{i=0}^n f(t_i) \\left[ \\int_a^b H_i(x)\\right]\\, dx.\n",
"```\n",
"\n",
"Thus we can identify the weights as $w_i = h^{-1} \\int_a^b H_i(x)\\, dx$. Using areas of triangles, it's trivial to derive that\n",
"\n",
"```{math}\n",
":label: hatintegral\n",
"w_i = \\begin{cases}\n",
"1, & i=1,\\ldots,n-1,\\\\\n",
"\\frac{1}{2}, & i=0,n.\n",
"\\end{cases}\n",
"```\n",
"\n",
"Putting everything together, the resulting formula is\n",
"\n",
"```{math}\n",
":label: trapezoid\n",
"\\begin{split}\n",
" \\int_a^b f(x)\\, dx \\approx T_f(n) &= h\\left[\n",
" \\frac{1}{2}f(t_0) + f(t_1) + f(t_2) + \\cdots + f(t_{n-1}) +\n",
" \\frac{1}{2}f(t_n) \\right].\n",
"\\end{split}\n",
"```\n",
"\n",
"```{index} ! trapezoid formula; for integration\n",
"```\n",
"\n",
"::::{proof:definition} Trapezoid formula\n",
"The **trapezoid formula** is a numerical integration formula in the form {eq}`quadrature`, with\n",
"\n",
"$$\n",
"w_i = \\begin{cases}\n",
" \\frac{1}{2},& i=0 \\text{ or } i=n, \\\\ \n",
" 1, & 0 < i < n.\n",
" \\end{cases}\n",
"$$\n",
"::::\n",
"\n",
"Geometrically, as illustrated in {numref}`fig-trapezoid`, the trapezoid formula sums of the areas of trapezoids approximating the region under the curve $y=f(x)$.[^comp]\n",
"\n",
"[^comp]: Some texts distinguish between a formula for a single subinterval $[t_{k-1},t_k]$ and a *composite* formula that adds them up over the whole interval to get {eq}`trapezoid`.\n",
"\n",
"The trapezoid formula is the Swiss Army knife of integration formulas. A short implementation is given as {numref}`Function {number} `.\n",
"\n",
"\n",
"```{figure} figures/trapezoid.svg\n",
":name: fig-trapezoid\n",
"Trapezoid formula for integration. The piecewise linear interpolant defines trapezoids that approximate the region under the curve.\n",
"```\n",
"\n",
"(function-trapezoid)=\n",
"````{proof:function} trapezoid\n",
"**Trapezoid formula for numerical integration**\n",
"\n",
"```{code-block} julia1\n",
":lineno-start: 1\n",
"\"\"\"\n",
" trapezoid(f,a,b,n)\n",
"\n",
"Apply the trapezoid integration formula for integrand `f` over\n",
"interval [`a`,`b`], broken up into `n` equal pieces. Returns\n",
"the estimate, a vector of nodes, and a vector of integrand values at the\n",
"nodes.\n",
"\"\"\"\n",
"function trapezoid(f,a,b,n)\n",
" h = (b-a)/n\n",
" t = range(a,b,length=n+1)\n",
" y = f.(t)\n",
" T = h * ( sum(y[2:n]) + 0.5*(y[1] + y[n+1]) )\n",
" return T,t,y\n",
"end\n",
"```\n",
"````\n",
"\n",
"Like finite-difference formulas, numerical integration formulas have a truncation error.\n",
"\n",
"```{index} ! truncation error; of a numerical integration formula, ! order of accuracy; of numerical integration\n",
"```\n",
"\n",
"::::{proof:definition} Truncation error of a numerical integration formula\n",
"For the numerical integration formula {eq}`quadrature`, the **truncation error** is\n",
"\n",
"```{math}\n",
":label: truncint\n",
"\\tau_f(h) = \\int_a^b f(x) \\, dx - h \\sum_{i=0}^{n} w_i f(t_i).\n",
"```\n",
"\n",
"The **order of accuracy** is as defined in {numref}`Definition {number} `.\n",
"::::\n",
"\n",
"In {numref}`Theorem %s ` we stated that the pointwise error in a piecewise linear interpolant with equal node spacing $h$ is bounded by $O(h^2)$ as $h\\rightarrow 0$. Using $I$ to stand for the exact integral of $f$ and $p$ to stand for the piecewise linear interpolant, we obtain\n",
"\n",
"```{math}\n",
"\\begin{split}\n",
" I - T_f(n) = I - \\int_a^b p(x)\\, dx &= \\int_a^b \\bigl[f(x)-p(x)\\bigr] \\, dx \\\\\n",
" &\\le (b-a) \\max_{x\\in[a,b]} |f(x)-p(x)| = O(h^2).\n",
"\\end{split}\n",
"```\n",
"\n",
"```{index} ! Euler–Maclaurin formula\n",
"```\n",
"\n",
"A more thorough statement of the truncation error is known as the **Euler–Maclaurin formula**,\n",
"\n",
"```{math}\n",
":label: eulermaclaurin\n",
"\\int_a^b f(x)\\, dx &= T_f(n) - \\frac{h^2}{12} \\left[ f'(b)-f'(a) \\right] + \\frac{h^4}{740} \\left[ f'''(b)-f'''(a) \\right] + O(h^6) \\\\\n",
" &= T_f(n) - \\sum_{k=1}^\\infty \\frac{B_{2k}h^{2k}}{(2k)!} \\left[ f^{(2k-1)}(b)-f^{(2k-1)}(a) \\right],\n",
"```\n",
"\n",
"where the $B_{2k}$ are constants known as *Bernoulli numbers*. Unless we happen to be fortunate enough to have a function with $f'(b)=f'(a)$, we should expect truncation error at second order and no better.\n",
"\n",
"::::{proof:observation}\n",
"The trapezoid integration formula is second-order accurate.\n",
"::::\n",
"\n",
"(demo-int-trap)=\n",
"```{proof:demo}\n",
"```\n",
"\n",
"```{raw} html\n",
"
\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"We will approximate the integral of the function $f(x)=e^{\\sin 7x}$ over the interval $[0,2]$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "56a6354d",
"metadata": {},
"outputs": [],
"source": [
"f = x -> exp(sin(7*x));\n",
"a = 0; b = 2;"
]
},
{
"cell_type": "markdown",
"id": "d18567c0",
"metadata": {},
"source": [
"::::{panels}\n",
":column: col-7 left-side\n",
":card: border-0 shadow-none\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.5\\textwidth}\n",
"```\n",
"In lieu of the exact value, we use the `QuadGK` package to find an accurate result.\n",
"\n",
"```{raw} latex\n",
"\\end{minipage}\\hfill\n",
"```\n",
"---\n",
":column: col-5 right-side\n",
":card: shadow-none comment\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.4\\textwidth}\\begin{mdframed}[default]\\small\n",
"```\n",
"If a function has multiple return values, you can use an underscore `_` to indicate a return value you want to ignore.\n",
"```{raw} latex\n",
"\\end{mdframed}\\end{minipage}\n",
"```\n",
"::::"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "1115641e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Integral = 2.6632197827615394\n"
]
}
],
"source": [
"Q,_ = quadgk(f,a,b,atol=1e-14,rtol=1e-14);\n",
"println(\"Integral = $Q\")"
]
},
{
"cell_type": "markdown",
"id": "ff238def",
"metadata": {},
"source": [
"Here is the trapezoid result at $n=40$, and its error."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "e9f7f718",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(T, Q - T) = (2.662302935602287, 0.0009168471592522209)\n"
]
}
],
"source": [
"T,t,y = FNC.trapezoid(f,a,b,40)\n",
"@show (T,Q-T);"
]
},
{
"cell_type": "markdown",
"id": "555f873a",
"metadata": {},
"source": [
"In order to check the order of accuracy, we increase $n$ by orders of magnitude and observe how the error decreases."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "b8a9ae1a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"┌────────┬─────────────┐\n",
"│\u001b[1m n \u001b[0m│\u001b[1m error \u001b[0m│\n",
"├────────┼─────────────┤\n",
"│ 10 │ 0.0120254 │\n",
"│ 100 │ 0.000147305 │\n",
"│ 1000 │ 1.47415e-6 │\n",
"│ 10000 │ 1.47416e-8 │\n",
"│ 100000 │ 1.47417e-10 │\n",
"└────────┴─────────────┘\n"
]
}
],
"source": [
"n = [ 10^n for n in 1:5 ]\n",
"err = []\n",
"for n in n\n",
" T,t,y = FNC.trapezoid(f,a,b,n)\n",
" push!(err,Q-T)\n",
"end\n",
"\n",
"pretty_table([n err],[\"n\",\"error\"])"
]
},
{
"cell_type": "markdown",
"id": "64585d10",
"metadata": {},
"source": [
"Each increase by a factor of 10 in $n$ cuts the error by a factor of about 100, which is consistent with second-order convergence. Another check is that a log-log graph should give a line of slope $-2$ as $n\\to\\infty$."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "ebab3311",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAFACAIAAADF7xm+AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ0BT198H8N+9NyHsMARBlCGKOAFFcICogAP3AK1WbR3YYR31sWq1jrbOtu7aioq1jio4ClqcOHEgqKAyRBkOENkbktzxvIh/SjEg6CWB8Pu8Cje555wMvhzOPTmH4DgOEEIIqRdS1Q1ACCHEPwx3hBBSQxjuTUlCQsKqVaumTZs2fPjwGTNm7N27t6SkRNWNahry8/MXLFjg4+Pz4MGDt++Ni4u7c+eO8ltVR+PGjTty5IjCu5KTk729vZOTk2s5/aOPPvrzzz/rVWwlZb4yVeuqS9tQ7QSqbgCqq+PHjwcEBHTs2HHQoEG6urrJycnHjx+PjIzcunWrlpaWqlvX2IWFhT19+nTFihW2trZv33vu3LkXL164uLgov2F1MXLkSHt7e5UUq8xXpmpdDfSUmxUM96YhNTV19+7d48aNmzVrFkEQ8oPDhg2bO3ducHDw1KlTVdu8xq+0tNTMzKxXr16qbsj7mDZtWhMqlheNuW1NBYZ70xAQEGBiYjJ9+vTKZAeADh06+Pj4ZGZmVh4JDQ0NCwtLT09v3br14MGDR48eDQAbNmzo27dvdnb2+fPnX758aWdn99VXX+3fv//JkycHDx6sPHf58uXp6en79u0DgJs3bx49ejQlJcXAwMDb23vy5MkURa1bt87DwyM9Pf348eMLFy7s2bMny7IHDx68dOlSSUmJk5PT4MGDly5devz4cX19/ZoKUdgYa2treRuOHz9+7ty5169ft23bdvLkyc7OzvLjCotS+EIpfAWWL18eGRkJAN7e3hs2bOjevXvVU+bOnZuQkCC/9/fffw8KCqr2NKVS6eHDhyMiIl69eqWnp+fq6urv76+jo1NYWPj999/Pnz//2LFjkZGRAoHAxcXF399fU1OzlmZnZGS8nVyrV6/u06dPTe338/MbO3bsxIkTAUAqle7bt+/GjRtlZWUODg7Dhg2rLKSmdtb0oaostqY3pdorI/+np6b3oqYPw6+//lr19XRwcFDYyGp1LV26tPIp1+uDXflZQtSqVatU3Qb0DhzHbdu2bfDgwZVhV6lXr159+/aV396/f/+ePXs8PDxGjRpFUdSBAwdomnZycgoMDIyLi8vMzJw4cWK7du0iIiKioqJGjBhx5syZPn36GBkZAUBpaenWrVvHjh3btWvXCxcurF27tmfPnuPGjTMwMAgKCkpPT3dzc9uzZ09qampUVFT//v2dnZ319fW3bt0aEhIyYsQILy+vJ0+eHD16lKbpCRMmiESimgpR2JiRI0cCQEBAwOHDh318fDw9PdPT0w8ePGhvb29hYVFTUW+/UDW9Ak5OTqWlpWVlZb/99puVlVW1Pwyurq7Z2dkikWjr1q3m5uaBgYHVnubmzZtDQ0MHDx48fPhwY2Pjv//+u6SkxNXVtaSkZOfOnXfv3iVJcuzYsaampiEhIY8ePRo0aBAA1NRskUjUq1evof+Tl5eXl5c3ZcoUPT29mtofHBzcsWPHLl26AMDq1avDw8N9fHy8vb2fPXsWFBQklUqHDx9uZGRUUzuPHz/evn17BweHai9XZbE1vSnVXhmSJGt5L2r6MBw4cKDq67l3716FjaxW1/Hjxyufcr0+2PLPEgLsuTcJWVlZ5eXl5ubmtTymsLAwKCho4sSJM2bMAID+/ftramoeO3ZszJgxAJCXl7d161aRSAQApaWlBw4c6N69u0gkunHjRrt27QDg5s2bDMN4enrKZLK9e/eOHDlyzpw5AODh4WFubv7LL7/I+1Cpqam7du1q2bIlAKSnp585c2bu3LnyzmP//v3nzJmTlJQEALUX8nZjpFJpQUHByZMn/f395Q0eMGDAtGnTzpw54+joWFNR1fpotbwChoaGOjo6QqFQ3vJqjIyMtLW1NTQ0Ku+t+jQBoKKiYsqUKZMnT5YXm56eLn+acrq6uj/++CNJkgBgZWW1Zs2a+/fvd+nSpZZmd+zYUX7utWvXIiMjV6xYYW5uXvs7KBcfH3/r1q158+YNHz4cADw9PZctW1Z5EbL2dtZO4ZtS7ZWp5W0VCoU1fRiqvZ41NfLtd+Gdb2tNzdbQ0Kjjs1ZvOFumCaBpGgCqDsi8LTk5WSqVenl5VR7x8vKSSqXyeRQ9evSQ/wIAgI2NDQAIhcJevXrduHFDfvDq1auOjo6mpqYvX77Mzc3t3r17zv/Ir2vJ/2V2dnau/N179OgRy7Kenp7yHwmC8PDwkN+uvZC3G8MwTHx8PE3TlaUJBIKAgID58+fXXlTdX4F6qfo0AWD58uWTJ0+WSqUvXryIiIiIj4+v+tU/b29vebIDgIeHh7a2dnx8fF2anZaW9tNPP02YMMHd3b2O7U9ISCBJcsiQIZVH5P8l1KWdtVP4plR7TC1PqpYPQ7XXs76NrO8H++1mN1vYc28CzM3NRSJRRkbG23fFxsZGRkZOmjQpJycHAFq0aFF5l/x2dnY2ABgYGLx9roeHx/fff5+RkSEWi+/evfv1118DgHwEf+XKldUeXFxcDABmZmaVR7KysnR0dCrHlwHA2NhYfqP2QhQ25vXr10KhUD5YLycfLI6Li6ulqKpqfwXqperTBIDHjx//+uuviYmJ+vr6bdq00dXVrXpv1RoJgmjRokVWVlbtrwAAlJaWrlq1yt7e/tNPP617+3Nzcw0NDQWCf39tTUxM6tjO2il8U6qp5UlJJJKaPgzw39ezvo18jw82ksNwbwJIkrSzs7t69eq0adOEQmHVu0JCQh4+fOjv7y//xOfm5lZeQMvLy4P//Y4p7PW7urpqampGRETI80I+cir/bdm3b1/r1q2rPf6ff/6pOlptbGxcWlpa9b/g/Px8+Y1aCjl79qzCxhgZGclksrKyMm1tbfmRjIyM4uLiWoqqpvZXoF6qPs2ysrIFCxZ4eHh899138iTdvn17YmJi5QNyc3OrnpuXl2dkZFR7szmOW79+vUQiWbZsWWVddWm/iYlJQUEBTdOV+V5Z+zvbWbva/y+Uq+VJhYWF1fRhgCqv53s08j0+2EgOh2Wahs8++ywrK2v37t0sy1YejIqKunnzZr9+/QDAxsZGKBSGh4dX3nvx4kWhUKhwWrechoZGnz59bty4cfXq1X79+skny1tZWWlqal65cqXyYZcvX542bVpWVla10+UDx5cuXao8cu3aNfmNuhdSyd7eniCIylNYlv32228PHTpU96Le4xWoiydPnshkstGjR1f2kdPS0qo+4OLFi5UDC9evXy8pKenYsWPtzT5w4EB0dPTKlSurdjzr0v5OnToxDHPhwoXKI5W339nOD1fLk6rlw1DVezSygd7W5gB77k2DnZ3dp59+unfv3vj4+F69eonF4sTExEuXLrVs2VJ+ocnQ0HD8+PF//fWXTCbr1KlTfHz8sWPHJkyYIJ8MUxMPD4+VK1cKBIJ169bJj2hra0+ZMmX37t0FBQUODg7JycknT57s3r27qalptXNtbGy8vb23b9+en5/funXrS5cuFRYWAgBJknUvpFKbNm0GDx68Y8eO/Px8c3PzK1eupKenz5kzp+5Fvd8rICcUCl+9ehUbG2tnZ1ftrtatWwuFwr17944dO5ZhmNOnT6ekpAgEgpSUFPkgUkpKyooVKwYOHJiRkfHXX3916dJF/jWcmpodHR194MCBESNGkCRZecnR0NDQxMTkne3v0KFD7969t23blpmZaWNjc+PGjdu3b7+zne98+nV8ZWp/L2r6MNTxxWzbtq3Cd+FD3tZmDsO9yZg4cWLHjh2PHDkSFhZWXFzcqlWr8ePHT5o0qfLrqdOnTzc0NDx79uypU6csLCw+//xz+XTgWjg7O+vo6Ojo6FSdJ+fn5ycWi0NCQs6dO2dgYDB69OhJkyYpPP3rr782NjY+ffo0x3EDBgxwcXHZtGmTfOC17oVUmj9/vpmZ2aVLl7Kyslq3br1y5Ur51M+6F/Uer4DcoEGD7t+/v2zZsh07dlS7y9jYeMWKFYGBgT/++GPr1q1Hjhw5Y8aMlStXbtmyZcWKFQCwcOHC27dv79y5UygUent7+/v71/4yxsfHA8CpU6dOnTpVWYuvr6+/v39d2r98+fLAwMDLly8fPnxYKBR+/vnnO3furL2ddXkF6vLKWFtb1/Je1PJhqMuLuW3btprehfd+W5s5Apf8Re+nqKjo2rVrffv2NTQ0lB/ZvXt3RETE/v37VdswpcnJyfnoo4/k876VX3tpaSlJko1k5Qn8MDRCOOaO3pOmpmZgYOC2bdseP35cWlp6/fr1f/75B79CojQ6OjqNJNkBPwyNEvbc0fuLj4/ftGnTs2fPAEBLS2vMmDFTpkypOlFPvam2597YNPMPQyOE4Y4+VEFBgVQqNTExaW7z0mQy2ePHj62tres1o1y9NdsPQyOE4Y4QQmoIx9wRQkgNYbgjhJAawnBHCCE1hOGOEEJqCMMdIYTUEIY7QgipITUJd5lM9uDBA1W3Qn3gjgeIXxzH4axrJVOTcM/NzR08eLCqW6E+ysrKVN0EpFZompZKpapuRfOiJuGOEEKoKgx3hBBSQ80o3B88eLDt9z0FBQWqbghCCDW4ZhHu2dnZPYeM9/r+yNcJLRz95n65dLWqW4QQQg2rWSzIuXbrb/d6LmStnAHgWddhf4fMX5SWZm1trep2IYRQQ2kWPfd/rt+VJ7tcRvfpG3/bp8L2IIRQQ2sWPXcdEQW0BASiNz/npN7TsUwr5qz1cMlphBqL8PDw5ORkVbdClRwdHeW7q/NCTdZzz8zMdHJyevXqlcJ79x068sWZVxUeXwJBQHkRsW86N+uAUKT1qR25sjvZShsjvrri4mI9PT1VtwKpD5lMxrKsSCSq5TE9e/Zs06aNiYmJ0lrVqKSkpBgZGR09epSvAptFz33KhPEZmb8dPOxbomFoo01/tX7RRW2dvY/ZgET24FN2ZgdymSNl2lh2o0So+fr222+dnZ3f/Th1FBQUdPz4cR4LbBbhLhAIli386tuv5+Tl5RkbGwOAL8D8LuTKu+yxVHZbHBuYxH7ZifzWkdIXqrqtCCHEh2ZxQVWOIAh5sst1NCCCPKnYsQJfG7JEBhtiWdujsg2xbAWuqoIQavqaUbgr1NWICPKkbowQeJgTORWwJIppH0QHJLI0q+qWIYTQB2ju4S7XpyVxZZjgwlBB9xbEy1JudgTTPpgOSGQZdbjYjJCaY1nsiymA4f4vLwsierQgyJOyExNpxdzsCMbhBB2cip8bhBojhmFCQkK8Bw/R0taZPPnjmzdvqrpFjQuG+38QAL42ZMJ4QZAnZaNHxOVzfuFM71D6Ugb24RFqLHJyctauXWtp3Xb06NFXnpdLhy4Juna/b9++nbs5BgQElJeXq7qBjQKGuwIkAb42ZKKvYJcbZaYFt7M4zzDa+wwdnYMRj5CKcRzn5jFg5ZoNGe2HweoYeuFFGPYtvTIGFl1M1LL7/Ms5M2bOUn6r4uLili9frvx6a4HhXiMNEvztyacThOt7UgYacDGd6/k37X2Gjs3DiEdIZcLDwx/HP6IXnIVJW6FVp3/vsOvHzjrEfrQ1KCiopu8zNpzMzMwzZ87Ib5uYmOTk5Ci5AW/DcH8HHQEsdiCTJwgXO5DaAriYznU/SfuFM8lFGPEIqcCOX3cKbXuCdQ3fdeo1idDU3bNnT12Kev78uVQqTUlJuXHjhvxIRkbGjRs3srOzKx8jk8mioqKuXr368uXLqmfJb+fl5eXl5VUtMyMjIz8/Py0tLSUlRT5AlJ2dfe3atcjISCVvcIbhXidGIljfk0qbKFzsQApJCE5lOx6jZ0cwGWUY8Qgpz4sXL06fOiXrN7vGR2ho066Tt+/8XSaTvbO0Pn36fPzxx5MmTdq6dSsAzJs3z8PDY+vWrc7OzvI/D7m5ud26dVu+fPmePXsGDx786NEjAHBzc0tMTJSX8Msvv6xfv75qmRs3bmRZduHChbNnz46Liztx4oSTk9OuXbvWrl07atSoD3jq9dYsvqHKFxNNWN+T+qIjuSaGxdULEGpQYS+4h28NgV74bTenpQ89fWs708M/+9KvM3493XlA9TD1siB6tPjPWlIWFhZBQUEAcOzYsZs3b8bFxWloaKSnpzs5OY0dO/by5ctWVlZnz56VP7gucy63bNmyY8eOY8eOyRfJWb169c8//zxx4sQ6ns4jDPd6s9QldrlRuHoBQg2qjObypdXDvbRcApQQCKq2MwUiAKKwrOLt0yVM9VUCx48fL79x/vx5sVi8efNm+Y9CofDhw4cODg6RkZHTpk0bNWqUt7f3e6ym16tXryVLliQkJPj4+PC44mNdNK9wl6TEASMTtXf88KLkqxdEZZPf32dOP+c2xLJ7H7P/15Wa14XUrPWDhxCqi/E25Hib6gdTVs5ud/BnuHcSXCbUeObVAEMj46CFY0Wid/8q6uvry2+Ulpbq6+sbGhrKf1y5cqW1tbWVlVVsbGxwcPDmzZs/++yz8PDwrl27AkDlYro0Tdde/rJly/r06RMaGjp+/PguXbr8888/JKmkwfDmNObOsgXB27N/XZK9c6ksPYWXInuaEKcGCSKqrF5gh6sXINRg2rZt6z1osODarhofQUuEt//84jP/2pcXfpuzs3NxcbF/FVZWVgBgaWm5cOHC69evDx48ODQ0FABMTU0rL65GR0e/XZSWlpZEIqn8ccCAAZs3b05KSrp48aIyp/E0o547x7HaPQYUXwySJN1//cscHRdv/aFTKbHxu898l74tiSvDBBfTucVRzL0cbnYEsy6WXepAzuhAUrhWPEK8mvPF5+dHjoT0R2DRRcHdUceYkrxZs+o91f3zzz8/evToiBEjBg8eXFxcfPbs2dOnTx87diwiIqJ79+7l5eWXLl368ssvAcDPz2/hwoWPHz+Ojo7OzMx8uyg3N7dZs2Y5ODhMnz591apVtra2lpaW0dHRnTt3Njc3r/8zfk/NYrOOqtiy4uLw4JKrJzlaRmiIdN1H6XlPIDV1eGkGB3AslV0ezSYVcgDQ2ZBY2Z30tWl6/x7hZh2IX3XcrOO3335753ruDMPYtLPLLKqQuc0A95lg8L+4zEiAK78LIg8N9R4Y+vfJurTq9OnT7u7uYrG4suSQkJC4uDg9PT03N7cePXrk5eWdPXv26dOnenp6w4cPt7OzAwCO406cOBEXF9e/f38jIyOWZbt165aZmXn//v2hQ4cCgEQiuXz5cnFxsbu7e1ZW1rVr13JycmxsbMaPH6+jU2PUyNdz53GzjmYX7nJMflZh2P6y6EvAcaSOvt7A8boeYwgBP9dDaRYOJ7Or7rGpxRwA9DIl1jhTA1s1pT48hjviF4/hDgDPnj3bsmXL3n1/lJaWEo4jGHtPQXQQnXjVvLXll5/5f/nllwYGBvy1XUl4D/em16nkBWVoajR5kemCraJ23djSosJTga83fFYecx34+FMnIGFqe1y9AKGGYmVltXnz5syM9IDff+tEpxGH5ni0EZ08efJFWsqyZcuaYrI3hGY05v42DUs7kzkbJUn3C07ukr1Ky/1jjYZlB/GomSLbrjwUToK/PTm5Hbkjjl0fy1xM5y6m014WxC+uVDejptSLR6hx0tbWnjFjxowZM3Jyclq0aKHq5jQ6zbTnXpXIzqnlN78ZTf4/St9I+vxx9vZF2TuXyl6l8VL426sXOOHqBQjxqkGT/T3WDEhLSzt37lx+fv57l8ALDHcAACAI7Z5eZsv26g/5mBBpSZLul946w2PxNa1e8Eo1bzpC6D9iYmKCgoJCQ0MPHz7MMEzV4/VdJj4nJ+fWrVv9+vWbOnWqvKiIiAj5ugVKhuH+L0KkpT/kY7NlgbruI/UHTeK9fPnqBUm+An97kuUgIJFtFySbd4vJwtWnEVIRmqYXL1787NkzPz+/kSNH2tvbT5kyRX6XVCoNCwvz8vKqb4EURWlpaTEMU1RUBACDBg0KDQ2ty1o3/MJwr47SNzQY9wWpK26g8uWrFzwcJ/C1Ictp2BbHtguSLYliipT91iOE4PPPP+/SpUvlkl7du3d/9uzZ3bt3AeCvv/4aMmRI7aeXlZXlV1FcXGxmZubn5/fkyRMnJ6fK77sOHDjw+PHjDfpE3tasL6iqEK5egFDtKuKjZK9S3z4uaNFKy8Gt8ke2pLA08pzCEnR6+5DaupU/lt27IjCx0GjTvvLIlStXIiMjAwICqp4lFosfP37co0ePy5cvT506tfL4tWvXhEJhu3bt4uPjy8vL5bn/9OnTrKysysdoaWn17du3sLDwzp07y5cvl0gk8tmfPXv2DAgIkC8fpjQY7qokX73gxmvu2yjmWia3JIr5NZ5d7kROtyMF+D8Vat44WQVbVqLguOQ/45gcyyp8mPy+6ifS//kH+bfffvP19SWI/8xeS0hIkC88IJFIKu+Kj483MDCYNWvWpk2bPDw8PvnkEzs7u7Zt23br1q1anRKJZNWqVR07dvz22283btwoP0hRVEVFxbueMc8w3Out6Nxh2atU8fDpghb8fJO4b0vi6nDBxXTumzvM/VxudgSzPpZd4kDO7ECSOGcSNVdaDu5aDu7vfBilbygeMb0uBer0HlrtSExMzLRp06odKSsrc3Jy4jiuahyLRKJWrVppamr27dsXAF6+fKmrqwuK0DT98ccfA0CvXr2Ewn+/F4nh3thxtKzkeihbUlDx8JaO+0h974mkjj4vJXtZEHfHCCpXL5gdwWyPY1c0zdULEGoSxGKxhoZG1SNr165dtWqVtrY2AFRdvtHW1vbUqVP9+/cHgPz8/PLyclNTU4Vl6ujo9OjR4+3j1f4/UAIMjvohBMKWi37V6ePDcWzJlROvvp9WdPYgJ5PyUziArw0ZN06w34Oy0SMe5XN+4UzvUPryK5wUjxD/xo0bFxERUfnjzp07LSwsPv/8c/mPVfvdAHDlyhUPDw8ACAsLGz16dEhISL0mwFQrTQkw3OuNEhsb+s1tufh3LUd3TlJedPZg5prppTfDgKdtVt5evWDgP7h6AUL8W7hwIcuyAQEB58+f37lzZ6tWrSo36wAAMzOzqt8/ys7O7t27NwBYWFgIhUJDQ8O653V+fr58HF+ZmunCYXyRJMUUhO6RvXwKAEIzK/GI6ZqdXXksv5QG+eoFBVIgAIZZEmuclbF6AS4chvjF78JhvJNKpdXGZwAgMjIyNTWVlykugYGBTk5OTk5OtTwGFw5rXER2ji0Xbjf+ZJnA2FyW+Sxn98rsnUtlL5P5Kr/q6gWaAjj9HFcvQIh/byc7ALi6ur548aK8/EO/ZFhSUpKXl1d7sjcEDPcPRhBaju4tv91t6DeX1BXLdwLJ/WMNnatgFf/3I1+94BmuXoCQcs2ZM+fatWsfWMj169fnzp3LS3vqBcOdHwQl0OnjY/btHt3+YwlKUB5z/fV6/+IrJ3isQr56weP/rl6wJIrJl7z7XITQe9DS0ho8ePAHFjJ06FCF/xk0NAx3PpHaegaj/c2WB+r08eFoGcXTLMmqrP67esGGWNbqCK5egBCqDsOdf5SBiaHf3Jb/t0Pb2bOBqpCvXnB7lGC4JVEsgw2xrO1R2YZYtoJ597kIoeYAw72hCC1soYG/tuBiQpwaJIgYIehnRuRUwJIoxi6IDkhkaX7mZCKEmjD8hmqT9/bqBZsest86kh+3w9ULUFOiqam5Y8cOMzMzVTdENRITE2ufKlpfOM9dNZj8LKCElL4hj2VyAMdS2WXR7JNCDgC6GBLvvXoBznNH/KrLPPfk5ORjx44prUmNUK9eveRfguUFhrtq5OxeJXkSo+s+Us97IqmpzWPJMhb+SmZX3WNTizkA6GVKrO1JDTCvXx8ewx3xqy7hjviFY+4qwNEygiI5aUVxeNDrtTNLb4YBy9uVUGGV1Qta4uoFCDVX2HNXGemzxMKQPZKURwAgMG0t9pmm5fjuBU7r5b1XL8CeO+IX9tyVD8NdxSribhf8vZvOTgcADWt78chZorad+a0iTwIbHzDb4thyGkgCxlmT63qStvq1RTyGO+IXhrvyYbirHsfQZZHni84cYIrzAUCzs6vBmNmCFq34rSW7An55yGx5xEoYEJLwqR25qjtlXsNoP4Y74heGu/JhuDcWnLSi5Hpo0fm/OEk5QQm0XQeJfaaSugb81vKshFsbw+59zDIcaAvgq87k4m6U4Vu/cRjuiF8Y7sqH4d64MIW5RecOld4+CyxLiLT0BozT8/QjhDwvTJFQwK28yx5LZTkAPSF80Yn81pHSr7I2NYY74heGu/LhbJnG5c1OIIt2anbsKd8JJD9oG++11L56QUZGxrZde5OTeVu4GCGkfNhzb7wkSTGFpwMNJy4QtrJpuFpuvOa+jWKuZXIA0EaHMDw2J7ecy+gwsmXKBTsy+/yRQOxtoQ+HPXflw3BHAAChz9jld9mHt6/Cs3vgPV9+UPNe8PZeMHPqJNW2DakBDHflw2EZBAAw0oqMGSPoknAQ+kytPFjRbdRvh/hckh4hpDQY7ugNkgD3jm2I3Gf/HirNe0kYxeWrw/92CDU3GO5NWPnDmyVXTnAMzVeBS7+aZRqxDWQVAACMjDj3S5bzp91O0H7hzLMSjHiEmhIcc2+qOIZ+vXYmnZspMDbXHzZN28mDl+XjT4T+s/73/elSTROi5NOPxr/qOkH+vSf5pPilDpRYBfuFoSYPx9yVD8O9CZMk3S8I2SNLTwYAjTZ24pEzRO0deCn5+fPnlpaW8ttPCrll0W8mxRuLYFE3an4XUkTxUg9qLjDclQ/DvYnjuPLYiMLQPXTeawAQ2TkZjPb/8KmTb3+J6U42tyjyzYxJOzHxozM53gb3AkF1heGufBju6oCTSkquhxRfOMpWlAJBaDsPFA+fTomN37vAmr6hejGdW3CbeZTPAYCLCfGTK9XPDBMevRuGu/JhuKsPtrSo+NKxkqsnOVpGaIh03UfpeU8gNXXeo6halh+gWQhMYlfeZTLLAQC8LKI/0moAACAASURBVIjNvaguhhjxqDYY7sqH4a5umLyswjP7y6IvAceROvp6A8fr9h9LUPXbLPeda8vIV4pfE8MUy96sMbm6B2Wm9WFNR+oLw135MNzVkzQtoSBkjzQ1DgAEpq0NRs3U7Nyr7qfXceGw7Ar48T6zM4GlWdARwJzO5DJHSk/4zvNQs4Phrnw4z109aVh3NJ33i8kX64Tm1nTWS+nzpIaoxUQTtvamHo4V+NqQpTRsiGXtg+mARJZmG6I2hFA9YM9d3bFM6e1z2j0GEKJ6DJq8x5K/t7K4RZHMjdccANgbEN/3IH1tsOuA3sCeu/Lhr5+6IymdPj71Svb309uUuD5CEORJ2eoTiQWcXzjTJ5S++Vodug4INUUY7og3BICvDZkwXrDLjTLVgltZnNsp2i+cSS7CiEdI2TDcmztOJgWWzzFyIQn+9mSyn3Bld1JTAMGpbMdj9OwIJqucx0oQQu+A4d7cFf3zR+aGz8pjrvNbrK4QVnWnknwF/vYky0FAImt/7N/NnhBCDQ3DvVnjGLo8IZp+/Tz3jzXZO5fIXj7lt/zWOsQuNyp2rGC4JZEvgSVRTPsgOiCRZXCcBqEGhrNlmjuOocsizxeG/cmWFABBaDm4iUdML9fQ5X2D7Ivp3KI7TEwuBwCdDYkNLtSwNvi91uYCZ8soH4Y7AgBgy0uKLx4tuRbCyaSEQKjRa6ixz1RSW5fnWjg4nsYuvsOmFnMA4GVB/ORCORpjxKs/DHflw3BH/2IKcorOHy69fRZYltTW1fP00+03mhDyvIK7lIXf4tlV95gCKZAEjLMmN7qQ1noY8eoMw135MNxRdbL05NwTu+jkBwAgMGqpP2yadvcBvOwEUlWeBDY+YOQ7gWgJYG5ncokDZYA7gagpDHflw3BHChQXF2u8eloQukf2MhkAhGZW4pEzNDu58F7R8xJueTR78CnLARiJ4BvcCURNYbgrH4Y7UuDN8gMcVxZ1sTDsT6YgGwB03UcajPuiIaqLyua+ucNcecUBgJUu8X0Pckp73AlErWC4Kx9OhUQ1IwhtF2+z7/YZjPmM1NIVdejeQPX0NCEuDxNcGCroakQ8K+GmXWVcQ+irr9Sh24GQqjSNnvuLFy9u374tv+3s7GxjU30bOey58+vthcPYshLeJ8+8jeXg4FN2SRTzqgwAwMuC2ORKdTXCTnyThz135WsaPfebN2/u3LkzJSUlJSWluLhY1c1pjpSQ7ABAEjC1PfnUT7i+J6UvhIvpnONJeuqVN1mPEKq7+m3Qo0K9e/devHixqluBlEFbAIsdyBkdyB/uMzsT2ANP2RNp7JzO5LeOlD7uBIJQ3TSNnjsAHDlypHv37kOHDn327Jmq24Kqk2Wkynd94lELTdjam3o0TuBrQ5bRsCGWtT0q2/oIdwJBqE4a15j7yZMnt2/fXvXIiBEjFixYUFZWpqmpSZJkQEDA6dOnQ0NDq52IY+78qu9mHdm/LpE8iRHZORmMmS00t+a9PbezuEV3mIhMDgA6iIkfnHEnkCYGx9yVT6nhnpOTU1FR0bp166oHc3NzHz16ZGtr27p1a6lUWlFRUfVeDQ0NTU3Nyh8zMzPd3d2fPHlSrWQMd37VL9xZpuj8X8WXj3OSciAp3d5D9YZMpvQMeW/Vqefs17fZp0UcAPQyJX5yodzM8Fpr04DhrnxK6v78/fffVlZWpqamw4cPr3o8JCSkQ4cOa9eudXJy2rx5s4aGhv5/yZM9Ojo6Nze3rKxs06ZNffv2VU6bUV2RlP6Qj81X7Nfz9CMIouTG6cwfPi08FchW8HwZdIQlGT9esMuNaqkFt7M499P0iPP0U9wJBCFFlNRzl89yefjw4c8//xwTEyM/yDCMra3tpk2bxo4dm5CQ4OzsnJaWZmJi8vbp+/fv379/v1Qq7dOnz7Jly8RicbUHYM+dX++xh6ocnfWyMGx/eWwEcBypo68/aJKu+wggef7KaYkMfn7IbHzAltMgJOFTO/L7HlTLBt9JEL0/7Lkrn1KHZY4cObJ+/frKcL99+7aPj092djZFUQDQp0+fmTNnTp8+/T1KzszMtLW11dD4d2mSGzdutGnThpdmN0MlJSW6uu8/95F58aTi3AH6WQIAkCYWmgMnCLv05q91b2SUwYZ44Z8pFMOBjoCba8/M70BrNZn5X80Lhju/tLW15bFZC1X+Kjx//rxNmzaVTbS2tn7x4sV7l6anp5eQkFD5o1gsJkm85vb+Pmg9907doVP3irjbBScD6Oz0sqObNCI7GoyaqWHTmb8GQgc9CGwJi5y4lXfZ4FR23SPBn6nCFU7kjA4khUPxjQyGu/KpMv4qKiqq9rVFIlFZ2fuP0hIEYVgFJrvKaXbu1XJpgKHfXErPUJqWkLXt/3L/WEPn8Dx01tGACPKkLvoIurcg0ku52RFMt+N0cCrOl0TNnSoT0MzMLDc3t/LHnJwcc3NzFbYH8Y6gBDp9fFou26PnNYEQCMtjrhed+bMhKvJsRUSPFgR5Um31iPgCzi+c8T5D38vBa62o+VJluDs6OmZkZMiHYmiavnXrlqurqwrbgxoIqakjHv6p2bJAHdfB+j5TG6gWAsDXhkzwFexyo0w04WI65/w37RfOyHd9Qqi5UdIF1efPn//111+xsbFXrlyZN2+ejY2Nn58fAMycOfPp06dff/310aNH09LSbty48X7l42wZfr33bJlGIl8CGx4wWx+xFQxokPBZR3J1D9wJRJVwzF35lNRzl8lk+fn5lpaWU6dOzc/PLykpkR//9ddfhw8ffujQIWtr69OnTyunMUjtGYpgfU8qyU/gb0/SHGyLY22PyjbEshWMqluGkLI0ruUH3hv23PnV1HvuVd3N4RZFMpdfcQBgqUsscyRndiBxKxAlw5678uGUEtRYlN4My/1jDZ2byW+xPVoQl4YJLgwVOBgRz0u42RGMawh9GXcCQeoOwx01DixbdPFIecz11+v9C0/vYytK+S3ey4K4N0YQ5ElZ6RLROdzAf2jvM/SDPIx4pLZwWAYpoJJhGaYgu+j8X6W3zgDHkdp6ep6+uh5jCAHPK7iX0bA9jl0XyxRKgSRgsi253oVspY3DNA0Lh2WUD8MdKaDCMXfpiyeFoXskT2IBgDI0FftM1Xb2BILn8M2VwE8PmC2PWAkD2gL4qjO51IES43SaBoPhrnwY7kgBlV9QlSTdLwjZLUtPAQANSzvxyJmidt14ryWpkFsezR5LZTkAYxEs6kYt6Epq4FBlA8BwVz4Md6SAysMdAIDjymMjCkJ3M3lZACCyczIY7S9sVX1v9A93J5tbFMlcy+QAwE5M/OhMjrfB2TQ8w3BXPgx3pECjCHcAAOCkkpLrIcUXjrIVpUAQ2s4DxSNmUPpGvFd0MZ1bcJt5lM8BgKspsdGF6oc7gfAHw135MNyRAo0n3OXYkoKiswdLb53lGJrQ0NQfOkVvwDjea5GxsC+JXXmXySwHAPCyILb0ojobYsTzAMNd+XB8ETUBpK6Bwfg5LZfs0nJ052QSaJgeiZAEf3vy6QTh+p6UnhAupnNOJ+nZEW+yHqGmBXvuSIHG1nOvSpqWILSwJYQNO7Ulo4xbfY8NTGJpFnQEMKczucyR0uN5WmYzgj135aux5x4VFZWenq7MpiBUFxrWHRs62QGglTaxy416OFbga0OW0rAhlrUPpgMSWUYd+kKoWagx3H/55RdcyQs1c/YGRJAndXOkoE9LIqOMmx3BdP3vTiBlZWXq8b8vUj81hvvYsWOvXLny+vVrZbYGoQ/ElhaxpUX8ltnblIgYIQjypGz1iYQCzi+c6XuKPnwjYdQnX3QaMaPjwNHf/7SloqKC30oR+kA1jrnfvHkzODg4NTXV2dnZzMys6masn376qbKaV1c45s6vxjzmXrv8I5vL7l/TGzBOz9OP99Eb+XSa7+4yWWUc/DoOPt4BBq2A4wR3Dn9nV7Bi0Xx+q1MnOOaufDX23K9evZqXlycWi588eXL9+vUrVSixeQjVB8cxxYWcpLzo7MHMNdNL71zgd16NfDpN4nihr+Qq0XEAGLQCACAI2nXy4bDLOD6DGhVBTXcsXboUAAoKCpKTk7Ozs42Njdu1a2doaKjEtiFUTwTRYtYqaWp8QegeaWp8/uFfisODxUOnaDm681iJoQhcuOQT5nZVd/7IYnTyi8uM9HV4rAihD1HbVMjAwMBjx47JZDL5j0KhcNy4cdOnTyf4XsXpw+GwDL+a7rDMGxxXHhtReDqQznkFACI7R4ORM4Wt2/FVfHFxcbfxc9J897z5ufAVnF5n/dn2Nc7kR7a4dIECOCyjfDX23P/+++8jR4589NFHXl5eLVq0yM3NDQ8PP3TokLGx8ejRo5XZRITqjSC0HN01u/YuizxfGPanJCnm9S9faTm4iUdMFxibf3jxenp6Ez26/RU0Pb3LRJ3cJ/rJl3QnrE0o5iZfZjY9ZDe6UANbYcIjFaux5/7ZZ585ODh8/vnnVQ/+/vvvMTExv//+u1LaVg/Yc+dXk++5V8GWlRSHB5Vc+5uTSQlKoO06SOwzjdQVf3jJhYWFv/9xyLFzh0GeAzkgjqexi++wqcUcAHhZED+7Ug5GGPFvYM9d+Wq8oJqent6pU6dqBzt37ozfbEJNC6mtKx4x3WxZoE4fH45jS2+GZa6dURwexMmkH1iyWCxePO+LwV6eBEGQBPjakIm+gi29KAMNuJjOdT9J+4UzacV4lRWpRo3h3qpVq/j4+GoH4+PjLSwsGrhJCPGPMmhh6De35dfbRXZObFlJ4anAvAMbeK9Fg4R5XcjkCcLFDqSQhOBUtvNxekkUU/ihf0cQqrcaw33o0KEnTpzYu3fvs2fPSktLnz9/vm/fvuPHjw8ZMkSZ7UOIR8LWtiZfrDP5Yp3QwlbXfUQD1WIkgvU9qSRfwZR2ZDkNG2JZ26OyDbGshHn3uQjxpbbZMvv37w8KCpJK3/Q6hELh+PHjp0+frqy21QOOufNLncbcFeM43rfuUygqm/vmDnPlFQcAVrrE9z3IKe2b43QaHHNXvhrDPSoqqlWrVrq6usnJyTk5OS1atGjbtq2BgYGS21dHGO78Uv9wV66L6dzCSOZBHgcALibET67NbicQDHflqzHcJ06cOGDAgNmzZyu5Qe8Hw51fGO68Yzk4+JRdEsW8KgMA8LIgNveiujSbnUAw3JUPFw5DqDpJ8sOSiNMcQ/NYJknA1PbkUz/h+p6UPu4EghoeLhyGFGjWPXeOy9o0V/riicDEQjzsEy0HN95H53Mq4If7zM6EZrQTCPbcla/GcF+3bl1iYqLCu/bv39+QTXofGO78atbhDiBJul9wcpfsVRoAaFh2EI+aKbLtynstjwu576LZY6ksB2CiCcscqS87kQI13fgSw1353nFBtanMasdw51czD3cAAI4riw4vPBXIFOUBgMjOyWDMbKG5Ne/13M7iFt1hIjI5ALA3IL7vQfraqGHAY7grH15QRQpguMtx0oqS66HFF46wFWXypQv0h0yh9PlfG/XUc/br2+zTIg4AepsSP7lSfVuq1bVWDHflwwuqCNWI0NDU8/Qz++4P3X6jOI4rvRmW+eOnhacCOQnPl0FHWJLx4wW73KiWWnAri3M/RfuFM8lFuHQBen94QRUpgD33t9FZLwvD9pfHXAcASmysP3iyTq/BQFLvPLFeSmTw80Nm4wO2nAYhCZ/akT/0oEy1+K1EBbDnrnx4QRUpgOFeE8nTh4Whe6TPHwOA0NxaPHKGZseevNeSXsp9f5/d+5hlONAVwsKu5OJulFaN63M3ARjuylfb8gMAkJaWlpCQUFxc7Ofnl5WVZWpqqrSW1QuGO78w3GtXEXe74GQAnZOh08fH0G9uA9WSUMCtvMsGp7IA0FqH+M6JnNGBpJrmUDyGu/LVGO4Mw2zcuPHSpUvyHy9cuDBt2jRbW9vFixc3wncIw51fGO7vxDF0acRpLSePhri+WlV4BvfNHeZeDgcAnQyIDS7UcMumF/AY7spX4wXVgwcPXr9+feHChXv2vNlLbNq0adHR0YcPH1ZW2xBqvAhKoOsxuqGTHQA8WxHRowVBnlRbPSK+gBtxnvY+Q9/PxWut6B1qDPcLFy74+fkNGTJER+fNnr8DBw4cM2ZMZV8eIaQcBICvDZngK9jlRplowsV0zvlv2i+cScWdQFDNagz3wsLCtm3bVjtoa2tbUFDQwE1CSB1wDA21XtCqLw0S/O3Jx77CxQ6kBgnBqax9MD3vFlOAO4EgRWoMdxsbm3v37lU7GBMTY2lp2cBNQkgdlFw58XrDZxXxd/gt1lC+E4ifwN+epDnYFoc7gSDFagz3iRMnhoWF/frrr0lJSQDw8uXLgwcPnjp1atSoUUpsHkJNE8eVRV+SZT7LCViRs+s7+TI1PGqjQ+xyo+6MEgwwJ/IksCSKsQum/3zC4jANqlTbVMgLFy7s3bs3NzdX/qOWltbHH3/s5+enrLbVA86W4RfOlvlw8uk0RWcPsuUlQBBaDm7ikTMFRi15r+hiOvd/kUxsHgcAPU2IjS5Uf/NGN50GZ8so3zvmuctkspcvX2ZmZhoYGFhZWWlrayutZfWC4c4vDHe+sGXFxeHBJdf+5mRSQqih22+0npcfqaXLcy0cHE9jF0Wyz0o4APCyIH5xpboZNaKIx3BXvneEe1OB4c4vDHd+MflZRReOlN46AxxHauvpefrqeowhBDyv4F5Gw/Y4dl0sUygFkoDJtuQGF8q8cfTHMNyVD8MdKYDh3hCkL5IKQ/dKnsQCAGVkKh46VdvZk/edQHIl8NMDZssjVsKAtgC+6kx+60jpq3onEAx35cNwRwpguDccSdL9gr8DZBmpAKBhaSceOUvUjv+dQJIKueX/2wmkhSYsV/VOIBjuyofhjhTAcG9YLFMaeb7ozAH5TiANt0DNnWxuUSRzLZMDgA5i4gdncrwNqZKReAx35VPDPV8QauxISqf3ULPl+8QjppOa2hqWdg1Uj4sJcXW44MJQQRdD4nEh5xfO9Aml5bs+IbWHPXekAPbclYYtKSC19YFs2G6WjIV9SezKu0xmOQDAcEticy+qnb7yOvHYc1c+7LkjpEqkrkFDJzsACEnwtyefThCu70npCeH0c67TMXp2BPOa5x2lUCOC4Y5Qc6EjgMUOZKKvwN+e5AACEtl2QbJV95hyWtUtQw0Awx2hxovOzZTv+sSjVtrELjfq4ViBrw1ZIoPV99j2wXRAIsuowwAt+heGO0KNV2Ho3qxN83J2r6RzMvgt2d6ACPKkbowQ9GlJpJdysyOYrsdp+a5PSD1guCPUWHGc0LQ1oSGqiIt8vc6/4OTvbGkRvzX0aUlEjBAEeVK2+kRCAecXzniF0fJdn1BTh7NlkAI4W6bxYApzi84dKr19FliWEGnpDRin5+lHCDX4rUU+nea7u0xWORAA423I9S5kWz3eptPgbBnlw3BHCmC4Nzay18+Lzhwoj7kOAJRBC/1Bk3R6DeF9mk2+BDY8YLY+YisY0CDhEzvyR2fKRJOHkjHclQ/DHSmA4d44SZ7EFITskb18CgBCM0v9IVO0HN15r+VlKffDfXbvY5bhwEgE33Sj5nUhNakPKhPDXfkw3JECGO6NF8eVx0YUngqkc18BgMjOyWDkTGFrW97ricvnFt9h/nnBAUAbHWK5EzmzA/neaxdguCsfhjtSAMO9keMYuizyfGHYfrak8M1OICNmCIzNeK/oYjq36A4Tk8sBQI8WxE+u1ID32gkEw135cLYMQk0PQQl0+viYfbtHb8A4ghKUx1wvOLajISrysiDujhYEeVLWesTdHG7gP7T3GVq+6xNq5LDnjhTAnnsTQue9LvrnD93+YzXatG+4Wspp2BbHro9lCqRAEjDOmvzJlbTSrWsvHnvuyofhjhTAcEcK5Ulg4393AlnqQInrMC0Tw135cFgGIVRXRiJY35N67CuY0o4sp2FDLGt7VLYhlpUwqm4ZeguGO0Kofqx0iT/7U5GjBB7mRK4ElkQx3U7QwamsOgwCqBEMd4TUXPnDW/lB2+S7PvGopwlxZZjgwlBBVyMiqZDzC2d6hdDXcCeQRgPDHSG1xnFFYftLb4Zl/ji96MwBTsLzCu5eFsS90YJdbpS5NtzJ5jxO095n6Ef5GPGqhxdUkQJ4QVWd0Nnphf/8UR4bARxH6ujrDRyv238sQQn4raWUhh1x7NoYpkgGQhI+tSNX96DMtN7cixdUlQ/DHSmA4a5+pM8fF4bskSQ/BACBiYV42CdaDm5A8LzTXk4F/HCf2ZnA0izoCGBOZ3KZI5WbnrZ682/lFZLVX3/eoUMHfmtENcFwRwpguKsrSdL9gpO7ZK/SAEDDyl48aqaobRfea3lcyH0XzcpXh9ePDIS0u0Xuc0Aoahm5a2r3lhtXLOa9RvQ2DHekAIa7OmOZ0tvnis4eYIryAUCzs6vBaH+BiQXv9Vx9xX0TSd9ZPQ6+OFb5L0LbgxPiz/2F4zNKgBdUEWpmSEqnj4/Z8n3iEdMJkVZFXOTr9bPzg7Yxxfn81uNhTuy0fKhr273q4E+G1YArV6/yWxFSCMMdoeaI0NDU8/Qz+3aPTq8hHMeW3gx7vWZG8eXj/NbSrWtXs/z4qkcqkm7vlTqnFavDgEEjh+GOUPNFiY0NJ843W7xLy9GdrShjSwr4LV8oFHr1sNe5/htIy4CWCCIPUToGwZl6dsH07Agmi+dpmeg/cMwdKYBj7s2Q5OkDYau2pLYuv8UyDPP36bCNAQdlNDP/E79+Q0eve0TJdwIxFMFiPnYCQQphuCMFMNwRv6rNc48v4FbdfTOdprUO8Z0TOaMDSfE8LbO5w2EZhJCydTIggjypC0MFTsbEy1JudgTjcIKW7/qE+ILhjhB6N05SzpaV8FumlwVxd4wgyJOy0SPi8rnh52jvM7R81yf04TDcEULvVnThSOaPnxSHB3EyKY/FEgC+NmSir2BLL8pAAy6mcz3+pv3CGZxO8+Ew3BFC78JxslfP2LKSwlOBr9fNKrt3BXi9VqdBwrwuZPIE4WIHUoOE4FS203F6SRRTyOffkWYHL6giBfCCKnqbJOl+QcgeWXoyAGi0aS8eOUPU3rGO59Z94bDnJdyaGHbPY5blwFgEi7pR87uQIpxOU38Y7kgBDHekGMeVx0YUhu6h814DgMjOyWDULKFF23eeV99VIaNzuG8imcuvOACw0iW+70FOaU/ibJp6wXBHCmC4o1pwDF0acbro7EG2vAQIQsvBzWDkLMrItJZT3m/J34vp3MJI5kEeBwA9TYifXCgPc0z4usJwRwpguKN3YsuKi8ODS66e5GgZoSHSdR+l5z2B1NRR+OD3Xs+d5eDgU3ZJFPOqDADAy4LY5Ep1NcKIfzcMd6QAhjuqIyY/qzBsf1n0pX93AvEYQwiE1R72gZt1lNGw/X87gQhImG5HrupOmWt/cOvVGoY7UgDDHdWL9FliYcgeScojANDs6Nxi9o/VHsDLTkw5FfDzQ2bzQ1b6v51AvnWk9Kv/HUFv4FRIhNCH0rCyN5n7s8kX64Tm1jqugxuolhaasL4n9XCcwNeGLKNhQyxre1S29RFLsw1UYdOGPXekAPbc0XtiGSDIt3fv430P1cgsbtEd5nomBwAdxMQPzqSvDXZV/wNfDoQQf0iK931ZFXI1Ja4NF4QOotqLiceFnF840zuUjshUh64qXzDcEUJN1QhLMm6cYJcb1VILbmdx/U7TfuHM0yKMeAAMd4SQ0sgyUktvnwWWzzFyIQn+9uRTP+H6npSuEIJT2U7HcCcQAAx3hJDSFIbszj+yJXPD7PKY6/yWrCuExQ5kwniBvz3JchCQyNoGyVbdY8ppfutpSjDcEUJKotNriKCFOf36Re4fa7J3fCN9kcRv+RY6xC63N9NpSmSw+h5rF0wHJLJMsxynwdkySAGcLYP4VTlbhmPossjzRWcOMMX58qULxMOnC1qY815jeAb3zR3mXg4HAJ0MiFU9mt10Ggx3pACGO+JXtamQnKS8+PJx+erwBCXQdh0k9plK6hrwWykHcCyVXXKHTSnmAMDLgtjoQjkZN5elCzDckQIY7ohfCue5M4W5RecOyS+xEiItvQHj9Dz9CKEGv1VLWfgjiV0ezWRXAAEw3obc4ELa6Kl/xGO4IwUw3BG/avkSkyzzedHZA/JLrJRBC/1Bk3R6DQGS5yGUfAlseMBsfcRWMKBBwmcdydU9KAOe/440LhjuSAEMd8Svd35DVZIUUxC6R/byKQAIzazEI2dodnLhvRkvSrkf77/ZCcRIBN90o+Z1ITXVdCeQ5nWFASHUOInsHFsu3G78yTKBsZks81np7XMNUUsbHWKXG3VnlGBgKyJPAkuimA7BdEAiy6pDF7c67LkjBbDnjvhV97VlOJm05FqIVtfeAtPWDdqki+nc/0UysXkcADi3IH5ypfqr104gGO5IAQx3xC/eFw7jBcvB8TR2UST7rOTNdJpfXKlu6rITCA7LIISaKZIAXxsyfrxgfU9KrAEX0zmnk/TUK292fWrqMNwRQk0KxwGv4w3aAljsQCZPEC52IIUkHHjKtguSLYliimQ8VqICGO4Ioaak9M75rM3zJE8f8lussQjW96QejhX42pDl/9sJZEMsK22yO4FguCOEmpLSG/9Inydl71iUu3c1/foFv4W3FxNBntTtUYJ+ZkROBSyJYroep4NTm+RsGrygihTAC6qIXzxeUOWkkpLrIcUXjrIVpUAQ2s4DxSNmUPpGH15yNRfTuQW3mUf5HAC4mhI/uVDuZk3pWiuGO1IAwx3xi/fZMmxpUfGlYyVXTnAMTWho6rqP1POeSGpq81W+HM1CYBK78i6TWQ4AMNyS2ORKtRc3jYjHcEcKYLgjfjXQVEg6O73wnz/KYyOA40gdff1Bk3TdRwDJ81dOS2nYEceuiWGKZSAk4VM78vseVEstfivhH4Y7UgDDHfGrQee500rVTwAADatJREFUS58lFobskaQ8AgCBaWuxzzQtR3fea8ko41bfYwOTWJoFXSF82Ylc7kjpCnmvhzcY7kgBDHfELyV8iaki7nbB37vp7HQA0LC2F4+cJWrbmfdaEgu4FXfZ4FQWACx0iBVO5IwOJNUox2kw3JECGO6IX8r5hirH0KU3w4rOHWJLCoEgdFwHG06c3xAV3criFkUyN15zANDRgFjdKHcCaXQNQgih90NQAl33keYr9otHTCc0NAXGZg1UUW9T4voIQZAnZatPJBRwfuGMVxgt3/Wp8cCeO1IAe+6IX8pfW4YpzCW1dAmNhq1RxsK+JPa7u0xW+ZudQNa7kG0bx04g2HNHCKkhSmzc0MkOAEIS/O3JZD/hyu6kpgCCU9mOwfTsCCa7oqFrfjcMd4QQ+iC6QljVnUryFfjbkwwHAYlsh2DZhli2glFlqzDcEULNDltSKHuZzG+ZrXWIXW5U7FjBcEsiXwJLohi7IFXuBILhjhBqdorOHnj9y5zcP9bQuZn8ltzZkDg1SHBhqMDRmHhRys2OYFxC6EsZKgh4DHeEULNDaOoQlKA85vrr9f6Fp/ay5SX8lu9lQdwdLQjypKz1iLs5nGcY7X2Glu/6pDQ4WwYpgLNlEL8a4U5MTEF20fm/Sm+fBZYltfX0PH11+40mhBr81iJl4bd4dtU9pkAKJAHjrMmfXEkrXWVMp8FwRwpguCN+NcJwl5NlPisM3VsRfwcAKENTfe+JOr2HAsFz+OZJYOMDZssjVsKAtgC+6kwudaDEPP8dqQ7DHSmA4Y741WjDXU6SdL8gZI8sPRkANNrYiUfOELV34L2W5yXc8mj24FOWAzAWwaJu1PwupIjnVc7+heGOFMBwR/xq5OEOAMBx5bERhaF76LzXACCyczIY7S9sZcN7PVHZ3KI7zNVXHAC0FxNrnMnxNmRDDNM0xnDPz8+/e/duWlrazJkzKw/eunUrMjLSzc3N2dn57VMw3PmF4Y741QTCHQAU7gQyfDolNua9oovp3NeRzMM8DgBcTIifXKl+ZgQA5OXlxcXFubm5ER88NNQYZ8ts3bo1ODh4/vx/V/w5derU0qVL27RpM2/evEuXLqmwbQghNUZoiPQ8/cyWB+q6jyRIqiwqPHffjw1RkZcFETNGsN+DMteGO9mcx2l6wLFcjwn+Th99PeL3Wx0Hjvk18OAHVtEYe+5yBgYGBQUF8tteXl5r1651cXG5dOnS9u3bT548We3B2HPnF/bcEb+aSs+9Kjo7vfD0Hzq9h2ja92i4WuQ7gayNYYr+2QyWTmDfX37c8sgnD47vFIvF711yY+y5v+3x48edOnUCgE6dOj1+/FjVzUEIqT+BiYXxp8saNNkBQEcAix3IhPEC8fOblckOAOmOU7bv2f8hJQs+tGkf4Msvv4yJial6ZMuWLT179nz7kVKplKIoABAKhRKJREntQwghpTDXBnMtKOS4ylmYBMsIhR+0z1NDhTtN0y9fvjQwMDAwMKh6/MGDBzk5OT169BCLxVu2bGFZtuq9GhqKZ362bt365cuX7du3f/HiRZs2bRqozQghpBIEQXw0tP8PiRfpjt4AABxrEXvgi9W7P6TMBhmWmT17tlgsbteu3Y4dO6oe/+STT0aPHv3zzz/b2dlFR0cLhULRf8kvEOfl5aWkpLAsm5KSkpeXBwB+fn7bt28vLS3dsWOHn59fQ7QZIYTqRZISV3DiN7akkJfSFnw+Y0R+mNXRT40uruv014QVM8Z94HWvBrmgevfuXQsLi6+++srBwWH58uXygzdu3Bg3blxCQoKhoeHGjRsvXLhw4cIFhacfOXLkxIkT8tujR4+eNGkSwzA//PDDzZs3BwwYsHjxYpKs/jcJL6jyCy+oIn41xQuq75S1daE0NY7U1Nbz9NP1GMPL8vFFRUWJiYkuLi4fXlQDzpbx9fWtGu5ff/11UVHRnj17AODVq1cWFha5ubmGhoa81JWZmWlra1t1VOfGjRuWlpa8FN4MlZSU6OrqqroVSH2oZbiz2enl4Udkj24BAKlvpDnAV6OHJ7zV9WwI2trab/dxq1HeBdWXL1927dpVftvc3FwoFKanp/MV7gCgp6eXkJBQ+aNYLH7nk0c14TgOwx3xSC3DHXQ76M9cKU1LKAjZI02NKwvZJb19Rjx0ipaju6pbBqDMcK+oqKh68VdDQ6O8vJzH8gmC4PFPBUII1YWGdUfTuT+Xx0YUng6kXz/P/WONyM7RYORMYet2qm2Y8vq2ZmZmubm58tvl5eUlJSXm5uZKqx0hhBoKQWg5urdcutvQby6payBJinn9y1e5f6yhc1V5FVB54d67d+8rV67Ib1+7ds3KyqpVq1ZKqx0hhBoUQQl0+viYf7dPf8jHhEBYHnP99dpZ+UHb+JpOU+/2NMQF1QsXLty7d+/QoUPm5uYDBw4cMmSIg4NDWVmZvb39mDFjevfu/d1333311Vdz587lq0acLcMvnC2D+KWeY+41Y/KzCsP2l0VfAo4jtfX0B0/S9Rij5DY0SM+9rKwsPz/fx8fHyckpPz9f/p1SbW3tGzduiESisLCw1atX85jsCCHUqFCGpkaTF7Vc/LtmZ1e2rFiWkar8NjTehcPqBXvu/MKeO+JXc+u5V1WREC00t6YMWii5XlWuLYMQQmpPs6OCLSiUAGeCI4SQGsJwRwghleEYmq0oa4iSMdwRQkhlSq+HZv7wSXF4EEfL+C0Zwx0hhFRGkvKILS0qPBX4esNn5bERwN8MF5wtgxTA2TKIX815tsw7lT+6XXgqkH79HAA0rDsajPbXsO744cXibBmEEFIlrS69tDq7lkWHF54KlKYlSJ89xnBHCCG1QBDaPb20uvUtuRmm4zaclyIx3BFCqFEgRFp6A8bxVRpeUEUKdOrUqbi4WNWtQOpjx44dP/74o6pb0bxguCMFCgsL1eNKO2okKioq+N2/Ab0ThjtCCKkhDHeEEFJDajLPPTc3t1+/fp07d1Z1Q9TE5cuX3d3dBQK83o74kZaWJpVK7ezsVN0QNbFx40Zra+vaH6Mm4Q4A9+7dS05OVnUrEEKowXl7exsYGNT+GPUJd4QQQpVwzB0hhNQQhjtCCKkhDHf0DhzH3blzZ9u2bYGBgYWFqtnHHaml3Nzc4ODgiooKVTdEPWG4o3d48eLFzz//rKurm56e7uLigl9FQXxZsGDBrFmzsMfQQHCuG3oHS0vLoKAg+e1z584lJiY6OTmptklIDYSGhtrZ2Zmamqq6IWoLe+6orp4/f56ZmYlTldGHKygo2Llz5zfffKPqhqgzDHdUJwUFBX5+fjt37tTR0VF1W1CTN3/+/BUrVmhoaKi6IeoM57k3dyUlJVlZWdbW1iT571/6ioqK6OhosVjctWtXACgqKho2bNiCBQvGjh2rupaipkEikaSnp1tYWFTdd4lhmLt37xIE0aNHD5qmDQwMOnXqBACPHj3q0KHDpUuXjI2NVddk9YRj7s3Xs2fPfHx8kpKSaJrOy8szNDSUH09OTh4wYEDbtm1fvnzp4ODwxx9/DBs2zMfHx9HRMSUl5f/buWOXZMIAjuOPERkVyTXkNSUnNEsUNkcEjYEK0R9gW0GR0tRY0BI4aC6RU7hEEEEcVCYNYUVES9GSk2BeRJQgh+9wINk7yUve6+P3Mz33nMNvePjx+HB3qqr29PTYmxz/J9M0x8fHHx4eKpXK5eXlxMSENf/29jY5OelwOKrVamdnp67rn5+f1q2RkZGTkxOa/TdwLNO+FEVJJpOPj48/5tfX12dnZ8/Ozu7u7u7v7/f394eGhm5vb6PRaDQa/fv3gKWjo2Nra6tQKLhcru/zsVjM7Xbncjnr72A8Hq/dmpmZ6e7ubnrStsCxTLsrFAqqqtZ27tVqtbe3N5vNjo6OCiHW1tby+XwqlbI7JlqJoijHx8e1nbvP54tEInNzc0KIvb29WCx2dXVla8C2wM4ddYrF4tfXV+2Dcx6PJ5/P25oILe/l5WV4eNgas6KahnJHHet1wdpjDE6nk7eW8I/K5TIrqvkod9QZHBx0OByvr6/WZbFYVFXV3khodaqqsqKaj3JHHafT6fP5MpmMdZnJZPx+v72R0Or8fv/FxYU1ZkU1DY9CtrXNzc2Pjw8hxPb2dl9f38rKihBieXl5dXW1q6vr6ekpm83u7OzYHRMtI5lMlkqlcrmcSqXOz88XFhZcLtfS0tL09LTH4zFNMx6Pn56e2h2zLVDubc0wDCFEJBIpl8uVSsWanJ+fdzqdBwcH/f392WzW7XbbmhGt5P393TCMxcVFIYRhGKZpCiH8fv/h4eHu7q7D4Tg6OrIexMJv41FIAJAQZ+4AICHKHQAkRLkDgIQodwCQEOUOABKi3AFAQpQ7AEiIcgcACVHuACAhyh0AJES5A4CEKHcAkBDlDjQskUiEw2HDMBKJRDqdtsZ2hwLqUO5AY9LpdCgUEkKEQqFwOBwMBjVN46v3+N9Q7kBjNE1TFEXX9Y2NDWsml8tpmmZvKuAHvucONMwwDE3TakcxiqLkcjmv12tvKuA7du5Aw3RdHxsbs8bX19cDAwNer/f5+dneVMB3lDvQMF3Xp6amauNAIGAYxs3Njb2pgO8od6BhpVIpEAhYY6vldV0PBoO2hgLqcOYOABJi5w4AEqLcAUBClDsASIhyBwAJUe4AIKE/zEMiRfv3fsoAAAAASUVORK5CYII=",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {
"filenames": {
"image/png": "/Users/driscoll/repos/fnc-julia/_build/jupyter_execute/localapprox/integration_18_0.png",
"image/svg+xml": "/Users/driscoll/repos/fnc-julia/_build/jupyter_execute/localapprox/integration_18_0.svg"
}
},
"output_type": "execute_result"
}
],
"source": [
"plot(n,abs.(err),m=:o,label=\"results\",\n",
" xaxis=(:log10,L\"n\"),yaxis=(:log10,\"error\"),\n",
" title=\"Convergence of trapezoidal integration\")\n",
"\n",
"# Add line for perfect 2nd order.\n",
"plot!(n,3e-3*(n/n[1]).^(-2),l=:dash,label=L\"O(n^{-2})\")"
]
},
{
"cell_type": "markdown",
"id": "dfa99523",
"metadata": {},
"source": [
"```{raw} html\n",
"
\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"\n",
"\n",
"## Extrapolation\n",
"\n",
"If evaluations of $f$ are computationally expensive, we want to get as much accuracy as possible from them by using a higher-order formula. There are many routes for doing so; for example, we could integrate a not-a-knot cubic spline interpolant. However, splines are difficult to compute by hand, and as a result different methods were developed before computers came on the scene.\n",
"\n",
"```{index} ! extrapolation\n",
"```\n",
"\n",
"Knowing the structure of the error allows the use of **extrapolation** to improve accuracy. Suppose a quantity $A_0$ is approximated by an algorithm $A(h)$ with an\n",
"error expansion\n",
"\n",
"```{math}\n",
" :label: extraperror\n",
" A_0 = A(h) + c_1 h + c_2 h^2 + c_3 h^3 + \\cdots.\n",
"```\n",
"\n",
"Crucially, it is not necessary to know the values of the error constants $c_k$, merely that they exist and are independent of $h$. \n",
"\n",
"Using $I$ for the exact integral of $f$, the trapezoid formula has\n",
"\n",
"```{math}\n",
" I = T_f(n) + c_2 h^2 + c_4 h^{4} + \\cdots,\n",
"```\n",
"\n",
"as proved by the Euler–Maclaurin formula {eq}`eulermaclaurin`. The error constants depend on $f$ and can't be evaluated in general, but we know that this expansion holds. For convenience we recast the error expansion in terms of $n=O(h^{-1})$:\n",
"\n",
"```{math}\n",
" :label: traperrorexpansion\n",
" I = T_f(n) + c_2 n^{-2} + c_4 n^{-4} + \\cdots.\n",
"```\n",
"\n",
"We now make the simple observation that\n",
"\n",
"```{math}\n",
" :label: traperrorexpansion2n\n",
" I = T_f(2n) + \\tfrac{1}{4} c_2 n^{-2} + \\tfrac{1}{16} c_4 n^{-4} + \\cdots.\n",
"```\n",
"\n",
"It follows that if we combine {eq}`traperrorexpansion` and {eq}`traperrorexpansion2n` correctly, we can cancel out the second-order term in the error. Specifically, define\n",
"\n",
"```{math}\n",
" :label: nc-simpson\n",
" S_f(2n) = \\frac{1}{3} \\Bigl[ 4 T_f(2n) - T_f(n) \\Bigr].\n",
"```\n",
"\n",
"(We associate $2n$ rather than $n$ with the extrapolated result because of the total number of nodes needed.) Then\n",
"\n",
"```{math}\n",
" :label: extraplevel1\n",
" I = S_f(2n) + O(n^{-4}) = b_4 n^{-4} + b_6 n^{-6} + \\cdots.\n",
"```\n",
"\n",
"```{index} ! Simpson's formula\n",
"```\n",
"\n",
"The formula {eq}`nc-simpson` is called **Simpson's formula**, or *Simpson's rule*. A different presentation and derivation are considered in [Exercise 4](problem-simpson).\n",
"\n",
"Equation {eq}`extraplevel1` is another particular error expansion in the form {eq}`extraperror`, so we can extrapolate again! The details change only a little. Considering that\n",
"\n",
"```{math}\n",
" I = S_f(4n) = \\tfrac{1}{16} b_4 n^{-4} + \\tfrac{1}{64} b_6 n^{-6} + \\cdots,\n",
"```\n",
"\n",
"the proper combination this time is\n",
"\n",
"```{math}\n",
" :label: nc-sixth\n",
" R_f(4n) = \\frac{1}{15} \\Bigl[ 16 S_f(4n) - S_f(2n) \\Bigr],\n",
"```\n",
"\n",
"which is sixth-order accurate. Clearly the process can be repeated to get eighth-order accuracy and beyond. Doing so goes by the name of **Romberg integration**, which we will not present in full generality.\n",
"\n",
"## Node doubling\n",
"\n",
"Note in {eq}`nc-sixth` that $R_f(4n)$ depends on $S_f(2n)$ and $S_f(4n)$, which in turn depend on $T_f(n)$, $T_f(2n)$, and $T_f(4n)$. There is a useful benefit realized by doubling of the nodes in each application of the trapezoid formula. As shown in {numref}`figure-node-doubling`, when doubling $n$, only about half of the nodes are new ones, and previously computed function values at the other nodes can be reused.\n",
"\n",
"\n",
"```{figure} figures/node-doubling.svg\n",
":name: figure-node-doubling\n",
"Dividing the node spacing by half introduces new nodes only at midpoints, allowing the function values at existing nodes to be reused for extrapolation.\n",
"```\n",
"\n",
"Specifically, we have \n",
"\n",
"```{math}\n",
":label: nc-doubling\n",
"\\begin{split}\n",
" T_f(2m) & = \\frac{1}{2m} \\left[ \\frac{1}{2} f(a) + \\frac{1}{2} f(b) + \\sum_{i=1}^{2m-1} f\\Bigl( a + \\frac{i}{2m} \\Bigr) \\right]\\\\[1mm]\n",
" & = \\frac{1}{2m} \\left[ \\frac{1}{2} f(a) + \\frac{1}{2} f(b)\\right] + \\frac{1}{2m} \\sum_{k=1}^{m-1} f\\Bigl( a+\\frac{2k}{2m} \\Bigr) + \\frac{1}{2m} \\sum_{k=1}^{m} f\\Bigl( a+\\frac{2k-1}{2m} \\Bigr) \\\\[1mm]\n",
" &= \\frac{1}{2m} \\left[ \\frac{1}{2} f(a) + \\frac{1}{2} f(b) + \\sum_{k=1}^{m-1} f\\Bigl( a+\\frac{k}{m} \\Bigr) \\right] + \\frac{1}{2m} \\sum_{k=1}^{m} f\\Bigl( a+\\frac{2k-1}{2m} \\Bigr) \\\\[1mm]\n",
" &= \\frac{1}{2} T_f(m) + \\frac{1}{2m} \\sum_{k=1}^{m-1} f\\left(t_{2k-1} \\right),\n",
"\\end{split}\n",
"```\n",
"\n",
"where the nodes referenced in the last line are relative to $n=2m$. Hence in passing from $n=m$ to $n=2m$, new integrand evaluations are needed only at the odd-numbered nodes of the finer grid. \n",
"\n",
"(demo-int-extrap)=\n",
"```{proof:demo}\n",
"```\n",
"\n",
"```{raw} html\n",
"
\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"We estimate $\\displaystyle\\int_0^2 x^2 e^{-2x}\\, dx$ using extrapolation. First we use `quadgk` to get an accurate value."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "8d4b2704",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Q = 0.1904741736116139\n"
]
}
],
"source": [
"f = x -> x^2*exp(-2*x);\n",
"a = 0; b = 2; \n",
"Q,_ = quadgk(f,a,b,atol=1e-14,rtol=1e-14)\n",
"@show Q;"
]
},
{
"cell_type": "markdown",
"id": "eb803355",
"metadata": {},
"source": [
"We start with the trapezoid formula on $n=N$ nodes."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "cb1776b8",
"metadata": {},
"outputs": [],
"source": [
"N = 20; # the coarsest formula\n",
"n = N; h = (b-a)/n;\n",
"t = h*(0:n); y = f.(t);"
]
},
{
"cell_type": "markdown",
"id": "f5f35686",
"metadata": {},
"source": [
"We can now apply weights to get the estimate $T_f(N)$."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "a8bded71",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1-element Vector{Float64}:\n",
" 0.19041144993926787"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"T = [ h*(sum(y[2:n]) + y[1]/2 + y[n+1]/2) ]"
]
},
{
"cell_type": "markdown",
"id": "1a0b8fec",
"metadata": {},
"source": [
"Now we double to $n=2N$, but we only need to evaluate $f$ at every other interior node and apply {eq}`nc-doubling`."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "64873727",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" 0.19041144993926787\n",
" 0.19045880585951175"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 2n; h = h/2; t = h*(0:n);\n",
"T = [ T; T[end]/2 + h*sum( f.(t[2:2:n]) ) ]"
]
},
{
"cell_type": "markdown",
"id": "e52432c2",
"metadata": {},
"source": [
"We can repeat the same code to double $n$ again."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "c619a079",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.19041144993926787\n",
" 0.19045880585951175\n",
" 0.1904703513046443"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 2n; h = h/2; t = h*(0:n);\n",
"T = [ T; T[end]/2 + h*sum( f.(t[2:2:n]) ) ]"
]
},
{
"cell_type": "markdown",
"id": "01236779",
"metadata": {},
"source": [
"Let us now do the first level of extrapolation to get results from Simpson's formula. We combine the elements `T[i]` and `T[i+1]` the same way for $i=1$ and $i=2$."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "fa4b9546",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" 0.19047459116625973\n",
" 0.19047419978635513"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S = [ (4T[i+1]-T[i])/3 for i in 1:2 ]"
]
},
{
"cell_type": "markdown",
"id": "ea51ec98",
"metadata": {},
"source": [
"With the two Simpson values $S_f(N)$ and $S_f(2N)$ in hand, we can do one more level of extrapolation to get a sixth-order accurate result."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "ed064f33",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.1904741736943615"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R = (16S[2] - S[1]) / 15"
]
},
{
"cell_type": "markdown",
"id": "bb574bd9",
"metadata": {},
"source": [
"::::{panels}\n",
":column: col-7 left-side\n",
":card: border-0 shadow-none\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.5\\textwidth}\n",
"```\n",
"We can make a triangular table of the errors:\n",
"\n",
"```{raw} latex\n",
"\\end{minipage}\\hfill\n",
"```\n",
"---\n",
":column: col-5 right-side\n",
":card: shadow-none comment\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.4\\textwidth}\\begin{mdframed}[default]\\small\n",
"```\n",
"The value `nothing` equals nothing except `nothing`.\n",
"```{raw} latex\n",
"\\end{mdframed}\\end{minipage}\n",
"```\n",
"::::"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "71b32942",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"┌─────────────┬────────────┬─────────────┐\n",
"│\u001b[1m order 2 \u001b[0m│\u001b[1m order 4 \u001b[0m│\u001b[1m order 6 \u001b[0m│\n",
"├─────────────┼────────────┼─────────────┤\n",
"│ -6.27237e-5 │ nothing │ nothing │\n",
"│ -1.53678e-5 │ 4.17555e-7 │ nothing │\n",
"│ -3.82231e-6 │ 2.61747e-8 │ 8.27476e-11 │\n",
"└─────────────┴────────────┴─────────────┘\n"
]
}
],
"source": [
"err = [ T.-Q [nothing;S.-Q] [nothing;nothing;R-Q] ]\n",
"pretty_table(err,[\"order 2\",\"order 4\",\"order 6\"])"
]
},
{
"cell_type": "markdown",
"id": "501b4823",
"metadata": {},
"source": [
"If we consider the computational time to be dominated by evaluations of $f$, then we have obtained a result with about twice as many accurate digits as the best trapezoid result, at virtually no extra cost.\n",
"\n",
"```{raw} html\n",
"
\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"## Exercises\n",
"\n",
"(problem-integration-tests)=\n",
"% must be kept as #1\n",
"\n",
"1. ⌨ For each integral below, use {numref}`Function {number} ` to estimate the integral for $n=10\\cdot 2^k$ nodes for $k=1,2,\\ldots,10$. Make a log-log plot of the errors and confirm or refute second-order accuracy. (These integrals were taken from {cite}`baileyComparisonThree2005`.)\n",
"\n",
" **(a)** $\\displaystyle \\int_0^1 x\\log(1+x)\\, dx = \\frac{1}{4}$\n",
"\n",
" **(b)** $\\displaystyle \\int_0^1 x^2 \\tan^{-1}x\\, dx = \\frac{\\pi-2+2\\log 2}{12}$\n",
"\n",
" **(c)** $\\displaystyle \\int_0^{\\pi/2}e^x \\cos x\\, dx = \\frac{e^{\\pi/2}-1}{2}$\n",
"\n",
" **(d)** $\\displaystyle \\int_0^1 \\sqrt{x} \\log(x) \\, dx = -\\frac{4}{9}$ (Note: Although the integrand has the limiting value zero as $x\\to 0$, it cannot be evaluated naively at $x=0$. You can start the integral at $x=\\macheps$ instead.)\n",
"\n",
" **(e)** $\\displaystyle \\int_0^1 \\sqrt{1-x^2}\\,\\, dx = \\frac{\\pi}{4}$\n",
" \n",
"2. ✍ The Euler–Maclaurin error expansion {eq}`eulermaclaurin` for the trapezoid formula implies that if we could cancel out the term due to $f'(b)-f'(a)$, we would obtain fourth-order accuracy. We should not assume that $f'$ is available, but approximating it with finite differences can achieve the same goal. Suppose the forward difference formula {eq}`forwardFD21` is used for $f'(a)$, and its reflected backward difference is used for $f'(b)$. Show that the resulting modified trapezoid formula is\n",
" \n",
" ```{index} ! Gregory integration formula\n",
" ```\n",
"\n",
" ```{math}\n",
" :label: gregory\n",
" G_f(h) = T_f(h) - \\frac{h}{24} \\left[ 3\\Bigl( f(t_n)+f(t_0) \\Bigr) -4\\Bigr( f(t_{n-1}) + f(t_1) \\Bigr) + \\Bigl( f(t_{n-2})+f(t_2) \\Bigr) \\right],\n",
" ```\n",
"\n",
" which is known as a **Gregory integration formula**.\n",
"\n",
"\n",
"3. ⌨ Repeat each integral in Exercise 1 above using Gregory integration {eq}`gregory` instead of the trapezoid formula. Compare the observed errors to fourth-order convergence.\n",
"\n",
" (problem-simpson)=\n",
"4. ✍ Simpson's formula can be derived without appealing to extrapolation.\n",
"\n",
" **(a)** Show that\n",
"\n",
" ```{math}\n",
" p(x) = \\beta + \\frac{\\gamma-\\alpha}{2h}\\, x + \\frac{\\alpha-2\\beta+\\gamma}{2h^2}\\, x^2\n",
" ```\n",
"\n",
" interpolates the three points $(-h,\\alpha)$, $(0,\\beta)$, and $(h,\\gamma)$.\n",
"\n",
" **(b)** Find\n",
"\n",
" ```{math}\n",
" \\int_{-h}^h p(s)\\, ds,\n",
" ```\n",
"\n",
" where $p$ is the quadratic polynomial from part (a), in terms of $h$, $\\alpha$, $\\beta$, and $\\gamma$.\n",
" \n",
" **(c)** Assume equally spaced nodes in the form $t_i=a+ih$, for $h=(b-a)/n$ and $i=0,\\ldots,n$. Suppose $f$ is approximated by $p(x)$ over the subinterval $[t_{i-1},t_{i+1}]$. Apply the result from part (b) to find\n",
"\n",
" ```{math}\n",
" \\int_{t_{i-1}}^{t_{i+1}} f(x)\\, dx \\approx \\frac{h}{3} \\bigl[ f(t_{i-1}) + 4f(t_i) + f(t_{i+1}) \\bigr].\n",
" ```\n",
"\n",
" (Use the change of variable $s=x-t_i$.)\n",
"\n",
" **(d)** Now also assume that $n=2m$ for an integer $m$. Derive Simpson's formula,\n",
"\n",
" ```{math}\n",
" :label: simpson\n",
" \\begin{split}\n",
" \\int_a^b f(x)\\, dx \\approx \\frac{h}{3}\\bigl[ &f(t_0) + 4f(t_1) + 2f(t_2) + 4f(t_3) + 2f(t_4) + \\cdots\\\\\n",
" &+ 2f(t_{n-2}) + 4f(t_{n-1}) + f(t_n) \\bigr].\n",
" \\end{split}\n",
" ```\n",
"\n",
"5. ✍ Show that the Simpson formula {eq}`simpson` is equivalent to $S_f(n/2)$, given the definition of $S_f$ in {eq}`nc-simpson`.\n",
"\n",
"6. ⌨ For each integral in Exercise 1 above, apply the Simpson formula {eq}`simpson` and compare the errors to fourth-order convergence.\n",
"\n",
"7. ⌨ For $n=10,20,30,\\ldots,200$, compute the trapezoidal approximation to\n",
"\n",
" ```{math}\n",
" \\int_{0}^1 \\frac{1}{2.01+\\sin (6\\pi x)-\\cos(2\\pi x)} \\,d x \\approx 0.9300357672424684.\n",
" ```\n",
"\n",
" Make two separate plots of the absolute error as a function of $n$, one using a log-log scale and the other using log-linear. The graphs suggest that the error asymptotically behaves as $C \\alpha^n$ for some $C>0$ and some $0<\\alpha<1$. How does this result relate to {eq}`eulermaclaurin`?\n",
"\n",
"\n",
"8. ⌨ For each integral in Exercise 1 above, extrapolate the trapezoidal results two levels to get sixth-order accurate results, and compute the errors for each value.\n",
"\n",
"9. ✍ Find a formula like {eq}`nc-sixth` that extrapolates two values of $R_f$ to obtain an eighth-order accurate one."
]
}
],
"metadata": {
"jupytext": {
"cell_metadata_filter": "-all",
"formats": "md:myst",
"text_representation": {
"extension": ".md",
"format_name": "myst",
"format_version": 0.13,
"jupytext_version": "1.10.3"
}
},
"kernelspec": {
"display_name": "Julia 1.7.1",
"language": "julia",
"name": "julia-fast"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.1"
},
"source_map": [
16,
20,
41,
43,
50,
53,
57,
60,
64,
67,
244,
247,
272,
275,
279,
282,
286,
295,
299,
306,
426,
431,
435,
439,
443,
445,
449,
452,
456,
459,
463,
465,
469,
471,
496,
499
]
},
"nbformat": 4,
"nbformat_minor": 5
}