diff --git a/rabbits-and-foxes-kmc-question.ipynb b/rabbits-and-foxes-kmc-question.ipynb index 657efb4..48b5de7 100644 --- a/rabbits-and-foxes-kmc-question.ipynb +++ b/rabbits-and-foxes-kmc-question.ipynb @@ -23,7 +23,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -79,13 +79,141 @@ } }, "outputs": [], + "source": [ + "def monte_carlo():\n", + " \"\"\"Perform a KMC simulation of rabbits and foxes population dynamics.\"\"\"\n", + " times = []\n", + " rabbits = []\n", + " foxes = []\n", + "\n", + " time = 0\n", + " rabbit = 400\n", + " fox = 200\n", + "\n", + " times = [time] # list of times\n", + " rabbits = [rabbit] # list of rabbits\n", + " foxes = [fox] # list of foxes\n", + "\n", + " end_time = 600\n", + "\n", + " k1= 0.015\n", + " k2 = 0.00004\n", + " k3 = 0.0004\n", + " k4 = 0.04\n", + "\n", + " while time < end_time:\n", + " rate_rabbit_birth = k1*rabbit\n", + " rate_rabbit_death = k2*rabbit*fox\n", + " rate_fox_birth = k3*rabbit*fox\n", + " rate_fox_death = k4*fox\n", + "\n", + " rates = np.array([rate_rabbit_birth, rate_rabbit_death, rate_fox_birth, rate_fox_death])\n", + " cumulative_rates = rates.cumsum()\n", + " total_rate = cumulative_rates[-1]\n", + " ut = np.random.uniform(0,1)\n", + " \n", + " if total_rate == 0:\n", + " break\n", + " \n", + " event = -1\n", + " for n in range (10):\n", + " u = np.random.uniform(0,1)\n", + " for i, cum_rate in enumerate(cumulative_rates):\n", + " if cum_rate / total_rate > u:\n", + " event = i\n", + " break\n", + " #carry out event\n", + " if event == 0:\n", + " rabbit += 1\n", + " elif event == 1:\n", + " rabbit -= 1\n", + " elif event == 2:\n", + " fox += 1\n", + " elif event == 3:\n", + " fox -= 1\n", + " else: \n", + " raise ValueError(\"event not in 0,1,2,3\")\n", + " time = time + ut\n", + " times.append(time)\n", + " rabbits.append(rabbit)\n", + " foxes.append(fox)\n", + " max_rabbit = max(rabbits)\n", + " rabbit_max_time_index = rabbits.index(max_rabbit)\n", + " max_rabbit_time = times[rabbit_max_time_index]\n", + "\n", + " max_fox = max(foxes)\n", + " fox_max_time_index = foxes.index(max_fox)\n", + " max_fox_time = times[fox_max_time_index]\n", + "\n", + " #plt.plot(times, rabbits , label=\"Rabbits\")\n", + " #plt.plot(times, foxes , label=\"Foxes\")\n", + " #plt.xlabel(\"Time (days)\")\n", + " #plt.ylabel(\"Population\")\n", + " #plt.xlim(0,600)\n", + " #plt.legend()\n", + " #plt.show()\n", + "\n", + " return [max_rabbit, max_rabbit_time, max_fox, max_fox_time]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[401, 1.2889238279960686, 891, 600.6679963468558]\n" + ] + } + ], + "source": [ + "values = monte_carlo()\n", + "print(values)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "#code to determine how many times we should run the code for enough precision\n", + "\n", + "num_runs = 100\n", + "max_fox_list = []\n", + "max_rabbit_list = []\n", + "time_max_fox_list = []\n", + "time_max_rabbit_list = []\n", + "\n", + "for run in range(num_runs):\n", + " values = monte_carlo()\n", + " max_rabbit = values[0]\n", + " max_rabbit_time = values[1]\n", + " max_fox = values[2]\n", + " max_fox_time = [3]\n", + " max_rabbit_list.append(max_rabbit)\n", + " time_max_rabbit_list.append(max_rabbit_time)\n", + " max_fox_list.append(max_fox)\n", + " time_max_fox_list.append(max_fox_time)\n", + "\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "5137", + "display_name": "comp_mod_class", "language": "python", "name": "python3" }, @@ -99,7 +227,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.13.7" } }, "nbformat": 4,