Skip to content

Instantly share code, notes, and snippets.

@terakun
Created July 19, 2019 06:32
Show Gist options
  • Save terakun/a5075d6f2281fda33fbaae9de29fe169 to your computer and use it in GitHub Desktop.
Save terakun/a5075d6f2281fda33fbaae9de29fe169 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
{
"cells": [
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import rc\n",
"rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})\n",
"rc('text', usetex=True)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/meip-users/anaconda3/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
" warnings.warn(message, mplDeprecation, stacklevel=1)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAETCAYAAAAoO4PeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXd4VNXWh989Jb1MCoQmJfQmEBKaioqhqJ9XvUZQURBpKnKxASI2wCvtKlaQUFSwE2yICIQqSAtBQEBFQkdawlDSy/n+ODMhPZNk2hn2+zzznLbP2SvJnF92WXstoSgKEolEUlV0rjZAIpFoEykeEomkWkjxkEgk1UKKh0QiqRZSPCQSSbWQ4iGRSKqFFA8PQAihCCFM1bnPso0TQiypQf3jhBCryzi/RAgxvYL7ooQQh0raI9EGUjwk9iABiC1DwOKAuS6wR+IEpHh4OJZWxSFL62SJ9QW3thSEEBdKlI8SQlwQQkRZjkdYjhUhxE4hRGTJOhRFSQGSgf5FnhMLpCiKklKeDSXqLWaPECLWcs+FEnZHCiFWCyGmW+0pcVypvRL7IMXDg7G8OEuAkUCI5fR0AEVRelu2ISXKrwHuVxQl2XJ6LnCb5f4Uy7PKYi5wf5Hj+4G5FdlQlKL2WITCek8TIA2YV6R4LGAChpdxbKu9khpicLUBEocSB8QripIIIIQYD+yk7BfKBKwGvraWtxCiKIrZcn+apVxZfI0qFiZL+f5A5yraYKU/kFjinsNFCyiKMtJyLbLEsa32SmqIbHl4NmFA4YCkpXtR3ssUizp2MaJEt2KCpfm/Gii3C2B5YROB/pYuT5qlvqrYYKUpEGfpflxAFY6i96SUKF/02CZ7JTVHiodnk4r6IgJgEQVzOWUTFUUZjyog0y3l41BF5TZLt6KyGZklqN2VAVwdKK2KDVYOAQmKooRYP0WfUR7VsFdSA6R4eBBCCJP1YzllbUlYZ0LmoXYvysL6Qg+33BMJhKK2IMyW+0dazpXH16gv7whL3VW1odhzrPcIIeZi26xNVe2V1AApHp7DhaIfIUScpYtwP+qLZ51VGV/knoSSvhWW7scMYK6iKPFQOAOyxnJvrGUmpRRFui7WLgs22FCUBCGEYnlO0XsiKT4YWyZVtVdSM4Sj4nkIIaKKjNiXvBaH+p8uSlGUGQ4xQCKROBSHtDwsSl9mf9PqP2AZSTdbjyUSibZwiHhYhKHkiLiVAVztX6eg9pElEonGcMWYhwnV6cdKmAtskEgkNUQOmEokkmrhCg9TM1enz0yofgDFEEKMQJ3uw9/fv3OrVq2cZ51Ecg2yc+fO84qi1KrKPU4TjyJuy18B0ZbTkahTe8WwTLnFA0RHRytJSUnOMlMiuSYRQhyt6j2Omm2JA6ItWytrAKzTt5YZGXN507kSicS9cUjLQ1GUBK56GFrPdS6yH++IeiUSifOQA6YSiaRayCX5EoeSm5vLiRMnyMrKcrUpEsDHx4cGDRpgNBpr/CwpHhKHcuLECQIDA2ncuDFCCFebc02jKAqpqamcOHGCJk2a1Ph5stsicShZWVmEhYVJ4XAAKSnlOXGXjRCCsLAwu7UCpXhIHI4UjuqTnFz2ZGRycjJms7rKIyEhgcTEROLjK5+HsOffQoqHROKmJCYmcv/9ZUciSEpKIioqiuTkZCIjI4mNjSUyMrJcsXEEUjwkEjfFKgiVMX68Gh4lJSWFqCjnLVKXA6YSjyc5OZmvvvqKmJgYTCYTKSkpjBgxwm7PT0hIwGQyFXYj4uLiKrmjZiQnJxMdrTppR0VFERkZSUhICPPmzavkTvsixUPiNCYt28f+U5fs+sw29YJ49a62FZYxmUyEhYURGRlJVFQUvXv3pn///iQlJWE2m4mNjcVkuhpfOSUlhYSEhDKfNW7cuGLHKSkprF69mrlz5zJ+/Hhuv/12EhMTa/zcikhKSioUP7PZjMlkYsKECQwfPrxQTJyBFA+JxxMZGcmOHTsYN25cYetg6tSpTJ8+HbPZTHx8fLGXNzIy0uaXOSEhgd69ewOqOKxYscIuz7WV+Ph4JkyYgMlkIjIykoSEBLvXUR5SPCROo7IWgiOxikZiYiIjR45k9Wo1ta7JZOLQoUPFylalhZCamlqsm2Ktp6bPLY+UlJRyWxZxcXE2zbjYCykeEo/H6g+RmJhIWloaI0aMYMeOHYUvetOmxbM6VKWFMHLkSJKTk0lISCh8qe3xXFBbNUlJSSQkJBQKVGJiYrHxmnHjxjFjxgwiIyMLfzZn4bAAyPZCLsnXNgcOHKB169YutSE+Pr5wOtNKSkpK4bRmybGJmtRRdLrUHs8tq56aCkRZfxMhxE5FUaLLuaVMZMtD4tGYzWbmzp3LhAkTip23vuj2ICUlpbAOez63JGazmdBQ90lDI8VD4tGYTCZ27tzp0DoiIyMdXgeoXZairSdXI8VDItEIjvYfqSrSw1QikVQLKR4SiaRaSPGQSDRKVZfk2xspHhKJG2PLknwrM2Y4N+2zFA+JxE2xZUl+0bJWr1lnIcVDInFTbF2S7yo8Zqo2vyCfdcfX8cvJXziTfgYfgw+RwZF0r9edqNpR6HV6V5socRGevCTfehwbG8v06dMdWm9JPEI8Tl05xbPrn2Vf6j6CvYNpGNiQf9L/YcPxDczbO48IvwgGtBzAA60eINAr0NXmXruseAFO77XvM+u0h9unVVjEk5fkA6SlpVVQ2nFoXjzSstJ4bOVjXMq+xNSbpnJ749sLWxnpuelsOrmJbw5+w7u73uWjfR8xquMoBrQcgEGn+R9dYiOevCTf2upwBZp/gyZvmczZjLN83O9jrq91fbFr/kZ/+jbuS9/Gfdmfup9ZO2cxbfs0vvv7OyZ2nUjH2h1dZPU1SiUtBEfiqUvyU1JSSElJIS0tjbS0NJKTk50WilDT4rHj9A7WHFvDmKgxpYSjJG3C2hDfO56VR1cyc/tMHlnxCANaDuD56OfxMfg4yWKJK/DkJfnW8/Hx8aWmbh2Nppfkj1ozin3n9/HzfT9XSQDSc9N5f9f7fHrgU5oGN2V6z+m0DG1pL5MlRZBL8j13Sb5mp2pPp5/mlxO/cF+L+6rccvA3+jO+y3jmxs7FnG3moeUP8dmBz3B3IZVUHeuS/JL/lSMjI4mLiyMuLq7GL7h1Sb7ZbLbrc0sil+TbiTXH1qCgcFfkXdV+Ro/6PVj6r6W88usrTNs+ja3/bGXqjVMJ8Aqwo6USVyKX5DsOzbY8fjn5C02Cm9A4uHGNnhPmG8b7vd5nfMx4fjnxCw8uf5CUi65dMyCRlIUjWjM1QZPikV+Qz+6zu4mOqFIXrVyEEDzc5mHm9ZnHpZxLPLT8IdYfX2+XZ0sknoomxeNv899cyb1Cp9qd7PrcmDoxfHnnlzQKasTotaOZ89scCpQCu9YhkdgLuaq2Guw6uwuAqAj7z2fXDajLJ/0+4V9N/8Xs3bMZt3EcWXn2ySoukVQVW1fVlixnnT4uz6/EHmhSPA6kHSDUJ5R6/vUc8nwfgw+v3/A6z3V+jlVHVjF01VBSM1MdUpdEUh62rqotq9zUqVOJi4srNiVtbzQpHocvHqZJcBOEEA6rQwjBo+0e5a1b3uKvtL8Y+NNAUsxyIFXiPGxdVVuyXEJCAjExMYDqueooj1OHTNUKIeIAMxClKEqpCCVFrkcqilLlFFeHLx4mtpFzpqxiG8US4RfBU2uf4uGfHmbWrbPoWrerU+qW2AdPX1Vbkh07dhSWS0xMdFj6SbuLhxAiCkBRlEQhRKQQIkpRlOQS11MURUkWQsSWvF4ZF7IuYM420ySoib1NL5f2tdrz+Z2fMypxFI+vfpxXe7zKPc3ucVr9kprh6atqyyIsLIyoqCgSExOLubfbE0e0PAYA1pBGKUAsUFIcpgO9UVseiVV5+OGLhwFq7N9RVeoH1GfxHYt5dv2zvLz5Zc5lnGNY+2EO7Tp5GtO3T+ePtD/s+sxWoa0Y32V8hWU8eVVtWViFElTh3LFjh2bEwwQUDTAQVvSipcWRIoS4AAwv6wFCiBHACICGDRsWu/ZP+j8ANAhsYD+LbSTQK5DZt83mpc0v8e6udzmbcZYXurwgAw1pAE9dVVsWcXFxhfWYzebC8Q9743T3dCGECXW8YyowTwiRrChKsZFIyzhIPKgL44peO5NxBoAIvwin2FsSo97I1JumUsu3Fp/s/4TUrFSm3jQVb723S+zREpW1EByFJ6+qLatcZGQkJpOJhIQEUlNTHdfaURTFrh/ULkmsZT8OGFfi+jjAVN71kp/OnTsrRZm6barS7bNuijvw8e8fK+0+bqcMXjFYuZh90dXmuCX79+93tQnK3LlzldWrVxc7d+jQIWXJkiXKkiVLlAsXLtitDns/t6x6akpZfxMgSaniu+6IqdqvAGu7KhJIhMIWR0nhSkBthdjMmfQz1ParXVMb7cLgtoOZftN0dp/bzeAVgzmdftrVJklKIFfVOg67d1sUdUwjWggRC5iVqzMpa4DOiqLMEEKME0KkAKFKFadqz2acdVmXpSzuiLyDEJ8Qnl73NINWDGJen3k0CmrkarMkFuSqWsfhECcxRVHiFUVJLCoMiqJ0LrI/Q1GUhKoKB8DZzLPU8qtlL1PtQvd63fmo30dk5WUxaMUgu88oSCQgV9XWGHOWmRDvEFebUYo2YW34+PaPMeqMPPbzY4XrbyQST0VT4pGVl0VWfhYmH/dR36JEBkey6PZFhPqGMmLVCDad3ORqkyQSh6Ep8TBnq4Newd7BLrakfOoF1OPjfh/TOLgxo9eO5ucjP7vaJJejyPCOboM9/xaaEo+L2RcBMHm7Z8vDSrhvOAv6LuD68OsZt2EcCX85blm0u+Pj40NqaqoUEDdAURRSU1Px8bFPtgBNxTC9kH0BsK947D5uZuHmw2w/nMblrDzqBPvQLTKU/tHXcX2D6tcT5BXEh70/5Jn1zzBpyyQu51xmSLshdrNbKzRo0IATJ05w7tw5V5siQRXzBg3s452tKfGwZ7dl0rJ9fLT5SKnzf5+9wt9nr/Dp1mOF5964tz0PdrmuyutYfA2+vHfre7y46UXe2vkWF7MvMiZqzDW1HsZoNNKkifMWMUqch6bEIz0nHVD/q1eX81eyiX796lq8AG8D4/u15OYWtQn2NXLCnMGGv86xeMtR/rmoRhB78du9vPjtXro2CWXhozH4e9v+azPqjUy7aRoBXgEs+H0Bl3Mu82LXF+V6GInm0ZZ45Kri4Wvwrdb9py9m0W3qGgAMOsHuV/uUEoJgv2Da1gvmyVuaAXA8LYPhi5L44/Rlth1Oo+2rK2lVJ5Bvn7wBXy/bBECv0/NKt1cI8gpi4e8LuZxzmf/e9F+MOmO1fg6JxB3Q1IBpRl4GAH5Gvyrfm5tfUCgc93aqz99v3GFTC+K6UD9+fronKW/cwWM3qM3vP05fpvUrP/PMV7/ZPBAohOCZzs/wdNTTrDiygmfXPUt2fnaVfw6JxF3QnHh46byq9R+76xuqcLStF8SsAVVPcK3TCV65qw2Hp97BoO6q+/m3u07SZMJPrNpn+5qWoe2HMrHrRNafWM/oNaPJzMussi0SiTugLfHIzahWq+PXv8+Tlp4DwI+jb6yRDUIIJt/djj+m9KNpLX8ARizeSeMXlnMlO8+mZzzQ6gGm3DCFbae38fjqx7mSc6VGNkkkrkBz4uFv9K/yfQ/N3wbAt0/2sNtMh49Rz5rnbuHnp28qPNfu1ZV8uf1YBXdd5Z5m9zD9punsObeHEatHFPqwSCRaQVvikZdR5cHSDX9d9S/o1ND+a2Ja1QniyLQ7eaSb2pV54Zu9tHp5BTl5lSeL6tekH2/e8iZ/pP3B0JUyvYNEW2hLPKrRbRm8cDtQ8+5KZUy5px0bx94KQFZuAS1eWkHysQuV3terYS/e6/UeRy4dYcjKIZzNOOtQOyUSe6Et8cjLwM9gu3hczsot3G9X3/HrYRqG+XF46h10j1TDtv579q88+9Vvld53Q/0bmBM7hzPpZ3j050c5deWUo02VSGqMpsQjPTe9SmMer36/D4BhNzrPw1EIwRcjuvHREDXo7De7TtL4heVczMit8L6YOjHE94nHnG1m8M+DOXrpqDPMlUiqjabEIzMvs0otj292nQTg+b4tHWVSudzasjb7J/ctPO4weVWx8Zey6FCrAwv6LCA7L5tHf36UQ+ZDFZaXSFyJ5sTDx2DbisCs3PzCfR+ja1zB/bwMxQZTBy/czivf/17hPa3DWvNRv48AGPLzEA6kHnC4nRJJddCUeOTk59ic4mDxFrXZf2f7uo40ySam3NOOz4erKSoXbTlKi4kryC8o3zO1qakpH/f7GG+DN0NXDWXPuT3OMlUisRnNiYdRb5t36f9W/Qm4pstSFj2ahrPrZTWzWE5+AU1f/IlT5vK9SxsFNeKTfp9g8jYxfNVwdpze4SxTJRKb0Ix4KIpCTkEOXjovm8pnW/wsmoRX3anMUYT4e3F46h00ClPHbXpMW8vKClzbrVHJIvwjeDLxSbb+s9VZpkoklaIZ8cgrUF2/vfSVi4c7R60SQrBh7K0817sFACMX7+S1H/aVW762X20W9l1Ig8AGPLXmKTaf3OwsUyWSCtGMeOQUqGtTbGl5bD+spsptVSfQoTbVhNG3NWfJ490B+PjXI9w4fW25ohfuG87CvgtpHKTGRd14YqMzTZVIykQ74pGvioctYx6fW9aXPNS1YSUlXUtM41B2vqQm8TlxIZMmE34qNktUlBCfEBb0XUDzkOaMWTeGtcfWOtNUiaQUmhMPW7ot3/+memje06m+Q22yB2EB3vz939sLj1u9/DMnyxlIDfYOZl6febQJbcNz659j1ZFVzjJTIimFdsSjCt0WK0E+2ojUZdDrODLtTjo3Uhfu3TBtLVsOlb1ILsgriLm959IuvB3jNo5jxeEVzjRVIilEM+KRm6+6d9vS8tAqS5/owX96qeEPH5y3lQWbDpdZLsArgA97f0jH2h154ZcXWHZomTPNlEgADYmHNWRfZS2PoovhtMizfVoyf1A0AFN+3M/TX5adttLf6M/s22YTExHDxE0T+fbgt840UyLRjnhYuy2VDZj+amnu92zhXsmwq0JsmwgSn+0JwHe/neKmGWUPjvoZ/Xj/tvfpXq87r/z6Cl//+bUzzZRc42hHPGwcMN1oWXzWs3m4w21yJM1qB7L7lT4AHE/LpPELy8t0afcx+PBur3e5qf5NTNk6hc8PfO5sUyXXKJoRD6uTWGXBjzceVMXjZg23PKwE+xk5WGQmpumLP5UZJ9Vb783bt77NrdfdytTtU1m0b5EzzZRco2hGPPIV1f9BLypeIXs8TZ3mbFY7wOE2OQOjXsfhqXdg8lNFs92rK8ucyvXSe/HmLW/Su1FvZibNZMHeBc42VXKNoRnxKFDUtSqViYcVT0rpKITgt1f6FI7j3DBtLbuPm0uVM+qMzOg5g9sb387byW8zd/dcZ5squYbQjHhYuy3XcprGRY91YWTPSADu/mAzvxwsHVzIoDPwxk1vcFfkXbz/2/u8v+t9t17rI9EuDhEPIUScECJWCDGunOtRljJxtj6zqi0PT2XCHa2Zcd/1AIz6LJmdR0sHWTboDEy5YQr3NLuHuXvm8k7yO1JAJHbH7uIhhIgCUBQlETBbj0swQVGUBCCynOulsI556ET5JuflV57uwBPoH3MdG8beQoi/FwPnb2XdH6Ujrut1eib1mERcizgW/L6AN5PelAIisSuOaHkMAKwd8hQgtuhFS2tjB4CiKDMURUm25aGFLY8Kui1HUtVE2O4Uw8NRNArzZ+kTPWhaK4Dhi5L4/reTpcrohI5Xur3CAy0f4JP9nzBt+zQpIBK74QjxMAFpRY7DSlyPAcIsXZcyuzVlUTjmUUG3Zd+pSwC0qRtk62M1TXiAN1+M6EZUoxCe/uo3Fm85UqqMEIIXu77Iw60f5vM/Puf1ra8XCrFEUhNcNWCaam1xlDXuIYQYIYRIEkIknTunDgpav/AVdVv2W8Wj3rUhHqAu/lv0WBd6tazNy9/v4701B0u1LoQQjIsZx5B2Q/j6r6+ZvGWyFBBJjXGEeJiBUMu+CSi5PDQVtTtjLRtT8gGKosQrihKtKEp0rVrq9KT1y24QhnIrPnhWTRjd3EN8PGzFx6jnw0c6c2+n+ry5+i+m/HiAgoLSAvJM1DMMbz+cpQeXSgGR1Jjy38Tq8xUQbdmPBBIBhBAmRVHMQAJgbW2YsIx/VEaeonZbKmp5pJxTxSOy1rUlHqA6k715fweCfY0s3HyYi5m5TL+vPQb91d+XEILRnUYDMG/vPBQUXu3+aoW/U4mkPOwuHoqiJAshooUQsYC5yIDoGqCzoigpQgizpbsSpijKDFueW1BQ+YDp8Quq52WDkKolw/YUdDrBq3e1weRn5O3Eg1zKyuW9BzsVy1tjFRCd0DF3z1wKlAIm9ZgkBURSZRzR8kBRlPgyznUu43qCrc+0ZarWunDMVUme3AEhBE/HtsDka+S1Zft59KPtzB8cQ4C3oViZUR1HIYTgw90foigKk3pMuqYd8CRVxyHi4QhsXdsiUXn0hiYE+xl5fskeHp6/jU+GdCHY7+qiwkIBQTBn9xwUFCb3mCwFRGIzmhEP6WFade7t1ABfo4HRXyTz4LytLB7ahbCA4hn3nuz4JALB7N2zURSFKTdMkQIisQnNdHRt6bZIStOvXR3mDYrm0LkrPBC/lbOXskqVeaLjE4zqOIplKct4afNL5BeUHcFdIimKZt7EwqlanWYaS27DLS1r8/GQLpw0Z9J/7pYyl/Q/3uFxnur4FD+m/MjEzROlgEgqRTPiYfUwlS2P6tG9aRiLh3YlNT2H/h9u4ajFlb8oIzuM5D+d/sPylOW8uOnFwt+5RFIWmnkT5ZhHzencKIQvhncjIyeP+z/cwt9nL5cqM/z64YyJGsNPh3/ixV+kgEjKRzPika/kIxAeFeTHFbSrH8yXI7pToMCAuVsLXfqLMqz9MJ6OepoVR1Yw4ZcJUkAkZaIZ8VAURXZZ7ETLOoF8PbIbXgYdD8Rv4bcyopINbT+UZzo/w89HfuaFX16QAiIphce8jWVFFpeUT2StAL4e2R2TnxcPz99WmBy8KI+1e4znOj/HyiMrGb9xPLkF2s6JI7EvHiMelzLVL3awrzZSTLoD14X68fXI7kQEeTNo4TY2HTxfqsyj7R7l+ejnWXV0lRQQSTE8RjzMFvEw+UnxqAp1gn34amR3Gof589gnO1hz4EypMoPbDmZs9FhWH13NuA3jpIBIAE8Sjww1KZRseVSd8ABvvhzRjVZ1Ahm5eCfL9/xTqsygtoMYFzOOxGOJjN0wtjB3sOTaxWPEw5oMKdBHOpFVB5OfF58O60rH60yM/iK5zLCGj7R5hPEx41lzbA3Pb3heCsg1jseIR3q26hHp5yXFo7oE+RhZNLQLXZqE8sxXv/HtrhOlyjzc5mFe6PICa4+v5bkNz0kBuYbRlHgIyvfxyM5TxcP3Gl6Obw/8vAwsfDSGrk3CePbr3SzdWVpABrYeyIQuE1h3fB3Prn+2MI+w5NrCJvEQQvxbCNHR0cbUhJw81QPVqHdzPVQUuHAU/tkNF45AgfuFArQKSPfIMJ5P2E1CGQLyUOuHeLHri6w/sV4KyDWKrW18ATwghHgRUFBjkK4GTIqifOMo46pCbr7q5+FlcEMP1HN/wuJ74VLpcYRiNO8Ld70NQfWcY1cF+HrpWTA4huGLkhibsJsCRaF/9HXFyjzY6kF06Hh92+s8s/4ZZt0yCy+9l4ssljgbW/9NX1AU5QVFUforijIANbVCU2CC40yrGrn5btjyyMmAyWHwQZfiwhFYFyLaQ1CD4uUProS3WsNrwfDdKMgtvXzemfh66Zk/OJobm4UzfukevtpxrFSZAa0G8HK3l9l4YiNPr3ua7PxsF1gqcQW2tjxChBBzUFsbyUCKoihLhRCJjjOtarhdt+XsAZjd7epx3EJod1/55fPzYNuHsGqievzbp+onoj08+iP4mhxrbzn4GPXMGxTNiMU7Gb90LwUKPNilYbEy/Vv2RwjB5C2TeXrd07x969t4673LeaLEU7DpTVMUZSkwA7W1MRLYaTl/2HGmVY0cS8vDy+AG4pF66KpwtLkbXrtYsXAA6A3Q4ym17MTT0OEh9fyZvTC9EcxqDzmll9E7Ax+jnvhHOnNzi1pM+GYvn28r3QK5v8X9vNr9VTad3MSYdWNkC+QawOY3TVGUw4qizFQUZYKiKEccaFPZ9VPx2pXCbovOxWMe+XnwniX97s3jof+iqj/D6Av3zoFXzdBzrHru4jF4ox6sGG8/W6uAj1HP3Ec6c2vLWrz47V4+3Xq0VJm4FnG81v01Np/czJi1UkA8HTf4N10FKtAFa5I0ly/Zt7Y46naEW1+s2bOEgF4vqSISM0w9t+1DdUwkZX3Nnl0NrMmlerWqzUvf/V5mesv7WtzHpB6T2HxqM0+ve1rOwngw2hIPd+fCUUg9qO6PWG+/5woBd74JE06Cf2313KK7YXYPp0/1ehv0zHk4itjWanrLT349UqrMv5v/u7AL88z6Z6SAeCgeIx56S3fFpUvz429Rt/fMUV94e+MdAGMPwmOr1OOz+2ByCBzdYv+6KjLDoGf2wM70bhPBqz/s46PNpYe+4lrEFc7CSD8Qz8TjxCPPVeKRkwGZlpgYHR9ybF0Nu8IrF6BuB/X4o37w5UDH1lkCL4OODx6Kok+bCCYt28+CTaUFpH/L/rzU9SU2nNggXdk9EI8RD0Nhy8NFHpsrLS4vUYOdU59OByM3wuBl6vEfP6pjIVkXnVM/FgEZGEW/tnWY8uN+Pi6jBTKg1QDVE/X4enUxnVzO7zFoRjz+Sf+nwlB41oTOLmt57PxY3fZ9w7n1NukJLxZZQj+tIRxa67TqjXod7z3Uib5tI3ht2f4yZ2EebPVg4WI6GQ/Ec9CMeCxPWV7hdWvLIy/fBeJRtLXjHeD8+r38VP+QTo+ox4vvhZ+d5/xr1Ot478GowlmYr3ccL1VmYOuBjI8ZT+KxRBmRzEPQjHhUhkvHPA5aBjDrRTm/7qLc/f7VbszW2fBOx6tz2A7Gy6Bj9sAoeraoxfhv9vBNctnL+a0RyV6mWl8TAAAgAElEQVTYKIMqax3NiYdSzstg9Sy1uqk7ld1fqNvoIc6vuyRNesJzf6n7Fw7DJJM6mOsErJ6o3SPDeH7Jbn7YfapUmUFtBxXGRJV5YbSN5sTDmrO2JNY4Hll5LkiTuP87ddvmbufXXRaBEfBy6tXjN+qCubRLuSPwMaqL6aIbqwGFVuwtHdJwcNvBPNv5WVYcWSEz02kYzYiHQahr+KyZ40ri56WKR0a2C7+IPsGuq7skeoM6DtK0l3r8dns4tcspVVvjgaghDXexen/poMpD2g1RE0sdXiGTa2sU7YiHJcF1ef+lAn3UwMeXs+R/sWI88i3cMEbdj78F/vzZKdUGeBv4aEgMbesH8+RnO1n3x9lSZYa2H8qYqDEsT1nOy5tflgKiMTQnHuWN0gf5qteleJRB78nwf2+r+18MgO3znFJtkI+RRUO60LJOICM/3ckvB8+VKjOs/TCe6vgUy1KW8cqvr0gB0RAOEQ8hRJwQIlYIMa6SchVeL4pRp7YsKm95OHkK0EmzGTUmeggMTFD3f3oeNs50SrXBfkYWP9aVyHB/hn2SxK+HSieWGtlhJE92fJIfDv3Aa1teK7drKnEv7C4eQogoAEVREgGz9biMcrFAb1ufW7l4uKjlkXNF3Rr9nFtvdWjeG0ZsUPfXvg5r/+uUakP8vfhsWFcahvox9OOkMlNbPtHhCZ7o8ATf/f0dk7ZMkgKiARzR8hgAWDMnpwCx9nioUa+KR3ndlgBLyoXLzh4wtYqZXiPJpup1hMc3q/sbZ8DqV5xSbViAN58N70pdkw9DPtrOzqMXSpV5osMTjLx+JN8c/IbJWyZLAXFzHCEeJtQYp1bCShYQQkRZWiY2Y215lCceOlcFAbJ6lwrNDB9BnXYwaru6v/kdWPGCU6qtHejDF8O7USvQm0cXbmfPCXOx60IIRnUcxfD2w1l6cCmvb31dCogb46pvfGhFF4UQI4QQSUKIpHPn1EG2ygZMXYb1yy00li+mVksYnazub5sD66Y6pdqIIB8+H94Nk7+RRxZs54/Tl4pdF0IwutNohrYbypK/lvDGtjfKdQyUuBZHiIeZq+JgAlKLXrSl1aEoSryiKNGKokTXqlULqLzl4TKsa1myzBWXc0fCmsJTSer+hmmwLd4p1dYz+fL5sG74GvU8PH87KeeuFLsuhGBM1BiGtBvCV39+JQXETXGEeHwFRFr2I4FEACGENfx3pGU2ZgQQWt6AakkKxcPdYkIYfdWtVr0kw5tfjXq2YizsWeKUaq8L9ePTYV0pUBQenr+Nk+bMYteFEDwT9QyPtn2UL//8kuk7pksBcTPsLh6KoiRD4WyK2XoMrLFcT1AUxTJniM35BCpzEpPUgHqdri6o+2YY/O2cjBrNagew6LEuXM7OY+C8rZy9XDxPjRCCZzs/y6A2g/jswGfM2DFDCogb4ZAxD0u3I1FRlPgi5zqXUaZpEXGpkMpmWyQ1pElPGPCZuv/pfXDuL6dU265+MB8P6cLZy9k8Mn875ozi4QqFEDwf/TwPt36YTw98yv+S/icFxE3QzBSB2455eBKt/w/6vK7ufxADmaWnUx1B50YhzBsUzeHUdAYv3F7K0U8IwbiYcQxsPZBF+xfx1s63pIC4AVI87IkbJq2uMj1GQ7s4dX96YzUPjRO4oVk4sx+KYt+pSwz9JInMnOJu6kIIxseM54GWD/Dxvo+ZlTxLCoiL0Yx4uO1ULUB4S3V7eo9r7bAXcQsgqL66P6WUm47DiG0TwVsDOrLjSBqPf7qzVGwWIQQvdn2RAS0H8NHvH/HervekgLgQzYhHVWZbCpwdTay5xcv+r5XOrdeRPLPv6v5n/Z1W7b861GPav9uz4a9zjPlyV6lUGlYBiWsRx7y98/hwz4dOs01SHM2IhzVxckX5P8ID1DJnLzs5zWGHB9StNQiyJyDE1cDKB1fCrs+cVvWAmIa8/H9tWPH7aV76bm+p1oVO6Hi528vc3fRuZv82m/l75zvNNslVNCMe/kZ/ANJzy0/23DhMXZx2+LyTE0LXaa9uL5cOu6dpvPzgiV/V/e+fhPMHnVb10Bub8NStzfhi+3Fmrvyz1HWd0DGpxyTujLyTd5LfYdG+auQEltQIzYiHn2XVakZe+fE4m9VWvT3/LuGxKKkBEW3hjv+p++9HO20AFeC5Pi14qGtDZq8/xPxfUkpd1+v0vH7D6/Rp1IeZSTP54o8vnGabREPiYW15ZORWLh6HzrpAPIKvU7cnkpxft6PpMhzqW9x03mzhtGqFEEy5ux13tK/D68sPsHRn6YjsBp2BaT2ncet1t/LGtjdI+CuhjCdJHIF2xMNQebelsOXhCvGIfU3drhjv/LqdwbA16jYjFXYscFq1ep1g1oCO3NgsnHFL95BYRjxUo87I/27+HzfVv4nJWybz/d/fO82+axnNiIct3ZbmEYGAi8Sj3X3q9qQHtjxAHUD9z2/q/vJn4XLpl9hReBv0zH2kM+3qBTHq82S2paSWKuOl92LWrbPoVrcbL29+udIkYZKaoznxqKjlUTfIB4DTl7LKLeMwRJF4Ik7MF+tUQptAr5fUfSd2XwD8vQ18NKQLDUJ8GfZJEvtOlf4de+u9eafXO0TXiWbipomsOrLKqTZea2hGPGwZ83BZQCArXZ9QtysnutYOR9Jz7NX9DTOcWnWovxeLh3Yl0MfA4IU7OJpa+h+Jr8GX93u9z/W1rmf8xvGsPea8vL3XGtoRDxvGPFxO7Kvqdtdi19rhaMZaZj7W/dfprax6Jl8WDe1KfkEBgxduJ/VKaZ8eP6Mfs2+bTZuwNjy34Tk2ntjoVBuvFbQjHlY/jzw3Fg9rbA+A03tdZ4ej8Q+D7k+p+9MaOr36ZrUDmD84hn8uZpW5DgYgwCuAOb3n0NzUnGfWPcOvp351up2ejmbEo3DAtIJui1vwwOfqdpGbpJ50FH2LRF4/7Pz/7J0bhfDug53YfcLM6C+SycsvvSgxyCuI+N7xNA5uzJi1Y9hxeofT7fRkNCMetox5wNUUDOfLaM46hVZ3qtuMVMj2cGe1YZbxhE/uckn1fdvWYdK/2pJ44Cyv/LCvzEVyJh8T8/rMo35AfUatGcWus85JuXktoBnx8NGrMymVdVuiG4UAlBna32lYB04X3+M6G5xBgyLxnXZ/5RITBnVvzOM3N+Xzbcf4YN3fZZYJ9Qllft/5RPhF8ETiE+w55yGrn12MZsTDuiS/slD80Y3V2MsuFY9+lkjkJ3ZArgumjZ2J1ffj2xEuM2Fc35bc07Ee/1v1FwlleKEChPuGM7/PfEJ9Qnl89ePsT93vZCs9D82IhxC2TcN2trQ8ko6UzkrmNISAjg+r+4v+5To7nEFoEzCorUL2/+ASE3Q6wYy4DtzQLIwXlu5h41+lc+ICRPhHsKDPAgK9AhmxegR/ppVecCexHc2Ih610vE6NqZx8zMWpEP71rro9vg3SDrvWFkczYr26/foRl5ngZdAx5+HONKsdwBOf7uTAP5fKLFc3oC7z+87HR+/D8FXDOWQ+5GRLPQePEw8fo5skX9Lp4e7Z6v67HV1ri6Op3frq/iXXhSUI8jHy8ZAuBPgYGPZJEufKietyXeB1LOi7AL1Oz7BVwzh66aiTLfUMPE483IpOA6/u7/DwgDX9pqvbpcNcakadYB/mD4ohNT2bEYuTyMot7QMC0CioEQv6LCC/IJ9hq4Zx6oqHxWJxAh4tHiVD2LmEpy3OYsufg0v/uNYWR9LFMmB6dLNr7QDaNwjm7QEd2XXMzNiEPeXGOY00RRLfJ5703HSGrxrOuYyyx0okZeOR4tG+fjDg4kFTK6aGVxeTvdUKPDVgr67IVyk3s/xyTqJfu7qM79eKZbtP8c6a8iOgtQptxZzYOZzLPMeI1SO4kOXCWTqN4ZHi0bdtBAAr9zlv2XiF9Bx7dUZiXi/X2uJIrCkbnBjvoyIevzmSuM4NeDvxIN//drLcch1qdeD9Xu9z/PJxHk98nMs5l51opXbxUPGoA8DKfaddbEkRJlj8D04lw88TXGuLo4gZqm73fOlaOywIIXjj3vZ0aRLK2IQ9JB8rv1XRpW4X3rrlLf668Bej1oxy/2UQboBHioc1oljJ5MkuRW+EZyyOSVtnw7a5rrXHEVzXTd260aJAL4OODx/uTN1gH0YsSuKfi+V/J3o26Mn0m6az+9xuxqwbQ3a+i5Y4aASPFA9bHcqcTnB9eHyTur9iHOz71rX22Bude36dQv29mD8omsycfJ78LJnsvLJnYAD6NO7D5B6T2frPVp7f8Lx7JhlzE9zzr10JlbmoF+Vippv98eu0h4eXqvtLHoUjm1xqzrVC84hAZt7fgV3HzEz5sWLX9Lub3c3ErhNZf3w9E3+ZSH5B+WJzLaMp8fAzVB6K0Mp9UQ0AWJJ03KE2VYtmseoUbngLWPxv2O8hAXvzyk/I5Q7c0b4uI2+O5NOtxyr9XjzQ6gGe7fwsK46sYPLWyVX6h3WtoCnxiPBXZ1HOZpyttOyQGxoD8NHmIw60qAaYGsJjK6FuB/h6MKx93ak5URzCAcvalshbXGlFhYzt05IeTcOY+N3v/H6y4ihoQ9oNYeT1I/nm4DfM2DFD5sUtgabEo36Amnz55JXyp92stLP4erjVoGlJ/EJh0PfQcSBsnAmf/B9cLHtVqCZY/py67TnOtXZUgEGv470HOxHu78XIxTu5kF5xa2lUx1E80uYRPjvwGe/tes9JVmoDTYrHictVe8Hc+j+Glx/c8wH8e546SzHnBtj1KRRorJmckwFZlsWIjW9wrS2VEBbgzZyHO3Pucjb/KSOZdlGEEIyNHluYWPuj3z9yoqXujabEo0GAOo5hS8sDoGGoOkZS0fy+23B9fxi5EWq1hO9HwcK+8M9uV1tlO7PaqtuY4a61w0Y6XGdi0t1t+eXgeeaVkcqyKEIIXur6Ev0a9+OtnW/x7UEPmyWrJtoSj0BVPGxteYy5rTkAb676y2E22ZWwpjDkZ3U1bloKxN+iCkmqmy8bX/cGZFqWAtwx07W2VIEHYq7jjvZ1eHPVn+w9UfH4h16n540b36BHvR68tuU11hxb4yQr3ReHiIcQIk4IESuEKLPzK4QYYflMr8pzC7stV2wTj3s6qeV/PVQ6w5jbotOpq3FH74QuI2HPEjXBdMJQOOOG0a9+GgcbLH/GMbuLJ79yc6weqGH+3oz5chcZORUPWBv1RmbdMot2Ye0Yt2HcNR9Q2e7iIYSIAlAUJREwW4+LXI8FEhVFiQciLcc2US+gHmB7t0VfJAlUTp7GxhB8TXD7NHVKt/tT8OcKmNMdFvZTx0RcHVw5PRVeC4btFk/Zkb9ASGOXmlQdTH5evDWgA4dT05ny44FKy/sZ/fjgtg+4LvA6Rq8dfU2HM3REy2MAYA3jlQKUFIfIIudSLMc2EeQVBFQt8dPt7dR1LpX1a92WwAjoMwWe+V1Npp1+Tu3K/K8FfDNS9VJ1ZuKliydgZjOYWeTPNvYQ1L3eeTbYmR5NwxnZsylfbD9m03ook4+Jub3nEuwVzBOJT3Dk4hHHG+mGOEI8TEDRtfBhRS8qihJvaXUARAE2Z4aujtv5q3epA3kzV2o8XqVfKNz4DDyVBI+tgnb/hoMrVS/VGZHw8f+p072H1tlfTLKvwIFlaktjVltVwEANNfDaRfAPt299LuDZ3i1oVz+IF5bu4ezlyoNWR/hHMLe32uoauXokZ9LdZAW3EzG4qmJLdyZZUZTkMq6NAEYANGxYs4xkdYJ9Cvfz8gsw6DU1RlwaIaBhV/Xzf2+rEdoProS/VqmOZlbCW0DtNhAaafk0gYAI8DGpXSK9sfhzCwrUQc8rZyH9LJw/CGd+h9O/w+k9kF/EH+Kud6HzYOf8vE7Cy6Dj7QGduOPdX5j0w34+GBhV6T2NgxszJ3YOj618jJGrR/LJ7Z8Q7B3sBGvdA0eIhxkIteybgPJGK2MVRRlf1gVLyyQeIDo6usZOGr3bRLB6/xneWXOQ5/q0rOnj3Ae9ARp1Vz+xr0GmWV3yf2InnExSX/o/foSCMgYCjX4gdFCQD0oBFOSq26L4mNS1OF1HQvO+0LBbadHxIJrVDmDMbc2ZufJP7t53mj6W0A4V0SasDe/1eo/HVz/Ok4lPMq/PvMLshp6OsLcDlaVFEa0oSrxltiVRUZRkIYRJURSzpcwIa9dFCBFrGVwtk+joaCUp6WrPpv0n7QHYO9j2Zd/mjBw6Tl4NwJFpd1b9h9Iy+Xlw8ThcOKwOcmaZVZGxOnQJAUKvioJ/LbUL4l9Lba0E1dfU7Ik9yM0v4K73NmHOyGXNczfj723b/9c1x9bw7Ppn6Va3G+/3eh+jxkRWCLFTUZToqtxj95aHRSiiLbMo5iLdkjVAZ8v56UKI8agtlPurU092fjbeem+bypr8vAr3U85dIbJWQHWq1CZ6g9plCW3iaks0gVGv47/3tuO+OVuYs/4Qz/e1raV6W8PbeK37a7zy6yu8/OvLTL1xqvuGhrATDhkAsAyKJhYZGEVRlM6WbaKiKCGKojS1bMttdZRF8xDV8etvc9mpBcvj3Qc7AXDvbJktXVIxnRuFck/HesT/ksLxNNsjit3b/F7+0+k/LE9Zzru73nWghe6B5kYPo2qrA1k7T++s0n3/6qD6iFzMzOVylpvF+JC4HS/c3hqdgFmrq+adPKz9MOJaxDF/73y+/vNrB1nnHmhOPGLqxABUy7tvZE/VN6H/3K12tUniedQJ9mFQ98Z8+9tJDp6xPSCyEIKJXSfSs0FP/rvtv6w/vt5xRroYzYlH5wg1M/v209urfO8Lt7cC4MA/l7iYIVsfkop5/Oam+Bn1fLCual1kg87AzJ4zaRXainEbx/H7+d8dZKFr0Zx4hPmoPmcZeVWPbi2E4OlYdczkhulr7WqXxPMI9fdiQExDftzzD2cuVe44VhSrG3uoTyij1ozi+GU3jGhXQzQnHjUdwX46tgUAV7Lz2HnUDZJCSdyaR3s0Jl9R+Gxr1fPZhvuGMyd2DnkFeTy15imPywejOfGwB58P7wrAfXO2uHegIInLaRjmxw1Nw/l+96lqfVeaBDdh1i2zOHbpGOM2jvOoYMqaFI8AoyUvi42ra0vSo2k49Sxu6/3nbrGbXRLP5F8d6nE0NYO9lcQ8LY8udbswoesENp3cxJs737Szda5Dk+JxZ6TqJfrz4Z+r/YwN424FYMeRCyzf48EJqCU1plfr2kDN4sL0b9mfh1o9xOL9i/nm4Df2Ms2laFo8lh9eXu1nGPU6lj11IwCjPk/maKrty/wl1xbhAd5E1vKvceL0sTFj6VGvB1O2TiHptM2Lyd0WTYrH9eFq7IiDF8rPfm4L7RsEM76fOn1788z1cvpWUi4tagdy+HzN/sEYdAZm3jyTBgENeH7D85zLOGcn61yDJsVDr9Pb7VlP3NKU3m3UfDAdJq8iPVvjuVMkDqFOsA9nLtU8d22QVxCzbplFRl4GYzeOJa+sFc8aQZPiAWAQ6pq+Y5eO1fhZ8wZF066+GqXsrvc3VRrLUnLtIQTYa5lbs5BmvNztZXae2anpXDCaFY8nOj4BwPy98+3yvB9H38TEO1pz5Hw6/eduqbJTkMSzuZCeQ6CP/Rah39X0Lu5vcT8Lf1/IumPr7PZcZ6JZ8Xig1QMAfPu3/XJoDO8ZyfzB0Rw+l85d723i10Pn7fZsibbZe/IireoG2fWZ47uMp3VoayZunmhTClV3Q7PiYQ2GDNg1CXGvVhEsfbIHgT4GBs7fxsyVf5CV6zmOPZKqs+/URQ6dS+fmFrXs+lxvvTczb55JTn4Ok7dM1pzDombFA6B9uBpVbOnBpXZ9bqs6Qfzw1I3ERTXgg3WHuOPdX9iipdwvErsya/Vf+HnpubtjPbs/u1FQI8ZEjWHDiQ0sS1lm9+c7Ek2Lx7SbpgEwectkuz/b39vAzPs7sOixLuTmF/DgvK0MX5TE32ddnC9F4lQ+33aMxANneTq2ebGIdPZkYOuBRNWOYtr2aZrqvmhaPBoGXY2snpXnmAHOni1qserpmxnbtyVbDqXSZ9YGRn+xi32nnJgrReISvtt1kpe//51bWtbisRscF8ZRJ3RMuWEKOfk5vJP8jsPqsTeaFg+Afzf/NwDjN5YZiN0u+HrpGXVrMzaMvYXhPSNZ98dZ7nx3Ew/N28ryPf9oLxudpEKycvP57/L9PP3Vb8Q0DuH9h6IcnrKjYVBDHmr1EMsOLaux86OzsHv0dHtTMnp6SXILcolarIYmTH442SlRqy9m5vLp1qN8vu0YJ82ZhAd48X/X1+OuDvWIamjy+MC3noqiKKzcd5oZP/9Jyvl0BnZtyKt3tcXL4Jz/sRezL3L70tvpXKcz7/Vyrv9HdaKna148AMasHcPa42uJqRPDwr4LnWQZ5BcobPzrHF8nHWfNH2fJySugXrAPt7WOILZNBF2bhOJjtJ83rMQxZOXm88PuU3y8+Qj7/7lE01r+vPavttzU3L6zK7YQvyee93a9R8JdCbQMdV6OoWtWPPIK8ui0WI2O/t3d39HU1NQZphXjclYuK/edYeW+02w6eJ7M3Hy89Do6XmcipkkIXZqE0blRCAE25gGROJas3Hy2pqTy455/WLnvNJez8mgREcDwmyK5t1N9l2UWNGeZ6bWkF/1b9ueFLi84rd5rVjwAVhxewbiN4wBIfiQZo851SXeycvPZciiVXw+dZ/uRC/x+8iL5BQo6AW3rBRPV0ESrukG0qhNIi4hAmxMLSapPTl4B+05dZOfRC2z++zxbUlLJyi0g0NtA33Z1uC+qAd0iQ92iyzl67WgOXjjIin+vcJo9bpH0yVXc3uR25u6ey6GLh4haHMWeQXtc9kXwMeq5tVVtbm2lxoFIz85j1zEz2w+nsv1IGkt2niAj56rjWcNQP1rVCaRVnUCa1g6gQYgfDUP9CA/wcosvs9Y4fyWbP09f5o/Tl/nz9CXL9jLZloHtJuH+PBDTkJ4twunRNNztupZd6nRh/fH1pGWlEeYbVvkNLsJjxAPgm7u/ocOiDgDctuQ21vZ3jyDH/t4Gbmwezo3N1WzyBQUKJy5k8keRL/aB05dIPHCGgiINQV+jngYhvlwX6sd1Ib40CPGjdpA34QHe1ApUtyZfIzrdtSUw6dl5nDRncvJCJicsW/U4g2NpGZy/cjUpd5i/Fy3rBDKwayOiG4cQ3SiE2kE+FTzd9dT1rwvA2YyzUjychU7o2PnwTjp/2plzmed4K+kt/hP1Hww69/oxdTpBwzA/Gob5FUumnJWbz7G0DI5bPxcyC7c7DqdxuYxwAQadIDzAm/BAL8L8vQn2NRLkayDIx0igz9X9IF8jQT4G/LwM+Bh1+Br1+Hjp8TXqMTq5f68oCpm5+VzJziM9O5/07DzLfl7hOXNmDmlXckhLzyE1Xd2q+9lk5RafGjfqBXWDfalv8qVXq9q0iAikVZ0gWtYJpFagbSlJ3YnMvEwAvA3ubbt7vVV2wEvvxa5HdjFt+zQ+2vcRu8/tZlKPSTQObuxq0yrFx6inRYQ6DlISRVG4lJXHucvZnL+SzbnL2cX2z1/JJjU9h6Op6VzKyuNSZi55BbaNZ+l1QhUTox5vgw6DXqDXCQw6gV6nw6ATGPTWY2GxBxQABQoUBcViY4HlfE5eATl5+eTkF5CTV0BuvmI5V0BOvm1+Mb5GPaH+XoQFqJ/mEQGE+XsR6u9NPZMPDUJ8qW/yo1agd6FdnsCvp37F3+hPw8CGlRd2IR4nHqBGbHqp20t0rN2RN7a9wX0/3MfA1gMZ2n4owd7BrjavWgghCPY1EuxrpFntyhN1W/+7X8rM41KWmmLzYmYumTkFZOXmk5mbT5blk5mbr57PU48LChTyChTyCxRy8xXyCwoKj/PyLYJkiW+h0wkMQiAE6CzjM0IIvPQ6vA06vAw6jHqBl0GHl16vbg1qyyfAW4+/twF/bwOBlq2/t4EAbwNBvmor6Vpjx+kd/HT4Jwa2Huh2LeaSeMxsS3mcyzjH28lvs+zQMrz13twReQcDWg6gTVgbO1opkdScVUdW8dLml4jwi+DL//sSf6O/0+q+pqdqK+PghYN8duAzfjr8E5l5mbQNa8ttDW+jZ4OetAhpIWc1JC5BURR2ntnJ/L3z2XxqM+3D2/POre9Qy8+5DmpSPGzgUs4llh1axg+HfmB/6n4AIvwiuLH+jVxf63rahrWlqamp2zcZJdolNz+XXWd3sfHERjac2MCRS0cweZsY1n4YD7V+yCU+SlI8qsi5jHNsOrmJjSc2su30tsJ0gL4GX1qFtqJFSAsaBTWiUVAjGgc1pl5APSkqkiqRm5/L0UtHOZB2gD/S/uBA2gH2p+4nPTcdo85IdEQ0/Zr04/Ymt+Nr8HWZnVI8aoCiKBy7fIzfz/9e+Dl08VCx/KIGYaBeQD1q+9Umwj9C3fqp21q+tTB5mzB5mwjyDkInNL9gWVIJuQW5XMi6QGpmKqlZqaRmpnI6/TQnrpzg5JWTnLh8gjMZZwoj3XnrvWkR0oI2YW3oUa8H3ep2w8/o5+KfQkWKh51RFAVztpmjl45y5NIRjl46yonLJzibcZYzGWc4m3GW3ILSuV50QkeQV1ChmAR6BeJv9Mff6I+vwbdw39/oj5/RDz+DX7FjH70PXnov9aNTt1KM7EOBUkBWXhZZ+Vlk5WWRmZdJVl4WGXkZV48t19Jz07mUc4nLOZe5lHOpcP9yzmUuZF3AnG0us45avrWoH1Cf+oH1aRDQgEZBjWgV2oomwU3ctuV6TbunOwIhBCE+IYT4hNCxdsdS1xVF4UL2Bc5mnOV85nnM2WbMWWZ1a/1kmTmfeZ5jl4+RkZtBem46GXkZVbbFoDPgrffGW++NUWfEW+9dTGCMeiM6ocMgDOh1+mL7emH5VLZf5FzJAWTB1WnYYsclzpdXvkApQFEU8pV8FBTyC/IpoIACpeJPvpJfeF/R83kFeeQW5MvF12EAAAS5SURBVBb/5KvbUtcs53Pyc8jKr1rQKIMwEOgVSKBXIEFeQQR6BVLHrw4mbxNhvmGE+YSpW98wwn3CCfcLd2n3w5k4RDyEEHGAGYhSFGVGVa9rBSEEoT6hhPqEVum+AqWAzLxMVUhyM0jPS78qLLkZZOdnk5Ofo24LcsjJzyk8zi3ILbxeWCY/hzwlj9z8XDKVTPKUPPIL8slXLJ/K9ouccwU6oUOHTt3a8NELPUadEYPOgFFnxKg3YhAGvA3eBOgC1HOW84X7lo+v0Rdfgy8+eh98Db6FHx+DT7Gtr8EXP4MfvgZfORNXDnYXDyFEFICiKIlCiEghRJSiKMm2Xr8W0AldYTfFnSj6Hx5A9RulMKp3yePC+8orV+S89aUXQhTfIuTLqVEc0fIYAKy27KcAsUByFa5LXIQQojATn0RSGY74ppiAounESy4LrOw6QogRwAjLYbYQ4ne7WuhYwgGtZIvSkq2gLXu1ZCtAlcOWueW/GUVR4oF4ACFEUlVHgV2JluzVkq2gLXu1ZCuo9lb1HkfM/5kB6wiiCSiZLamy6xKJRAM4Qjy+AiIt+5FAIoAQwlTRdYlEoi3sLh7WmRMhRCxgLjKTsqaS6+URb28bHYyW7NWSraAte7VkK1TDXrf3MJVIJO6J9Hm+RhBCxAkhYoUQ4yopV+F1ifax+lqVc82m7wm4mXhUZnhVfjBHY4OtIyyf6c62rQxbCh3zAHN5Xx5LV7K3M20rx47KfrdRljJxzratLKrwvR1R1nVnYvkbLynnmk3fEytuIx6VGV7VH8yR2GBrLJBomXKOtBy7kgGos1xw1THPLbHx7zxBUZQE1N+ty74HYPP3NsVyPcXV9lrtKOdylb4nbiMeVG64O70AldkSWeRcCldnl1yFLY55UZYvlqup8HdraW3sAFAUZYYbLG2w5XtpbX1GuoG9FVHp96Qo7iQeNfZMdSIV2qIoSryl1QEQBbgmpkDVqNrqPsdR2d85BgizdF1c3n2l8u9CMmqL40KJcprHncTD47A0UZPd4L9NhY55btTqsJXUIlP+bjHuUR4W/yYzMBWYJ4RwdSu0IqrkwOlO4qElz1RbbYlVFGW8c0yqkMoc9yItg3ojgFAX98sr+92mcrXPbkZtibiSyuwdAUy1hJ4YDrid2FXXgdOdxENLnqmV2YoQYoQ1VomrB0xtcNxLsAxAgvoCuJLKfrcJRa6bsIx/uJBKvwtWLL/jssOPOQlLSy26RIutWg6cbuUkZvnPl4I6sGRdGLdTUZTO5V13FRXZWmQ6LA31v9L9GusWuBQbvwdpQIw7tOxssHec5Xqoq7+39sStxEMikWgHd+q2SCQSDSHFQyKRVAspHhKJpFpI8ZBIJNVCiodEIqkWUjwkEkm1kOIhkUiqhVtGT5doE4tb+wBUr08zbuDMJ3EcsuUhsSdmLGtPLB6197vYHokDkeIhsRuKoqSguownl7W2Q+JZSPGQ2BuraMQCc11piMSxSPGQ2A1rrArLwsDQIit1JR6IXBgnsRvW1aVyBfG1gWx5SOyCZYxjJK6PByJxErLlIZFIqoVseUgkkmohxUMikVQLKR4SiaRaSPGQSCTVQoqHRCKpFlI8JBJJtfh/urdroAUK8qwAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f36767d89b0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def f(q,q_old,p,p_old):\n",
" f1 = -q + q_old - dt * (( q + q_old ) * (p + p_old + q + q_old - 2.0 ) + 2.0*p*q+2.0*p_old*q_old) / 4.0\n",
" f2 = -p + p_old + dt * (( p + p_old ) * (p + p_old + q + q_old - 2.0 ) + 2.0*p*q+2.0*p_old*q_old) / 4.0\n",
" return np.array([[f1],[f2]])\n",
"\n",
"def J(q,q_old,p,p_old):\n",
" df1dq = -1.0-dt*(2.0*q+3.0*p+2.0*q_old+p_old-2.0)/4.0\n",
" df1dp = -(3.0*p+p_old)*dt/4.0\n",
" df2dq = (3.0*q+q_old)*dt/4.0\n",
" df2dp = -1.0+dt*(3.0*q+2.0*p+q_old+2.0*p_old-2.0)/4.0\n",
" return np.array([[df1dq,df1dp],[df2dq,df2dp]])\n",
"\n",
"def H(q,p):\n",
" return q*p*(1.0-p-q)\n",
"\n",
"def LotkaVolterra(q0,p0,N,dt): \n",
" q = q0\n",
" p = p0\n",
" q_old = q\n",
" p_old = p\n",
" dt = 0.1\n",
" qhist = [q]\n",
" phist = [p]\n",
" thist = [0]\n",
" Hhist = [H(q,p)]\n",
" N = 200\n",
" for n in range(N):\n",
" t = dt*n\n",
" while True: \n",
" new_z = np.array([[q],[p]]) - np.dot(np.linalg.inv(J(q,q_old,p,p_old)),f(q,q_old,p,p_old))\n",
" new_q = float(new_z[0])\n",
" new_p = float(new_z[1])\n",
" eps = (new_p-p)**2 + (new_q-q)**2\n",
" q = new_q\n",
" p = new_p\n",
" if eps < 1.0e-5:\n",
" break\n",
" qhist.append(q)\n",
" phist.append(p)\n",
" thist.append(t)\n",
" Hhist.append(H(q,p))\n",
" q_old = q\n",
" p_old = p\n",
" return [qhist,phist,thist,Hhist]\n",
"\n",
"dt = 0.1\n",
"N = 1000\n",
"[qhist1,phist1,thist,Hhist1] = LotkaVolterra(1.0/8.0,1.0/8.0,N,dt)\n",
"plt.plot(phist1,qhist1, label=r\"$p_0 = q_0 = 1/8$\")\n",
"[qhist2,phist2,thist,Hhist2] = LotkaVolterra(1.0/4.0,1.0/4.0,N,dt)\n",
"plt.plot(phist2,qhist2, label=r\"$p_0 = q_0 = 1/4$\")\n",
"[qhist3,phist3,thist,Hhist3] = LotkaVolterra(1.0/16.0,1.0/16.0,N,dt)\n",
"plt.plot(phist3,qhist3, label=r\"$p_0 = q_0 = 1/16$\")\n",
"\n",
"plt.legend()\n",
"plt.xlabel(r'$p$')\n",
"plt.xlim(0,1)\n",
"plt.ylabel(r'$q$')\n",
"plt.ylim(0,1)\n",
"plt.title(r'Lotka Volterra')\n",
"\n",
"plt.axes().set_aspect('equal')\n",
"plt.savefig('orbit.pdf')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAETCAYAAADH1SqlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFmFJREFUeJzt3c+PI+ldx/HPd7UTRtFu4uneidBmtJmpgSibgIS83iCurOewQogD7pkTyoVxi3vUnUEIxGnp5i9oL/cwO4YLiD2Mh0gcGbeFEkFyYBxFjAhSb/c6cGDRBr4c/FRvtcf2U+6usqu73y9pNK56nnrqO56u+nT9cNncXQAAzPPSqgsAAFQfYQEAiCIsAABRhAUAIIqwAABEERYAgCjCAheSmbmZ1U6zXPi7ZWaPzrD+LTN7PGX+IzPbmbNc3cyeTdYDrBphAZSjK6k5JbBakvZWUA9wJoQFLp1w1PAsHH08Snfo6ZGAmX080b9uZh+bWT1Mt8O0m9m+mSWT63D3oaSBpLuZcZqShu4+nFXDxHpP1GNmzbDMxxN1J2b22Mx20nompqP1AjGEBS6VsKN8JGlT0rUwe0eS3P1O+PvaRP8nkjbcfRBm70l6Jyw/DGNNsydpIzO9IWlvXg1Z2XpCMKTL3JJ0JOn9TPempJqk+1Om89YLzPTyqgsAlqwlqePuPUkys21J+5q+A61Jeizpg7R/cM3dR2H5o9Bvmg80Doda6H9X0lsL1pC6K6k3scyPsx3cfTO0JRPTeesFZuLIApfNuqTjC8jhdNGsnWdT42sP7YnTRA/C6ZzHkmae0gk76J6ku+EU1lFY3yI1pG5LaoXTSR9rHBTZZYYT/bPTueoF5iEscNkcarzjlSSFEBjN6Ntz922NA2Mn9G9pHCLvhNNEsTumHml8+umePruwvUgNqWeSuu5+Lf2THWOWU9QLTEVY4EIzs1r6J8xKjxTSO5Xe1/h00TTpDvx+WCaRtKbxEcIoLL8Z5s3ygcY763ZY96I1nBgnXcbM9pTvrqpF6wWmIixwkX2c/WNmrXDKZ0PjHW1619N2Zpnu5GcbwumkXUl77t6Rju9QehKWbYY7nV6QORWVnoJSjhqyumbmYZzsMolOXjyfatF6gVmM77MAAMRwZAEAiCIsAABRhAUAIIqwAABEXZhPcL/22mt+8+bNVZcBAOfK/v7+R+5+PdbvwoTFzZs31e/3V10GAJwrZvaTPP04DQUAiCIsAABRhAUAIIqwAABEERYAgCjCAgAQRVgAAKIuzOcsTu3H/yA9/D3p5/8jvcTbAeAceuePpV9vl7oK9o4vXZFe+2XpS1+XPvfKqqsBgMV96c3SV0FYfOU3pN/vrboKAKg0rlkAAKIICwBAFGEBAIgq5ZqFmbUkjSTV3X03T3vmC+TvuPt2nnEAAMtR+JGFmdUlyd17kkbp9Lz2EBQbYV49zJs7DgBgeco4DXVP46MBSRpKasba3b3n7pthXuLugxzjAACWpIywqEk6ykyv5203sy1Jm7F+AIDlqtQF7nBdYtPMann6m1nbzPpm1j84OCi5OgC4vMoIi5GktfC6Jukw1p69RqHxKad2jnHk7h13b7h74/r16FfIAgBOqYyweCgpCa8TST1JyhwtTGtv6mQwDGeNAwBYvsLDIlycTm+FHaXTkp7Mae9ISsysHfp054wDAFgyc/dV11CIRqPh/X5/1WUAwLliZvvu3oj1q9QFbgBANREWAIAowgIAEEVYAACiCAsAQBRhAQCIIiwAAFGEBQAgirAAAEQRFgCAKMICABBFWAAAoggLAEAUYQEAiCIsAABRhAUAIIqwAABEERYAgCjCAgAQRVgAAKIICwBAFGEBAIgiLAAAUYQFACCKsAAARBEWAICoUsLCzFpm1jSzrbztZtYOf3Yy83bStjLqBADkU3hYmFldkty9J2mUTs9rN7OmpJ67dyQlYVqS2mb2TNKw6DoBAPmVcWRxT9IovB5KauZoTzL9hmFaku67++0QLACAFXm5hDFrko4y0+uxdnffzUzXJT0Mr9OjjPpEHwDAElXqAnc4RTVw94EkuftuOKpYz5yayvZvm1nfzPoHBwfLLhcALo0ywmIkaS28rkk6XKC96e7b0nEQtML8Q312auqYu3fcveHujevXrxdVPwBgQhlh8VCf7dgTST1JMrNapL2dnmoKRxH9tE3S7TANAFiBwsMiPYUUdvijdFrSk1nt4fWOmT0zs48z/e6Go4tnmXEAAEtm7r7qGgrRaDS83+fgAwAWYWb77t6I9avUBW4AQDURFgCAKMICABBFWAAAoggLAEAUYQEAiCIsAABRhAUAIIqwAABEERYAgCjCAgAQRVgAAKIICwBAFGEBAIgiLAAAUYQFACCKsAAARBEWAIAowgIAEEVYAACiCAsAQBRhAQCIIiwAAFGEBQAgirAAAEQRFgCAqJdXXQCAi+nTTz/V8+fP9cknn6y6FEi6evWqbty4oStXrpxq+VLCwsxakkaS6u6+m6fdzNqh+ba7b+cZB0B1PX/+XK+++qpu3rwpM1t1ORfKcDhUkiS5+7u7Dg8P9fz5c926detU6yz8NJSZ1SXJ3XuSRun0vHYza0rquXtHUmJmzdg4AKrtk08+0fr6OkFxSoPBYOb80WgkSep2u+r1eup0OnPHMjOtr6+f6SivjGsW9zQ+GpCkoaRmjvYk028YpmPjAKg4guJ0er2eNjY2prb1+33V63UNBgMlSaJms6kkSWaGS+qs/xdlhEVN0lFmej3W7u6dcFQhSXVJ/RzjAMCFlAZAzPb2tqTxaal6vdyTL5W6wB1ONQ3cfZAnBcN1jrYkvfHGGyVXB+C8GQwGevjwod5++23VajUNh0O12+34gjl1u13VarXj00KtVquwsacZDAZqNBqSpHq9riRJdO3aNb3//vulrlcqJyxGktbC65qkwwXam+nF7RzjKByNdCSp0Wj4mSsHcKHUajWtr68rSRLV63XduXNHd+/eVb/f12g0UrPZVK1WO+4/HA7V7XanjrW1tXViejgc6vHjx9rb29P29rbeffdd9Xq9M487T7/fPw670WikWq2mBw8e6P79+8fhUZYywuKhpEZ4nUjqSZKZ1dx9NKe9nbkzqjmrH4Dz50//5p/1L//+n4WO+fXXv6A/+e1vzO2TJImePn2qra2t49/+33vvPe3s7Gg0GqnT6ZzYWSdJknvn3e12defOHUnjMPjwww8LGTevTqejBw8eqFarKUkSdbvdwteRVXhYhFNIjbDDH7l7etXliaS3prWH1ztmtq3x0cTGnHEAILc0JHq9njY3N/X48WNJ46OOZ8+enei7yBHA4eHhidNO6XrOOu4s826XbbVa0TuizqqUaxaZi9XZeW/Nag+3x17LMw6A8yd2BFCW4XAoaRwUR0dHarfbevr06fGO/fbt2yf6L3IEsLm5qcFgoG63e7wTL2JcaXzU0u/31e12jwOp1+uduN6ytbWl3d1dJUly/G8rk7lfjFP9jUbD+/3+qssAEPzwhz/Um2++udIaOp3O8e2lqeFweHyb6eS1hbOsI3v7ahHjTlvPWQNh2v+Jme27e2PGIscqdTcUABRlNBppb29PDx48ODE/3bEXYTgcHq+jyHEnjUYjra2txTuWiLAAcCHVajXt7++Xuo4kSUpfhzQ+BZU9OloFwgIAKq7sz2/kwSPKAQBRhAUAIIqwAIBzJr0leJkICwCooDyPKE/t7pb/dT+EBQBUTJ5HlGf7pp9KLxNhAQAVk/cR5cs089ZZM/tNd//7ZRYDAEW6yI8oT6ebzaZ2dnZKXa80/3MWu2b2lxp/3ek/lV4JABTsIj+iXJKOjo7m9C7WvLDoSfqZpD80s3c0/ta6nqR9Sbfd/cGcZQHgMx9+R/qPHxQ75i/+qvTun83tcpEfUZ4eVSzLzLBw9++El++b2e9qHBQNjb/2tCmJsABQeRf1EeXD4VDD4VBHR0c6OjrSYDAo9atV816zcHf/mcbfSfHEzPgiIgD5RY4AynKRH1Gezu90Oi/cSluGmY8oN7N/lfRI0lNJa+7+F6VXcwY8ohyoFh5RfnkeUb4jaajxaac7Zrap8XWLgaSau//B6UsGgHLxiPJizbtm8X54+UTSn0uSmX1R4+sVm+WXBgCnxyPKi7XQI8rDdYu/MjO+DxsAluTcPqLc3X9cdCEAgOricR8AgCjCAgDOGR5RDgCQlP8R5ZP90s9+zPoQ4GkRFgBQMXkfUT6t33vvvadWq3Xi8yRFWOhuKABA+fI+onyyX7fb1dtvvy1psQcU5lFKWJhZS9JIUt3dX/gKp1ntZlZ390Fmesfdt82s7e6dMmoFcHFd9EeUT3r69Olxv16vV2hgFB4WZlaXJHfvmVkyJQCmtptZU9KepOxDVdohWPgQIICFXfRHlE+zvr6uer2uXq934tlSZ1XGkcU9Sel3/A01/sT3INYewmPyEv99dy/2Kg2Apdv5xx396OhHhY75tbWvafub23P7XORHlE+TBqM0DsqnT59WOixqGj9DKrW+YHtWEo44pp7OAoCYi/qI8mlardbxekaj0fH1iyJU+gJ3GhBmdsfMmu7Oo9GBcyh2BFCWi/yI8mn9kiRRrVZTt9vV4eFhta9ZaHzhOn08Yk3S4YLtkiQza0s6CqehDiW9EKmhT1uS3njjjTMXDuBi6fV62t7ePvEQvs3NTfV64987z3KxO33KbKfTUb1eV5IkhYwrjY8Q8pw+mtYvXXfRF9vLCIuHGn+jnjTewfckycxq7j6a1T5FX+NrGtL4ovfeZIdwh1RHGn+fRRHFA7gYeER5sQoPi3BnUyNcaxhl7oR6IumtWe3hrqeGmbXcvRv6tc3sSNKz7B1VABDDI8qLVco1i2mfiXD3tyLtXUndiXl8tgLApXduH1EOALhcCAsAQBRhAaA07tx3UhVn/b8gLACU4urVqzo8PCQwKsDddXh4qKtXr556jEp/KA/A+XXjxg09f/5cBwcHqy4FGof3jRs3Tr08YQGgFFeuXNGtW7dWXQYKwmkoAEAUYQEAiCIsAABRhAUAIIqwAABEERYAgCjCAgAQRVgAAKIICwBAFGEBAIgiLAAAUYQFACDq0j9I8L8++VR/94Of6tP/LecxymYljKkSBlU5tUoqqdoy6y3lP60U5b23Jf2MlTJqSdvZOfr5+pUvf0G/9KVXCx8369KHxd9+/6d68Nc/WHUZAHBqf/RbbxIWZfudX3td37y1pld+4eXif5Mo4WClrK+RKev7abykisurt4QxSyr2vH2n0Hn6GTtPP1+StPb5z5U08mcufVh8/nMv6/b1V1ZdBgBUGhe4AQBRhAUAIIqwAABEERYAgCjCAgAQVUpYmFnLzJpmtrVIu5nVFxkHALAchYdFusN3956k0ZQAmNpuZk1Jj/KOAwBYnjKOLO5JGoXXQ0nNPO0hFIYLjAMAWJIywqIm6Sgzvb5g+6L9AAAl4wI3ACCqjLAYSVoLr2uSDhdsz93PzNpm1jez/sHBwZmKBgDMVkZYPJSUhNeJpJ4kmVltXnvecbLcvePuDXdvXL9+vYDSAQDTFB4W7j6Qju9uGqXTkp7MazezlqRG+HveOACAJbOyHp+8bI1Gw/v9/qrLAIBzxcz23b0R68cFbgBAFGEBAIgiLAAAUYQFACCKsAAARBEWAIAowgIAEEVYAACiCAsAQBRhAQCIIiwAAFGEBQAgirAAAEQRFgCAKMICABBFWAAAoggLAEAUYQEAiCIsAABRhAUAIIqwAABEERYAgCjCAgAQRVgAAKIICwBAFGEBAIgiLAAAUaWEhZm1zKxpZlt522fM2wl/t8uoEwCQT+FhYWZ1SXL3nqRROj2vfc4ybTN7JmlYdJ0AgPzKOLK4J2kUXg8lNXO0z1rmvrvfDiECAFiRMsKiJukoM72eo33WMsm801kAgOWo9AVud98NRxXrZjZ5hCIza5tZ38z6BwcHK6gQAC6HMsJiJGktvK5JOszR/sK8EAStMO9QUjK5InfvuHvD3RvXr18v8J8AAMh6uYQxH0pqhNeJpJ4kmVnN3Uez2mfMSy9s35a0V0KtAIAcCj+ycPeBJIXTRqN0WtKTWe1z5t0NRxfPMuMAAJbM3H3VNRSi0Wh4v99fdRkAcK6Y2b67N2L9Kn2BGwBQDYQFACCKsAAARBEWAIAowgIAEFXG5yzOlY/++yN990ff1c//7+eSJJON/zY7MZ012fbCdLpMZtHJtlljTPafN/6sZbLzZ/WN9Y8ue4pl5sl7V56r4H4F3w1Y9X9H3vGypv3MzvwZnNNnclsyW7DPrL6WY/sad3px3mTNC/7czjNt3/FCn5zri41184s39eVXvpxrrNO69GHxvX/7njrf70iSrrx0RVJmgzr+67MNLH2dbpyn2fgAoEjfbnxb3/rGt0pdx6UPi42vbmjjqxuFjTctRCbnzQqjySCaNW/a/BfGnrbMtHFn9J8WgvP+bfOWKfK3NSnfb2xScb+1LTpeXiv7d+TsJ03/2Zz8f37hZzBPH/mLP09Ttoc8286sPvO2o2njFXkUl2esItf3+iuv5xrrLC59WBRt6uF0sfsYAFg6LnADAKIICwBAFGEBAIgiLAAAUYQFACCKsAAARBEWAICoC/PlR2Z2IOknp1z8NUkfFVhOkapaG3Utpqp1SdWtjboWc9q6vuLu12OdLkxYnIWZ9fN8U9QqVLU26lpMVeuSqlsbdS2m7Lo4DQUAiCIsAABRhMVYZ9UFzFHV2qhrMVWtS6pubdS1mFLr4poFACCKIwugIGZWn5humVnTzLZm9J/bXnJt7fBnZ0b/nbTfkuuau95lvWfZusysbmZuZs/Cn70p/Zfyfq3SpQuLKm3AE+utxMa76HpX+H5VagM2s6akR9n6JMnde5JGU3aKc9tLrq0pqefuHUlJmJ7UNrNnkobLqiu23mW9Z1PqWnN3c/fbkjYkTdtGl/F+vbCPWOb+7FKFRZU24In1VmLjnWHlG+8MldiAU+E9yK7nnqRReD2UNPl/Gmsvs7Yks75hmJ50391vh2WXVVdsvUt5zybrmqil4e7Tfp5Kfb+m7SOWvT+7VGGhCm3AEyqx8c6w8o13mipswBE1SUeZ6fUF20vj7p2w05GkuqT+lG5JUb+RLmjeelf2nknHO+wPZjSX/X5N20csdX922cKikhswG+/prXgDPtfCb5oDdx9Mtrn7bgjZ9RlHuqVY1XpzuuPuo2kNZdc9Yx+x1P3ZZQuLSmPjPZWVbcARI0lr4XVN0uGC7cvQdPftyZnhvHgrTB5q+pFu4XKsd9Xv2dTTOMt8v+btI8p22cKi6hswG+/iVr4Bz/Aws85EUi/UVZvXvixm1nb33fC6OVFbP1PPbU0/0i3D1PVW4T0zsxd+flb0fmX3EUvdn122sKjsBszGu7gKbcAKwdRIAyr9zS/8X44yvwk+ibSXXltY5064i+zjTNdsbXdD/2dl1TbjPZu23qW+Z5N1ZUxeF1v2+zW5j1jq/uzSfSgv3EY5lJSk5wDNbN/d35rVvoSa0lv1jjT+TWDD3XtT6joKde0uo65Z6131+5WpLZG07e6bmXkrf8+Aos3ZRyxtf3bpwgIAsLjLdhoKAHAKhAUAIIqwAABEERYAgCjCAgAQRVgAJTOzZBlPvgXKRFgA5WtqeR+kBEpBWAAlCs/y2dT4oYa1WH+gql5edQHARebuAzMbunt31bUAZ8GRBVCicDRxFO0IVBxhAZSrIenxkr9FECgcYQGUa6jxg9/WYh2BKuNBggCAKI4sAABRhAUAIIqwAABEERYAgCjCAgAQRVgAAKIICwBAFGEBAIj6f38oeAxkJ3lGAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f36767deda0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(thist,Hhist1, label=r\"$p_0 = q_0 = 1/8$\")\n",
"plt.plot(thist,Hhist2, label=r\"$p_0 = q_0 = 1/4$\")\n",
"plt.plot(thist,Hhist3, label=r\"$p_0 = q_0 = 1/16$\")\n",
"\n",
"plt.legend()\n",
"plt.xlabel(r'$t$')\n",
"plt.ylabel(r'$H$')\n",
"plt.title(r'Lotka Volterra')\n",
"\n",
"plt.savefig('hamiltonian.pdf')\n",
"plt.show()"
]
}
],
"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.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment