Skip to content

Instantly share code, notes, and snippets.

@Lirimy
Created May 14, 2019 01:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Lirimy/5ad7e50c4fae1ed0700629cf4877ad7e to your computer and use it in GitHub Desktop.
Save Lirimy/5ad7e50c4fae1ed0700629cf4877ad7e to your computer and use it in GitHub Desktop.
分岐図
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"gr(size=(400, 300), fmt=:png);"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"using NLsolve"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"using ForwardDiff"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f! (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = Ref{Float64}(2.)\n",
"\n",
"function f!(du, u)\n",
" x, y = u\n",
" du[1] = a[] / (1 + y^2) - x\n",
" du[2] = a[] / (1 + x^2) - y\n",
" nothing\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Results of Nonlinear Solver Algorithm\n",
" * Algorithm: Trust-region with dogleg and autoscaling\n",
" * Starting Point: [1.0, 4.0]\n",
" * Zero: [0.381966, 2.61803]\n",
" * Inf-norm of residuals: 0.000000\n",
" * Iterations: 5\n",
" * Convergence: true\n",
" * |x - x'| < 0.0e+00: false\n",
" * |f(x)| < 1.0e-08: true\n",
" * Function Calls (f): 6\n",
" * Jacobian Calls (df/dx): 6"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[] = 3.\n",
"r = nlsolve(f!, [1., 4.], autodiff=:forward)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.38196601125010465\n",
" 2.6180339887498962 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r.zero"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.338655 seconds (2.22 M allocations: 88.334 MiB, 8.83% gc time)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEsCAIAAABi1XKVAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deZxT1f0//vf7nJttkmF2mI1hH4Z9E5FFEBVRFHFp1Sq1ovb78YP2YxeLVn+tHytarV30Y/18Wq1La9WqrdadioJURFQE2WFYB2Zj9kkyM0nuPef9+2MQkSoMk2SSm3k//+BBJpOb9+SevHLuuSfnIhEBY4zZgUh0AYwx1lUcWIwx2+DAYozZBgcWY8w2OLAYY7bBgcUYsw0OLMaYbUQbWJZlLV68OC8vb/r06VVVVTGpiTHGvlK0gfXggw/6/f6Kiopp06bdeeedMamJMca+EkY5033ixIlPPvnkuHHjAoFAeXn5pEmTYlUZY4wdw4jy8RUVFc8999zs2bMHDx785JNPHvl5R0fHT37yk7179x75yQUXXLBo0aIony5+LMuSUiJiogvpEqUUIgphjyFIIlJKGUa0ja3HcGOIHyISQnS7MUTbhvx+PxFt3br1kUce+e53v7t27drOn+/cufO111677777jvzmoEGDlFJRPl38NDY2ZmZm2uVNFQgEENHn8yW6kC4xTdPv9+fk5CS6kK7ixhA/pmlGE1jRHhIWFhZ++umnBQUFtbW1Q4cODQaDnT//7LPPrr766k2bNkWz8Z5UX1+flZXFbTQeTNNsbW3Nzc1NdCFdxY0hfkzTVEq53e7uPTzabuTcuXOfeuqpcDj86KOPnnLKKVFujTHGjiPawPrFL36xYsWKfv36vfvuu3/84x9jUhNjjH2laDu9+fn5y5cvj0kpjDF2fPY4s8AYYxB9D4sxxrouaBIQdHPInXtYjLEes6FBP7YLdwe6P8GNe1iMsbjrMOF32/WBIOV7yGd0fyoV97AYY/HVFIHHdqq3K8kiun4oDopixhj3sBhjcbS+kV7YS6bGKwbjpUPQC2Y0X3jhwGKMxcueVvifLcpStGiE48wCQADTjGqDHFiMsdhrCcNf9ygQOLWfGJsjpubFZrMcWIyxGNvth598YgYscUEJLC4TInbLXnBgMcZi6YV96rUK0IAz+sL1pbFMK+DAYozFCgHcstZ6pYKGZYr/niin5MV+QTEOLMZYDOwL0K83600NNNAHd08Sp+TGZflDDizGWLTW19P/levNjXRusbh9nHTKeD0RBxZjrPsUwJIPrTeqYZgPvj9KfHOIiFtYAXBgMca6rTEEP99A71SBRLimzLh0QNyfkQOLMdYd/gjcv0ltaKIJOXDXJGNIn554Ug4sxthJe6OS9vgp2wWXDaTryhyeuB4HHoUDizF2EiyCG1brNfV6fCbeOVEMz+yprAIADizGWNdVBGFHK1W2k0fQ90aJ4Zk9felGDizGWJfsaKXvrlJuB/50vBiVLXNcCaiBA4sxdgJaw6uV+t2DOsuD47Nxel+UPXog+AUOLMbY8YQUPFWu3zqg871w3yliZFZPHwYejQOLMfa12ix4YZ/e2EQTcnFRqRiQnsi0Ag4sxtjXefMAvXZQD0qn8/vLc4rRmQQLqnNgMcaO1WbCvRvNj+tAIl491JjaL9EFfY4DizH2JU1h+P5qc3UjDsug+04R43MSXdBROLAYY194s5J+u0UpLSbnqqWTHMMyEjxodQwOLMbYYY/u1A9uJqX1tWXyx2OTYczqWBxYjDHQBB/W0T8PUo5LX1tmLBqWhGEFwIHFGGsz4bNm2tKkLx2AMwsdxd5EF/T1OLAY69X2++mxcp3phPP6y7IMMJK0a3UYBxZjvZQi+OAQ+U1IkzAhV4zOSnRBXcCBxVhvFFLw00+sjxrhykHilnHSldwdqyM4sBjrdRpCtKJaBxSWpesLBwi7pBVwYDHW26yspj/v0sU+uGYontrXEdsLncZbtNE6ffp0/NwNN9wQk5oYY3Hy7G714BbdpvT0PDytX4wvy9wDouphEdGOHTsqKyszMjIAwDC4v8ZYktraSHdtVNkOVZZpXD/cMSwj0QV1S1Q9rEOHDkUikQULFhQUFCxcuNDv98eqLMZYDP3zoL5ylfV+DYAw7p0sbZpWEGUPq7a2dvLkyb/5zW9KSkp+8IMf3Hzzzc8999yRe+vq6ubPn3/k5ty5c6+66qponi6uWlpawD6dxGAwiIimaSa6kC4xTdPv98tELVJ58lKsMTyzz/n4fqMjQvML1H8Pj/hberK6Y5mm6fF43G539x4e1S4ZP378ihUrOv9/3333jRo16uh7vV7vwoULj9wcPHiwx+OJ5uniyuPxeDweu7RRy7IQMZlfz6MZhtHZTBNdSFelTGPoUPDoLvxXLWQ78bZxsKBIAiR4LxiGIUT3D+yi2iXr168PhULTpk0DAKfT6XJ9aVV6r9d7+eWXR7P9nuRyudxut13aqGmaiNjtj6keJqUMh8N2qRZSpTE0h+m/N8L6Jj0xm3481ij2JaS6Y0kplVLdfnhUY1htbW0XX3zx9u3bI5HI3XfffdFFF0WzNcZYrHzaAHdtsBoiek6RXjo5WdIqelF9hsyYMePOO++cP39+a2vrvHnzHn744ViVxRjrHgJ4ZAutOKTTDbymVMwptM3QYVdEFViIuHjx4sWLF8eqGsZYNFoj8NgOWl6tfQ76yVgsS+gVbuLBHkfpjLETCprwwj71WRNMzdP/OdLRzzYnOU4CBxZjqWC3H5ZXabfAbw2BswsdrpQ6EPwCBxZj9tZmwq92SFOqIV6YVSBGZqbaYeDROLAYs7GWCP3Hx46PG3FSLv1/44w820wd6SYOLMbsanMT/M8WrTVMzNQPTnGmfFoBBxZjNrWjlR7fqZoieH5/dcUAyuyTykeCR3BgMWYzJsHvt1N1UJVmiiuHwgi3RuwVaQUcWIzZy7oaeLNG7/ZTvhuvGioyHBAIJLqmHsSBxZhtbGiEK1dbIQvumoiXDZHe3vf27X1/MWM2pAl+s1X9fR9lOOiMfLp6mFPaZyH2GOLAYizZtZlw4wfmh/WY6aCfTxLn9U/RWaFdwIHFWFIr99OvN6otrWJgBjx8qijN7JU9q89xYDGWvJ7apX+7VblQXDoQvzdC+JyJLijROLAYS0Ya4NFt6sHtEDLhP0fAkrGi10xdOB4OLMaSTlMEfvqx2tICQ9L1tcONSwdwVh3GgcVYcqntgHs3qA1NNLAP/GqyUejltPoCBxZjSWR1DT2wSTkFXDoIbhhh9MKZVsfHrwdjScHUcM8Gva6BpIArS/HiAb137sJxcGAxlnh7/PCTjyPbW0WBR//vFMfQlF7TKhocWIwl2Ls1dN0qq0Ph6f3oN6c5S1LlCjfxwIHFWCK9W0V3brAiGiZm45/PMNL4HXlc/PIwlhgdJvxhh/V2FeY64buTHd8ZluiC7IADi7EEqAvDt941d7filL7qjvHOcTmJLsgmOLAY62kb6/Xt63VtGwzy6QemOAfwoFWXcWAx1qNe22fds5EIxDeH4u1jpZNnL5wMDizGekjIhH8c1OsbIcuDVw6Bbw/lrDppHFiM9YTdLXTDh2ZzSFwyEJ+eJXN7wRVu4qFXr63DWM+oCMBl7+nNTZjn0deWcVp1H/ewGIuvxjCsqNH5aTQpj3431eniTkIUOLAYixdFsK5ehbUszcT/PU0O7MNvt2jxK8hYXFS1w8PbVHMI5vXXF5QIyd8OjAUOLMZib1MjPbpDlQdgal+anW9wWsUKBxZjMXbvBvP3u2BQGi4qxSsGG25+k8UOv5aMxdJze9Q/9gEoOK/QuqaUTwfGGAcWY7HREIKHt+oDAZqWL64vk6OzHYmuKAXF4BTrli1bvF5v9NthzL52tcBdG9SHh/SwTPjlaXJ0dqILSlHR9rBaW1uvueaa9vb2mFTDmB19XAfP7VXBCHxzsFw0HA0eYo+bqHpYRHTNNdfcdtttsaqGMXuxNNy8Vn1vrRVScMNI+d0yTqv4iqqHdf/99w8ZMuQb3/jGV95bV1e3YMGCIzfPOeecK6+8Mpqni6uWlhYAMAx7DOoFg0FENE0z0YV0iWmafr9fStt817eLjaHSr7/zaVp1hzHYq64sCpcaurm5R+r7Mts1Bo/H43Z383RE99+fK1euXLZs2fLly7/uF9LS0i6//PIjN4cOHerxeLr9dPHm8Xg8Ho9dAsuyLERM5tfzaIZhdDbTRBfSVV1pDC9UqNs/dQRNGJKun50J+WmuHivvGLZrDEJ0/8Cu++/Pd999d9WqVU6ns/MmIr7//vszZsw48gs+ny+Zu1THcLlcbrfbLoFlmiYidvtjqodJKcPhsF2qhS40hjer6L4tVpspZhTAi7OlJ6HnA23XGJRS3X5496Nu6dKl9DkAIKKj04qxlGQquHGNuvkD1dcNP50Ir5+T4LTqbezRoWAsGTSE4bGt1tpDlGbgzyYaZxbwAHtPi01gdXayGEthNe1w/0ZrSzOeOwD+s9Qo5oXYE4F7WIyd2MoqdddG8kn8xmC8ephMs80Jz1TDgcXY8RDAx4f0qwfBIeC8YrihjNffSyQOLMa+VnMYPqhVzSYOTsebR8uBfBiYaBxYjH21pjAsXGU2hsSS0XhxKfr4bGAS4MBi7Cssr1bP7AVBMKIPze3PaZUsOLAY+xIiWF5rPFNj+MP0X6PkpYOE4NkLSYMDi7EvtFtw32fWykr34ExaMkbO4plWSYYDi7HDGsLwjXf0gSCM7kO3jaEROZxWSYcDizEAgPJW+L/tOhihoZnwy1FtQzOyEl0R+wocWIzB29X0j/3UYeol4/HiEqOlKdEFsa/BgcV6NaXh+X366XLdNw1uGikn5/FhYFLjwGK9V0MHzF9uNYbo7EK8oUyMzeW0SnYcWKyX2tOqr31f72zFAT5cOllmJ2wBPnYSOLBYb/S7bfqxckUaLy6h3003PPw+sAneUazXeXYP/GqziihcMkZ8fwx/mdlOOLBYLxLRsKZWb2rQo7PEwsHiiqE8aGUzHFist9jZpH+3EyTQtH7GkvHAg1Z2xP1h1ivUh+jHn+p3qyjTCef357SyK+5hsdS3skqvqKHSDJici7eNEw7+mLYtDiyWyojgb/v1s7t1tgtvGStGZHJW2RsHFktZEQVPl+vXDqgin7h5tCzNSHRBLGocWCw1bWlS925S/oic1Q+vHyGzeNAqJXBgsRT0r2q9aDWFNV43FH441pB8IJgqOLBYqnl6l3p2j3ZJmF1A/z1JIs+1SiEcWCx1KID/WKNWHIShGeL30+XM/EQXxGKN+8osRZgEN36g3qqgNNR3ThCcVimJe1gsFRwMwh+267oOmpYHv57mKElLdEEsPjiwmO39Yx/9Za9FWpxXTNcOd/AQewrjwGL29uOP9F/2UIYL7xwLlw9xCE6rlMaBxexKEdyzUT2zWwPQDaXiW8NkoiticceBxWypMQwPbFGb6mhMFv5wnDG3iCcv9AocWMx+DgTo7/usLQ1iUAbcPcngpRd6Dw4sZjObGvUP1mqXxP9XJs8qhnRuwr0J721mJ49sgz/s0D4HnVVECwYAz2LvbTiwmD1YGh7bSU/tUqbGG8vEVcP4dGBvFO1eX7Zs2ciRIzMzM0eOHPn222/HpCbGjtEagm+tMP93hxrkxadmSU6rXiuqHpbW+qqrrnrhhRdmz5790ksvLVq0qKqqKlaVMdapth0uf9fa5Ydhfej+yTCIl7XqxaIKLMuynn766TPPPDMYDLpcrszMzFiVxVin8lY49y2rVdPp+eKpmTLTkeiCWEJFFVhOp3PevHnBYLBPnz6IuHr16qPvrauru/jii4/cnDNnzre+9a1oni6uWlpaAMAw7DGoFwwGEdE0zUQX0iWmafr9filPemLnM/uN5yukJjktSz0xMUxBaI5Hff+GG0P8mKbp8Xjcbnf3Hh6DXeLz+YLB4EMPPXTzzTd/8sknR36elpZ26aWXHrlZWlrq8Xiif7o48Xg8Ho/HLm3UsixETObX82iGYXQ205N61At78Yl9os2ipZP0N0oEQM/9sdwY4scwDBHF96ei2iX79+9/5JFHHnjgAa/Xe9111917771H3+vz+RYuXBjN9nuSy+Vyu912aaOmaSJitz+mepiUMhwOd73asAWnv2FVBGFyHnx/lDi7yBnX8v4dN4b4kVIqpbr98KjOthQWFj7++OOrVq0ioueff37ChAnRbI0xAGiJwLdXWdtaUWu6Z6I8u4i/Ici+EO0Y1ssvv/zDH/5w3759ZWVlTzzxRKzKYr1Tc5ie3q1qgjAjjx6aZgzP4Imh7Eui7fTOmjXr008/jUkprJdbc4ie2W1lO8XiMWLBAEzjrhX7N/Y4Smcp7+X9+v6NJIW4bYy4YADyd27YV+LAYon3ww+t5/dRjpN+OM4xfyBnFftaHFgskSyC2z9Rf6sAQ+CSiY6FgxNdEEtuHFgsYTosWLpef9RAw9Nh6Sk4pV+iC2JJjwOLJUaFH366QVUEaEwm3TXZkcOL8LEu4MBiCfBZPd2+zqoKibML4J7JDjefEGRdw4HFetpnjXDzR1bQxMsG0a1jDYOXimFdxoHFetTP1qm19ZTpEpcNlDeOTnQ1zG44sFgPURquWKH+Wa1yPOIvp4vpfCl5dvK4O856gkXw/U/UihpyC/Gz8WJ6Pk+2Yt3BPSwWd/4IXLfOu6OVRmTR/ZMc0wsSXRCzLQ4sFl+17fqez6AqAMMy4M9nOfr29FIxLKVwYLE4eusg3PaJdgCeV6R+OsWTxs2NRYfHsFi8PL5DXfmeuTeoR+fQfw0NcVqx6HFgsbjY3Kh+tUlbBOcWGX84DRw8NZTFAn/qsdj7615reSWNzaHr+zl+NBpNs/tL4jJ2NA4sFktEcPsn+q2DUNJH/HKyKMvk6QssljiwWMy0mXDFCmtDI+V74K4JRhlfppLFGgcWi42WCJz1amRHUBR44MnT5ZjcRBfEUhEPurMYqG6D85dF9rVBP7d68WzHmFxuVywuuIfForWpEa5dFa6PiDHZ+vmzXfn2uKAnsyUOLBaVT+usq1dDfQdOydXPnenyORJdEEtp3HVn3beyVt+6jkIWnpUPf5/DacXijntYrJte3KP+6yNtIN42mm4cw18RZD2Be1isO96o0D/fQATivP76xjH8scd6CDc1dtLeOAh3rCfDgF+MwkUjuG/Feg4HFjs5Sz5UKw9RrhuuH25cwZcRZD2LA4udhKtWmm8eFBlO+J/TcCqvccx6HAcW6xKt4cpV6s0D4JP08FQ5ldc4ZonAgcVOzFLw44/Mf9VghgOfnCnOLua0YonBgcVOoN2Cc9+MbGuVw/voR093jMpOdEGsF+NpDex4COCmD62trcLrpN/NMDitWGJxYLGvFbLgmXLllnBqX3ztbJyQw0eCLMH4kJB9td0BmLfM1IQ3jpS/nYIuXuOYJQHuYbGv0BqB2W9Y+/wYtvTCIZxWLFlEG1ivvPLK6NGjMzMzZ86cWV5eHpOaWGK1RmDqa2ZzmPLc6vXznHm8XAxLGlEF1oEDBxYuXPjYY4/V1NRceOGFixYtilVZLFH2tsKYv0UqArq/Fzdf4hqXleiCGDtKVIG1d+/eK664YurUqR6P5zvf+c7OnTtjVRZLiNp2mPe2VRsSuS7x8QIjh/tWLMlENeh+xhlnnHHGGQCglPrZz352+eWXH31vfX39pZdeeuTmnDlzjvmFpNLS0gIAhmGPsxDBYBARTdOM4TZrOmDee+7GiByURi9Na7eC0ByjLZum6ff7pbTNSBg3hvgxTdPj8bjd7u49PAa75J133lmyZMk555yzdOnSo3/u8XgWLFhw5GZpaanHk7wf2R6Px+Px2KWNWpaFiDF8PZWGP5eLVgv6uunduZQb0z1lGEZnM43hNuOqlzeGuDIMQ4juH9hFtUuI6Pbbb//ggw/++te/lpaWHnOvz+e7+uqro9l+T3K5XG632y5t1DRNROz2x9QxGkJwzjJLAX1/tOOmkZATm61+QUoZDodjVW0P6M2NId6klEp1/8K6Ue2SNWvWvPzyy2vXrjUMIxgMAoDP54tmg6znHQrB5e+G9wRksRdvHQdu2xy3sd4oqsB67733du7cmZX1xZkkIoq6JNZzQhZcszKyrVkOy6B/zHFwWrEkF9VZwjvuuIO+LFZlsR7QFoFZr1nrmnFQH3jmDEexN9EFMXYiPNO9lyKA77wf2REgl8BnZhnD+bLyzA7sMazIYm7yS+HtfjnAq/8y0xiakehqGOsa7mH1Rt9eqbb7DSHgt1PkxH68BgOzDe5h9TqXvmO+WQnpDvzzTDG3P39iMTvh9tq73LMhsqwSgehnE+S5nFbMbriH1Yssr4SHd6BT4i2j8aaRfCTI7IcDq7f4sFovXkMOxP9XBndMcCS6HMa6gw8KeoWKIFz9gWruUPP6088n8/RQZlccWKnvQBtcuCwSCItpBfB/U7lvxWyMAyvFKYIFy63yAPZ167+d5Yzie/KMJR6PYaW4C9+2tjdDplP/9WyXk48FWZyFw7CuEcr99E612tWqdrdiOwEROgwRsgiUJhRjs2HDJd3cPgdWKnt4s/7gEHkM+P1U50j+8g2L2v56+KhVf1Bv/auG6oPQgRSISNIEEkEToQAiQA0kAAAQgCQggdYWaSABiEAEovtnqDmwUtarFfq+zZQm6GcT5UWDeRIDOwGtYWcAlh2ItJnGR/Whjc1GfTtamhBAE6LQmgQAAiAQAAigzpxCAA1KgCAgAiREwxAWonAAGmi6BY3INeYWwsxC2VfK/ummUt1fJYEDKzUFTHhuryWEuHCg4waecsUATBM21MOWNlpdZa1v1v4O9GtoDQsCAETSGhCBNACCUAAO0AhIQAhAgEgkABCBBGiUAgW4gHId2uWhIR4xu58+v8RZmgN4uK197ehDlCs5c2ClIK1h4XvmoXb88Vi8qYzTqrdoDcH7NbCqmja0qoqAVdshvAL8FkVIEgAoAgGgCYQAEKB153EcIABQZ1QhIqBIM1QfQQrRZ2i3gIsHu+bk4aQiSEuCMVAOrBR0+ydqbR0WeeGaoZJPC6YMpeD9GtgTpA+raV1j5FBYBiJkobC0GzQRRg4PDxEBCUAJBO2oQSMIAtCAgCjAEA5BgnSmh0qzSWjHBf1xUq44NRuOWmM5CZLpa3BgpZp/7NV/2kNOAf9zmuzjTHQ17GSEw7DBD+8d1MtqIgfa0VC6ok0oNICAOg/WUAFIIAIQQPD5ARgBCgBwGBoUgoEuVOmGNISeV+wYn4eD3Di1AHyuxP5xscGBlVLqQ3DbOh2xaPFomJHPB4NJx7Lg43p4u8L6yG/tbsQGC4OmkGApcGgiODyQRCCMwx0l0IAaQAAIAEKSgpTLiU6BXkH9vTCur5ye017mEZNKjr6cQspOD+bASimXrTArO+C0PPzpRN6zCdPUAW9UmO/ViY0Nke0BYSkhBFkWagAgOHzUBsbhkSPUGiWQAkBAgUCIKA1yEIzMonTpzHWqBUXywmHgPdxfPnbPBgKIveaziZt16lhTp6uD0M+Nj87g3Rp3G+vhnSq1uZk+bFRNHdhhCSEpotC0CFADCtAKpAM0AGpQAKABBKBGLUGAIckjdJ80zCC8dJA8M09MKQLnV+833ptf4NciRbRGYNkBPb2fuH6EGNwn0dWkBKVgSwusqtH/OGg2tENVB/ojAjQREBGCQNAIAkALQALUoASQAgGIhlBaGOgUWkqd6xQXDXScmY+nFsT+mo+9DQdWijj7Lau6He6aJKb17TWHBzGyownePKjW1YkNrVaFn0xNCH00EmkTBBwe4UYCAkAN2DkVQCBpQEwzyCshy41D063vDHaO62cMPfxpkbwn2myNAysVvLCbtrdoJ4q5hZxWX60+BH/fDW9Wmp81qmZTkIaQlkCaEA8frEEEEAAkIABZQAgIgsDpACLKdVGhl/o5cf5A4/wiR0H6vz9Dyo5zJxUOLNvTAL/dptIM8YNR2L/XX3h7RzO8VmG9VaV2+6E5jCGFREBIoDtnbMPhCQFAIKjzP0gopHYJ6XMqn9NKB3HtwMAFQ7MGZvK7I+nwLrG9J8qhzYJZ+XTbuF60N/cF4dmd5pp62NpCNe2oAIGAiIAIBABIAA2EoAEkAAlAQkBEdKBOM2hSrriwRIzLEjOKjt7q4Rewvh6yen30J6de1MRTUmMInt5lOSX9ZJwjJc9t17XDW9Xq9f30SaM2kBo7RMAC6vxCSefaAIcHmDSgACCUAjQ4hXI6cWg6DUqHy/s7Lxz0dSfgmM3wbrS3J3ZG6jtgVpGYmJvoUqLjj8CyfeHXq+Rb1SoYQg2gEKhz5qQAsAAEAuKRL9cKiQ4gt0FpkvI88pRca+EQY2bBkfu5Yacm3q82FjThxf2GRH3rKDvtx01N8K+D6uVK/VEdhQmJEBG1hsPrBJAE6BxdUoBoILokuJzgRDqzSMwrFiP6iPF9/32rfFauV7BTQ2fH+OdBTajOKsKByTrxatUBeKFKr6yxQiFfddiyFBEAAMLhfhACASCQVogIGrxSaJOy09S8IuPiYvcZ/cHJX4dkR+HAsisCWFlL+S7xvZFJ0bmobIPfb4u8VKEjSgRM0dxBCgCkBi0+b2YEhCAIERyofZIQobgPXDFYfmeos5/n6I3xFAH21Tiw7GprM7aEaVaBHNLj3auqAPxlj15ZZ+5pxdo2sAAjCgEBCUk7QCgEIkQgjSSkJIegsT6YVmhMypXfHAQGNzrWXdx27Oq5CtwdgMsGqrjuRE3wbiX9sdzaH1AHOiRp0WGpNguoc8UlhM6v7yKhITDTaXUoOT5XX1TiurAEh2Qc3ohpmq2trbm5Nj8vwJIAB5YtKYKP6pEEzSqM5dFTUwd81qCXVelPGow9gYgCoQiaQ9oi7FzeRKDlktIpUQJlunWJR59R6Jw3QEw/PArOh3IsvjiwbGlNvWiMiLIsyIhiVbaqdlh+UH/cSFsarTrTaAjqgEIkQgRLWYjCYVC2A3NcGCLon6ZP72ucP9A4twhScsIXswUOLFsKWVTiwwtKTmK4vSEE79fq1w5YH9dTTQdKkEFTm4QIhEIYaCEKJ1JBOgxOE4OzYEquOLsAiwG2Q2kAABUxSURBVL6Y8M29J5Z4MQsspdSoUaN27NgRqw2y4yjLguuG6PMGfe2C7R0adjeq9U345F4FGioC2BQGE8hSQhMgokfodAM8Tir1wuR8ObtAzuwrnBxKLLnFJrAeeuihZ599dufOnTHZGjuhtgj0T1PpR81Rqg7Cx3X09/1WGHCPX+9uRafENKnrI+AEnebAPg7IcUGhh0Zm47R+eEGJgy8EzWwnNoE1duzYIUOGzJ8/PyZbYycUjGBdWL9TSQVe/N+t1r9qwSLdQaKlA3wO8rnQY8AArz5/gAyZel5/16l5/GU6lgpi04pnz5797z+sr6+/7LLLjtw866yzjr6ZbFpaWgDASNY5Qp/Uwvc2O/f7DQLQJEi4HWBOzg1NydJ7/NJUkO+GsgzlBJjdLzI1T3u//He0BaAtQZUDgGmafr9fStv06JK8MRwjGAwiohnlFUp7immaHo/H7e7m0qtx3CVut/u88847crOsrMzj8Rzn9xPL4/F4PJ7kaaMr98PPt6mtQaM9LEB0XuXp83UvgYDI48L5xeLCEsjzgBtBSoDDX3hJuss5GYbR2UwTXUhXJVtjOD7LshDRLi+vYRgiiotlxnGXpKenL1q0KH7bjy2Xy+V2uxPYRj88ADeva9vkNyxLAAChAOrcrxpIgkDUJARMzFC3jPZUBoNBZfxookvaYYaBlDIcDnf7Q7XnJbwxnBTTNBHRLi+vlFIp1e2H22OXpCSl4NrV5j/2qXZLUGfPCQxAeXgdcSJESJMwIQduGWfML/nSYx/fjLuC0BCifh47JBZjMcKB1aPe2A13bAxta0INSNi5VoEBYIKQQICEUqgxOeKuMY7zBx9vO00R2NKK71Xry4fYZmCIsejFMrCoc+0QdhQieGSr/tUWs6oNEUFR54Fe5/ri1Lkg3QAP3DnW+e1RJ7HZqbn0URMcCnH3ivUu3MOKvaY2WPiv0MpaYWqkzmV8SQAQACCAEOCWalaheH66My2t8xEnPV9zbDYJoOf20NxCPTyr+0OYjNkLB1ZsVIXg1g/NF3aTBkGooXMxA1CgNSJmOahvuvXklLRTi068qa5AgKHpuLpe/+gjeP1cDizWW3BgdV9rB8xfEfm4DiwlCVXnlew6LwWMgiSqmX3puTnO3MPTDGK8dOZ/DDcf3eVeXoP3b9G3jubMYr0CB9ZJu/XD0BPl0GwahASEgASgQSMK8Eg9tR+8MNWZmRn3MnIccH5/eGEv/W6z+tYAoySdM4ulPg6sLnlmM/zXxg5/yNAEhxeuAwUgEbSBcG5/8dzp0uOBHl7S4E9nGE0Rc1WtnP66XnuBKPqKyxEzllI4sL6WZcGZ/wyvrUFNSKAPXwUPEUG5BcwuwsdmGPlpCS7yudmOEX83a4M09iXr/QXGyPj37BhLIA6sY62shB9+0rGlwdCIAAgoADQACsKCNHh6tiO2i3xGyeeA9893THo50mrqyf+wnjodvjmE9ylLWdy4D/vrXvjpjnCbpTQYQA4QBICohGFYVw/BR89Iui/oHTEwHTZ9wznxb5EmC65drf9VG354evJWy1g0entg/XID3bM50mYCUSagBoWASjhotE//7lTP9BKwxUtUlAaVVzkvWWG9W4V/2AlvV0f+MNU4o5iH4VmqscG7MR7uWku/3B4Ja0kEgAiAiCQAZxfT3053pdtw9Noh4bU5xsOb1L2baX9AzX9HXVBCD02Rfb2Jroyx2OldgbXkE+ux7SoQFoQAIIAABXldeFF//auR/qysLLt8Qf/rfG+svH4UXLZCra7Bl/fTywesqbnqqRmuATwYz1KCvd+fXbS6Gha8HQoQKo2gRedcdJdQC4fhH2YdXpSjvj6xNcaMR8Jrc1zbm/RlK9WOVlxdJ0e+Yp6aRz8d6zizmL97yOwtlQNLa7j2X5Fn96ACAJCABERC0s3D8VfTYzzvPNmMyBabLxVrKvWPNphbGuXaQzD3bcvrwItL6MGpjowU/+tZykrNwFq8Bv64M6yUBEIQAAqE1NMyadUl9ljkLFamFYsPi12VQbjl4/CrFRgw8dk98OK+SLoDLywWv5kuvam5/1nKSqkGGwrBqW+a25tAA3R+qQ8BZ+bRozOdQ3vxIE6xD/56pgsAfvmZ+ss+vbtV1rXrJ/bo5/errDScVyT+Y5gcm5foKhnrghQJrNf2wMIPwsGIBNCgCYSQAm4aIX4zNUX+wJhYMl4uGS9rg7BkvbW+Th4IQlUAH9uu/7BD5bvFjHy8ZABcOEi6eToES1a2fz//53uRP+5FrQWAANIocEQOvH+2MzMj0ZUlq3wf/HmmCwCCCp7cpv5vJx1sEy0mvVVJb1fqH67VDgEDfXTuQMc1Q7Bfor97xNjRbBxYo//esaPVSZ1r46FyCrxlrLh7ko3/oh7mk/C9MfJ7YyQAbGnWz+yiTxtgZysdCmN1SKxtUPeu13keKs2QQ9PF3P7ijH7g5dF6llD2e3tHIjD8bx0HQwZpA1ADUR+X/tts91nFia7MzkZniV+cCgASAFZWqr/tw3WNZnUYg0quPYSravSfdmmLwBA0rS+OzMIJ2TA2S47NTXTdrJexU2B1hKH0xVBNCAkMII0ocjx6/RxXEQ8Yx9TsYjm7GDqvb9gchtcPmCurjc0tVrkfwxasb4KP6vTvtQhbEYcUg3xQnK5n5xsTM6HQhyOz0D7XS2X2Y4/ACoVgyuvhra2SNAIiAo3MMTZdwu+MuMtywbeHOb49DDqX+qoMQGUb/WmP2t4KHx9CEFgdoup22tmq2yIUUZDrhoI0QoEZTlo0RAzoI4dnoj2u8MnswAaBNe0V8+MGIBIAGqWcnI0fXsRRlRjF6VCcjqfldzYbqTWUt+p1jUZ9SK+s0TuatRTiUBjrO4AIK4PklnpUNqw9BIFIn74ea2A6TO4L+S5Z4oWSDOzvBRefkWQnI6kD65J/hl49IAgQiFDA3P74xtykLri3EQLKskRZFgCIH3y+rnwgDO/XUXmL6uOUe1qhX5r6qA7bLdgf0NXtYn0jKa0jpEMKXQK9hs5wwuQ8imiclW9YBA4F4/PUiGzDxZ9K7N8k6fv/oS3mjz4CIgkCkGh8jlh3SZKWyo6R7oJ5/XFe/yP7S/zncLOuqbXJmVsdoKoQlfut8hZR3kIhBe0WRjSuriN/GDc0mIdCwh9RGQ4p0EKAAi91KMx0wJxi7OcRwbAqSRcD03FgOma7wOBvRvY+SZcCmw/BjDcjQUUAAoj6ptHui1xeng1kcy4JY7JgTBYC4NEXEFIKmiNQ30Gbm8Djwtf308ZmkeHQtR3Cb0FLBBs6oFZgh0URZTVEEMhKd0BJuvRK2uknJ5IhRI4HLIVZLn1GPtabYkZfyHFiSEOuC4t94Eu6Ns66L7l25py3QisqDRACANwG/XOOc0aMLuTHkpOUkOuBXA+OyAYAnN8fOqdWdApbUBuCqjaSAAfa6aM6aFOiw6J8D7aEaU8Q2i1sC2N9u/YrEAK3NlGLqf+6GxCxNqAjBGkO7OfRLkNKpCE+OLNItkWsugg2tOtBfeQgL+Z7cGgfbGyH/lmQLsElAbnjlsSSJbAe26ZvXGMpkICEBHdMsO6axCeXejuXAQN8MMCHADAF8JuDjr5TAIAGaA6DPyw+a1JtEQlAHzfqAjcdDOL7EhraKc1BTilCEWrXoAkyGvSuVjwYpJawEIK8hs51o0S9L4A+J/bzQFOYEHBkus8w6NMGc0gGTs3DAwFqNinbDTP6GiXpoBWGFE3MBpMg3YV9DBAIBo+49YjEB1YkAoUvhJpDAkCAgJn5uHKe7OHrZTGbEgA5Lshx4aA+nS0ZF8KR845fRAgRtFsQ0WASNbbD7lZa34wCdMhCtwMiJN6r1OlO9DgwZEKHpoAJDW1UFzZCzbolpCtDFIwIBPrwkOrrpqp2tAj6uXWEJGmd74E2CySIfh4q90OOW5dlyqAFTR1U0Mdwk1XVQeMzcW5/fLNSDM/Ug9MQBaY70SNJacr2CJcAlwF8vrQrEhxY96wL37mRCA1A8Aq9/uJevawCixNE8DrACwCAfd0wIhvnD4KjEw0mfun36+pa0tKz9nUIBwppwJYGvTMAgZAqzZQKYFeAqgI0yCf3B1W7ErkerPJDGKjdUn5TmoSE1BiixhBubTY1YlMYNzfhR02wvl65BBR5IaS1AExzQIVfu6TKScOSNDzURnuCYEhlKWNkphqTI/s5UQgoDxCSEtIYkqaL02WrCWlSIYnakDqv2OGRdMgPYQVlDmgNk1ciIOS6gQQYlIKHtwkLrLY2yHs+FFYSCNCA80tCr8zxJaoYxo6GCG4HjPl8TGLY4atqH+kDHYmBY44DjZAJQoBFUBeG+jYCgKYIbWrAAem6JF08s8fq63bkufXBdkQEFaE0iS4JDoQCLwUiKJDCSoYUVLVhmgOaDWq06EAbBkLCAiryYobDagwhIjWHVZuF21qtNIkb610RTVkeVRGgoAk+J0zOw9oOynFRQwcaSP3cwq/BVFDXrjssLE6HdAM7LJIC5w+QTgkVfnNynoEIW1rAobXHgaWZwlQq0yP2tFJ5qxjTB8qyKd2J/jAUeMlADFngcYCpwClBQOewc09ITGA9tI1+tMYijSAh00F7rnBkOjmtmO25HQAAToCBBgz0duYazi2Gzmg7rW/nQMcJ3tztYXA6pKUhpACADgYwrKhNgVOi05BVbWRpXdcOe/36ghJHh6Yil1XXYRV78VO3o6nDMiQOTRdEKsspAhGFQnZo8ocwYFF1B3aYEFTKJWXAJEHkNEgCVbdhdZsKWrCxQUdAOFH392mPARbp1bW6OWxluODUXOEyqDUsJuaQ1xA7WrVAqgsJSwEJvc+PLkED09EizHGCApqcTbURbOzAXa2UJlVfrzE+C9sJ8510el8c3d2VNBMQWHl/CTeFBCChgDOKI+/M5eu6MPaFNBcAgCHAbQAAZrrgqD4dTMrBz3t2nf/itIwIIvp8nZNFjgz+Gkf9DigFIMFU0GKC0IZJ1BRWpmkUplOHBQfbYFA6+kNqa9ARjlgKRHGaaAmTQ1KeCz+px9JMPTVPCNAV7XhKHpJGS5MhUaI2CSMKayUZiKQhpKE1AhGF1R1Y2QEHA1TZTiiM2jD4TQIil4Ah6TC6uy9OjwbWthaY8LJpaQTQUoiKS50FffiSn4zFXec30qWE/MMJhkWHl8dGABjUBwAAfHJULhwVeQgACwYevZkv+obfHHy802IRAgdCyAJEqG0DKSCiKdOFDSFKA+0Rutt/SM8F1sJ3w8/t67z4O5ZlWlu/wVHFWGpyIgCAxwAAGNiZhoAAkONC0wSlur/leI2VKaXU53WFw5D7p8hz+yUQOSSsX+DY+o2km7re1NSkonkhe1ZbW1t7e3uiq+gqy7KampoSXcVJ4MYQP1E2hmgDq7m5ef78+dnZ2RdeeGFzc/ORn+/evXvv3r0A8O5B8D5jNkcQNOZ7MbTIMS4pV30777zz9u3bl+gquurXv/71I488kugqumrr1q2XXHJJoqs4CfPmzetsvbZgr8awbdu2c889t9sPjzaw7r///gEDBtTU1JSUlPzyl7885t6Zr3ecs8wkpUHA3AGhqit5hV3GWPdFG1gvv/zyTTfd5HK5brrpppdeeunou/TSLR/UGiC1kLjuIseb5/DZQMZYVKIddK+qqhowYAAAdPazjvzc5/NpTyEQGS3VhQ+fecndUFBQMHDgwCifLn6CweCSJUt8PntMB9u2bZuUcteuXYkupEtaWloaGxuvvPLKRBfSVYFA4NZbb+XGEA8tLS3RvLBIRNE8vdfrbWxsdLvd7e3teXl5bW1tR+568MEHt2/ffuRmYWHhoEGDvmobjLFeZPbs2f379+/eY6PtYRUWFh48eHDYsGFVVVVFRV9aC+b73/9+lBtnjLGjRTuGNX/+/CeeeIKInnjiiQULFsSkJsYY+0rRHhK2tLRcddVVGzdunDhx4tNPP52RwRdcZozFS7SBxRhjPSYuM92/bjZpklNKlZWVJbqKE3vllVdGjx6dmZk5c+bM8vLyRJdzAsuWLRs5cmRmZubIkSPffvvtRJfTJVu2bPF67TELZ/r06fi5G264IdHlnIBlWYsXL87Ly5s+fXpVVVV3NkFxcOutt954442hUOjGG2+87bbb4vEUMffggw+eeuqpcXpBYqiiosLn861Zs6a9vf2BBx6YNm1aois6HqVUdnb2O++8o5R68cUXCwsLE13RibW0tEyaNCn5WwIRaa2zs7MrKysDgUAgEOjo6Eh0RSfwwAMPXHXVVW1tbbfccst1113XjS3EZa+UlpZu376diLZv315aWhqPp4i5FStWvPbaa8nfTFeuXHn99dd3/r+uri4nJyex9RxfOBx+4403tNZ+v//VV18dOXJkois6Aa31RRdd9OKLLyZ/SyCimpoan883adIkn8+3YMGCQ4cOJbqiE5gwYcJnn31GRH6/f926dd3YQlz2itfrbW9vJ6L29vb09PR4PEWc2KKZdrIs64Ybbli8eHGiCzmxQCAAAIj4wQcfJLqWE/jFL37xox/9iGzSEjZs2DB79uwNGzY0NjZeffXVV1xxRaIrOoHs7Oxbb701Kytr0qRJmzZt6sYW4rJX0tLSOnunbW1taWlp8XiKOLFFMyWi5cuXT5gw4dZbbzVNM9G1dEkwGLznnntOOeWURBdyPCtWrJg1a1YkEiH7tIQjqqurs7KyEl3FCRiGsWTJkurq6jvuuGPKlCnd2EJc9srQoUPLy8uJqLy8fNiwYfF4ijhJ/maqtb7ttttOP/30nTt3JrqWE9u3b98tt9zS+f/a2lqv15vYeo7vjjvuOGaE9/333090Ucfz6aefHum0NjQ05OfnJ7aeEyooKKiuriaimpqa7jWGuJwl5Nmk8bNmzZqXX3751VdfLSwsDAaDwWAw0RUdT2Fh4eOPP75q1Soiev755ydMmJDoio5n6dKlR94YAEBEM2bMSHRRx9PW1nbxxRdv3749EoncfffdF110UaIrOoG5c+c+9dRT4XD40UcfPeWUU7qzidgk55c1NzfPmzevqKho/vz5LS0t8XiKOInTCxJDS5cu7YE9GEPvvffexIkTs7Kypk6d2nkqxhaS/4UlIq31I488MmTIkNzc3Kuvvrq1tTXRFZ1ATU3N2WefnZGRMXPmzF27dnVjCzxxlDFmG3y5WcaYbXBgMcZsgwOLMWYbHFiMMdvgwGKM2QYHFmPMNjiwGGO2wYHFGLMNDizGmG1wYDHGbOP/B7YWrfTsNXA3AAAAAElFTkSuQmCC"
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 10000\n",
"\n",
"as = rand(n) * 6 # list of a\n",
"result = rand(n, 2) * 6 # use as initial u\n",
"\n",
"@time for i in 1:n\n",
" \n",
" a[] = as[i]\n",
" r = nlsolve(f!, result[i, :], autodiff=:forward)\n",
" \n",
" if r.f_converged\n",
" result[i, :] = r.zero\n",
" end\n",
"end\n",
"\n",
"scatter(as, result[:, 2],\n",
" markeralpha=0.2, markersize=2, markerstrokewidth=0,\n",
" legend=:none, xlims=(0, 6), ylims=(0, 6))\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEsCAIAAABi1XKVAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO2deZwU1dX3z71V3bPvMBvIIoKsSgQV0BgejXGLS8zjkxg1BkXcwYiiGEWMYMSERVl8UJ88PmjQoL5JcIkRNYqoYMQF0AiyjsIswDBrz/RU1b3vHxeuRVV3T0939VLV5/vHfLpruqtvVd361Tnnnnsu4ZwDgiCIG6CpbgCCIEi0oGAhCOIaULAQBHENKFgIgrgGFCwEQVwDChaCIK4BBQtBENeAgoUgiGtAwUIQxDWgYCEI4hrUVDcASRmMMcYYIYQQAgDyBYKkLShYGQfnnHNuGAbnnBBCKRXzSRlj5omldhWzv0CQJENw8nNGYRgGY4xSSilljAGAoij2j5l7hXgtZM78L7NpJj9s34IgDoIWVqYgHUBVVbsVFLsGWTCLl1ndhAiKLWbxMqubsOziPiAkE0HB8jiccyFVlFJFUZxSinCWlJQqi6KJZojX8r8WRYuwWwQRoGB5GekA+ny+JP90NNJjljOLyxnZ2URFy1hQsDyI2aqKxgFMFRHUxyxeuq4rimIfGQg5uIljnd4GBctTmKUq+VaVg5i1TAwRWGTIHkSTQmZXNAgzRIC4DhQsjyAyFRhjiqK4WqqipFv1iRBEgyOhNAyiuQ4ULNcjk6qEVYU3myCy9AgjFINorgMFy8WIgI6wquTth0RPlEE085YIQTRMrE0CKFiuRASqOOcZ4gAmnxhcTsMwRBZuuCAa4JhA3KBguQxxYwCAg0lVSAxYFE0EyFT1qBvK7GaGdDntu0KXMzIoWK5BWFVwZNQs1c1BuieCPYWTn2IDBSvdke6GmKiMUuUNIk9+MkuYnBggCDn5KXOCaChYaY20qtI5/xNxHLP0dGujyXEA0VVC7soz/iYKVprCGBPJCpipgIRE9orI9TbMbqaQNlc//FCw0g45AZAQgpF1JDbCuYeyDloqGuUAKFhphJQq8QwUo4EIgkhQsFKPSP4Uzz1Xm+sIkmhQsFKM2apKdVsQJN3BmyQ1SKsqzSvAIEhagYKVbBJUAhRBMgEUrORhKauQ6uYgiPtAwUoG5mJVGKtCkJjBmyexSKtKUZSQCX4IgkQPClaisM9VxnpVCBInKFjOI6UKY+oI4iwoWE6CFWAQJKGgYDmAeY0DlCoESRwoWPEil4BHBxBBEg0KVuyI4T/AWBWCJAsUrB4jU9UZY6qqYrICgiQNFKyeYZ6rLDzBVLcoLjDTAnEXKFjRImNV3pirLOxEQogYK5CFwzOhLjjiXlCwukdaVd6IVcmlDQRyTNOy2IF9JVH5MfPnPXBCEBeBghUJOVfZM1aVkBghUvYjirAl5PIt0kCTH7a7yR44b0j6gIIVGm9bVTHsIZoF8iwGl3hrUTTLZ9BMQ3oECtZReK8CjMWqSiiW+FdIRaOUmhd0sS8dat8njgwgEhSsw3DOdV330mKlyZSqKImw3qfd5RSIuqxmTxODaJkMCtZRJUDF4lqpblFcSLMlraSqW8Kt8Wl+hEQTRLNLmNsvKGImowXLPFfZ5/O53fWIP1CVzjgbRDPLOpppLiJDBctjZRX4EdxlVTlOSNvKkSAaKlqakHGCJcLq4C2pAgBPWlWOIyXMbodGCKLZzTT5ArUsybj+jjWzZ8+eG2+88Sc/+UlNTY3lXyJQpeu6YRgir8rtaiUL2girCm+bOLFrmUystfyVXxHnn5mQq7dJAOc/OYp3LKzdu3ffc889S5cuJYTk5OSY/yVn1XjDqoIjj320qpJJNEE0AECXM6F4R7AWL178yCOPlJSUyC2tra2tra3l5eXglQow5n7vgcPxHhHyNgQhXU4Io2iW2Z0IeMklPHToUGVl5auvvvr66683Nzc/8MADAwYM2LBhg1hZy+1X3ewAolq5l5Aup8ynCWn+MxuZ7G96xMJqbGzMz8+/5pprzjvvvNra2q+++mrfvn1DhgzxQHBH9ktvOLNIOCJkookXUU7n9PbkJ48I1o4dO1avXv3Xv/519OjRnPPrrrvuqaeeuuSSS1LdrrgQT1FvpLMi8RNDJhqYXE6ZfebqvuQRwdI07dRTTx09ejQAtLa25ufnp7pFcSENfs+MEiDJIbKBRik1m2NuxCM3g6qq48aNE69XrVp1/vnnp7Y9MWOJVaW6OYh3COdyuguPCNaIESN27twJAJ9++un69et/9KMfpbpFPQbD6gjSLR5xCfPy8ioqKiZPnlxcXLx06dJUN6dn4KwaBIkSjwgWANx7772pbkKPQalCkB7hHcFyETJTAaUKQXoEClZScWmxKgRJE1CwkoSMqWNSFYLEDApWwjE7gDj8hyDxgIKVQCwOYKZN+0IQx0HBShSyoE2qG4Ig3gEFy2G8XVgdQVILCpZjoFQhSKJBwXIA81zlVLcFQbwMClZcyFJqKFUIkgRQsGKEH1mOGB1ABEkaaBf0GMtyNaluDhKCxiDsaQOGaSSeA++3HmCOVbndsBKyaxxBFgtPdbscoKYd7lyvr/jaCOipbgriNOgSRoV5BDDVbXEAwzAsRySkyqxZ5qVfwhXnTUPadLj2PX1fG5xSCfm+VLcGcRoUrG7wWAUYuUSjXPcs3HGZF2WRcgZHL4VgFrIIa1sljaABt76v1wfIyb3Z5CHYtz0IXtTQeG8CoPD7KKVRLtEYTUVdeZYsWM6YXdcScT4Zg+VbYXsrjCjhj5/uU1x/xZAQoGBZ8V4FGOEAUkodX58x5KotcGRpKctG+1J64ZzN2Bp53yfsxV1sbCksOU3NUWLYAeICULC+w3up6j21qhyk24CX3eUU2ONohmGAqSxPyL395l/GU9uMEj+5ZYRS4nfwOJD0AgULwOTUeMOq4pxLq8rnS9PIc7e+odQyGTuzrOouv/7at/DkVs4YmXWSMq7CU+uGIhYyXbDMDqAHermQKvDEgoZS0SJYiJzzt77Vpn2o+BR27RD+X/2JYXznlkIoQ88DVzmTyVzBsptUrs5CElLFGFNV1e1SFT37OshVa0mrxv+zv/rbU446aruzaU7dEJ6ylLD0GehEIpOJguWxWJXIVAAAQojP5/PAEUWJzuDejRrjdFihsfwH1p4c2Z4SuR0QKm8Djh4WsP/NnDOchmSWYNnHp1yNlCrhALraQuwpmgbnv6HvbocfD4QFJ/uzem5TRpm3YUndADjK5ZS7koYbOp4JJVMEy57A7Wq8FKuKjQc+Mza18BygvztJKUrYsGCEB5vF2RSaZR4WiBxE80Y/TD7eFyz5VPTAjS0e8marKtUtSg1/2sme2c56qfzPZ5LynNS0IWTehqp+d0NFDqKZv4hBtOjxrGA1NTXt3bvXSyOAclZNJksVAKyr43/YxHL95N4T6ciy9D0P0edtQJggmiXSikE08LBgtbW17dq1yxtXV/oayc//TDd2tPBFXxhZCtwyXLn8OHefimgUTdpi3QbRQv71Hp4VrL59+55++umpbkW8iNovxDRXOZP5ugkufUvnQKeNINcOyYizEfmJa5/zZBa4kG6m27NqPStYbscwDE3TFEVxfAKgS+k04Kdvde1pI2PL+a+Ow7mCAN0ZU5YgmvQ6XR1PQMFKO2SsSiRVoVoBAGPwH6927WqFvvnw/H+oPuy2UeBJ3xCvfLoghv/Mc5VF4gICAH+p4d+2k0IfvHSWrzw71a1BUgcKVlog5ioLqyrVbUk73q3jr3/DzulLbxlOh5ekujVISkHBSiV2qyrVLUo7Pj7An9wGCpBZo+nxpXh+Mh0UrNTgigowKacxCFPfZ+06LJgAqFYIoGAlH5xVEyWHumDKe12Mk3P78DMr8UQhALjMVzLhnOu6ruu6qFaMahWBLga3fGh82khHlZA5Y/BEIYdBCysZyFX/RKwq1c1xAeNXa9tb4OTeMH+8jxBcXxA5DApWYjFXgMGYepTcvUHb2gwq8Lkn+Qv9oKNeIUdAwUoUMlYlh/8yql5VzDy3nT+3i5T5YcF49dSKVLcGSTNQsBwGK8DEwyf7Yf4Xuo/Ab05SfzoADVLECgqWk2BZhXjY2w6T1+lNGkwdqUzKjLnNSE9BwXIGlKo40Tk8+IneosHEKjJtOJqlSGhQsOJFrFWDFWDiweDw+03G/i5+aX94cKwPzyISDhSs2BGp6qqqYgWYOLn0LW3dXjijmtz3PTULsz6Q8KBgxYKwqgBAVVXMq4qTFdvg4wai+MiU40kBrjKPRASDBT1AZCpomgZHpMrthlXKMy0+3s+f3tHVJxeePp2e1w97I9INaGFFi7CqPDNXWS6iJxc+ENgr6kLCisC9vc+4eT0vVpR7T4Tz+7lb+pHkgILVDeYKMN6IVUmrSuSImTPFLGtSiRJd5krh4q1l+YPYyoTv64Bfvs1aDH5SP/WigWhbIVGBghUWs1R5xqqSK0eF/IBlu73GrpgOaVnKBSK6liEVTWNw3t/1Jh0GFZInJ6JaIdGCghUC7xWrklIVZ+Z9t+t9WhxMucWsaISQy97S9wRIRS78/Ww11/U2K5I8ULCOQkqVoigemFUTcj3OhBJufXbZDM75rR9ob+xTChS++iylKv8ojYvZwUQyBBSsw7AjeKNSlXTW4reqnELI0PO7YNUunq3w345VRpSFcDntb8XggJRCp4JoiBtBwTqqWJVYWSvVLYoLs1Sl27FsOsjnfGL4qHLHKPX6oSE+EDKIJuaQy1HLyEE0+yGn20lA4iGjBcsyATDlSUlxYo6pp+FdWtMCP3/baNbYLwaRGaNi3EmEIJp0OS1fsYfV5NAnWmeuI0MFy1xXL008pnjgR0gfB9BCqwbjX9EbO/lplfDwKQkZxwi3bqg9iCb/Zc9ECylhqGjpQ2YJlrirDcMQN3Z63ts9wqxT6WwhXv6OdqADcv3kxR/6lBTd/mZFM9uh0QTRzDsJGURLQvsR8J5gPfnkkzk5OVdeeaX9X9Kq8kz+Zzo7gGauXdf1YS3tk8dXTqTF6T1bMOZMNDloE9LEQ8fTKVxvYlh44YUXurq6LBsZY7quM8YURfGAWlnC6ml+OH/cylZtpzqHRROUce5frctioFET4fxxMaTDTEgXHtJgOqe78JSF9d57740aNcrcaRhjmqZ5Zl3lCGNh6cnOZjbnU64qcP1wuMjrc5sdD6LJxxIaaBJP9aH//d//Pf/88wGgsbHxzjvvfPPNNwHA5/N5QK3MYXW3HMvXrXD+G4bG2KRB8PAYTz0aY8NinQmLzGKmgW12p8VA4zYgk8w073Sjr7/+um/fvjk5OXDkMo8ePdozYXVw20JhjMPlb2p72mFUMZk/3jvdLNHY8zYsFz1c3oZl8pN84THrzDs9admyZTNmzPjyyy/feOONmpqa+fPnX3LJJaluVFwIqRLmoev63Af7oTIPANgLP8xyW9vTmsiZsRY5kwMC5i2u60tmPCJYjY2Nzc3Nq1atev3118eMGTNr1qxUtyguZLdzl1UleXSz8dIefv4xytQRvlyPdDF3ELK3mDfa02jdhUd607Jly9ra2s4444yTTz5527ZtqW5O7MhAFbgnsm5h4RfGvM95SRafWEVRrdINl3YqiUc61HXXXVdRUQEA69evd2mR9TRPVY+SPS0w71NoM+h1w9i4chcfCJKeeESwhFoBwLhx48aNG5faxvQUF6WARqa5Cya/z3J87Kxj6JwxXqgjhqQbHhEsN2LO/3S1VSUwDLhrI+xuYyf3Vp49w5VGLpL+oGClAEuqeqqb4wzn/6Nr40FyYilffpoH5BdJU1Cwkgo/ugRoqpvjGJesMd6qhRwVZo9WSrJS3RrEu6BgJQmRywpecQDNvPwN+8e3DIBfM5icXu2pQ0PSDRSshGOJqXvJsAKAuk644V1GOLm4n/LoBDVz5oggKQGfh4lFWlXekyoA6NDhqre7NCDjKvlLZ+PDD0k42MkSgrmsgvd0SvLrDfqXTWRYEV/9I0xiQJIBCpbDeCBVPUpW7uB72vngQrJ0PORjP0KSAnY0x5AjgB6LqYfk1Ro+Y4NRlkVe/ZHatyDVrUEyBhQsB3BpBZiY2dkMd3ykd3G47FhAtUKSifdtgYQikxW8saRFNDR0wflrOg8GyYXHwL3fw9AVklQy4h5LBGapygSrSvKrd4xv29Vyn77sNFQrJNmgYPUYcwpoRkkVAPx+M9veZBybT144KzsL5wsiSQdjWD3AUgHG1Woly4TLLVJ/wx3XI5/BH7YY5Vnk2TPosLIktRNBzKBgRYU3ilUJ5LoGYkBTLtAit8tPmhdN+KjeeOBzwhi5ZRgf1StljUcyHBSsSJjnKntAqgDAMAzGGKVUVVXDMODoNVrMmJXrQAe/+C2mG3Rsb37PiVTXdbHdbpRlSA4akipQsELjsWJV0gEUUhWNmpg/8+iXhsagOof/8wK/euRkmBeY4qaVkMVae/LswdGWGnhxKRckaaBgWfFYsSqzVPl8sYzrTf/I+OsufnZfZfF41W+S7pCWlFmnLIomByu4aT1Ri6KF2y2CCFCwvsMzpYoFnHPDMESgKjapAoAP6/mru7nOyc1DlfKcnn03Gukxy5llWVD5deG6WtxPD1wgJAZQsAC86AAKqVIUJZ7D2d0Kv95g5Pr4g99Tz6hOiEBEUB9LsUNplJlHBuzOpnk74j1QsEDEXCilHqjlJHwuxlicUgUAzV3wq7XawQ7y435w2YAU3P9Sg8QsArsG2YNoMkvDrmgQRd4Gkv5krmB5rAIMY0xYVaqqxuwAmjnn79rnh+B7JfDIqWma0R6lyxkyiCaCeoBBNLeRiYJljlWlui0OIBxAAHBwUfub1ulbmvV8n2/ROOpzs4sc7oRIwYomiGbeW8jtSNLILMHyWAUYmaounCanXNp39vFndgAD39Jx9JQKL5yoCEQTRLNsCRlEs6TdopwliEwRLC9VgBHHYhiGjO84uPPaNrjrEz2PsosGkP8alCndIyTdWlKW8Jk0dSF8EA28EoJIFd7vkV6aVQMmqyrK/M8eEWRw8mq9sQsuHUCf+H6ahq7SB/vQpKpabyhLEM3uctp3hS5nBDwrWAcOHNizZ4+0qlLdHAeQUiViVYn4ievXag2dPFul80/xbMdIMhHsKbuzaVE0HmamfSYrmmf7paZptbW1HjC/zcGRxEkVAKyuMdbV84os8tT3lYrcBP0I8h12DbJgrqVhVjc51glhnE3i3clPnhWsqqqqcePGpboV8SLnJydUqgDgs/0w9QPODLJoAj2nnwc7uhuJMogGJn/TonF2H9PtQuZZwXI1oufpuh79XOV46NDgv97R9gf4uf3JTwdiXT7X0G1mf8jcWldHSFCw0g5zBZjkuLQ/XKPvbYd+hfTpM7A/eArvRbtcrLUeQwyK67ouBpsUJUmWzo3vaR/VgUrI82fQAhwYRNIbfKKmHm6qAJPoWJWFt/cZK3cQhbAbhtATy/HphaQ7KFipxJEKMDGzuxWufJcYnJ/Xl84bjz0BcQHYTVODUxVg4uHCt/QDQVqeQ1adhd0AcQfYU5NNOkgVALxawwNBWuRnz5yuunp6M5JRoGAlDxGoSsQEwJ6ydh9c/76Rp8Dq/yCnV2EeA+IaULCSgZCqlFtVgjYNfvp2sKVLmdCPjS/HDoC4CeyvicVcASZN0mGuX2e0a0qZn/0Rs64Qt4FdNiGY50lIBzAdSjA//iWs2csrc/nLZ/lyFNb9FxAknUDBch4Zq0pyUlW3bGtiD37KdAa/PoGO6EWO1G5CENeAguUkSagAEzMdOpz3htGm8bOPIbcOw+uOuBLsuM4gJgCmoVUlOec1fXcrKVThTz/ACTiIW0HBihfLXOVUNyc0/28P33iIqVR58GSWjdcccS3YeWNHZiqks1QBQH0H3P6+pgK98jh+y3B/qpuDILGDghULwqoSapW0sgox872/aPUBPqQY/vt0vNyIu8Ee3DPkXGVVVUXQKtUt6oa/1vBDXaBQuuQUsDc2HTItECR6ULCiRcaq0jasbufrZrh1nVHs53PHqGf1OyrDXqSJEUJkUV1RCByiqGOJIKkCBasbUlsBJk4mvKw1dcH43uSa479TK8u6Z+apQlKz5CchlBVmVjS31whH3AUKVljkupgpn6scG1eu1RuDRAX23JmHA+3RLNFoVx+5RXzXsnyL2UaTH7PsHxUNcQoUrBDIEqDpMFc5Nt6v4/9vN1OA3HGCUp3LGTvs68VzOPYC4fYtlrie0DizotmNMjTTkOhBwToKuQKgGx1AM5PeNXSNDi3lD445LAfJUQT7YshwtKIJ+8uyPhXY1na3gIqGCFCwDsMYEw6ge60qyRX/1HcFoEBlb5ybpHV3oiekognsLqfA7nhiEC1jQcGyxqrc3vu/amSvf8MpwG0jlcocNylvuDWp5HWJJogWUsLcfk0RSeYKlr0CjNuTksQRXfSWEWDkrAo+a0y6Z7T2lGiCaBa1ChdEk64ooJnmKjJUsNKnWrFTiHvy6rX6zlbip/DMmRk6BSekbWUPohlHautEE0RDRUsfMk6w0rkCTAyYb7MdzfCXPQoAu3wQL8tObbvSFylh4awzDKKlM54SrDfffHPFihU5OTmLFi3Kycmx/Df9K8D0CJlUJY/lojeCQYP0yYH/+X6GmldxEi6IJrdHCKJJMG8joXjBGxI899xzH3zwwVNPPXXRRRetXr3a/C/DMDRNAwBVVdO8skI0yBCMeYhgXT0/qJE8H/zvaV4LXaUPZutMIC4BNQFHTx4QcVKZ2SdeWwCc1Bk1HrGwNE177bXXnnnmGQDo6OjIy8sT2xljuq4TQjygU2AKuFhGM4M6nPWaYTD64Mn8zP4oWKkkQt6GIKTLadkuviUNOg90XafwiIW1ffv2iy66SLz++9//Pnr06BtvvPH1118HAFVVPeAD2pe0MHPhmk7d4ITw20Z45AnkYYgNaZ1Jk01+WFx0s8spTbbMtM480r+HDRs2bNgwANixY0dFRUV+fn55efm4ceM8MAJojlWFlN1NB+CtWgoE7h7Fs11/uBlNyLwNy3W3GFzRZKJ5yUzzWgefO3fu1KlTi4uLH3jggeLi4lQ3Jy6kVRVBrQDgJ291AVcKFfjtWBfPJUKixG6d2YNoFoEz22WWmequwyMWluCf//zniSeeWFlZmeqGxEu4WJWdtbXGvnaqEvbcRE9dSiRmQmaiSdzuc3inlxuGsXz5chF3dy/mCjDRxCYue1PXASZW0nMHhEggQhCP4W65NbN8+fJJkybV1NRMnTp1zZo1qW5Oz7DME4oy3PD4v3mjRlUgT3wfnUEkI/CIhdXY2Lhy5cpRo0YVFxf/5je/qaioSHWLokV6f6Tnxaru/lhnwE+t4AMLEtM4BEkzPCJYCxcuPPPMM6dNm1ZWVpbqtkSLWapiGMG5dZ3e1sUpkBfPwrx2JFPwiGDNnDkzNzc31a2IFpk+E5tUCf5nOwGg3+tt9MJpg0jG4BHBcotaxeMAmpm/2QjqBiH8tbOynGsdgqQ7HhGs9MecqRD/3hZu5sDphN60V178O0MQ1+CdUcK0xZyt50i28Zv7jPoOpqrw5x/h5UMyC7SwEoi5WJWDEyN+9U+dgTIwh1VZK+ggiMdBwUoIEcpXxkmbDnVBFYDP/B7mXiEZBwqWw8gRwATNgZj1scEZz1eNScfjtXOSz+rhld3wVQf/uK6jppMDEMOgHJT/+4Hv8uNS3TjkCNjpHSP6CYDx8E6t4VPozwclaPceZH0NvFILr+9t2xugXQZp0ygHyoFxoEAIMMYJAeDAATgBwoD4gHPgAJQC5TM/1i8/Dm+TdAGvhAN0WwHGKZqC8O9DVKF83qmYzQAAUNsOz30Fr37b8UUjM4jSrAEHCsA5qEA5ZxwIAHAgFLgfgAJnQIS3ToEBUAKEAjAAhVAO3ACiECFlBifU8AH87Yd4qtMIFKy4cDZZoVse3sS7OC+hpDhjkttf+Rqe2NG8vS17dzPVOeWcAKOcMqCEcwAOQAlwBYACUCAcQJhLDADEWwIKcA6MEMoJAUIZAShWjGFFfFRJ9mll6sRKqOglfg1vh3QHr1CMiGQFe/mhhPK33QYADC9xd0kjCzWt8N9bYVO98f5BLcAoMwjnDIByoAAGAAFeBIQBEAAAwoEaQI6oElEJY0CAgEJUTjXw+/T+2eS4EjqmRPnlYOhXEvpH5bVD3AUKVo8xV4BJtA9owUf0XEWdNtJ9V625DR7eBq/u7vqmxQjohBHKgQOnnHMgHAgFoAAqcAMIHE4P5AyoOLmcMCAUCBiKjw/I5RWF/hmD1PMwFp55uK/rp5CkxapC0s5gW5uiEv6TAc78NKVULCYkMC8JE/M+D7XCo1tg+c5AU1BhQDnnnCgcODADKAFOgKtACXAGoAAwoBQYJxyAMAJAKaiGUZ5Hrzo269LjYHQvAABd1xVFrCGCmRyZDgpW9zg1ATBOPqozdB0UFahDUkkI8fl8MmVMvBCFdOUHxAt51GYt+6yWzNvMXq/tCugK5yBUh4uBNlAPR5QIAdCBU6BAGAXCCSWKohf5oH9O8IbjCv5zABQWOnM4SCaAghWJOCvAOEtdG/GrZECuk2VFIx+XzCn7pA7u/Kzl04asAKOcEQ6Ecw6EAmHAKBBiinZzQhiAQoERQ++Vq57b13fHSDqsl3nHotdhlQmkx6BghcaRCjDOUhNgBFhZdsLrIC/5DFbsDnzZrAS1I9pEFeD5wAEIA0KAi5VYAAglTKcKFPm7zqj2/3oYPaWcE6JyLj4hxjIZADOM79aDsa8NgyBRgoJlJU0cQDtdGlDgqs/hdVL/vgvu+ri9JuBrNzjnCgcOQIArh4PfIp2ScwKMUEIJjCyCH1Ybs0/JzjncEBUAOPfJhaSkEtlXzbMvpSefB/a1rZw9TMQboGB9h7idkj/2FyVDS5WiLKjOjiunYX8QHt2kP7/TqGkh7PA6FxyYHygAocAZACOcAoAC4FO0E0vV537g62ddL+1w8DuyuHdrSZnFy6xr5jgaAGialpWVBWHUDckoULAA0ixWFY5BxZwC7Ovs2SX7qgme3Gb8dTi+zmUAAB5VSURBVLdxoIO0a4RzA4ACVYBwkQJOOAPC/AqUZBkX9cv+3Ym0OIoK8Y6cscjqwzk3DIMx5vf7xVKgslCPWc7sbiaKmofJdMFKXFkFxxmQR9sN9mWjYRhUCe8X1rTCQ5/pb9ax2jYIcuCMAjeA0sNBcaoQgBzVKPSzS/vxe8fkVmT3zMdMQnRPSJWwdv3+bpL6ZXskcmNIZxOOrISciJYjiSZzBcucVCW3pLZJkSnNhv65pEkjNe0w0JQKwDn8o4Yt28o/Pqgf6KAGJ8A5EEVEuwlnPoWW5fMzSsn0Eb4xVeJLsVz3JET3hFQBgFzEuFsi2FMhg2gAINRQrqoNoZzN9H+AZSaZKFiJrgCTOJ6cSN/ey7oY1w3y9NfGqp3sk4PQonPGFS5SMRkQwgiQijz9jCrlrmH+E3oDfHfTfrdQuf2GFHZHuDs/0Rmz0qpSFMWp62IfBBDHqKoiDTVEEE2eIIvpbfmLcpYqMkuwzHOVXdfn/m8z3Lyxs5OpMzdox+QqjRoPGJQBVzjJocYxBWRcBbl2sHJaVYjv2g/Wbk6GCw8lesokY0yuIJs4wy2kbRhNtMvsZoZ0OeVRyEdgugXRDAOCHPwKqOnSotjJFMFKgo2QCNbXw3n/CHRoKiOEEwBQgDEgapPOT+5N833svL7qVYNobs+LN4S0sMQLcaLkUJ20vL4zz2xfiWCdRSDJUhXbdY/wRSle4lypqhp93gY4qmhBHfYHoTHA9wSgNkjyFL6hXv+wnjTpvFc28VMo8JFL+sG1xzucE5N8vC9YUqpc5ABOWN3ycX0WUxTODAA/UAYcgHMCVKHGsGLj2CLlxDK4b7TfqWk6Fhhj9hvV/NZioInPmxXNLmHmLTJWpShK4nzMJIz8ysEBkXhhaYD5hSVvA44207oNonUxeL+OA+MqpdtbeUUO7GiDPc1saws7aJBdLSTYxXJV8FGSrfJ8Fdo0qOvgnICfGP0KlGNy+KAixT1P6rB4VrDq6+t37twpHcBUN6d71uyEn7zfEQxSRgjwnMNpB0AJ4wpho4v1DZflAwCA2tIFV7+rr94NQY3PGO1kbazoo3shtcyiaPRwntfhtzLgLW5XGasyGyMxmGkhjyLRUhVNxC2a+L084TqDgx18azMryiK7W4w2RgNBtqaONHdBoQKcKJ8fYt+2cb9ChxdrBif9ClTC+cEgb9MJM3iBAsU5fHiJryJLr86HCeUKVWguhb65UJELigu6f7R4VrAopY2NjenvAN6zHuZv6TS4yoEDqEABGCcAFHiur/ODCfnDB1FLlYJCP/z+FPU3H+srd2rP7eS3j/TdOjLeLul4dM9uL0gH0OfzhUwsCBdEg55ombD1EveIEt4fY6yngwP7g7C7mVPCd7dDc5BvbQOti+1ph12tfG87KJQa3GBcGZDPCSjZCslT4atWIECrc1hVDhuYT9o16F/Af9YfclQ4oSTYK1vJ99E8FXIVQikQogAwQhQ4fN4ckP40xLOC1bt377Fjx6a6FWG54LWuNXsZA8oJABCgBDgnnKhcv21w9sMTAQAMQxVKYu95xxXBglPVn72rfVRPpm8wVu7Ql57mP6mX7WeiIAkus5QqswPYbRDN4kNFdjnNJlVCBweEVeXzHfUIYRzqO2BvO+wLMADY2cL2d8InDVBVyKtz6Ht1XDegLghNnaQ0hxX7SLvGDSA+yg1G2xhlwPJUKM1SsjicfQz0zeVFOfTEEl+bBr2yoDpPUcNfGfOAgPkFHDmNcLTub9mypbW1tbq6urKy0i3rpZvxrGClJ//1WuffaikjhDMAqoIwagj0ye94e2L+sRVgviKUUl3XzV3QHOmoyoV3zlF++qb2dh359CCMX60VZcPNw9X7R0d7uyZhIEI4gISQHsWqoolSm8Nh5sEBMIlXhCBaDw6BQXsX291sBBh8coi0GUrQIIc69XfqoL6DZVE6pIjnqGR7C9kfAEpYWTbhHLoYPxSEiiAd1ws0BgV+1l+lvbPZ2DI6sJBUZIOfsMHFakk27Z0FeT4lZjc2cvzefELWr18/d+7czz77bMCAAXV1dXV1dZ988snxxx/fo7ORclCwksELX8NV7wR1onCuAOHAKRBGGe9XzHdcJkJQIQJRhBDzk9zc+Y5IGH/ph0p9B/xqLVtbDwc74LcfG498CsNLjSmDya8Gq4oS+gZImlRRSmXSk7OQI5N1AEAJlfgfLohmVv+ABu067G6D5k6+voEV5ZCaVviyiWmMbm9mpTlgMNgboGVZBgOFM3ZIIwRYrxzQDajrIIZBslUwdFqeB7SA983lJdlwarlSSHmvbJZFlUElpH9BuHuMyoYxlqiIm9jhunXr5s2bt2XLlttuu+2ll17Kzg5d2Ketre2uu+7q7OysqKh46KGHnG2JU6BgJZC2Nuj3QnsLyzpyjzCghBI6oJh9fWmPQ+XhnqV9CmDNBVDfDlet7fyg3tcF8OkBcvNB5bYNRnUOL8+Fiyvh6qGkLPc7e8QeqHLwVkm0VEF3a38EDWjXYVsTEMJ3tpO39gQHFSubDtEvDhq5Pt4vl37VwjUABSDAwK/woEY0AwI6z/MRg9OAzv2Utetc64KKIpKrwqAiJVdh5dk8oPPqPGVwPh9eQlqCvH+RrzIXCkKUQSVSj7o9ioQ64++9997999+/e/fuu+666y9/+YvFk7Vwww03zJo1a8iQIfb8lfQBBSsh3P5+cMm/qUEIcD9wBoRSoMfk6zt+lp2ge7giD944LxsANjSw+zfCpmZ+qAP2tJOd7Wx9A7l3C+SqnHAoy4FTe/ETimBCGTmliitKVDPvIApFE6HoREhVlw4HNejSYW8bf2uf1qrR3GzY3841qh6fb+xqhQ8bWJNOclWeT+CbAOnSwa9CmwaFfgho0KqRAh/TOQ8YJFthTV3QqoFfob1zWB6DfgWEMajIg6ZOcmIp7VvAdEYH57ESlfYpoQVHUi2FPcoPL9QDIqTNOSOEmAcPojzqRFu4nPOXX355zpw5gUDgzjvvvPLKK0MaoWbWrVs3YcKEIUOGQHqPqqNgOUzlysCBNh9XFAACAIQbPpU9f4r/4hHAeTLO9qnl9PXz/AAQ0OGZ7WzpFtYQhLYu0mkYGqctLWRXK32eM6BMBVCAUEJ8FAYWssocUp1LunSe5zOOK1Iqc2FgHsmmpDCLMZ33ygdVA7+PqGqIpCEpVT6fr7MTiAKGCqBDwICABi0a6AAtnczgoAOvC8DmJij1G7vbFB8hu9tJMTWKsuHLJq4zXuKnuSr5rJG1amAY0M6JSoAwTilhwJqDCqF6NqWEgI/o/Qt4kJFvA8RgvNDHwa/4CctSaXU+tOp8SBEv9pFtTXBeXzWP8gDQkaXqKb0gVwU/AVW13sOMMcYIACiKL4a8DTD5mxAmiCYiepCwlAvG2CuvvDJ79mzDMG6//farrroqSvV55ZVX7rjjjhdffPFf//rXL3/5yxEjRjjeNkdAwXKGlz6DX2wMGkThhh8IA+AKsDP76a+fnZ+qJuUo/LohcN0Qn7g36jrV/97CP2nUv2mHna2gGRSAaYwzBh0GbGki/26iQAydU2YolIICXKFgcOCcG5z7KPFR3mVwjROV8iyAgmwS6OIMeJ4fWoOgUtIvn/XO5l8dYk0aIYTnqcQA0qkxHShw5qMkWyGE8KDGA4yqlBLOCAWN82xCSrKgqYtwzkuzoTgLDgbB4JwqiqqzPAUKcniRDwqz6I4WnkvIqFLWO0/p5WdnVvtyVOCM5GbBMXmQQ4HSHidz8yOp6tCTRFb7AJx9t2ZFi8bPikfCNE1buXLl7373u7KystmzZ1944YU92tvBgweXLFnSt2/f6dOn33zzzatWrUrPrAgUrHh5+FOYtbHLAA6MAOGEcl8W3/qTrH55KWtSyCGnymyYPZZYUro0Bl8cgs0HjVaN1XdAfYB+2cwPdRJGxLw70tLJOw0eMGge5QqFVg46IyoxVEpUDoSCQhQfBYVwnwJ5KivIUvKzoJPxLMrLsgkn0OmDLp37/VBASHmOkedXGOPfdLAhBUwzaN9c6CBkSCEdUUTbuqAsS6nOh6pcwriSp4JCv/OcEnSuZMTN8Zx7qWjyWkhjx563IYht8lMwGHz++ecffPDBioqKRx555KKLLoqhtV9//fXgwYMnT54MAMOHD9+3b1+fPn1i2E+iQcGKnUlr2XM7CAMDRHVzykcUdX5+WaTyd+bB+AQFL+TOo9m/j8LoMhhdpgBEa5hwzg3jsPkQJjIS2gfhnB/5FcWWPcQADhd7kc3mnIsJQpCA0yWtKkJIcgYHImdpmLdb9iBehCwj0d7evnz58gULFhxzzDHPPvvsuHHjYm5qIBC45ZZbxOu9e/dWVYWaQ58GoGDFwux/dc3dTEQ9JTA4VeDqY4JPnZMfMjtBQo4UP4Aj4Qzzdvsne4TZqkposSoeawUY+x1r37/EvDHkzDs4OiutR4pmHsfs6VFEiVNhdbvLKV60tLQsXrx48eLFp59++s9//vMf//jH8agVAAwdOlTkkdbU1OTn56dt3B0Fq2dc8y48s1VjBAA4EFCY8cfx/itPgMhSJTF3X0sE1+wjRHANQm5J2/ko0RPOGLF8xvzXrGVmzyvcX2lVJTrlQrxI0IlqaWlZunTpggULJk6c+Pbbbw8fPtyR3ZaVldXW1hJC7rjjjiVLljiyz0SAghUtL++Cn/6zy+AEKAAHlbLnz/BdPIg40i0jRHAjxDvsG+XXnfKhIsxHcYqQEbeQhPOkLHsDm7Emj0LUsTEP1UW52+gPxKm92dm/f/+SJUsef/zxc845Z926dc4mqd95551z5szx+XyLFy8uLy93cM/OgoLVPfUd0P+5Lo0ZABQ4p5TMGOmfOw5EgZREE/IGsKRNkh7OvDPvORwhJwA6S/RSFT2W0yXcWEKIz+ezFIeAo4JocS07JqUqQVZVfX39ggULnnjiicsuu+zjjz/u16+f4z9RXV29dOlSx3frOChY3VC1IrA/qHBOgfooGFNHkvnjnavn0nPMrpN5e+SYrkWtuK0gn1nRhFSJu7TbhMM4j8JBqbL/RMiIW4Sfs7iccrvZmLUbZTyRS8Pt3r17wYIFzz333OWXX75ly5b0HLlLJihYYZn1Mfzuc42BHygnHMb21tdfnMrV1eWNFGfAG8IrmrSqpLEQeWQgBscz0cYIxDE40K03Zw+fic/ruh5yJ/G4nDt27Jg3b95LL700efLkf//73716xVSLw3OgYIVA06DoT4Gg4QNOALiPa53X5VoymJKJU0NOkX9CWFX2aHSEIBr0pHxVEo4i0TWXzaH9cGH7kEE0+4CA3KF9++bNmx955JE1a9Zcf/3127dvLykpcfxA3AsKlpXL3mr6y+48DlnAOQX25Hj/r0aElSqLs+B4wDWZN3m4WFW4gzLHiSIH0eQH7F9x6qCSUx4eoshxi9AH5CmSu7JkV23YsOHhhx9+7bXXCgsLr7/++hkzZuTlpS7/OC1BwTqKvD8GgyQXOBDGBpYEv74s0sQaRVHMYSAeKmMIjg42mW2QyC1J9Og4j2k+SkgiGA7mU9FtsRf7bqNRtBQuutNTIlz9d999d+7cudu2bbv99tsnT568d+/ehoYG+7As53z16tUXX3xxzG1wOyhYh3lkI9zzmc45BQYKhbXn+8b17d4HjNB97RFcywMWwkc6un2Mx0ni5qNIQkbcItyxPXU5xX+TI1UJtXDfe++92bNn79y5c9q0aS+//LJ9MQsz//jHPxYuXIiClekc+3xwTysF4EBhRAnb9FMHguvRuIdmQ8Ns8lhuj/gjuPLnUjgfJTI9cjlFIVbzT8ixAvndOP3NREsV5/zll1+eO3duW1vbjBkzrrjiimjS7p999tlBgwY53hgXkemC1aRB5QpNAwAClPO/nek/f2Ci7Bo74maQN6TostLIkh8zixocbaBFjuOaSY5VlYib3HxE0gG0CK7FQOvpsmOWo5AfSMSJEhVgHnjgAV3Xb7/99miKVQk2bdo0cuTInTt3Ot4kF5HRgrVoM7/jX4Z4VJdks/1XJDvBKkJwxK5B9u+CzdmEUElDMlVdSFWC7AXxInHJCpHLw1u22FWbR1Ex2Sy44ISZZoEx9uKLL86aNau0tPT+++/vaQWYpUuXPvTQQzNnznSqPW4kcwXrxBcDW5r8wBkhfNpwOn9Clj2MkjjMUhXbLdGtb8g5NwzDMAwR4pEWhz0qFI/LKW/yHn2rRzhSc9nSQst++BHMguvIsmMCUQFmzpw55eXlsVWA2bt3b35+fllZGWNsxYoVl156aX5+ykqtpZAMFayS/+to0VXgBiXw0YVZ36tI3k/HL1XR/ESEFYktnwTT7Wo2OiI7m2Z7JNHrgyUt4mb5V4QgWvSTn4LB4NNPPz1nzpzjjz/+6aefPu2002Jr55IlSyZPnrxixYq1a9defPHFGZvukHGCVd8Bff8UZMwHlOdn8eZfHnU/yykpUcY7ekTSpAqiHjiLYBnZnU1+9ERicjQQ3n6JgbSKuEUfKJR7BoC2tralS5cuWrSooaHh0UcfnTp1asxNbW1tfeGFF7799ttJkyZNnDjxwgsvjHlXbiezBOupf8MN72ucEKD8uOKurf951GOKHJk6F2W8I3pFS4Ixwo+sAeFUBZiQYSBxFEJE7EG0kEXm7H8jq4O0DRNtVTl+OWRrm5qaFi1a9Pjjj5999tlr1qwpLy8vKIhU1rFbvvrqqxUrVkyYMAEAnn/+eQfa6loySLAe2Ai//UwHIISTOSfwu08Ja1R3ay9YtCmkopn/Fc1jPGbMmZPJrAATTdzKfGbsZ8myE7PgujSv6sCBA4sXLxYVYNauXTt06FBHdnvyySfL15MmTXJkny4lUwTr+y/DB/U6EE4A/vYfvgviy2WxhyrkaxnvMIuXZYvlWzH7m+kzHyUcEb5lNmBFXpW0cMVYQbdBtJ4eSEKN3IaGhvnz5z/55JMXXHDBhx9+mLhsqfHjxydoz64gIwSr+I+drUwBAn6FH7ran53IRCtimo9iv9kg1giuZYuL5qOEQxyRiLipqmr/iQhBNLvLGcEoTrRUmSvAbNq0qW/fvon4FUTgfcHKe7qzkxEAXqBC09WJzbSS91iEeyOaCG64IJosVSxuP3NFupjNtJBHkWg3NpqIW7eWVOQgmvkoxMmMMogWPTt37ly0aJGQqs2bN1dWVjqyWyQCnhas3Lz/bDifAQOqlGXpDVcmsJqV1JQezUcJR0h7QY7NhQxUOZU0JMdJ0z/iFk7RzLYh2FI3wCTxIXcVjcu5ZcuWefPmiQowW7duLS0tjedAkOgJUefIM6j/3caVbAB+Qqn+6U9zEvQr8kkOCcuc7FG1YrvLad9ilzDz2yRIVRLc2GiOwqxl5i1gOiGW07Jx48a5c+d++OGH06dPv+mmmzIzezOFeFOwNm3adMkll+y54S0o6ZP3z8f0v86qrKysqqoqLy+vrq4uLy8Xb3v37l1dXV1RUZGdHYvxZYmkJILI81FixiJnIuPJ/IFug2g9/bkeZYfF9hM9kqoo9wmmc7Vu3bqHHnro888/z83Nraury8rKOnjwYOIuPRISbwpWbW3tddddt/jx5VXVfbIV6OjoqKurq62t3b9//759++rr6+vr62traxsaGsTbrKwsKWcVFRUVFRXmt+Xl5ZbpqUmwqmTmpCM+ZkgiRNwsRpm9k9ibFNJq40dqhCUhUyFkqxxhzZo1c+fO/eabb+6+++6rr77a7/cDQCAQEAv5STZv3vzYY491dnbee++9zi5pg0i8KVg9pbm52axfFjk7cOBAaWmpWb+qqqqEqAl1c7betpyPkqAFICDifJRovmtxMOV281tzNryiKE4NCFhgRy/Y5Thr1qy577779u/fP2PGjGuvvTZCBZjNmzcvWrRo4cKFnPP77rvvscceS0R7EE8H3aOmqKioqKgoXJof57yhoaGhoeHDDz/cu3dvQUHBt99++8knn9TW1gp1a2lpkSpWUVEhvE7z22hmfsmBszSZjxKOcAIht4hUdc65PIroZ95FfxTyKwnSwVdeeeXBBx8MBoPTp0+PpgLMokWLFi9enJubW1tbm7ET/ZIAClb3EEKE9IwaNSrkB7q6uhoaGqR+7du3b9u2be+88458yzmX1lmfPn169+5dVVVVWVkpY2qqqkYYAXQEeZMnvwJMyLwN+9toMtHkhxM3MvDiiy/ef//92dnZM2fOvOyyy6L8of79+wsP8ZlnnvnZz36WiLYhgC5hcggEArW1tXV1dcLNFOpWV1f30Ucf1dfXA0BpaalZvyorK8XbPn36lJeX9+7dOx730I0RN3sQzRJxiyaI1iM0TVu5cuXcuXN79+591113xVABRuxk0qRJzz77bGxtQLoFBSuVLFu2rKOj44YbbggGg2Y5q6urE2/FYgQHDx4sKysLN8RZUVFRVlYWcv9mBzBBh5D8iFuEIJo9e0O+iKBlogLM3LlzjznmmN/+9rdnnXVWzE2dP3/+mDFjJk6cGPMekMigYLkAxpgIokk5k2MC4m1bW5tliFN6nULdHI+qmCNuCR3HjEFwZX6DZZTTLmdtbW2PP/74ggULTjjhhIceemjs2LHxtLahoeG+++5bvnx5PDtBIoOClWxaWlr+/Oc/a5o2evRoUTAkfoLBoDlj489//vPAgQMPHjwoRzwBIMKYQHl5uRiqjxIRVk+CVZU4N/bQoUOPPfbY0qVLJ0yYUFhYeNddd40YMSLOfd56663Tp08fMGCAEw1EQoOClVR0Xb/88sunTZtWWFj49ttvNzQ0PPTQQ0n43Y6OjkOHDtXW1u7bt0++kG9ramp0Xa+urq6qqiopKREvzG/79esnhgWSaVVBYqTKXAHmnnvuGTZsmCO73bhx46uvvnrTTTctXry4V69et956qyO7RSygYCWVL7/8cseOHbJi5OzZs6+44orBgwentlUAIM0xOSYg3+7fv//AgQO9e/eWUTPLEGdVVVVxcXH8bUh0xK2hoWHp0qXLli0799xzZ82a5eBp55yfe+65Q4cO1TRt6tSpTpXBQuxgWkNSGT58+PDhw+XbysrKAwcOpINglZWVlZWVmdsmCQQC48ePv+66604//XShX7W1tbt27dqwYYNUt87OznBDnOKtJSncgpSqBKVc7NmzZ/78+aKswqeffup4BZiXXnqpf//+06ZNO/bYY53dM2IBLayU0draOmXKlJUrV3pgPlpnZ6clY8MyRKAoikhAs0zqrKioEIqWoAS0Xbt2LVy4UEjVzJkzq6qqEvErHR0dOTmJml2PmEHBSg1dXV1Tpky5++67M8R9aGtrEykaMmNj69atH3zwQWVlpXhbUlJits7MCRwiE62nsi4qwLz22muTJ0+eMWNGuMwPxF2gYKWA9vb2KVOm3HrrrePGjUt1W7rhb3/726pVq/70pz85vmfOuaZpcnRS5G1IObMkcBw6dEhaZ+YENGmsFRYWyj1//vnnf/jDH958880pU6bcfvvtRUVFjjceSRUoWMmmsbHx5ptvnj17dvpP6K+pqZk1a5bf73/iiSdS2xJd10MmoMm3uq5XVFTk5+fv37+/vr6+b9++M2fOvOmmm1LbbMRxMOieVGpqaqZPn/7oo49WV1enui3doGnajBkzli1bdvfdd6e6LaCqanV1dYST1tHRMW/evCeffPLSSy8dOnRoY2NjXV1dMluIJAe0sJLKwoULv/rqK/nW5/Ndc801J510UgqbFI7p06f/4he/GDNmzJQpU1JuYUVDR0eHz+eLUAEG8QB4dZPKr3/9a/NbTdPS84GxatWqwYMHjxkzJtUNiQpN03w+H47TZQJoYSFWGhsbzz777LFjx/bq1auiomL16tWzZs0aMWJEGg60tbe333jjjXl5eYFAoLy8/OGHH07cbCEkHUALC7FSWlq6ceNGOJL+/sYbb3zzzTd9+vRJQ8GaM2fOjBkzRo4cCQCvvvrqF198ccIJJ6S6UUgCQcFCwiLS36urq6+44opUtyU0+/fvF2oFABdccEFqG4MkARQspBvSdinjQCCQlZX1+9//fuvWrcXFxTfffPPAgQNT3SgksWAMC3ErBw4cGD58+KpVqyZOnLhv374777xz9uzZ6TAxE0kciSrvjSCJhhBy9tlni/Ke1dXVS5YsWbp0aaobhSQWFCzErZSWlppXwC0pKWlpaUlhe5AkgIKFuBVCSF5e3tatW8Xbmpqa0tLS1DYJSTQYw0JcTFtb22233VZRUZGdnb1ly5bFixeXl5enulFIAkHBQlzPnj17OOdYTD0TQMFCEMQ1YAwLQRDXgIKFIIhrQMFCEMQ1oGAhCOIaULAQBHENKFgIgrgGFCwEQVwDChaCIK4BBQtBENeAgoUgiGtAwUIQxDX8f9TKFPYun4dxAAAAAElFTkSuQmCC"
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter(as, result[:, 1], result[:, 2],\n",
" xlims=(0, 6), ylims=(0, 6), zlims=(0, 6),\n",
" markeralpha=0.2, markersize=2, markerstrokewidth=0,\n",
" legend=:none)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"stable (generic function with 1 method)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"negative(x) = x < 0\n",
"stable(J) = reduce(&, negative.(eigvals(J)))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r = nlsolve(f!, [3.0, 3.0], autodiff=:forward) = Results of Nonlinear Solver Algorithm\n",
" * Algorithm: Trust-region with dogleg and autoscaling\n",
" * Starting Point: [3.0, 3.0]\n",
" * Zero: [1.3788, 1.3788]\n",
" * Inf-norm of residuals: 0.000000\n",
" * Iterations: 5\n",
" * Convergence: true\n",
" * |x - x'| < 0.0e+00: false\n",
" * |f(x)| < 1.0e-08: true\n",
" * Function Calls (f): 6\n",
" * Jacobian Calls (df/dx): 6\n",
"J = ForwardDiff.jacobian(f!, similar(r.zero), r.zero) = [-1.0 -1.3106; -1.3106 -1.0]\n",
"eigvals(J) = [0.310602, -2.3106]\n",
"stable(J) = false\n"
]
},
{
"data": {
"text/plain": [
"false"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[] = 4.\n",
"@show r = nlsolve(f!, [3., 3.], autodiff=:forward)\n",
"\n",
"@show J = ForwardDiff.jacobian(f!, similar(r.zero), r.zero)\n",
"@show eigvals(J)\n",
"@show stable(J)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.143750 seconds (407.72 k allocations: 67.396 MiB, 24.38% gc time)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEsCAIAAABi1XKVAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd4BU1dk/8Oc5907bmS2zhW0sbalLbyKgCMQSCwpGE/KqRKMxBjGxBQu/xCR2TYy+xuQniT2v9ReJmhgUFVBqpInUpSzL9l6mz73nPL8/FgnyCiw7Mzszu8/nH3ZnmHufnXvmO+eee+69SETAGGPJQMS7AMYY6ywOLMZY0uDAYowlDQ4sxljS4MBijCUNDizGWNLgwGKMJY1IA8s0zYULF+bk5EyfPr2qqioqNTHG2DeKNLCefPLJ9vb28vLyadOm3XfffVGpiTHGvhFGONN9woQJL7zwwtixYz0eT2lp6cSJE6NVGWOMHUeP8PXl5eWvvfbarFmzBg0a9MILLxx9PBAI3HPPPQcPHjz6yCWXXHLddddFuLrYMU1T0zREjHchnSKlREQhkmMIkoiklLoeaWPrNtwYYoeIhBBdbgyRtqH29nYi2rlz5zPPPPOjH/1ow4YNHY/v3bv3vffee+SRR47+z4EDB0opI1xd7DQ1NWVkZCTLh8rj8SCiy+WKdyGdYhhGe3t7VlZWvAvpLG4MsWMYRiSBFekuYUFBwebNm/Pz82trawcPHuz1ejse37Zt24IFC7Zv3x7JwrtTQ0OD2+3mNhoLhmG0tbVlZ2fHu5DO4sYQO4ZhSCntdnvXXh5pN/KCCy548cUXQ6HQ0qVLJ02aFOHSGGPsJCINrIcffviTTz7Jzc39+OOP//KXv0SlJsYY+0aRdnrz8vJWrFgRlVIYY+zkkuPIAmOMQeQ9LMYY6zyPAaCgi0Pu3MNijHWbzY2wdC/s93R9ghv3sBhjMecz4MkdsjoAfVMg3dr1qVTcw2KMxVZ9EP68R62qUUB0w2Do7+z6oriHxRiLoXU18u3DQATXDbPMKQI7GpGc8MKBxRiLlT0t6tlSIEU3Dten5wEiGEZEC+TAYoxFX11AvXYALEizC7UxGTA+SudlcWAxxqJsdyst2SR9Eub2w2uKQUTvshccWIyxaHphL31QrQTg+QV4wzAtimkFHFiMsWiRALesMT+ogWGp+PBEHJujRX0VHFiMsSjY1UKPbZdlPhyRDr+dJIZnxuTyhxxYjLFIfVwtX9kPpe142QBx5yjUYnatVg4sxljXhSQsWm+uqsYSNy0ZhxcWYnQHrY7DgcUY66IqH/1yi9pYT04LLhymXVAU86vgc2AxxrqiMUi/2yEPtMNZufjLCVpBSneslAOLMXZ6TAl/PSjrfFDowAlDxfeK0dJdJyVzYDHGTkNIwoJPjS+acFIOPDhB65/arTdD48BijHXWlibY56FAWOTY6J7R3Z1WwIHFGOukzY1w81ozxYK/n6INywB79KeFnhoHFmPsFAwFL5fS5w2qn4um5sLY+N0SlwOLMXYyrWF4bq/8tEb2TxGPn6F3/27gsTiwGGMnVB+gdyqotI1m5OEPBovslHimFXBgMca+ERG8eECtrKLxWWLBEDE1R4gEuKA6BxZj7Hh1Qbj3c3nQQ2lWnF0AYzMTIKsAgG9CwRg7ziEP/GSt+UkNZFjgmalibGyuu9A13MNijP3H0j3yhb2UosNl/fCXE7RMa7wL+joOLMbYEQ9tMV/YDzrSjSO064Ym4u4XBxZjDEISVlfTthYakoGLx+ozcxNoN/BYHFiM9XZ1ftrWTOVeunaIdk6BcCZwKiRwaYyx2NvbTs/uUnkp+N2Bor8LMUG7VkdwYDHWS7WF4b3D5NRpQBpMy8EBcZ3C3kkcWIz1Ri1huPvf5o5WvGYw3DJCS/CO1VEcWIz1OrtbaGMjpVrF7Dy6alDSpBVwYDHWq4QkvHZQLa+gwelw/RAcnpFMaQWRz3SfPn06fuWmm26KSk2MsViQBC+Uyr8fIqdQc/vjCLdIrrSCCHtYRLRnz57Kysr09HQA0HXurzGWoN47TI9+oYZmqLPztP8qtuR3yz0joi6iHlZdXV04HL7sssvy8/Ovvvrq9vb2aJXFGIuiZ3bTnevlQS/lOcQdo0WSphVE2MOqra2dPHnyE0880a9fv9tuu+1nP/vZa6+9dvTZysrKQYMGHf11/vz5t956aySri6mmpibDMJKlk+j1ehHR7/fHu5BOMQzD4/EopeJdSGf1pMagAJ7a63ivxmpFuGNY4Jq+4fr6bi7wawzDcLlcdru9ay9HIopKHTU1NSNHjmxubu74ddu2bd///veXLVt29D+43e7s7OyorCsWGhoa3G53srRRj8eDiC6XK96FdIphGG1tbYm89Y/TYxpDhQ+e2qH2tmOqFRaPpHHZ8T890DAMInI4HF17eUSbZMuWLcFgcNq0aQBgtVptNtuxz1osluHDh0ey/O6kfSXehXSKpmmImCzVKqWS6L2FntIYDrTTA9vUvlZxUT+4dZSWkhjxq5SSUnb55RElrs/nmzdv3u7du8Ph8P333z937txIlsYYiwoCeLecHt4uNYRbRml3j0mUtIpcRH/HWWeddd99982ZM6etre2iiy56+umno1UWY6xrAibct8nc0U4DneLHJWJMZrwLiqqIAgsRFy5cuHDhwmhVwxiLRLmH/rDL/LIVC1Px15O0bNupX5JcekpPkbFer9IL71er+gDOGwDXDdatPfHD3RP/JsZ6GUmwpk59Xg+Zdlw8TozMSLYJ7J3GgcVYctvrgV9+Ycu2m7OLLOfnY74z3gXFUvznZTDGumx/G9y91bK9mQwSVw7o4WkF3MNiLEkRwfsV8oVSyLKK7/eXvzrDmsiXNo6WXvAnMtbjmASf1agPq8hhgav7GtNzyZlg9+OKEQ4sxpJMa4ge3iYDkqbnabPyhcMMAPTYUfbjcGAxlkz+exetrZUODUe6xeUDhEWAxxPvmroRBxZjyUEpePWg+v12ZSp8aQbM7tsbj5hxYDGWBBpCcOc6c1szlKTDFYPF7L69ZR/wOBxYjCW6g+1w0zqz3AtD0uCJM/ShGfEuKH44sBhLaG+WyZf3kd+EWfn42ymaq3d/ZHv3X89YAgsruGWt+VE19nHB4lFi3oDeOGh1HA4sxhJRfRAe3Gx+WksOCz0wUf9WXi8dtDoOBxZjCWdLE/xqs9kcxtmFuHis1r+nn3DTeRxYjCWWT6rVy/upOQTfKlT3jrbYLPEuKJFwYDGWKAwFL+5Tbx1Q+U58aCLOKOCP5/H4HWEsIVR54Za1RpPEYam4aJQ2xh3vghISBxZj8ffKPvX7naZhinE59MRUvZfPXTgJfmMYi7NHvlC/204C4cYR+IuxPfPSxtHC7w1jcROS8N87zb+XU7ZdXD1EWzIuOe6EGEccWIzFx+cN6rEvVVsIx2biPWO1AanxLigZcGAxFgcfV8pbN1K7Sd/pR/dPtvSGi4VGBb9PjHUrqeBPe9Sr+5VV4Hf6waNTLBY+5abTOLAY6z6mgqd3qVcOqHQL/Xq8dl4RZ9Xp4cBirJtsb5AvHUQLymnZYuFwS0kWnx542jiwGOsOz+6Wv92hrIC3jNJuGCZ07lp1Cb9tjMXcmjr59C7lN3DuIPzxcE6rruMeFmMxRADlHjrYLs4txBl5dPkAnmkVEQ4sxmKl2kfvV6rCFDGlD3x3kLBzWEWM+6aMxcTHleqWdeb6erLoOCwdOa2igntYjEWZSXD/FrW7lcIkLiiEc/PjXVAPwoHFWDT5TViw0vy0Hka58YVz+GKhUcaBxVjUNITgTztlVQCGpsEjE4nTKuo4sBiLjo316smdyibg5hE4f7Cm87TQGIjCoPuOHTucTv4qYb0XEbxzSD3xpQqYNKcIrx7CaRUrkfaw2trarr32Wr/fH5VqGEs6dQH4yy6zPIhFTrh2uGVUL74tczeIqIdFRNdee+3dd999omeNYyilIlkXYwloYx19Z4X5Vjnl2ugXEzROq1iLqIf16KOPFhcXX3HFFd/47OHDh1NSUo7+ev311993332RrC6mmpqagsGgrifHoJ7X60XEZNkTNwzD4/EYhhHvQjqrk43hpXL77/elAOK3+gQX5Pr8TRCXHY2kawxpaWl2u71rL+/653PlypXLly9fsWLFif5D//79t2/f3uXldzNd191ud7IElsfjQUSXyxXvQjrFMIyUlJTs7Ox4F9JZp2wMRPDjtfTmAdOi4feK8ekz0xDTurPCYyVdY5BSdvnlXd8l/Pjjj1evXm21WhERABBxzZo1XV4aY8miIQAL15qrq1WqBR+YhH+YqiEPsXeXrgfWAw88QF8BACI666yzolcYY4mo1kdXfKTer6RJfdTmy/UfD+czbrpVcuwBMZYIPqqi50slIkzPFc9O11x8E/luF53A6uhkMdZTSYDnS2FllWwL421jxJy+KPi6AfHAPSzGTsEgeGan/J8DUOSCJ6doQ3juQvxwYDF2MiFJ7x2mSh+cky9+OAQ4reKLA4uxE/pbGX1So6bkwAV9cWY+8v244o4Di7FvQAQfVsknvlQ+Ka4ZrJ3ZJ94FMQDgwGLsfzMULNlKm5pwVJa4pB+nVQLhwGLsa6p88NDelAMBzLbRvWP0/qnxLogdgwOLsf945xA9tdP0BrVLB+EdYzQHTwtNMBxYjB3xSTXcs1kGJFzXN7x4lG7ltEo8HFiMgVTw6gH1QTUNT4PZheJ7OSGBKad+Get2HFist2sz4P4tZo0fi1y0cJLWz4kNDfGuiZ0ABxbr1T6opFf3q5oATMqCX4zXHHxt48TGgcV6r6W76L4tStfh/knimmKhcVglPA4s1huFJfxuh3z9AFks8L0BeO1gnsOeHDiwWK/TbsK8j8y9zTDMrX4/3To7N94FsU7jwGK9S0DCY1vkYQ/2ccH/nWrlk5mTCwcW60V2t8JH1dJmgVuGw49KeF5o8uHAYr3F49vpb2XmKDf+pESbmM0D7EmJA4v1fFLBujr5aQ1JwqsG6xOT5vY97HgcWKyH8xrw6Daj3ItXFuOsfK0oOW7fx74ZBxbryco99Lsv5UEPTu0D3x2g2bm9JznegKzH+qSK3i6XjQG6ciAuGMpz2HsCDizWAymChWvNNbU0JAPvHqtNyeF5oT0Eb0jW0/gNuOffxruHwVB4x2id06on4R4W61F2NsufbzT9pphZAD8frY3PindBLKo4sFjP8X6Fun09+KR2YZF6coolhe/M3ONwYLEeYnOD+tMeJQEu6Y9/ONOq8Y5gT8SBxZKeqeDl/WpFBfVzwI1DcU5/PuOmx+LAYsmt1g8/WG02BmhCH3HzSL3EHe+CWCxxYLEkVtoKcz82qwMwOgN+M0HL5+uw93QcWCxZ7WulezdLv4kT3HLZ+dZ0HmLvBTiwWFL6rFY+9AVoCPeOpuuH8xB7b8GBxZKMqeDB7aqsjQpT6NJ+2qX9Oat6EQ4slkzq/fCbL9S6OlngFG/M0p3cfnsZ3uAsaVT54Iefmntb4ex8+u0UjdOqF+JtzpLDBxX05C4ZkHB2AT491ZphjXdBLB4i3f9fvnx5SUlJRkZGSUnJhx9+GJWaGDvO73aoheuNslZ1eX948SyN06rXiqiHpZS66qqr3nzzzVmzZr399tvXXXddVVVVtCpjDACkhAWfmmvqQBfaHWPoRyN4n6BXi2jzm6b5yiuvzJ492+v12my2jIyv3TIpFApt3rz56K85OTn5+fmRrC6mDMMwDIOI4l1IpxiGgYiGYcS7kE4xvnK6L2w14KmdsLIadKSlZ9GMXNE9fzI3htiJsM6IAstqtV500UVerzctLQ0R16xZc+yzjY2NCxYsOPrrlVdeuWjRokhWF1Mej0cIoevJ8QXu9XoRUUoZ70I6xTAMj8djsZzezM7mMCza7Djk1SZnmAuHGWPtsq0tRgUejxtD7BiG4XR2/br6UdgkLpfL6/U+9dRTP/vZzz7//POjjxcWFm7fvj3y5XcPInK73cnSRm02GyK6XK54F9IphmFYLJbs7NO4Wc3r+9Wze6SfcFwf8fQ0e64jdtV9A24MsWMYRiTZGtGg+6FDh37+858DgNPpvP7663fv3h3J0hjr8LdDdNM6c1MTXFSIr84S3ZxWLJFFFFgFBQXPPffc6tWrieiNN94YP358tMpivdbKWnj4C+nQ6MIidd9ETed57OwYkY5hLVu27Pbbby8rKxs+fPjzzz8frbJYLyQVvF+pNjXCmExcMs4yj8+5Yf9LpHvp55xzzrGHAhnrGp8BP1un9nvVhX3F76Zobp5pxb5Jcgwrsp7tiwa4eYPZFKBR2eJHw5DTip0IBxaLsw318MPPzMYgnJ0Hr8zU7LwjyE6MWweLp9U1sGRTWEq4sEC9NtvCacVOjntYLG7eOKhe2q/sQrtvovZfxfGuhiUDDiwWBwRw/1bj3XLMcdCScdr0vHgXxJIEBxbrbkEJt2+Ua2sx3wkPT9THZWO8K2JJgwOLdatqHzyzW21tpOI0+ON0PY9nsbPTwYHFus/Hh+GH66VFqB8OwZ+O1F18nxt2mvioDOsm/zqsrlpjNvlVfyfeOZrTinUF97BYd3i10vrCAakDXD6QXjrHgjxsxbqEA4vF3JvlsHSfVQp8ZLK4egh36lnXcWCxGAoruHO9sb2NBrjUwtH6RX3jXRBLchxYLFY8Bpz9jjzkEyPS5WPjg2f0TYl3RSzpcWCxmAiacO6/zFIvZFpo6TTMFyreFbGegAcUWPTVB+DpPVKaMDYdNs21DHfHuyDWU3APi0XZlka6b7MZVNpdY8XlA1ETkCT3c2FJgHtYLJr+XkYLVpt72uHCvurKYtS4fbGo4h4Wi5p7Pw//cRfaNPhpiX77aJ5qxaKPA4tFx72bwkt3C0K4dhgsmcBpxWKCA4tFyjThtn+H3z0MTissGStuHMH7gSxWOLBYRAImzF9p/rseMx349BRtdiH3rVgMcWCxrgtLuPNzc1sTpFvFs9PFWbmcViy2OLBYF9WF4JJ/Go1hnOiGx6dicRqnFYs5DizWFfV+mL9KHvRifxe8NFtP5WvFsG7B46PstB1sh/9aadb5YEYerZrDacW6D/ew2OlZV0u/3ylbwnhOAT091aLzVx7rRhxY7DS8eVDevp4sOt1agotG6RoPW7Huxd+PrLNe2kc3r1WeME3PFYtKNE4r1v04sFinvF+hlmwyDYBrhuh/nalpWrwLYr0S7xKyU3toi3q2VDpteP1g+vVk7lmxuOHAYqfw+FbjsR1CR7xjDP20hI8IsnjiwGIns6uFtrRjroN+WAycVizuOLDYCS1aZ2xrgG8Xivsv0Aanxbsaxjiw2InM+Ed4fT1mWNQTZ1o4rViC4KOE7BvMXxHe2gg60E9LLGfkxrsaxr4SaWC98847o0aNysjImDFjRmlpaVRqYvF1+cfmv2rAruPb5+q/mMhfaSyBRNQcDx8+fPXVV//5z3+uqam59NJLr7vuumiVxeLl4n+F/1mOAuDlGdqF/TitWGKJaAzr4MGD8+fPnzp1KgD84Ac/eOSRR4591uv1vvzyy0d/HTp06Lhx4yJZXUyFQqFgMKjryTGoFwwGETHq1X5nNa5pAIcufzVezeoTDgajs1jDMDre3ugsLva4McSOYRhCdP2LMKI/cubMmTNnzgQAKeUvf/nL733ve8c+GwgE3nnnnaO/nnfeecOGDYtkdTEVCARsNluybPWONhrd+eZ/KrWuqbPoRH+eGpyVIwOBqC3ZMIxAIBCI4hJjjBtD7BiG4XA4uvzyKGySjz76aPHixeeff/4DDzxw7OM5OTl/+9vfIl9+9zBN0+12J0sb1XUdEV0uV7QW+JO18qCXJmfRTSWWywfaorXYDh1fqm530txPtZc3hpgyDENK2eWXR7RJiOjee+9du3bt66+/PnTo0EgWxeJo8Ubz1f2QYcMNc/V8e7yrYezEIgqsdevWLVu2bMOGDbque71eAEiWmGdH3b7B/J8DkG5Rj0yycFqxBBdRYK1atWrv3r3HdvWJKOKSWPe5bQO9WaYcunh4ov79wXxWM0t0ER23XrJkCX1dtMpi3eDRL4znS1UojM9Mwe8P5hkMLAkkx7Aii7qfrA69Wi6cGv18rH7xAO5bseTA36u90Uv74a/lIiThthJx2yhOK5Y0OLB6nad3qVvWmkh04wjx8/HcAFgy4fbau6yppl98rkKSLi6C/z4zOaYaMnYUj2H1ItU+uuYz0wCYmQ+vzbbGuxzGThv3sHqLhhDM+8j0h2lGLvzzPL52KEtKHFi9QlDCZSuMA+0wOku8e74lSc45Yex4HFg9n2nC1SuNsnYsTIHXZukW3uYsaXHj7fmu+CT8r8PoEPD38y05XT9PnrH4432DHm7xRvVhjdAF/WqCPjA13tWw3qPF9B+u9myt8lc2yqDSUIXbAypoEEDKuKIBP5ndtaVyYPVkyyvkc6WEQLeM1BYM5QmiLEoMae5raSmrCpU3hcpbIRhSYSVJEgFIJRARkYgAhEIlAAmkiYJIgRIgQB5s7fKaObB6rN0t8NONSgf8rwH6g3y7ZtZ5AWlWNNdtPWQerA3X+1AIMlGBhKBJigAUACEIhQQEAEoAEgmFJAAUEgESgoYARESgNLJYHZCRYkmzK1Pa8zNdF3X9Qp4cWD0TATy23fCExfg+tHQGpxX7D9Mf9u+uDlS1hMuaQ1WtMmxQ2EBSIIUCJQgVkRAAiKSIiASCAkQUQAoBCBEBgQgFgiZAamDVdLtVz3boWRnpYwusA3Jt2SccK43nBfxYwrpzo3nIg3P6wYOTdOS86lXaZaCyqWXf4fDeBk23mM1tZovfNJQAIEAiJVADIkICAkQgAEQkCQIVAkhQqAkAkIBCAJDQMxzWvAzbAHdKvz6OgjQt1wXxuxwzB1YP9GKpfO0AuG3i5fEih6/J1/M0BJt3VHj31xnVrSIYMjw+MDUJpqZpZBABIYECOjKSRASIgIJIEZIAJFQACAI0q0UQosshXFbHALcoysge1RczE/rQDAdWT7OjmX61TRGJW0dSkTPe1bAuMIxweWvrnjoZMHxbDkqPnyQiESkJAESIKBAkEBIAACICkSkAFCiAjnElBE1DIEumS9jtaBeW7JS0oQWO/jlaUXp8/7gIcWD1KAbBD9coTxAu7gc/Gs7nNicqKUM1nvYvK/w7aoyGNkANQqYMBkgRECAgACgEgYBEAAhKodAkoQZACACCiFBHFBo6dEuaTVmsGWcOShtbZEmzQo8+j6En/2290LUrw6Vtor8T/nw2p1WcUUugfUdV46oD0tuioy5DBgVMIkWEQKQQBSkgBARElIQCQBAqJEISIASCZtNBaIpIT7PYB+enjsxNLS5C9/FnrXs8niS6a06EOLB6jr1t8O96YRfq6WnCxnnVDcJhs8Hf9O+ycFlTqMmrAZj+oBmUIA0BKAkEAAkCQqIwEQECIWHHv4IABQFqAFqWK7VfujUzzTUg1z4sH1L5U3lC/Nb0EFLB+5XG+Fy8uMA6Iz/e1fQkAeXbVRGsavbtqwlX+2UoSFIIUCQVKSCgjoOwCgAFQcesSUJA0AQQgABE1LTMNFuOC7McaaPzU0f0gxQ+Ja6LOLB6iCtXmjub8OZRcM2QeJeSjMLkP1jnK63zbaoMNTYTag2moRMoQDhy3B8kgaCOISVTwpG5Iig0EERCWBzWlIGZtoIMS1GWLSfVPiArzn9RD8WB1RNsrKcVVQoJZ+Vqgr+8T8yobG7ccCC4q56CQdkWVqYBBAoAiQhBqI4hbSCUCCiJBByZGyAsutA0QLJmp1kH5aSNLUwZmI9Ofq+7GwdWT3DP59KKeNlAMTqLP0IA3rB3b51n52GztCHgDZChlCmBEIkUgkZKAiEJTQBRxz6dUEQaCrAigJAK9MI0V9/0tNFFqSX9efctoXBgJb13y1RdCEZliV53jXa/at5xMLCzLrC7WvlDylQCSZkEKIhICCCFBAoQBQIRIaDoOBEOBDjtWnaqwyFyzhupFeVYMr82v7ahocHtdus9en5AkuJNktyCJizdK22I/2es5uqp1z0OSF9phW9LRSgowzVtRnM7mIoIkZRUIBAEAQkUhASIIKDjSQS0CqvDactx2icNcA3IthVnx/svYZHiwEpu/7PPOOTHsZl0XmG8S4mcUv59ja0rtvuq28EfNgJhDQR0nGYCgECSSHR0llTHTh1adELUSQA4bCl9M91nDkkZWYBObtU9Fm/aJBZW8OJ+lApuH5lk25E8wZatFa0f75atXjMsBSKZBEBEQKA0EAQgUElQghCE0EAjMJw5GVpRpnBpzpF90wflQ2pP7VKyE0qyhs6O9Vkd+U0an4kTcxL3ggzBQ61t6/f46traDrcYAUmgNAkKFAEiKFAokAiRiBBRFyAJSKAlPdWSn5F6RlH6hH7CzsHEjuDASmLvlctsh1g8JlHSSpY116/fa7SFoa4t2OSloAECgIAIEFERISpUQglAEKiRJKHZhC0rI3XqgIwpA7UMPlebnQIHVrIq90KrgZOzYWw8pjKoVqP584OhPdWhqkYzQBg2lDIVIBIAAQpQSgEBKU0IBYDgtNidTsqwZJ8zIm10P7DzXAHWFRxYyerVMm1/m7qkKPbX5wupti1V/vLa0L5GpQzwmka7D0CpjoubHJkHjgpQ04VSoKXanUWZlsF9rNn29MnFgGgYRltbW3Y2H6RjkeLASkpEsL4ZgxKm94lyXMlA2L+n2ru7NlTerMKSQmGj2Q9KCRIKFSKCJhSAAE1YQHNaLVmpKQOz0qcNseal9uwLm7BEwC0sKW1qhgqvGJAK+SkRLSdY2RyqaPXvraTGUKDeo3wBAgAiUAoIQBdosWgWRCmE3YKZ6bZBOVlnF9tyM8CaKANnrFfhwEpKzSEtywrn9T29qe3hhnbPjtrgzupAdRNqQrYGSUlQCgBRIAEoSbrdgmk2a5rd2sftKMlLHVUg7MdfgImxeIlaYEkpR44cuWfPnmgtkJ3E0DT64RA5p/ikURI2feVN7Z8dECQDVc1Ga0CFFCkJCgBIEwI1RF0XaQ57Xqpe5E4b1dfRNxOsvez8HpZUohNYTz311Kuvvrp3796oLI2dUljRiDzTOD8AABSlSURBVDSVYzvmIZOCdW3tW8sxGA5UtoYrmlFoqIPRFtQsGggABejUrSmplgyH3ictY2qxfUAW8I4dSyrRCawxY8YUFxfPmTMnKktjp6QjWJUya5rJbmv9dJ9/b7UMmioQVv6QsOqkIZiKXJbUIX2k17CP7usqzrYVuIHnErAkF53AmjVr1v9+sKWlZfHixUd/nTJlyvnnnx+V1cWC1+vVdT2RT9C/bxu8cUhvCYswIYE+2d/yWnCbrb87dKjJbPGJFKtWkKahsg3Msw3PsWa6wPqffAoDhH2eeFVuGIbX67XZbKf+r4kh8RvDsbxe71e3hk8ChmFYrV0fFY3hJkFEu/0/V+2wWCyYwLf0xK/Eu5D/eOgL+Va55ZBfoyNnACNoOlDHqcBYZ3Pa+2SkjC1ynDkQNaFnnPBeu3GXgO/tySVXwUlXbSQvj2FgZWRk/OY3v4nd8qMrEAi4XK74fqk+uAWe3hFqMgWRIiXgyLRxDYgANQRAogwrnVsoinVpc9n6TJiqJUMr7bg7eRLd1iURGkPndZyGmSxvL9+qPom9fwju3Bwqb4UQaYQApKBjnOnIHTIBhZYpjCl91ZsznMeeAvzn7cGd7cJjQAZPOWC9CQdWtwqE4I4vjBd2K1MKUkiCQAoQHecHAyIiqDQB5/UTfzrD4j5yz/Bv2EYtBnzZDCurad6AZOhiMRYl0QysZBn262Z7vHDjKmNLE5iSDCUAEEiAJgDMjj16h4XOy9P+cLZW0OlhqPGZtLYBqgMEPDGB9Sbcw4qJp3bCw18EW8JCGkDYcXs6BaQBIpIUGvaxhZ6YnPLdoV1c/uRMsOvw1kG6tB8U8UVZWK/BgRU1f9il7t4QDildAQApEDooCUKAAquuXDr9YhwuGqkJ0TGVPKKL0mkC8uzyswb9ro3mq7N5I7Legtt6RO5eI/98MNwuLUoREB0ZMgdCAQjU1y4eO8Ny5bCYrHrRUPXXcnj7MD23y7y+hLcj6xW4oZ+2Zfvh7i2hgx6hFAISgN5xo3JUQtNUQYrxh2mOi/vHvIz8FDwrx/igyvLADrpsIGQn7jQsxqKGA6tTytphzkf+va0WUkQgAAV8NTCFCANT1cuz7WfmdHdVf5vtOO/D8OYGnPqusekyS7r91C9hLKlxYJ3MreuM5/aqgBSEAHRkyhMiIMGQdPXQKNvcER2PxecuCboO/+9c66i3jEMeGvG38ObLrfncz2I9GgfW8eo9sHAdvVdhSkQiBEAgAgBESNHkXeNsS8Yl0EyCLCv869uWmf8w64Mw7A3jH9/WZ+QlUHmMRRcH1hGb6+HHq1StJ2QqASQBdZASNQ2RpuZr73xLd1shXj2pkxuXBZ/O0We8b3hDMOcD866x4t5xfE0r1jP19sBaXgXXrgw1hpFkGqAAIEBChblO8+cj9VvH6knxFo3JhEPfscz5SG5qlL/ept4uUy+dYxmZGe+yGIu2JPg0xsLbZXDj6lCrqREpQAASIJQAGpZF78yyFqfHu77Tl+GAz+Zov9hEz+yiL1to2j/k9wfBE5O1lKS5pgtjp9a7AuvlHfB/tvmrQxYCAOoYnEKLhlMy6eVxbYWFbl1PxJ2+zrt/kn7baJj3QWhrC7xYql7aD+fnGy/PtKdzbLEeoVcE1v5WmLU82ODTDCCQOiAAoNBpeo5adcmRj3JDQ1xLjJ5MG6y+1Lauln6wmsp88l+VWtHr5swCeHCiPpp3ElmS6+GBtXij+u/dhmEKAA0QgDShycsGaP/vWz38asHT8nDf9/S/l8lfbpUH22lFBf2rgtKttKBYPHiG7uBBeZacemZg/e5LWLIhbHacIwMIoATi0HRj/aXOtN50V5i5A7W5A7U9bXTnelpRo1pC8MwuWLo3nOmgqwfir8dbLcm9B8x6nZ4WWDP/SWvrQkrpIASAQqHGuOGJM2wzCwCgl17sbng6/uPb6DPE/ZvMZZV02AdVXvH4l/BsqdnHTvP6iR+P1PrzJR9YMughgbW+Huas8LcGrUQEJACVILh4EP59NnchjnBa4JGp+iMA+zxwxwa5q4VqfLQ/DL/dAY9/aRSm4DlF8P1++K1CzdKL+qAsySR9YP1iPTyyx1AmAViAFGrQz0Wr5tj6cZfhBIakwrvnaQBQ5YEXD6jnS2WNDxpD6u8H4B8HyWUhm5WK0/CyfuKqQSKVA58lkiQOrGn/MP9dSyQICAFRR/rBEG3pDO4edFZhKiwZJ5aMEwCwsRZeKQt/2aIdbJWVHihvF6uq1d0bZL6TRmTow9LoggKcXiAsPfxYBUt0SRlYI94K7mvTiAhQgCKnlf440Xr1yHiXlcym5MGUPCsAEOlvl9G75caXLVgbhOaQWFUlP6igP+4WJkibwLNzabSbJmaJsVn6oCScYcuSWjIFlmHA0NcDFSGdQAAQIqbZ5eYrbAP5sirRgwjfGYTfGXTkAEWNH96rho/LzD1eOtgGfoPWN+Bn9RgyZchUFg2HuqGvE2b3wUlZWlE6DnRBktwfjyWl5AgsIjhjWWhrkyDQQRAqMdAd3ndFSrzr6vnyU+DGwXDj4CPtZK8XqlvVS2XGvnbL5gaFBIdaZXmbtqMJ/GY4LDE7BYsciDqlobp+mLWfE4rT0c4RxqIkCQJr5rvmmjoi1AAlgijJge2XJsc5yT3PMBcMc4lZfTtODxCmgu3N2qYmaPLLT+vFPg9ohBVeajQEKVHuN+2aGOmmjXXgNdw5DmNAGpyZLXKdon8K9HVg34yEvPwFS2AJ/bG/8hO1rIwICBAQ6Jw89fElfFJcAtEFTMjGCdkAoN/z1YOtBnxUqfa3yUybts9DhSm4vgbaTeVpF5U+tbmBTKXCREGD7BqmWihdh0m5wlR0di5KQJugcTna0FTg2RXsf0vQwPrjTvmzdUqhBNSQsMQtt3/HlpiXo2LHybDAFQPFV/fjAAC4eZhR3dBWp2XWBbXaoNrbDqVNcMAHQYleSf4wrqpT3rD4vFHWBUV7WKVblCYQAAocKmCS24bnFWAfl+YLy/5OvX8a9ndBphU0PmTZ+yRcYO2sg2nv+72mBohAItsmSy+3pzsTrk52WpxWOCO7I2C+FjMStKYA1Hppe4t0WS3vlps7WkS6jeoC1G5iS1g0BKnGjz4TwlI2hgHISLVg/1SRotO+drAgaQjZdjQlZdhhdoFeH1Zn9dGyLCqoRJYOfV3g4k55D5JYQXDe8tAnFQiogxB2YS6/wH52QbxrYrGkAfRxQB8HjsnRAWDugK91ov0GVHuh2qvQKmq8tK7B8BsYUiLXRm0m7m8HnwSfgQ1+apcokHY2m60GvrbPRIRajwwTpFgwzyGtmtARitNwVr7mM2R9WDX6YWCaGOgUeQ4YnCaa/FDkhlQNbBof5UxoiRJYS3eqRRtNKQEQUdA9I+H+KTxbobdLscBgNwx2CwCAXPxu8fFng0qC5iC0h9S2JvKZAjTaWE+FDlHhkZ8hNAUoxQoWTQQN9EuQQGlWta8NKrzYGkIhwKnLbDtqaJZ50GXFXAc0hwgBS1JTdZ02NxjFGXhmDlZ4qSUMmXY6K0frlwpSYcBUEzKFSZCWgqkCNASdR9y6RfwDK2hA39fDLWECpYEGZ/eBVZfwWBXrFA0hxwE5DlGcceSRBYM7/v3ajqdSEJAQJjAUNPlhf7va3AyCZEgKuwXCJFZVqlQrOiwYNCigyGtQgx/qw1qwRbYGoTKovGGBAOvrVB87VfpREuTazZBCAMizk89EDaGPHUu9lG2joanCJ6ElBLkplCJEhV+Nz8Rv98V/VMLwdByUigQy1SJSNKVAy7SDVYBDByTu3J1anAPrV5+rB7YbhAIInbraMs86OOPUr2LstAgBTgEdZ5f2scOITDFnAHwt1CYc/UkDgPr61pRU96Gg0EFoGuxsoj0e02uKoS5NKlXqpSqPGpiqHfYqH2G2FSt9EFYQNA1PSDclElFTiJqCaG0mhbI5hDtb4N9NtKWBbAIKnRhUJJRMsWG5x7BpmGWHIjvUhqjMCxqCKbHEjWMysY9dIeK+VgIEgVScpvV1QmuIXBZBRA0huKAQHBrWtFFY4TAN2vzktCEg5DpAAVgF9LwAjFtgeULQ57Vw2BSgBAo1r1i9NYtvqscSAiLYLTDqq/Y4JA2POULdEXPaMT8fpQWDICxgKqgPQ5NfAWJziLY14sBUs8hlee2AyrbKnBSs8BJqaBoqRRd2klaLKHIJb6vSEE2l/Apr/MqpY2sQmww67EVPiEzEglaVbqGmICKaLSHwmbirDR2Ctjfaw4rcNqPcJ7wGOa04MVPWhTDTAk1htCJl20W7QRKozg9BEwtSKNUKQYk60sVFmlXAYS9N7oNEuKNJ2qzCJnBoGpgEmXYobaPdbTguQw5362kWagtC31TUEEIEdgRJoCPoOoAC0S0HbeMTWA99bvzySyRAQMqwq8PX2Jw977uA9T52OwCAVYMBFhjg7PgE43mF0JF3U/qIr+cdfP3nI0+FQqDpmgIIERBBhYdCEvwSdZQWoVUHSCLU+eCQBy8spDBhkT1QH9CL0vRNjWabAQgwPF1oHsy0UogQFBlK+cPoUVgboKBBPhNsuvCGlQCwWZROUOXH2oDymvBFowoDWAUUOcGho0Gwpla2hkW6lc7INm06tIVgvJucVrGnHTSA2qBSSihUh9qFRciBqZpB5LaSUmJSDtUHoDGM+1vJoclcpz7WDQGJeQ6YkQ2juzpAHYfAynw52BYUoAECnV8g3r+QDzsz9h+2rz4QHYcYRmYe/TLXAGDi1/bzEACmpBGi6XLZTzBR8UggSgUSwRMGlBACaA+KMGKug4ISK7yqnwN9Ene0iZCBUlJfF7YbaNEp32F+3gTDUmFyDlqIyvwwOUcHIkXKqqO1XTOkCims10jTNEFkSuULg0FUH6Aqv6jyQaUPgERDCLyGAiVSNBjiwtFdfXO6NbA218GZ7xkKNNBQR3X4CmtuWneun7HeSxOgAWR9lYb5KR3BhwDQ70hnEIa7//Ngxw8XF31jfwIvH3h0D/AbdgXDYbBYISwBEJqDABJCBBl2rTkAVjLtgrr8V3RfYF3+YeidwwgKQODY9PAWPnWZsR7KagUAsGkAAPnHfNAzrGAYQkrZ5SXHaqDMMAzDMDp+bg1B+ouhdyoEINgssOtySwKmVU1NjWma8a6is9ra2trb2+NdRWeFw+Ha2tp4V3EauDHEjmEYNTU1XX55pIHV0tIyZ86czMzMSy+9tKWl5ejjBw8ePHToEAC8fQCy/2p6pQYK+jrJf511WELeHW/u3LkdBSeFp5566tlnn413FZ21e/fu+fPnx7uK0zBv3ryysrJ4V9FZSdcYLrvssi6/PNLAevTRR/v3719TU9OvX7/HHnvsuGfH/j145UqDCFDBlYOwfD6PrzPGui7SwFq2bNmiRYtsNtuiRYvefvvtY5+SD+3fUY8AKFDtnq+//q34z6pnjCU1JOr6iD0AuFyuhoYGh8MRCARyc3OP7ksvX758zqGpSk+xNJW5nzwHAPLy8oqKiqJQcmxs2LBhzJgxKSkJN7j2jQ4cOKBp2oABA+JdSKd4PJ49e/ZMnjw53oV01saNG0ePHs2NIRa8Xq9SatWqVV17eaSB5XQ6m5qa7Ha73+/Pycnx+XxHn/rLX/5SWlp69Nf8/Pz+/ftHsi7GWA8wderU/Pz8rr020t20goKCioqKIUOGVFVVFRYWHvvUDTfcEOHCGWPsWJGOYc2ZM+f5558noueffz6SwX/GGDulSHcJW1tbr7rqqi+++GLChAmvvPJKejrfqY4xFiuRBhZjjHWbmMx0P9Fs0gQnpRw+fHi8qzi1d955Z9SoURkZGTNmzDj2sEZiWr58eUlJSUZGRklJyYcffhjvcjplx44dTqcz3lV0yvTp0/ErN910U7zLOQXTNBcuXJiTkzN9+vSqqqquLIJi4K677rr55puDweDNN9989913x2IVUffkk0+eccYZMXpDoqi8vNzlcq1bt87v9z/++OPTpk2Ld0UnI6XMzMz86KOPpJRvvfVWQUFBvCs6tdbW1okTJyZ+SyAipVRmZmZlZaXH4/F4PIFAIN4VncLjjz9+1VVX+Xy+O++88/rrr+/CEmKyVYYOHbp7924i2r1799ChQ2Oxiqj75JNP3nvvvcRvpitXrrzhhhs6fq6vr8/KyopvPScXCoX++c9/KqXa29vffffdkpKSeFd0CkqpuXPnvvXWW4nfEoiopqbG5XJNnDjR5XJddtlldXV18a7oFMaPH79t2zYiam9v37RpUxeWEJOt4nQ6/X4/Efn9/tTU1FisIkaSopl2ME3zpptuWrhwYbwLOTWPxwMAiLh27dp413IKDz/88B133EFJ0hK2bt06a9asrVu3NjU1LViwYP78+fGu6BQyMzPvuusut9s9ceLE7du3d2EJMdkqKSkpHb1Tn8+XkpISi1XESFI0UyJasWLF+PHj77rrLsMw4l1Lp3i93gcffHDSpEnxLuRkPvnkk3POOSccDlPytISjqqur3W53vKs4BV3XFy9eXF1dvWTJkilTpnRhCTHZKoMHDy4tLSWi0tLSIUOGxGIVMZL4zVQpdffdd5999tl79+6Ndy2nVlZWduedd3b8XFtb63Q641vPyS1ZsuS4Ed7PPvss3kWdzObNm492WhsbG/Py8uJbzynl5+dXV1cTUU1NTdcaQ0yOEvJs0thZt27dsmXL3n333YKCAq/X6/V6413RyRQUFDz33HOrV68mojfeeGP8+PHxruhkHnjggaMfDAAgorPOOiveRZ2Mz+ebN2/e7t27w+Hw/fffP3fu3HhXdAoXXHDBiy++GAqFli5dOmnSpK4sIjrJ+XUtLS0XXXRRYWHhnDlzWltbY7GKGInRGxJFDzzwQDdswShatWrVhAkT3G731KlTOw7FJIXEf2OJSCn1zDPPFBcXZ2dnL1iwoK2tLd4VnUJNTc25556bnp4+Y8aMffv2dWEJPHGUMZY0uuVeYowxFg0cWIyxpMGBxRhLGhxYjLGkwYHFGEsaHFiMsaTBgcUYSxocWIyxpMGBxRhLGhxYjLGkwYHFGEsa/x/hAGGuVX0LRgAAAABJRU5ErkJggg=="
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"stability = zeros(Int, n)\n",
"dummy = zeros(2)\n",
"\n",
"@time for i in 1:n\n",
" \n",
" a[] = as[i]\n",
" J = ForwardDiff.jacobian(f!, dummy, result[i, :])\n",
" \n",
" if stable(J)\n",
" stability[i] = 1\n",
" end\n",
"end\n",
"\n",
"scatter(as, result[:, 2], color=stability,\n",
" markeralpha=0.2, markersize=2, markerstrokewidth=0,\n",
" legend=:none, xlims=(0, 6), ylims=(0, 6))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.1.0",
"language": "julia",
"name": "julia-1.1"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.1.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment