{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## The spruce budworm model\n", "\n", "We consider a dimensionless model for budworm dynamics. The budworm spread by diffusion on a one-dimensional domain, $0 \\leq x \\leq L$, and undergo logistic growth and predation by birds:\n", "\n", "\\begin{equation}\n", "\\label{equation:budwormPDE}\n", "\\frac{\\partial u}{\\partial t}=D\\frac{\\partial^2u}{\\partial{x}^2}+f(u),\n", "\\quad\\mbox{where}\\quad \n", "f(u)=ru\\left(1-\\frac{u}{q}\\right)-\\frac{u^2}{1+u^2}.\n", "\\end{equation}\n", "\n", "We suppose that exterior to the domain conditions are extremely hostile to budworm so that we have the boundary conditions\n", "\n", "\\begin{equation}\n", "\\label{equation:budwormBCs}\n", "u(0,t)=0, \\quad u(L,t)=0.\n", "\\end{equation}\n", "\n", "Note that $f'(0)=r>0$.\n", "\n", "We want to consider the minimum domain size for spatial structure. We do this by consider the system at steady state so that \n", "\n", "\\begin{equation}\n", "0=D\\frac{\\partial^2u}{\\partial{x}^2}+f(u).\n", "\\end{equation}\n", "\n", "Multiplying by $\\partial{u}/\\partial{x}$ and integrating with respect to $x$, we have \n", "\n", "\\begin{equation}\n", "0=\\int{}D\\frac{\\partial u}{\\partial x}\\frac{\\partial^2u}{\\partial{x}^2}\\,\\text{d}x+\\int{}\\frac{\\partial u}{\\partial x}f(u)\\,\\text{d}x.\n", "\\end{equation}\n", "\n", "Thus we have\n", "\n", "\\begin{equation}\n", "\\frac{1}{2}D\\left(\\frac{\\partial u}{\\partial x}\\right)^2+F(u)=\\mbox{constant}=F(u_\\text{max})\n", "\\quad\\mbox{where}\\quad\n", "F'(u)=f(u).\n", "\\end{equation}\n", "\n", "We can therefore find a relation between $L$, $D$, integrals of \n", "\n", "\\begin{equation}\n", "F(u)\\stackrel{\\text{def}}=\\int_0^uf(y)\\,\\text{d}y,\n", "\\end{equation} \n", "\n", "and $\\max(u)$, the maximum size of the outbreak, as follows:\n", "\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial x}=-\\left(\\frac{2}{D}\\right)^{\\frac{1}{2}}\\sqrt{F(u_\\text{max})-F(u)},\n", "\\end{equation}\n", "\n", "since $x>0$. Therefore $\\partial{u}/\\partial{x}<0$. Integrating, gives\n", "\n", "\\begin{equation}\n", "2\\int_0^{L/2}\\text{d}x=-(2D)^\\frac{1}{2}\\int_{u_\\text{max}}^0\\frac{1}{\\sqrt{F(u_\\text{max})-F(\\bar{u})}}\\,\\text{d}{\\bar{u}},\n", "\\end{equation}\n", "\n", "and hence\n", "\n", "\\begin{equation}\n", "\\label{equation:L}\n", "L=(2D)^\\frac{1}{2}\\int^{u_\\text{max}}_0\\frac{1}{\\sqrt{F(u_\\text{max})-F(\\bar{u})}}\\,\\text{d}{\\bar{u}}.\n", "\\end{equation}\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import integrate\n", "from scipy import optimize\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "\n", "#set the plotting parameters\n", "plt.rcParams.update({\n", " \"text.usetex\": True,\n", " \"font.family\": \"serif\",\n", " \"font.serif\": [\"Times New Roman\"],\n", " \"font.size\": 18,\n", " \"axes.linewidth\": 1.5,\n", "})" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# model parameters\n", "r=0.6;\n", "q=6.2;\n", "D=0.1;\n", "args=[r,q,D]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "f = lambda x, r, q : r*x*(1.0-x/q)-x**2.0/(1.0+x**2);" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "First root found 1.0366846996193442\n", "Second root found 1.7544014258756127\n", "Third root found 3.408913874504972\n" ] } ], "source": [ "root1=optimize.brentq(f,0.9,1.1,args=(r,q))\n", "print(\"First root found %s\" % root1)\n", "root2=optimize.brentq(f,1.1,3.0,args=(r,q))\n", "print(\"Second root found %s\" % root2)\n", "root3=optimize.brentq(f,3.0,5.0,args=(r,q))\n", "print(\"Third root found %s\" % root3)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEeCAYAAAB40PUWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXzV1Z3/8de52febnYSwXRL2NQFEQIUabN1qq7hWB1ttcPx12o7TunSctjPjjMV2Oq1jZwDF1n0BrbXVqkQrCqIQArLKkoSwBbIvJGQ/vz/uNxhjQm7u9r3fm8/z8eBRk3vPN58GeHPy+Z5zvkprjRBCCGuzmV2AEEIIz0mYCyFEEJAwF0KIICBhLoQQQUDCXAghgoCEuRBCBAEJcyGECAKhZhfgDqVUAVBrfOjQWj/iwhg7cAOwVGt9vTeuKYQQgcJyYd4Tulrr9cbHDqXUaq31ivOMyQUcOMPa4Y1rCiFEIFFW2wGqlNqutc7r87kSrfV4F8bmAo/3M97tawohRCCwVM/caJXk9vNSvVIqP1CuKYQQ/ma1NosDqO/n8/22T3x9TaXUDiAVOOzm1xZCiKHIBqq01rP7vmC1ME/i85uUvdUDdn9d0+ixFwBTo6KiwubNmzfSza8thBAu27lzJw0NDf2+ZrUwDwha6zXAGqXU+/Pmzbvk/fffN7skIcQwsHjxYjZu3NhvJ8BSPXNDUj+fc3dW7strCiGE31gtzIvoP2STgOIAuqYQQviVpcJca10PlBorUHqza60LA+WaQgjhb5YKc8NKnDcfgXNrxwt7fexQSq3rJ5yh/3bKoNcUQohAZ7kboFrrNUqpAmMNuB3n1vveOzUdQD7O4K4HZ8ADy4ClQK5SaiVQYtzIdOWaQggR0CwX5nBuNclArxUCiX0+Vwo8Yvwa8jWFECLQWbHNIoQQog8JcyGECAIS5kIIEQQkzIUQIghImAshRBCQMBdCiCAgYS6EEEFAwlwIIYKAhLkQQgQBCXMhhAgCEuZCCBEELHk2y3Cy72QjLxcd40T9WRwpMdw8bzRjU2LMLksIEWAkzAOU1ppVG0v55dufER5qY0xSDO8fqGTtpjLuv3wSdywah1LK7DKFEAFCwjxArfmglJVvfcZVMzL4j29MJyE6jMqmVv7ltT089MZ+zrR18sP8CWaXKYQIENIzD0BFR2pZ+dZnXDF9BI/eNJuE6DAA0uIi+b9v5XF9Xha/KTzEaztOmFypECJQSJgHmM6ubh58bQ8ZCVE8smwmNtsXWyk2m+Lha6czb2wS//zH3ZTXNJtUqRAikEiYB5iXi47z2akm/uWqycRG9N8FCw2x8ZubZmGzKR58bQ9aaz9XKYQINBLmAaSzq5tVG0uYNcrOV6eOOO97M+1R3LN0Ah8equbtvaf9VKEQIlBJmAeQt/ae4mhtC3ddMt6llSq3zR/DxPQ4fvHX/XR2dfuhQiFEoJIwDyBPbylnbHI0S6eku/T+0BAb/3TZBI7UtPDazpM+rk4IEcgkzANEeU0zW8tquX7OKEJsrq8fXzolnamZ8fzPe4dkdi7EMCZhHiBeKT6BUnBt7sghjVNK8YNLcyivaeHNPad8VJ0QItBJmAcArTWv7TjBouwUMhKihjw+f3I6Y5Oj+cPmMh9UJ4SwAgnzAHDgdBNHa1u4fFqGW+NtNsXyBWMpPlrPp8fqvVydEMIKLLmdXylVANQaHzq01o94OsZ4vYcdWKO19ksybjCWFuZPTnP7Gsvysvivdw7y+81l/Oam2d4qTQhhEZYL855Q1lqvNz52KKVWa61XuDtGKXUvfcJbKbUaGPCa3rRh/2lmjbKTFh/p9jXiIsO4LnckL2w7xs9b2rFHh3uxQiFEoLNim2VFTygDaK1LgXwPx8ztZxZer5Sye1ztIE43trLreIPLyxHP5/o5o2jv7Ob1T2WZohDDjaXC3AjX3H5eqldK9RvoLo5x9DPe7o82y+bD1QAsnpjq8bWmjUxgSkY864qOe3wtIYS1WCrMAQfQX8DWGq+5O+Y+YINSaiWca8us9qxU13xUUkNidBiTR8R75XrXz8li94kG9lc0euV6QghrsFqYJ/H5Tcze6nHetHRrjNa6EMgDCpRSGijVWhcPVIRSqkApVQTkVVRUDKH8L9Jas6WkhgvHJ3/pdER3XTNrJGEhSmbnQgwzVgtzn1BKOYA5wDjgEZyz9IKB3q+1XqO1ngNsz8hwbzkhwNHaFk7Un+XC8SluX6OvpJhwvjIpjT/vOklXt5ymKMRwYcUwT+rnc4PdqBxszH1GQNdrre/DOUtfOVAf3ls2H64BYMH4ZK9e96oZmVQ1tbG1rL8fSIQQwchqYV5E/8GdBAzUFjnvGCOwN/R+wWixXA8sdb/UwRUdqSUlNgKHlx/QfOnkNKLCQvjLLlnVIsRwYakwN1aXlPazZNBu9L29MsZQBNS4X+3gdhyrJ3e03esPZo4OD+Urk9N4a88pOXxLiGHCUmFuWAmc62crpXKBwl4fO5RS6/qE94BjjEC/sZ+vUwCs8W7pn6trbqesupnZoxN9cv2rZ2RQ09zOllKf/nskhAgQltsBqrVeY6wmycfZPnH02f3pwLkhKAljSaILY75rLEssMT62A+t9uc58p3GGyuzRvtmXtHhiGjHhIfzl0wouyvF8DbsQIrBZLszBGc7nea0Q+NJ0d5Ax9TjXmvvNjqN12BTMyErwyfUjw0LIn5LOO/tO8R9d0wgNseIPYUIIV8nfcJPsOFbPpBHxRIf77t/Tr04dQV1LB9vL63z2NYQQgUHC3ATd3ZqdR+t91mLpcfGEVMJDbGzYJw98FiLYSZib4GhtC01tnUwf6ZsWS4/YiFAWZCfzzr7TaC0biIQIZhLmJthnnJsyNdO3YQ7OZ4QerW3h4OkzPv9aQgjzSJibYN/JRkJsipz0WJ9/rfzJzqN1N+yT54MKEcwkzE2wr6KR7NRYIsNCfP610uMjmTXKLn1zIYKchLkJ9lc0MiXTO0feumLplHQ+Pd7AqYZWv31NIYR/SZj7WW1zOxUNrUzJ8G+YA/ztQKXfvqYQwr8kzP2s56ER/pyZ56TFkpkQyfsS5kIELQlzP9t30hnmk/04M1dKccnENDYfrqG9Uw7eEiIYSZj72WenmkiLiyApJtyvX3fJxFTOtHVSVC5nnAsRjCTM/exwZRMT0uP8/nUXZKcQFqLYeKDK719bCOF7EuZ+pLXmcOUZstN8v768r9iIUOaOTeJ9CXMhgpKEuR9VNLTS3N7FeBPCHGDxxFQOnG7iZP1ZU76+EMJ3JMz96HClc0t9jmlhngbAxoMyOxci2EiY+9EhI8zNaLOALFEUIphJmPvR4cozJEaHkeznlSw9epYobjpULUsUhQgyEuZ+dLiyiey0WK8/wHkoFk9Mpbm9S5YoChFkJMz9RGvNocozZKf5f1libwvGJxNiU2w6VG1qHUII75Iw95Oa5nbqWzpM65f3iIsMY/YoO5sOS5gLEUwkzP2ktKoZgPGpMSZXAotyUth9ooH6lnazSxFCeImEuZ8cqXGG+bgU88P8opwUtIaPSmrMLkUI4SUS5n5SXtNMqE0x0h5ldinMzLITGxHKh9I3FyJoSJj7yZGaFrISowgNMf9bHhpiY74jmU2HZfOQEMEi1OwC3KGUKgB61tY5tNaPeGOMUupeoL7nfVrr9d6p2DkzH5Nsfoulx0U5KRTuPx1wdQkh3GO5MO8J5Z6gVUo5lFKrtdYrPBmjlNoAXK+1rjc+rlNKFfZ87AmtNUeqW8gbnejppbxmUU4KAB8eqpYwFyIImP8z/9Ct6D1j1lqXAvmejDFm5Ov6BHeeN4IcnMsSz7R1BlRoOlJiyEyIlPXmQgQJS4W5UsoO5PbzUr1Sqt9Ad3HMA0Bh7xeNwPeK8gBaydJDKcWinBQ+Kqmmq1ubXY4QwkOWCnPAgbOn3Vet8dqQxxhhbwdQSi1TSuUrpe41Pu8VR6pbABiTHO2tS3rFopxUGls72XXcKz+ACCFMZLUwT+Lzm5i91WMEshtj5vT8t9Z6vda6EFgDrBuoCKVUgVKqCMirqKgYtOjymmZsCrISAyvMF45PBmCz7AYVwvKsFua+YgfOtVWMXnmSUqq/9gxa6zVa6znA9oyMjEEvfqSmhZGJUYSHBta3Ozk2gqmZ8bLeXIggEFjp4pqkfj43WEvkfGNK4VyA91bL4DdWXVJe08zYALr52duinBSKj9bR3NZpdilCCA9YLcyL6D+4k4Bid8YMcqPTK83k8toWRicFVoulx0XZqXR0abaWyZG4QliZpcLcmD2X9nNz0m70ut0dU6yU6nsD1YHzHwKPNLV2UN/SEXD98h5zxiYSHmqTUxSFsDhLhblhJVDQ84HR1y7s9bFDKbWuT3ifdwxwn/Gr9+ulWuuBZvsuO2E8PDkr0fwzWfoTGRbC3LGJchNUCIuz3A5QrfUaYzVJPs72iaPP7k8Hzl53EkabZLAxWutCpZTd2DwEkKy1XuqNek/UBXaYAyzMTuGRtw5Q2dRKWlyk2eUIIdxguTAHZzif57VC4Ev75s83xnjda+ew9Hb8XJgHZpsFYFF2Co9wgC0lNVwza6TZ5Qgh3GDFNoulHK9rISLURkqsOQ9xdsXUzAQSosJka78QFiZh7mPH684yMjHK1Ic4DybEplgwPpnNh6vRWrb2C2FFEuY+dqL+bEC3WHoszE7hZEMrZdXNZpcihHCDhLmPHa87G9A3P3ssynYeiSurWoSwJglzH2pp76S2ud0SYT4mOZqR9ihZby6ERUmY+9AJC6xk6aGU4qKcFD4qqZEjcYWwIAlzH+pZlhgID3F2xcLsFJpaO9l9osHsUoQQQyRh7kPH65znmI+yQJsFYIEciSuEZUmY+9DxurOEh9pIiY0wuxSXJMdGMCUjXtabC2FBEuY+dLz+LCPtUdhsgbvGvK9FOSlsL6/jbHuX2aUIIYZAwtyHKurPkmm31lknC7NTaO/qZtsRORJXCCuRMPehUw2tjIi3Rr+8x9yxiYSH2KRvLoTFSJj7SFe35nRTGxkJ1pqZR4eHkjvGLuvNhbAYCXMfqT7TRle3ZoTFwhycu0H3nmyktrnd7FKEEC4a0hG4SqmxQC4wF+e54Ek4n5VZD2wDCrXWjd4t0ZoqGloBLDczB2ff/FfvHGTz4WqunplpdjlCCBe4FOZKqUuBFUANzmdtFuJ8EHItzkB3GL8eUUolAqu11u/5pGKLONXg3DBkxZn59JEJxEWGSpgLYSHnDXOlVALOx61t11rfMMDbGoAy4F3gcWPcpUqpHwFrhutM/fOZubVugAKEhti40JHMh4ecR+IG8vG9QginwXrmN2itfznUWbbW+l2t9a+ApUqpePfLs65TDa2Eh9pIjA4zuxS3LMpJ4UT9WY7WtphdihDCBecNc631455cXGv9ynCemWckRFp2VrvQOBJXVrUIYQ1eWc1itFWG5Qx8IM415tbrl/dwpMSQkRAp682FsIghh7lS6kdKqWt7h7fW+l2cLZVrvVqdhVU0nrXkSpYeSikWZsuRuEJYhTszcwX8BKhXSh1SSv2fEeIlOFe0DHvd3ZrTDW2MsODNz94WZadQ39LBvpPDslMmhKUMOcyNG6JztNY24C6cq1l+Amz3dnFWVdvSTntXt6Vn5gALsp1H4krfXIjA51HP3Fi1cr/Weg6QjXMN+rB3yliWaMU15r2lxUUyMT1O+uZCWMCQdoD2UErF912lorUuU0rN9k5Zg379ApwblgAcWutHvDlGKbVOa329u/VZefdnXwuzU3j2k3JaO7qIDAsxuxwhxADcuQG6CjiilKox+uVf6XUz1Oc9855Q1lqv11qvB9YrpVZ7a4xSKhdY5kmNVt792deinGTaO7vZXl5ndilCiPNwp82yXWudBCzF2S9fj/NmaBfOM1p8bYURyABorUuBfC+O8fgfpIqGVkJtipQYazxh6HzmjUsm1Kakby5EgHMnzEuVUndqrYuNfnkSkKi1DtFaP+HtAntTStlxHvTVV71Sqt9wHsoYpdSy3qHvrlMNraTFRVjqCUMDiY0IZfZou/TNhQhw7qxmeRdYp5T6Sq/P+etx7g76n/3XMvCM2qUxSikHzsPDPFbZ1EaahTcM9bUwO4XdJxqob5EjcYUIVAOG+fl2dGqtG1w9r8XLO0N7jtztqx7nkbyejMnVWntlNU5lk3NmHiwWZaegNWwpqTG7FCHEAAYMc611o1Lqu8YZ5kOmlBpntGMCfseJ0W4pHML7C5RSRUBeRUXFl16vamojLT54wnzmKDsx4SHSNxcigLly0NZSpdTDSqlZrlxQKTVbKfULYLaPeuhJ/XxuoFn5oGOMnjpaa5dv3mqt1xhr67dnZGR84bX2zm7qWjpIiwueNktYiI35jmTpmwsRwAZdZ661ftw41/wGpdRPAI2zbVHC562K8UAykABsAB72UR+9iP6DO4mBNywNNqYAzi1JPEcpdS9Qr7VeM5QCq860AQRVmwWcR+K++1klx2pbGJUUbXY5Qog+XNo0ZATz40BPsDtwhqGdzx9MUerrG6Fa63qlVKlSyt5nJm3XWvfbJnFhzJfGKaVWurIRqT+Vjc4NQ6nBFubGkbibD1dz07zRJlcjhOjLndUsDVrrHcZW/leM/93hxxUtKzFm03BuRl3Y62OHUmpdT/vElTHeVNXUMzMPnjYLQHZaLGlxEdI3FyJAeeU8c38y2h71Sql8pdQyIF9rvaLXWxw4NwQlDWEM4LwR2rMzVCm1eqC16+dT2RPmQXQDFJxH4i4yjsTtliNxhQg4gz0D9DqcM9oVWusjfqnIBefrYxutk8ShjOkzthDnw6vdUtnUhlKQHBPu7iUC1sLsFF7dcYL9pxqZmplgdjlCiF4G65kn4Qy3cy0LpdQsrfVOn1ZlYVVNrSTHhBMaYrkfega1sFffXMJc+FpTawenG9to7eiivasbrSEhKpT4yDCSYyMICYId1t40WJiPB7byxZ2RNwIS5gOoamojNcj65T1GJESSnRbLpsM1FFw83uxyRJDQWlNS1cyW0hr2nmhg78lGyqqbOdPWOeCY8BAbY5KjyU6LJW9MInPHJjE1Mz4oJ1GuOm+Ya63vN9aMP6GU0hizdKXUEpwHbgX8hiB/q2xqC7plib0tyk7hxW1HaevsIiJUjsQV7unu1mw9Usvrn55k44EqTtQ7Txq1R4cxNTOeZXlZZCREkh4fSUxEKGEhCg2cae2k/mwHx+taKK1qZs/JBv665xTgbG1eNnUEV83I4EJHclCcjTQUrqwzvx+4Xyk1DueNxftwLlN0KKXqcK7jLga2AcWB1Fs3Q2VjGxPT48wuw2cWZqfwh4+OUFxez4Xjk80uR1jM6cZWnv24nFeLT3Ci/izR4SEsyk7h7iXjuSg7lVFJUSg1tBA+1dDKJ2U1FO6v5E87T/DC1qOMTY7m1vljuGHuKOIjw3z0/yawuPxwCq11Gc515uONgMd4GMUcIA/no+NylVIlOGfwK4dbsHd3a6rPtAXdGvPeLnAkEWJTbD5cLWEuXPbZqUbWfFDKnz89SWe35uKcVO792kSWTkknOtytZ+ScMyIhkmtmjeSaWSNp7eji7b2neGZLOQ+9sZ9H3z3EikvGc/uCscREePZ1Ap07/+8e7vkPrfUOYAfOmTrgPJMF55Gz9yulinx9LG4gqW1pp7NbB3WbJT4yjJlZCWw6XM2PvjrR7HJEgDtS3cyvNxzk9U9PEh0ewrcuGMO3F45lTHKMT75eZFjIuWDffbyB3757kF++fYAnN5Xx469O5IY5o4K2/TLkMB9sc5Dx+LhaoA4/PHkokJzbMBREx9/2Z1F2Co/97TB1ze0kBuESTOG5+pZ2/uudg7yw9SihIYr/t2Q8373IgT3af39epmcl8MTyuew4WsfDb37G/a/uZt324zz0jWlMzvDmYa6BwVe3fufQZ0njcHBuw1AQz8wBFk9Ko1vDB4eqzC5FBBitNa9sP86l/7WR57ce5aZ5o/jgx0v48Vcn+TXIe5s9OpGXVsznl8tmUFbdzNcf28SaD0qCbvObT5pIxgMsht0aoZ5zWYJtK39fM7PsJEaH8f6BKq6ZNdLsckSAOFrTwr2vfMrHpbXMHm3nmW9MZ0pmYMyAlVJcP2cUl05O54FXd/Gfb37GBwer+fUNM4PmJ+lhF7i+1DMzD+YboAAhNsUlE1LZeLCKriCb3Yih01rz8rZjXP7bD9h7spGHr53OK3ctCJgg7y0pJpxVt+bx8LXT2V5ex9WPbeLTY/54dLHvSZh7UVVTG3ERoUSFB//66yWT0qhtbmfX8eD4iyDcU9vczopntnPvK7uYkWXn7R9ezM3zRgf0TUalFDfPG82rdy8gLMTG9au38GrxcbPL8piEuRc5d38G96y8x8U5qdgU/O2A9M2Hq53H6rny0Q95/0AVD145mefuvIBMe5TZZblsckY8r39vEbmj7dzz8qf8z7uH0Nq6P2lKmHtR9Zk2UoZJmCfGhDN7dCLvH6g0uxThZ1prnv24nOtXfUSITfHq3Qu48yJHQM/GB5IUE84zd1zAtbNH8l8bDvLvf9lv2Rujwb2K3s9qmtuZkB5rdhl+s2RiKr965+Cw+olkuGvt6OKf/7iHV4qPs3hiKr+5cZZpq1S8JSzExq+un0lCdBhPbi6j4WwHjyybYbmDvGRm7kU1Z9pIjhk+obZkUhoAGw9Kq2U4qG1u57a1n/BK8XF+cGkOTy6fa/kg72GzKX561RT+MX8CrxQf54FXd1luhi4zcy/p7HI+yDk5Njj+cLtiSkY8aXER/O1AJcvysswuR/hQWXUz3/79Vk42tPK7W3K5ckbG4IMsRinFD/Jz6NKaR989RHiojX+/ZtqQz4oxi4S5l9S2tAOQHDt8ZuZKKZZMTOPNPRV0dHUTNoyPHw1mW8tqKXimCJtSvPDd+eSN+dKzX4LKP+bn0NbZxeqNpcSEh/LAFZPNLskl8rfPS2rOOMM8ZZhtb18yKZWm1k6Ky+vMLkX4QOG+09y69hOSYsJ57e6FQR/k4Jyk3P+1Sdw2fwyrPyjlqY+OmF2SSyTMvaQnzIfTzBycR+KGhShZohiE/rTzBHc9u51JI+J45a4FjE6ONrskv1FK8fOvTyV/cjr/+ue9bNh32uySBiVh7iXVZ5y7P1OGUc8cIC4yjLljk3h3f+D/YReue/bjcn740k7yxiTy3J0XDMsD1UJsikdvnsW0kQn8wwvFAb9TVMLcS3rCfLjNzAGWTknnUOUZyqqbzS5FeMHqjSU8+NoelkxM46nvzCNumDzcoT/R4aGsXT6XlNgIVjyz/dzJqIFIwtxLaprbCQtRxEcOv3vKS6ekA7Bh3ymTKxGeWrWxhIf/+hlXz8xk9W15RIYF/9EUg0mNi2D1bXnUtbTzveeL6ezqNrukfkmYe0nPGnOrLGPypqzEaKZkxFuirygG9sSHpfzir5/x9ZmZ/ObGWbI6qZepmQn84rrpfFJWy8N//czscvolv1teUnOmfVitMe9r6ZR0isrrzrWbhLU8uamMh97Yz5XTM/j1DTMtt/vRH745O4vbF4xl7aYy/vzpSbPL+RJLhrlSqkAptcz4da83xhivFyilViqlViulhvRgjerm9mHZL+9x2dR0tIb39stZLVbz9JYj/Ntf9vG1qSP4zU2zCJUZ+YD++crJ5I6285NXd3OstsXscr7Acr9rSqkCoFZrvV5rvR5Yr5Ra7ckYpVSB1nqN8es+YAOwfSh11ZxpG3ZrzHubkhHPSHsU70jf3FJe2naUn/5pL0unpPPozbOltTKIsBAbv71pNgA/fGlnQPXPrfg7t8IIZAC01qVAvrtjlFIOYHzvNxvvTVJKLXO1qOHeZlFKsXRKOh8eqqalvdPscoQL3tpTwQOv7ubiCak8dstswkOtGAf+Nyopmoe+OY3t5XU8+t5hs8s5x1K/e0brI7efl+qVUv0GuotjCvp5vRZIcqWulvZOznZ0Des2C8BlU9Jp6+zmg4PVZpciBvFRSTXff2EnM0fZWXVrLhGhsmplKK6ZNZLrcrN47L1DbC2rNbscwGJhDjiA/lbu1xqvDXmM1rpUa93fHmUHUORKUed2fw7jNgvA3HFJJESFSaslwO050UDB09sZmxLN72+fS3T48FtO6w3/es1UshKj+fH6Tznb3mV2OZYL8yScIdxXPTDQDcshjzF67IVa62JXivp89+fwnpmHhdj4yqQ03vusMqB6ieJzpVVnWP7kVhKiwnj6OxcEzRG2ZoiNCGXldTMor2nhV+8cMLscy4W5zxk99BVa66XneU+BUqoIyKuoqPj8kK1hHuYAX506gvqWDj4uDYwfPcXnTje2ctvarQA8c8c8RiQEx1PpzXTh+GRumz+GJzeXsb3c3D/zVgzz/vrYgy0jHMqYlcCl57uYseplDrA9IyOj11Z+meUsnphKTHgIb+wOvHW4w9nZ9i7ueGob9S3tPPWdeThSh88TsXztvssnkZkQxY/X76K1w7x2i9XCvIj+QzgJGKgl4vIYpdRK4D6t9ZBO1Klpds7Mk4Z5zxwgMiyE/CnpvLXnFB3SagkI3d2af1q3k30nG3nsllymjUwwu6SgEhsRyi+um05pVTOPvnvItDosFeZGyJb2s6HHrrUu9GSM0SdfbSxb7PncYEseAWfPPC4iVM6xMFwxPYO6lg62lNSYXYoA/rvwIG/uPsVPrph87lF/wrsuyknlutwsHv+wlMOVZ0ypwVJhblhJr6WESqlcoHcoO5RS6/qE92Bj8oGiniBXStldDXKQNeZ9XTLBaLXsqjC7lGHvTztP8D/vHeamuaO4Y9E4s8sJag9cMYmosBB++qc9aO3/54daLsy11msw1ogbm3rytdYrer3FgXNDUJIrY4wbnhuA7UoprZTSQJ3xOdeWJja3Dfs15r1FhoWwdEo6b+2VVouZio/W8eP1u7hgXBL/ZqFnWVpVSmwE935tEh+V1PC6CWe3WHKBqRHOA71WCHxp3fhAY4zZuEd/ymvOtDMqafg8hcUVV87I5LWdJ9l8uJrFE+VHe387UX+Wgqe3k5EQyapb82R3p5/cPG80Lxcd46E39rNkUhrxfjwLXn6HvaCupZ0kWa/7BRflpImbBwQAABL+SURBVBAXESqtFhM0t3Vy51NFtHV0sXb5nGH5lCCzhNgUD31jGtVn2vjNBv/eDJUw94K65g75C9NHT6vl7b2naO+UVou/dHdrfvjSTg6cauSxb+WSnRZndknDzowsOzfNHc3TW45QWuW/m6ES5h7q6ta0d3WTFDN8H601kKtmZtDY2snGg/KwZ3955O0DbNh3mp9eNYVLJqSaXc6wdc/SCUSE2viFHx9kIWHuoc5u513rRGmzfMlFOakkx4Tzxx3HzS5lWFi//TirNpZw6/zRLF8w1uxyhrXUuAjuXpLNO/tO83Gpf5boSph7qOcMEtkw9GVhITaunplJ4b5KGlo6zC4nqG07UssDr+5iYXYyP7t6qqxcCQB3LBpHZkIkD72xj+5u3y9VlDD3UEeX8zdJDizq37W5I2nv6uaN3XIj1FeO1baw4pntjEqM5n9vyZMHTASIyLAQ7rt8EntONPLHHSd8/vXkd91Dnd0yMz+f6SMTyE6LlVaLjzS1dnDHU9vo7OrmieVzSIiWezeB5OoZmcwcZedX7xzw+bktEuYe6pmZy9LE/iml+ObskWw7UsfRmsB6ZqLVdXVrvv/CDkqqmvm/W/Pk8KwAZLMp7vvqRCoaWnn+k6O+/Vo+vfow0NnVTYhNERdpyf1XfvGN2SNRCr/8qDmcPPzmfv52oIp//fpUFmanmF2OGMCC7BQWZifzu78dprnNd49UlDD3UGe3JjE6DJtNbjgNZKQ9ivnjkvnjjuOmnFkRjF7cepQnNpVx+4Kx3Dp/jNnliEH86LKJ1DS384ePjvjsa0iYe6ijq1tufrpgWV4WR2paAuZ5iVa2paSGB1/bw8UTUnnwyslmlyNcMHt0IvmT01m1scRnK7skzD3U2a2lX+6CK6ZnEBcZygtbfds3DHZHqpv5++e2MzYlhsdumU2orFyxjH+6bAJn2jpZ82GJT64vfxI81NmlSZTdn4OKCg/h2tkjeXPPKeqMh3mIoWk461y5ArB2+Ry/HuIkPDc5I56rZ2Ty5KYj1BhPJ/MmCXMPdXR1y7JEF918wWjaO7t5VW6EDllnVzffe76Y8poWVt2ax5jkGLNLEm74/qU5tHZ2sXZTmdevLWHuIecNUAlzV0waEc/s0XZe2HpUboQO0UNv7OfDQ9X8xzenMd+RbHY5wk3ZabFcMT2Dp7eUe713LmHuIa21zMyH4OZ5ozlceYai8jqzS7GMZz4u5w8fHeHOReO4ce5os8sRHvrekmzOtHV6fWWLhLkXyGoW1101I4O4iFCe+7jc7FIsYdOhan7++l6WTEzlgStk5UowmJwRT/7kdJ7cXMYZL647lzD3Ajn+1nXR4aFcl5fFG7srqGxsNbucgFZSdYa7n9tOdmosj948mxDZyxA0vveVbBrOdvCsFyc1EuZeID3zobl9wVg6u7VX/yAHm/qWdu58qoiwEBtPLJ9DnKxcCSqzRtm5KCeFJz4s5Wy7d85skTD3AumZD83YlBgunZTOc58c9fnhQ1bU0dXN3c8Vc6LuLKtuy5Pnywapf/hKDtVn2nlpm3f2XkiYe4E8Mm7ovrNoLDXN7aY8xTyQaa352et7+aikhv+8djpzxyaZXZLwkXnjksgdbWft5jK6vHDeuYS5h5RSxEXIIVtDdaEjmUkj4nhyU5ksU+xl7aYynv/kKHddMp5leVlmlyN8rOBiB8dqz/L23lMeX0vC3EOhNiVPdXGDUorvLBrHZ6ea2HzYP4/VCnQb9p3mP97cz9emjuDer040uxzhB0unjGBMcjRrPij1eFIjYe6hsBAJcnddMyuT9PgIfve3w2aXYrq9Jxv4wYs7mJaZwH/fOEtO4RwmQmyKOxaNY+exerZ7uPfCkmGulCpQSi0zft3rjTHuXBMg1GbJb2FAiAgNoeDi8WwpraHoyPA9TfF0Yyt3/KGIhKgwnlg+h6jwELNLEn60LC+LhKgwHv+w1KPrWC6JlFIFQK3Wer3Wej2wXim12pMx7lyzR6jMzD1y87xRJMeE89gwnZ23tHdy51NFNLZ28MTyOaTHR5pdkvCz6PBQbps/hnf2naasutnt61guzIEVRuACoLUuBfI9HOPONQHkCFIPRYeH8p1F43j/QBW7jzeYXY5fdXdr7nnpU/acbODRm2YzNTPB7JKESf5uwRjCbDbWbnJ/dm6pJFJK2YHcfl6qV0r1G76DjXHnmr2FSW/TY3934RjiI0N5+p2P4feXQ9Nps0vyqaqWKm5/63Z+9ubHvLX3FP98xWTyp6SbXZYwUVpcJNfMyuS9TfsoueVWOquqhnwNS4U54ADq+/l8rfGaO2PcueY5KbERg71FDCIuMow7FjmYUboGXb4FNq40uySfWrVrFdtPF/PyobXcNn8MdywaZ3ZJIgAsXzCWb+55m7YdxVT97n+HPF5ZaY2vMVNerbUe3+fz64BtWutHhjoGKB7qNXu9532l1CVhYbLV2hP1P44gKvTLP+Gc7dTYf+n9Q/zNkv1YNrbwL8+futu7Ofy94XnPQDh9MmYsEf0spmjr7uaC8iPnPu7o6EBrvVFrvbjve2W3ixuMG6YFwES73c6KFSvMLsnS/rurmcvVR0xWpUTaumjrDmGfdvBm10LuuSd4trK3lrdSZN9Nrb0SW5jG1m0joymDyZWT+fo9Xze7PGGi11tbydu9m6yTFYR3d9Nus3E8M5Oi6dO4J/Lzm+LPPfccx44d6/caVgzz/vY32z0cM6Rraq3XAGuUUu/PmDHjkocffniQLy8G9ed/RBeX0KbDCLN1MjPvImZf/d9mV+VV+ysaKVz/I2yhlYTbwumgg4XzFvIv8//F7NJEAKj4+c+pf+ll2kNCCe3qInfRIq78+c++8J4tW7YMGOZW65kX0X/IJuFsl7gzxp1rCm9rrkTlfZuPv/Iyz3ZeyonjwXWi4rHaFm7//VZCws5w1bhref7K57lh4g3UnJXdr8Kps7oG+003se1Hv+LNcRfSeHJoCwEsNTPXWtcrpUqVUnatde+blnatdaG7Y4Z6TeEDNz0HwMVac9P+KH59uon3mtuD4kTKUw2t3PLEx7R2dPPSN1YxaUQ8AA/Of9DkykQgGfXY/wBwdWMrCw910bpgLEP5mc1qM3OAlTj71QAopXKBwl4fO5RS64wlhy6NceF14SdKKf79G9Noau3kl29/ZnY5Hqs+08a3nviYuuYOnv7OvHNBLsRA0uMj+dq0EbxcdIyWdtefRGS5MDf61T1rxJcB+Vrr3ncgHTg3/CS5OsaFawo/mpAex3cWjuXFbcfYcdS6zwqtb2nn1ic+4UT9WZ68fS4zRw12a0cIp+ULxtLU2smfdrp+RLSl2iw9jPAd6LVCIHEoY1x5XfjXD/In8JddFfx4/S7+8g+LiAyz1nklja0dLP/9Nkqrmll7+xzmjZNzyYXr5oxJJCctlhe3HuXmea49xNtyM3MxPMRGhLLyuhkcrjzDrzccNLucIalrbudbj3/CvpMN/O5buVyUk2p2ScJilFLcNG80nx5vYN/JRpfGSJiLgHXxhFRuuWA0j39YyjaLnKpYfaaNmx//mAOnm1hz2xyWyjZ94aZrZ48kPNTGiy4+Vk7CXAS0n1wxmazEKH744k7qmtvNLue8Tje2cuPqLRypaebJ5XNZMinN7JKEhSXGhHP5tBH8cccJlx76LGEuAlpsRCi/uyWXqqY2fvDSTq88K9EXDlc2ce3/fsSphlae+vY8FuWkmF2SCAI3zxtNU2snb+yuGPS9EuYi4M3IsvOzr0/hg4NVPPruIbPL+ZJtR2q57v+20NbZzYsFF3KBI9nskkSQuGBcEo6UGF7cOnirRcJcWMIt80ZzXW4Wv333EK8WHze7nHPe2FXBt574hOTYcP549wKmZ8mZ5MJ7lFLcOHcUReV1HDrddN73SpgLS1BK8Z/XTmPB+GTuXb+LDw4O/bxnb+rq1vzir5/x/54vZvrIBF65awGjkoLnUDAROK7LyyIsRPHitv7PZOkhYS4sIyI0hFW35ZGdFstdz27no5JqU+qobW5n+ZNbWbWxhFsuGM3z372AxCA4dkAEppTYCC6bMoJXi4/TfZ4jyyXMhaXER4bx9B3zyEqM4tu/38b7Byr9+vU/PFTF5b/9gK1ltay8bjr/+c3pRIRaa0OTsJ7r52RR19JBfUvHgO+RMBeWkxYXyYsFFzI+NZY7nyri2Y99f8Li2fYufv76Xm5bu5XYiFBevXsBN851bWeeEJ66KCeVEfGRVDUN/LAWCXNhSUkx4bxQMJ9FOSk8+NoeHnh1l0trcYdKa81beyrI//VG/vDREW5fMJY3vn8R00bKjU7hPyE2xbW5I6lvGXivhYS5sKyEqDDWLp/L3y8ezwtbj51rf3jL9vI6blu7lbueLSYuMpSXCubz869Ptdw5MSI4LMvL4ny7LCx50JYQPUJsivu+NomLclK4d/0ubli9hcumpHPPZRPcOm62q1vzwcEqntxcxoeHqkmKCeenV03h7y4cQ2iIzH2EeRypsYxPjWXfAItaJMxFUFgwPoW3f3gxazeVseaDUt75zYfMG5vEsrwsLp6QyoiEyAHHtnV2UVxez98OVPKXT09ysqGVlNhw7r98ErfNH0NMhPw1EYEhNS5iwNfkT6kIGjERoXz/0hxunT+Gl7Yd46VtR7n3lV0AjEqKwpESy4j4SCLDbLR3aeqa2zlS08zhyjN0dmvCQhQLxqfw4FVTyJ+cTniozMSFdUiYi6CTFBPO3y8ez12XONhf0cSmw1XsOt5AaVUz+ysaaevsJtSmSIoJJ9MexZJJacwaZWdhdgqxMgsXFiV/ckXQUkoxJTOeKZnyqDYR/OTnSCGECAIS5kIIEQQkzIUQIghImAshRBCQMBdCiCAgYS6EEEFAwlwIIYKAhLkQQgQBCXMhhAgCltsBqpQqAHrOOXVorR/xdIzxOsB4wA7cp7Wu91LJQgjhc5YK855Q1lqvNz52KKVWa61XuDtGKVWgtV7T6/3LgO04g10IISzBam2WFT2hDKC1LgXy3R2jlHLQJ7SN9yYZoS6EEJZgmZm5UsoO5PbzUr1SKl9rXTjUMUApUADc1+f1WiDJhbKyd+7cyeLFi114qxBCeGbnzp0A2f29ZpkwBxxAf33sWuO1IY8x/gFIHGBc0UCFGK2bAiCxoaGhZePGjdvOV7gYVAZQYXYRFiffQ+8I9O9jNlDV3wtWCvMkPr+J2Vs9zpuWXhljBHWh1rp4oEKMHvsa4/1FWuvFA5ctBiPfQ8/J99A7rPx9tFKY+5zRQ1+htc4zuxYhhBgKU8Lc6Ff37VMPZIVx0xL672MPNCvvMZQxK4FLXaxLCCEChilhbvSqv3TDchBF9B/CScBALRGXxyilVuLe+vI1g79FDEK+h56T76F3WPb7qLTWZtfgMqVUCZDXO3CVUiVa6wHXhLsyplefvLTX5/pdISOEEIHIauvMV+JcRQKAUiqXXjN8Y0PQOmNJoqtj8oGiniBXStmNzwkhhGVYamYO52bRpTjbJ1/Ymm+E8DqcM/HSwcYYNzxLBvhSibKlXwhhFZYL80Dgzvkw4ouMn55uAJZqra83ux6rknOFPNfrzyI4v4fjgZW9J4RWIEsTh8id82HEFxmtLgfn3/AlBiHnCnnNFxY/GD/hb6f/DYUBy2o980DgzvkwohetdbHxPbTUzCeQyLlCXjXH+NWjFLD3ufcW8CTMh8CFs16E8KeCfj7n6rlCwqC1zuuzcs0B1FutXSVtlqFx53wYIbzO+IlwyOcKCZfcB3zX7CKGSsJ8aNw5H0YIv3DlXCExMKM9tRTnzU/L7TGRMBciCMi5Qp7TWq9XShUCK5VS9t73xqxAwnzo3DkfRghfk3OFvMDok69QStUppUqt9FOO3AAdGnfOhxHCpzw4V2jYM3Z893cjuRS40d/1eELCfAiMvyyl/SxZsluxxyaszwii1X3PFTKxJKvJx/lTTV92oMbPtXhEwnzoznvWixgSWULnATlXyCsK6XMct3H/IQmLnaAo2/ndcL7zYcTgjL8sPSsH8oFHgJLeuxnF+cm5Qt7T688jOFem5WHB7fwS5kIIEQSkzSKEEEFAwlwIIYKAhLkQQgQBCXMhhAgCEuZCCBEEJMyFECIISJgLIUQQkDAXQoggIGEuhBBBQMJcCCGCgIS5EEIEAXk4hRB+YJxm6ADytNYren3+XgA5rE14SmbmQviYcf69wzgV8oY+5+HfiPMETiE8IqcmCuFjSql8rXWhcfb9u1rrxF6vaeTIWuEFMjMXwsd6PYXqRuDlns8brZdSCXLhDRLmQvhPAbCu18dLkadUCS+RMBfCD4w+uR3nQ8F75AMbjNflcW/CI9IzF8JPlFJaa62M/+557FsezlUuhdJuEZ6QMBfCT4xnx9pxrl6pB3KNl0q11utNK0wEBQlzIYQIAtIzF0KIICBhLoQQQUDCXAghgoCEuRBCBAEJcyGECAIS5kIIEQQkzIUQIghImAshRBCQMBdCiCDw/wH6a43kRAU9xgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dx=0.001;\n", "u=np.arange(2*dx,4.0,dx)\n", "\n", "# plot the results\n", "fig,ax=plt.subplots(figsize=(6,5))\n", "fig.tight_layout(pad=3.0)\n", "ax.plot(u,f(u,r,q))\n", "ax.plot(u,np.zeros(len(u)),color='black')\n", "ax.plot(root1,0,'*')\n", "ax.plot(root2,0,'*')\n", "ax.plot(root3,0,'*')\n", "ax.set_xlim(0,3.5)\n", "ax.set_ylim(-0.02,0.1)\n", "ax.set_xlabel(\"$u$\")\n", "ax.set_ylabel(\"$f(u)$\")\n", "\n", "# to save the figure\n", "# outputpath='../Lecture Notes/Chapter_SpatialModelling/figures/spruce_budworm_f.pdf'\n", "# plt.savefig(outputpath,bbox_inches=\"tight\", pad_inches=0.1);" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEeCAYAAAB40PUWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU5dk/8O+dPSEhC5AQNsOAiCwiCZvYCkrQurR1R0URt0S7+NpfW6n1fdu+7ftWoYvV+mKJO1hFwC62tRWCglbCkrCIhC0EQoCQQMIkgSRkmfv3x5xgCDNZJpOcOSffz3VxXc6ceYY7Y/LlyX2e8xxRVRARkbUFmV0AERF1HcOciMgGGOZERDbAMCcisgGGORGRDTDMiYhsgGFORGQDIWYX4AsRyQBQYTx0qOqiro4xjjeLA5Clqk5/1EtE1N0sF+bNoayqq4zHDhFZoqqZvo4RkSfRKrxFZAkAr+9JRBRIrNhmyWwOZQBQ1UIA6V0cM9nDLNwpInFdrpaIqAdYKsyNcE31cMgpIh4DvYNjHB7Gx7HNQkRWYbU2iwOAp4CtMI75OmYBgDUiskhVFxhtmSXtFSMi2wAMAFDQ3muJiPxgJIATqjqx9QGrhXkCvjyJ2ZIT7pOWPo1R1WwRSQOw1uifz1bVrd6KMMI+A8DYyMjI0ClTpgzuxNdARBZQcaYehSfPwOVSpPTvg8SYcLNLwvbt21FZWenxmNXCvFuIiAPAJADDATwF9yw9U1WzPL3eeD5LRNZNmTJlxrp163quWCLqVmcbm/CLv+fjrY2HcfXgWDx/1+VwDIg2uywAwMyZM7F+/XqPnQArhnmCh+faO1HZ3pgFLVbDLBCRd+GepReqarYvRRKR9Rxz1uKxP27FjmInHvnqcPzwutEIC7HGqUWrhXkuPAd3AgBvbZE2xxgnPte0PKCqW0XkDgCzATDMiXqBzwpO4rvvbEN9owt/uDcVXxuXbHZJnWKNf3IMxuqSQg9LBuO8zaB9GWPIBVDue7VEZAWqisXrCnDfq5vQPzoMf/3OlZYLcsBiYW5YCPfJRwCAiKSixezZuCBoZavw9jrGCPQ5Hv6eDAAee+ZEZA91DU144t3tWPSvvbjxskH487euxIgA6Y93ltXaLFDVLBHJMNojcXBfmt/ySk0H3BcEJcBYktiBMY+IyEIAB4zHcQBWcZ05kX2dqD6LjGW52HbYiSe/dgkemzECImJ2WT6zXJgD51aTeDuWDSC+k2OccK81J6JeYM/xKjz0Ri7Kz5y1ZH/cE0uGORGRr9buLsXj72xDdEQIVmZOx/ghsWaX5BcMcyLqNZblHMJP3t+FsYP64pV5kzEwNsLskvyGYU5Etqeq+M3qfXjx4wLMGp2I398zEVFh9oo/e301REStNDa58OM/78SK3CO4a/JQ/M/N4xASbMWFfG1jmBORbdXUN+I7b2/DR3vK8Pisi/G99IstvWKlLQxzIrKlijP1eOCNLdh5xIn/vWUc5k69yOySuhXDnIhsp7iiBve/thlHnbV46d40XDd2oNkldTuGORHZyoETpzH35U2oqW/EWw9PxeQUT/vs2Q/DnIhsI/9YFea9tgkA8G7mFbg0ua/JFfUc+53SJaJeadvhU7grKwehwUG9LsgBzsyJyAZyDpTj4Te3oF90OP748FQMTYgyu6QexzAnIkv7eG8ZHl2Wh2EJUXjr4alI6mufqzo7g2FORJb1z50leHz5NoxKisGyh6YioU+Y2SWZhj1zIrKk9/KO4Ntvb8VlQ+Lw9iPTenWQA5yZE5EFLcs5hP/66y5cObIfXp43yXb7rPiCnwARWcof1h/As//cg/RLE/HiPamICA02u6SAwDAnIktQVTy3Zh9e+KgAN12WjOfmXI5QG26Y5SuGOREFPFXFL/6+G699dhBzJg3FL28dj+Age26Y5SuGOREFtCaX4uk/78TyLcWYPz0FP7lpDIIY5BdgmBNRwGpocuH7K3bg/R3H8J2rR+L7146y7Ra2XcUwJ6KAVNfQhO++sw1r8kux4Guj8djMEWaXFNAY5kQUcGrqG5GxNA//LjiJn39zLOZdkWJ2SQGPYU5EAaWqrgEPvr4FWw+fwq/vmIDb04aYXZIlMMyJKGBUnKnHvNc2Ye/xarx4TypuGJ9sdkmWwTAnooBQVlWHua9swuGKGmTdNwlXj040uyRLYZgTkemOnKrB3Fc24WT1WbzxwBRcMaKf2SVZDsOciExVeOI05r6yCWfOum/zNnFYvNklWRLDnIhMs7ukCve9ugmqwPKMKzBmUO+6O5A/McyJyBTbi524/7XNiAoLxlsPT8WIAdFml2RpDHMi6nEbCk7ikaW5vfo2b/7GLceIqEd9sLME81/fgiHxUVj56BUMcj+x5MxcRDIAVBgPHaq6yB9jRORJAM7m16nqKv9UTEQAsGxjEX7y1y+QNiwer94/GbFRoWaXZBuWC/PmUG4OWhFxiMgSVc3syhgRWQPgDlV1Go9PiUh282Mi8p2q4vm1+/G77P2YNdp9U4nIMN5Uwp8sF+YAMlU1rfmBqhaKSHpXxhgz8pWtgjuNQU7UdU0uxc/e34VlG4twe9oQPHvreITwphJ+Z6lPVETiAKR6OOT0FugdHPMUgOyWB1W1sCu1EhFwtrEJj7+zDcs2FiFzhgO/uv0yBnk3sdrM3AF3T7u1CuNYp8cYYR8HACJyu/HaVABZnJkT+a6ytgGPLstDTmE5nr7hUjxylbcfUfIHq4V5Ar48idmSE0Yg+zBmUvN/t+ip5wJYCWB2Vwsm6o2KK2rwwBtbUFR+Bs/NmYBbJnLnw+7G33fc4gCca6sYM/IEEfHUnoGIZBiBn1ZSUtJDJRJZw45iJ25ZvAFlVXVY+uBUBnkPsWKYJ3h4ztusvCNjCoFzAd5SBQCPfXhVzVLVSQDykpO5RSdRs9W7jmNOVg4iQoPwp29N54ZZPchqbZZceA7uBABbfRljrGzx9vexZ07UQa9/dhA//3s+LhsSh1fmTcKAmHCzS+pVLBXmquoUkUIRiWs1k45T1ewujNkqIo5WK1gccP9DQERtqG904Wd/24W3Nx3GdWOT8Ls5E7mG3ARWbLMsBJDR/MDoa2e3eOwQkZXGKpUOjQGwwPjT8nihqnqb7RMRgJOnz2LuKxvx9qbDeGzmCCyem8YgN4mlZuaAu19tnIBMh7t94mh19acD7l53Aow2SXtjVDVbROKMi4cAoJ+qciULURu+OFqJjKW5qKipxwt3T8Q3Jgwyu6RezXJhDrjDuY1j2QAu2N2+rTHGce7DQtRB7+84hidX7UBCVBhWPTod4wbHml1Sr2fJMCcic9Q3uvDsP/fgtc8OYnJKPF66Nw39o3miMxAwzImoQ46cqsF33t6G7cVOzJ+egh/fcCnCQqx42s2eGOZE1K61u0vx/1bsgMulWDw3FTeM5/UVgYZhTkRe1Te68JvVe7Hkk0KMSe6LxXNTkdK/j9llkQcMcyLyaF9pNZ5Yvh35JVW4Z+ow/OSmMYgI5bLDQMUwJ6LzuFyK1zccwsJ/7UFMeAhenjcJs8ckmV0WtYNhTkTnFFfU4Ed/+hyfFZQj/dJEPHPrZbws3yIY5tSjTp4+i51HKpFfUoWCstM4XlmHE6fPoqq2AU0uRaNLER4ShJiIEPSNDMXAvhEYmhCFoQlRGDGgD8YOikVsJO8b6W+NTS68/tkh/HbNPogAz9w6HndNHoo29i2iAMMwp26361gl/rajBJ/sO4H8kqpzzw+Oi0RybAQuToxGbGQogoMEwUGCsw0uVJ9tQGVtA/aWVmPtnjLUN7rOjRuaEIlxg2IxOSUB0xz9MHpgDIKCGDq++uJoJZ76007sPFqJWaMT8fObx2FwXKTZZVEnMcypW9TWN2FlXjHe3nQYe45XIzRYkHZRPH543SWYnJKA0ckx6BvRsRm2y6Uoqz6LfaXV+OJYJXYdrcKOI07884vjAIC4qFBMHZ6Aq0YNwKzRSRgYG9GdX5ptOGvq8bvs/Vi2sQjxUWF48Z6JuHF8MmfjFsUwJ7+qrW/C6xsO4tVPD6L8TD0mDInFL24eh69floy4qDCf3jMoSDAwNgIDYyNw1agB554/6qzFpsJybCwsx2cF5fhwVymexhcYO6gvZo1OxKxLkzB+cCxn7a00NLnwx41FeC57P6rrGnD3lGF48rrRiI1i+8rKGObkF6qKf+wswTMf7MFRZy1mXjIA35o5ElOGe7oviH8MjovEralDcGvqEKgqCspOI3t3GT7aU4oXPy7ACx8VIDEmHOljknDtmCRcMaIfwkN679I6VcXq/FL86sO9KCg7jStH9sN/3TQGowf2Nbs08gOGOXVZWVUdFrz3OT7eewKXJvfFb++cgKmOnr3DjIjg4qQYXJwUg8dmjsCpM/VYt68Ma/JL8ddtR/H2psOIDg/BjEsG4NoxSZh5SWKvOZGqqvhoTxmey96HL45WwdG/D16eNwnplyaypWIjDHPqkn99cRxP/elz1NQ34Sc3jcH901MQHABtjfg+Ybhl4hDcMnEI6hqakHOgHKvzS7EmvxT/+LwEIUGCaY5+mD0mCbPHJGGQDU/4NbkUq3cdxx8+KcSOYieGJUTh13dMwM2XD0JIMPdUsRuGOfnE5VL8ds0+vPhxAcYPjsVzcy7HyMRos8vyKCI0GFePTsTVoxPxvzePw/YjTqzeVYo1+cfx0/d34afv78K4wX1x7ZiBuHZsEi5JirH0jLW6rgErco/gjQ0HUVxRi6EJkXj21vG4LW0IQhnitsUwp06rqW/E4+9sQ/buMsyZNBQ/v3msZXrRQUGC1GHxSB0Wjx9dPxoHTpzGmvxSrN51HM9l78Nv1+zD0IRIXDtmIGaPScKki+ItMYtVVeQWncLK3GL84/MSnKlvwuSUeDx9wxjMHpMUEL8tUfdimFOnVNY04IE3NmN7sRP//Y2xmHfFRZaexY4YEI0RM6Lx6IwRKKuuw9rdZVi96ziWbSzCq/8+iPioUFwz2t2KuWpUf0SFBc6PjKoiv6QKq3eV4v0dx3Dw5Bn0CQvGjZclY+7UizBhqKf7mJNdBc53JgW8suo6zHt1MwpPnMHiuan42jh7bYOaGBOBu6cMw91ThuH02UZ8su8E1uSXInt3Kd7begThIUGYnJKAtIviMTklAZcPi0N0eM/+CFXXNWDLoQr8e3851uw+juKKWogAU4cn4NtXj8T14waiTw/XRIGB/9epQ06dqce9r2xCcUUtXps/GV+5uL/ZJXWr6PAQ3DA+GTeMT0ZDkwtbDlZgdX4pNh2swAsf7YcqECTA6IF9MWZQX4weGIPRA/vikoEx6B8d5pffVuobXSgoO438kirkH6tCXlEFdh6thEuBsOAgXDmyH749cyTSxyTxbj/EMKf2nT7biPmvb8ah8hq88cBkTB9h7yBvLTQ4CNNH9sf0ke6vu7quAdsOO5FbdArbDp/Cur0nsCrvyLnX9wkLxuD4SAyJj8LguEjE9wlD34gQ9I0IRUxEyHkXMTU2KarrGlBd14iqugYcr6zDUWctjpyqxTFnLRpdCgCICA3C+MGx+M7VIzFtRD+kDovndrR0HoY5tam+0YVH3szFF8eqsOTetF4X5J7ERITiqlEDzrsa9eTps9h7vBp7j1fjcEUNjjprcfRULbYePgVnTUOH3lcESIwJx+C4SFw+NA5fn5B8buaf0q8PT2JSmxjm5JWq4j//shM5heV4bs4EpHNPa6/6R4ej/8hwXDnywn/smlyK02cbz83AXarnjgUHybkZe5+wEG49QD5jmJNXr/77IFbkHsHj14zELROHmF2OZQUHCWIjQ3vNFadkjsBfQEumWL/vBH75wW5cP24gnkgfZXY5RNQOhjldoKSyFk8s34ZRSTH4zZ0T+Ks/kQUwzOk8jU0u/Mc723G20YXFc1MD6iIZIvKOP6l0nhfW7sfmQxV4bs4EOAYE5l4rRHQhzszpnM0HK/D7jwtwR9oQnvAkshiGOQFwb571w1U7MCQ+Ej/7xlizyyGiTmKbhQAAi/61F0XlNXjnkWnc24PIgjgzJ2wqLMcbGw7h/isuwhUjevYOQUTkHwzzXq6uoQkL3vscwxKisOD60WaXQ0Q+suTv0yKSAaDCeOhQ1UX+HCMiK1X1jq5XGviyPinEofIaLHtoCpchElmY5X56m0NZVVcZjx0iskRVM/0xRkRSAdzeTeUHlOKKGvzfxwW4cXwyvnrxgPYHEFHAslyYA8hU1bTmB6paKCLpfhzj8EeRVvCz93chOEjwnzddanYpRNRFluqZi0gcgFQPh5zewrkzY0Tk9ubZu92tyS/F2j1l+F76KCTH2u/O9ES9jaXCHO5Zs9PD8xXwPqPu0BgRcQAo7GqBVlDf6ML//CMfFydGY/6VKWaXQ0R+YLUwT8CXJzFbcgLwdvfajo5JVdWtHSlCRDJEJBdAWklJSUeGBJS3NhahqLwGT994KUItcOd5Imoff5IBGO2W7I6+XlWzVHUSgLzkZGvd1LiytgEvfLQfXxnZHzNG8aQnkV1YMcwTPDznbVbe7hijpw5V9dSKsZ3FHxegsrYBT90w2i83HSaiwGC11Sy58BzcCQC8tUjaG5MBnFuSeI6IPAnAqapZPlcbYI6cqsHrGw7h1olDMHZQrNnlEJEfWSrMVdUpIoUiEtdqJh2nqh7bJB0Yc8E4EVnYkQuRrOa3q/dBAPzgOt45iMhurNhmWQhjNg2cm1Fnt3jsEJGVze2TjozpDQrKqvGX7Ucxf3oKlyIS2VCnZuYikgL3mu3JcLcumleKOAFsAZCtqlX+LfF8qpplrCZJN2pwtLqS0wEg3ajN2cExAM6dCL3D+O8lAFZ6m/FbzfNrCxAZGozMGSPMLoWIukGHwlxEZgHIBFAOd585G+412RVwh6bD+LNIROIBLFHVj7qlYrjDuY1j2QDiOzOm1dhsuL9W29h7vBp///wYHpsxAgl9wswuh4i6QZthLiKxcLcn8lT1Ti8vqwRwEMBaAC8b42aJyA8AZHX3TJ3a9/zafegTFoJHvtprdiog6nXam5nfqaq/6uybqupaAGtF5DYRWcNAN0/+sSp8sPM4Hr9mJOI5KyeyrTbDXFVf7sqbq+p7XRlPXff82n2IiQjBQ1/hrJzIzvyymsVoq/T1x3uR/+w9Xo0Pd5XiwSuHIzYq1OxyiKgbdTrMReQHInJry/A22iqzReRWv1ZHXbJk/QFEhQVj/vQUs0shom7my8xcAPwY7i1k94vIS0aIH0Av2gs80B05VYO/7jiGu6cMY6+cqBfodJir6q9UdZKqBgF4FO7VLD8GkOfv4sh3r3x6EEECPPzV4WaXQkQ9oEuX8zevWgEAERkOgMkRAMpPn8XyLYdx8+WDebUnUS/h0wlQTyc7VfUg2t+9kHrAmxsO4WyjC5kz2PUi6i18OQH6BwCHRKTc6Jdf0yLcmR4mO322EW/mFOHaMUkYmRhjdjlE1EN8mZnnqWoCgNlw98tXwX0ytAmeb89GPWj55sOorG3Ao9yDhahX8aVnXigiD6vqK3Dv0/IjEYlV1Uo/10ad5HIpluYUYXJKPCYOu2B7GiKyMV9Ws6wFsFJErmnxHIM8AKzfdwKHK2ow74oUs0shoh7mNczbuqJTVSs7uisirwztOW/mHEJiTDiuGzvQ7FKIqId5DXNVrRKRR4w9zDtNRIYb7RhustUDDp08g3V7T+CeqcMQFmLFe44QUVe0+VNvbLQ1W0SeEZHLO/KGIjJRRJ4FMNHoq1MPeGtjEUKCBPdMGWZ2KURkgnZPgKrqy8a+5neKyI8BKNw3pTgA9+qVOAAjAPQDEAtgDYBn2EfvObX1TViRW4yvjRuIxL4RZpdDRCbo0GoWI5hfBtAc7A647zAUhy9vTFHIADfHX7cfRVVdI098EvVinV6aaAT2tm6ohXygqngzpwijB8ZgcgqXIxL1Vm32zEUkhatRAltu0SnsLqnC/dNTICJml0NEJmlvZl4IYI2IbIW7R76Cq1MCy9KcIsREhOCblw8yuxQiMlF7Yb5VVa/rkUqo08qq6vDPnSW4f3oKosK6tAEmEVlcewuSc3ukCvLJ25sPo9GluHfaRWaXQkQma286d6rlAxGZBWAi3Huy5LLlYp6GJhfe3nQYM0YNwPD+fcwuh4hM1l6Yn3dGrflmFCJSACBPRNbwwiBzfLjrOMqqz+LZ2zgrJ6L2w1y9PL8KvDDIVEtzijA0IRIzRiWaXQoRBYD2euaZIrJYRG5ptURRGeTm2V1Shc0HK3DftIsQHMTliETUsS1wrwXwHoBTxt2F3gWQLiITWr9QRJ7xd4F0oaU5RQgPCcKdk4aaXQoRBYj2wjxLVUcCiAdwHYBX4N6HJQ3AVhFpEpEPReT7xkZcvG1cN6usbcBfth3FNy8fhLioMLPLIaIA0V7P/Bng3CX82cYfAICIpAJIN/48DeBX8N5jJz9ZlXcEtQ1N3IeFiM7TZpi31RdX1a1wL1FcBAAi4gCwwq/V0XlcLsWynENIHRaHcYNjzS6HiAKI3y4bVNVCEclu/5VdJyIZcG/DCwAOVV3U1THGccDdRooDsEBVA+oG1Z8WnMSh8hp8b/Yos0shogDj12vAVfVH/nw/T5pDWVVXGY8dIrJEVTN9HSMiGaqa1eL1twPIgzvYA8bSDYfQPzoc149LNrsUIgowVry/WGZzKAPu3wjg7tv7NMZoD50X2sZrE4xQDwjFFTX4aG8Z7p4ylLeFI6ILWCoVRCQOQKqHQ04R8RjoHRyT4eF4Bdw34AgIb20sQpAI7pnK28IR0YUsFeZwL3301MeugPdlkW2OUdVCVfV0VwcHAmSjsbqGJrybW4xrxyQhOTbS7HKIKABZLcwT8OVJzJaa70XqlzFGjz3bWLFjuvd3HIOzpoHLEYnIK6uFebczeuiZqjq7jddkiEgugLSSkpJurUdVsTTnEEYlRWOaI2C6PkQUYKwY5p4Szdus3JcxCwHMauvNVDVLVScByEtO7t6VJduKnfjiaBXuu4K3hSMi76wW5rnwHMIJcF/A1KUxIrIQAba+fOmGQ4gJD8GtEwebXQoRBTBLhbkRsoXGCpWW4lTV4wVLHR1j9MmXGMsWm59rb8ljtzpRfRb/2FmC29KGoE84bwtHRN5ZKswNC9FiKaGxR0zLUHaIyMpW4d3emHS475xUaDyOMzvIAeDdLYfR0MTbwhFR+yw33VPVLOMEZDrc7RNHq6s/HXBfEJQAY0liW2OME55rjP9u/dd5WrLYIxqbXHhr42F89eL+GJkYbVYZRGQRlgtzwB3ObRzLhocQ9jbGmI0H3JnFNfmlOF5Vh59/c6zZpRCRBVixzdIrLM0pwuC4SMy6NMnsUojIAhjmAWhfaTVyCssxd9ow3haOiDqEYR6AluYcQlhIEObwtnBE1EEM8wBTWduA9/KO4hsTBqFfdLjZ5RCRRTDMA8zK3GLUNjRh/vQUs0shIgthmAeQJpdiaU4RJl0Uz9vCEVGnMMwDyLq9ZThcUYP7OSsnok5imAeQNzYcQlLfcHxt3ECzSyEii2GYB4iCsmp8uv8k7p16EUKD+b+FiDqHqREg3txQhLDgINzN28IRkQ8Y5gGgqq4B7209gpsmJKM/lyMSkQ8Y5gHg3c3FqKnnckQi8h3D3GQNTS689tlBTB2egMuGtHfDJCIizxjmJvvbjmMoqaxD5gyH2aUQkYUxzE2kqsj6pBAXJ0Zj5qhEs8shIgtjmJvok/0nsed4NR65yoEg7o5IRF3AMDdR1icHkBgTjm9ePsjsUojI4hjmJvniaCU+KyjHA1cOR3hIsNnlEJHFMcxN8vza/egbEYK503iREBF1HcPcBLuOVWJNfike+ooDfSNCzS6HiGyAYW6CF9buR0xECOZfmWJ2KURkEwzzHpZ/rAof7irFQ18ZjthIzsqJyD8Y5j3suex9iIkIwQNXDje7FCKyEYZ5D9p8sAJr8kuReZWDs3Ii8iuGeQ9RVfzyg91I6huOh77CS/eJyL8Y5j3kHztLsL3Yie/PvgSRYVxXTkT+xTDvAbX1TXj2n3swemAMbksbYnY5RGRDIWYX0Bu88NF+HDlVi+UZ0xDMPViIqBtwZt7N9pVW4+VPCnFb6hBMc/QzuxwisimGeTdqaHLhyVWfIzoiBE/feKnZ5RCRjbHN0o1+v3Y/thc78fu7JyKhT5jZ5RCRjXFm3k1yDpTjxY8LcFvqEHx9Are4JaLuZcmZuYhkAKgwHjpUdVFXx/jynt4cOnkGj/0xDyn9++C/vznW17chIuowy4V5c+iq6irjsUNElqhqpq9jfHlPb0qr6vDgG1sAAK/dPxnR4Zb7iInIgqzYZslsDl0AUNVCAOldHOPLe17gcHkN5izJQWlVHV6ZNwkp/ft09i16t+rjwOvXA9WlZlfSrU7UnMD8f83HydqTZpdCAaahrAyH7r0PjSdOdHqspcJcROIApHo45BQRj+Hb3hhf3rM1l0vxl21HcePvP0XFmXosfWgqJqUkdGQotbR+EXB4I7B+odmVdKs/fP4HbC3dipd2vGR2KRRgTi5+CbV5eTjxf4s7PdZqPQAHAKeH5yuMY76MqfDhPc/ZVVCEsXN+iLLqsxgSH4m7pwxD3odFyGtvIJ3z4LGnEILGL5/IfRXIfRWNCMFrg54xrzA/Wxy2GE3SdO7xir0rsGLvCgRrML5V/y0TKyOzTX/+BQQ3ffm94Vy+HM7ly9EUHIwN//H4uedLSkq8vofVwjwBX56kbMkJIM7HMZ1+T6PHngHgkpNHD+Hkqt8Axpt8vryN6smjn0YLfj07HDePDkWfMMGZesWf9zTgB6urUXqm06ctAlZIbAgG3jUQfVP7Iig8CK6zLlTmVaL03VJkVtrn66TO6x8cjCcTEzErOgaRQUGodbmQXV2NX50ow8nMjn1vWC3MA4KqZgHIEpF1U6dNm/Gn994zuyTLi/30Z4javQIaHIaosHrcPGcernnxp2aX5Xcv5L+AD45+gNCgUDSEN+C+O+/Dd3/6XbPLogBQ+9zvUP/3vwOhoYhsaMBt992Le5944rzX3H777cjJyfE43oph7qkZ7W1W3tExvrwnACAiPByDBnEdeZfpGWDSg8CkB4Dc1xF9uhTRNvxc6xbLghMAAAYdSURBVPbW4c5L7sQdo+7Ayn0rcbL2JL9/CABQXFuLqLvuQvycO3Hq3RVoPHHigu+NsDDvFx+KqnZ3jX5jnKw8parS6vk8AAtUNbuzYwDkdvY9W7xm3YwZM2asW7fO1y+JiKjDZs6cifXr169X1Zmtj1lqNYuqOgEUGgHdUpy30G1vjC/vSUQUaCwV5oaFcJ98BACISCqA7BaPHSKyslU4tzmmA8eJiAKa5XrmqpolIhnGGvA4uC+9b3m61wH3BT8JMJYctjemA+9JRBTQLBfmwLnVJN6OZQOI78yYjhwnIgpkVmyzEBFRKwxzIiIbYJgTEdkAw5yIyAYY5kRENsAwJyKyAYY5EZENMMyJiGyAYU5EZAMMcyIiG2CYExHZAMOciMgGGOZERDbAMCcisgGGORGRDTDMiYhsgGFORGQDDHMiIhtgmBMR2QDDnIjIBhjmREQ2wDAnIrIBhjkRkQ0wzImIbIBhTkRkAwxzIiIbYJgTEdkAw5yIyAYY5kRENsAwJyKyAYY5EZENMMyJiGyAYU5EZAMhZhfQWSKSAaDCeOhQ1UVdHWMcB4ARAOIALFBVp59KJiLqdpYK8+ZQVtVVxmOHiCxR1Uxfx4hIhqpmtXj97QDy4A52IiJLsFqbJbM5lAFAVQsBpPs6RkQcaBXaxmsTjFAnIrIEy8zMRSQOQKqHQ04RSVfV7M6OAVAIIAPAglbHKwAkdKCskdu3b8fMmTM78FIioq7Zvn07AIz0dMwyYQ7AAcBTH7vCONbpMcY/APFexuV6K8Ro3WQAiK+srKxZv379lrYKp3YlAygxuwiL42foH4H+OY4EcMLTASuFeQK+PInZkhPuk5Z+GWMEdbaqbvVWiNFjzzJen6uqM72XTe3hZ9h1/Az9w8qfo5XCvNsZPfRMVU0zuxYios4wJcyNfnXrPrU3mcZJS8BzH9vbrLxZZ8YsBDCrg3UREQUMU8Lc6FVfcMKyHbnwHMIJALy1RDo8RkQWwrf15Vntv4Tawc+w6/gZ+odlP0dRVbNr6DAROQAgrWXgisgBVfW6JrwjY1r0yQtbPOdxhQwRUSCy2jrzhXCvIgEAiEgqWszwjQuCVhpLEjs6Jh1AbnOQi0ic8RwRkWVYamYOnJtFF8LdPjnv0nwjhFfCPRMvbG+MccLzgJe/Kp6X9BORVVguzAOBL/vD0PmM357uBDBbVe8wux6r4r5CXdfiexFwf4YjACxsOSG0Ai5N7CRf9oeh8xmtLgfavuCL2sF9hfzmvMUPxm/4efB8QWHAslrPPBD4sj8MtaCqW43P0FIzn0DCfYX8apLxp1khgLhW594CHsO8Ezqw1wtRT8rw8FxH9xUig6qmtVq55gDgtFq7im2WzvFlfxgivzN+I+z0vkLUIQsAPGJ2EZ3FMO8cX/aHIeoRHdlXiLwz2lOz4T75ablrTBjmRDbAfYW6TlVXiUg2gIUiEtfy3JgVMMw7z5f9YYi6G/cV8gOjT54pIqdEpNBKv+XwBGjn+LI/DFG36sK+Qr2eccW3pxPJhQDm9HQ9XcEw7wTjh6XQw5KlOCv22Mj6jCBa0npfIRNLspp0uH+raS0OQHkP19IlDPPOa3OvF+oULqHrAu4r5BfZaLUdt3H+IQEW20GRl/P7oK39Yah9xg9L88qBdACLABxoeTUjtY37CvlPi+9HwL0yLQ0WvJyfYU5EZANssxAR2QDDnIjIBhjmREQ2wDAnIrIBhjkRkQ0wzImIbIBhTkRkAwxzIiIbYJgTEdkAw5yIyAYY5kRENsCbUxD1AGM3QweANFXNbPH8kwDAzdqoqzgzJ+pmxv73DmNXyDtb7Yc/B+4dOIm6hLsmEnUzEUlX1Wxj7/u1qhrf4piCW9aSH3BmTtTNWtyFag6AFc3PG62XQgY5+QPDnKjnZABY2eLxbPAuVeQnDHOiHmD0yePgvil4s3QAa4zjvN0bdQl75kQ9RERUVcX47+bbvqXBvcolm+0W6gqGOVEPMe4dGwf36hUngFTjUKGqrjKtMLIFhjkRkQ2wZ05EZAMMcyIiG2CYExHZAMOciMgGGOZERDbAMCcisgGGORGRDTDMiYhsgGFORGQD/x/ap66FjzyyjAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "F=np.zeros(len(u))\n", "err=np.zeros(len(u))\n", "for i in range(0,len(u)):\n", " F[i],err[i]=integrate.quad(f,0,u[i],args=(r,q));\n", " \n", "# plot the results\n", "fig,ax=plt.subplots(figsize=(6,5))\n", "fig.tight_layout(pad=3.0)\n", "ax.plot(u,F);\n", "ax.plot(u,np.zeros(len(u)),color='black')\n", "ax.plot(root1,0,'*')\n", "ax.plot(root2,0,'*')\n", "ax.plot(root3,0,'*')\n", "ax.set_xlim(0,3.5)\n", "ax.set_ylim(-0.02,0.08)\n", "ax.set_xlabel(\"$u$\")\n", "ax.set_ylabel(\"$F(u)$\")\n", "\n", "# to save the figure\n", "# outputpath='../Lecture Notes/Chapter_SpatialModelling/figures/spruce_budworm_intf.pdf'\n", "# plt.savefig(outputpath,bbox_inches=\"tight\", pad_inches=0.1);" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/baker/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in sqrt\n", " import sys\n" ] } ], "source": [ "# Evaluate the integral using simple quadrature\n", "# Note in this calculation the solution won't be defined anywhere that F(u)>F(u_{max}), hence the RunTimeWarning\n", "suminteg=np.zeros(len(u))\n", "for i in range(0,len(u)):\n", " integ=np.zeros(len(u))\n", " for j in range(0,i):\n", " integ[j]=dx/np.sqrt(F[i]-F[j]);\n", " suminteg[i]=np.sqrt(2/D)*np.sum(integ);" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.053197938273974885\n", "0.05319839977033716\n" ] } ], "source": [ "# Determine (semi-manually) the regions where the solution is defined\n", "x1=np.max(np.select([u" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot the results\n", "fig,ax=plt.subplots(figsize=(6,5))\n", "fig.tight_layout(pad=3.0)\n", "ax.plot(u,suminteg)\n", "ax.plot(root1*np.ones(101),np.arange(0,202,2),linestyle='--');\n", "ax.plot(2220*dx*np.ones(101),np.arange(0,202,2),linestyle='--');\n", "ax.plot(root3*np.ones(101),np.arange(0,202,2),linestyle='--');\n", "\n", "ax.set_xlim(0,3.5)\n", "ax.set_ylim(0,200)\n", "ax.set_xlabel(\"$u_{\\mbox{max}}$\")\n", "ax.set_ylabel(\"$L$\")\n", "ax.set_yticks(np.arange(0, 250, 50));\n", "\n", "# to save the figure\n", "# outputpath='../Lecture Notes/Chapter_SpatialModelling/figures/spruce_budworm_L.pdf'\n", "# plt.savefig(outputpath,bbox_inches=\"tight\", pad_inches=0.1);" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }