{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Resolvendo numericamente equações diferenciais com R\n", "\n", "*Esta é uma breve descrição do que é integração numérica e um tutorial prático de como fazê-la em R.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Software necessário\n", "### R\n", "*Para executar os códigos R deste notebook em seu próprio computador, você precisa instalar o seguinte software:*\n", "\n", "* [R](http://www.r-project.org/), juntamente com os pacotes do [CRAN](http://cran.r-project.org/):\n", " * [deSolve](http://www.vps.fmvz.usp.br/CRAN/web/packages/deSolve/index.html), uma biblioteca para resolver equações diferenciais\n", " * [ggplot2](http://www.vps.fmvz.usp.br/CRAN/web/packages/ggplot2/index.html), uma biblioteca para plotagem\n", " * [reshape2](http://cran.r-project.org/web/packages/reshape2/index.html), para manipular data.frames\n", "\n", "Para instalar o R, baixe-o de sua página inicial (Windows ou Mac): http://www.r-project.org/. No Linux, você pode instalá-lo usando a maneira preferida de sua distribuição, por exemplo:\n", "\n", "* Debian/Ubuntu: `sudo apt-get install r-base`\n", "* Fedora: `sudo yum install R`\n", "* Arch: `sudo pacman -S r`\n", "\n", "Para instalar os pacotes, tudo o que você precisa fazer é executar o seguinte no prompt `R`\n", "\n", " install.packages(c(\"deSolve\", \"ggplot2\", \"reshape2\"))\n", " \n", "O código R apresentado aqui e alguns exemplos adicionais estão disponíveis em https://github.com/piklprado/ode_examples (obrigado, [Diogro](https://github.com/diogro), e todos mais que contribuiram para estes tutoriais).\n", "\n", "### Executando comandos este notebook\n", "Você pode executar os comandos neste notebook diretamente no R de duas maneiras:\n", "* Copie os trechos de códigos R neste notebook e cole-os no prompt R. Alternativamente, você pode salvá-los em um arquivo de texto puro com extensão .r ou .R e executá-los diretamente em um shell R, ou do R Studio.\n", "* Baixe o [script R](https://raw.githubusercontent.com/piklprado/ode_examples/master/Integracao%20numerica%20em%20R.R) e execute os comandos em um shell R.\n", "\n", "### Executando comandos R do Jupyter\n", "Este tutotial está em um Jupyter notebook, uma interface que combina texto e comandos em várias liguagens, e executa os comandos, se você tiver as linguagens instaladas em seu computador. *Para executar os comandos R neste notebook a partir de IP[y], você também precisa:*\n", "\n", "* A interface [Jupyter notebook](http://jupyter.readthedocs.org/en/latest/install.html) e\n", "* o [R Kernel para Jupyter](http://irkernel.github.io/).\n", "\n", "Para instalar o notebook no Windows e Mac, recomendamos instalar a [distribuição Anaconda](https://store.continuum.io/cshop/anaconda/), disponível em http://continuum.io/downloads. No Linux, ele deve estar disponível nos repositórios da sua distribuição.\n", "\n", "Para instalar o Kernel do R no Jupyter notebook siga as instruções que estão em http://irkernel.github.io/\n", "\n", "### Executando da web\n", "Se por algum motivo você não quiser instalar nada em seu computador, você pode usar um serviço que executa notebooks na nuvem, por exemplo, [SageMathCloud](https://cloud.sagemath.com/) ou [wakari](https://www.wakari.io/). É possível visualizar notebooks disponíveis publicamente em http://nbviewer.ipython.org, mas nenhum cálculo pode ser realizado (apenas mostra os resultados pré-calculados salvos)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Como funciona a integração numérica\n", "\n", "Digamos que temos uma equação diferencial que não sabemos (ou não queremos) derivar sua solução (analítica). Ainda podemos descobrir quais são as soluções por meio da **integração numérica**. Então, como isso funciona?\n", "\n", "A ideia é aproximar a solução em pequenos intervalos de tempo sucessivos, extrapolando o valor da derivada em cada intervalo. Por exemplo, vamos tomar a equação diferencial\n", "\n", "$$ \\frac{dx}{dt} = f(x) = x (1 - x) $$\n", "\n", "com um valor inicial $x_0 = 0.1$ em um momento inicial $t=0$ (ou seja, $x(0) = 0.1$). Em $t=0$, a derivada $\\frac{dx}{dt}$ é igual a $f(0.1) = 0.1 \\times (1-0.1) = 0.09$. Escolhemos um pequeno passo de tempo, digamos, $\\Delta t = 0.5$, e assumimos que esse valor da derivada é uma boa aproximação ao longo de todo este pequeno intervalo de $t=0$ até $t=0.5$. Isso significa que neste tempo $x$ vai aumentar $\\frac{dx}{dt} \\times \\Delta t = 0.09 \\times 0.5 = 0.045$. Portanto, nossa solução aproximada para $x$ em $t=0.5$ é $x(0) + 0.045 = 0.145$. Podemos então usar este valor de $x(0.5)$ para calcular o próximo ponto no tempo, $t=1$. Em resumo, calculamos a derivada em cada passo, multiplicando o resultado da derivada pelo passo de tempo, e aí somamos o resultado ao valor anterior da solução, conforme tabela abaixo:\n", "\n", "\n", "| $t$ | $x$ | $\\frac{dx}{dt}$ |\n", "| ---:|---------:|----------:|\n", "| 0 | 0.1 | 0.09 |\n", "| 0.5 | 0.145 | 0.123975 |\n", "| 1.0 | 0.206987 | 0.164144 |\n", "| 1.5 | 0.289059 | 0.205504 |\n", "| 2.0 | 0.391811 | 0.238295 |\n", "\n", "Claro, isso é terrivelmente tedioso de fazer à mão, então podemos escrever um programa simples para fazer isso e fazer um gráfico da solução. Abaixo, comparamos com a solução que obtivemos na tabela acima com intervalos de tempo $0.1$ com a solução analítica conhecida desta equação diferencial (a *equação logística*). **Não se preocupe com o código ainda**: existem maneiras melhores e mais simples de fazer isso! Os pontos são os valores aproximados para alguns momentos, e a linha mostra a solução analítica." ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAAJYCAMAAABvmDbGAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deWATZf7H8acH5UZO14tFFBVl\nV1B01WUVcT3WXVLuUwRZTsVFEBQQlB8gorIIKEpXLYciAgooAiIgcssNAnIfpbRA+3ULlKul\nx/zmySRtkiZzPTPbJPN5/5HMTJ9+me72ZZM0TZiEEBKOlfYJIBQNARJCFgRICFkQICFkQYCE\nkAUBEkIWBEgIWRAgIWRB1kDK0upCwSXNNSa6cMWOqecK7Bmba8fUrIKrtozNs2Xq1Xx7xp61\nY2pOwTmtJecshkRaZUsXNdeY6PwVO6ZmSfaMvWrHVJLybBlbYMvU/EJbxub9146pOVKW1pIs\nQFIJkAiQCJCEAyQCJAIk4QCJAIkASThAIkAiQBIOkAiQCJCEAyQCJAIk4QCJAIkASThAIkAi\nQBIOkAiQCJCEAyQCJAIk4QCJAIkASThAIkAiQBIOkAiQCJCEAyQCJAIk4QCJAIkASThAIkAi\nQBIOkAiQyA5I2V94ty78u/MLK3yuAcnAWEByOqS5L3q3hr+8e0mLLcXXgGRgLCBZAyl1x5mA\nI4GQTiz/MU1jyKYvfsrQ+Gd+WLhD61QMQVr8mssLKcWVJklTRhRdA5KRsYD089Qpq9RXLGp2\n3d1vpKut2P2PGFZpuD+CAEiTqjJWK1ltyOGnGGONNqktmXe9vKTzafXTNQRp7cLRXkjf9ZIv\nNrYs9F4DkpGxUQ5p75svTDiuumJwgvy9+WymyorZjNdcZUV6I/eS4X4H/SHNc68ou1xlisu9\n5E6VH1s7r3EveVllCBm+abfQC2nmUPnikOu891q+XDNjxow5F7XKkXI115goJ8+OqZcle8bm\n2zH1olRgy9hCo5+woJL8fXfdJpUVyjc4ey/0iuwblSULQy/5TFlRgXwPFlzy3WumLGkbfMC5\ntLS0LcoKNv3wHr9+XldUL2VFlQsqX9DFi5dMQpoySr446Ur1XsuXwxs3bvy4niEokjs1oGmr\naQUqC7Jqur/x7lJZ00r53mwcekWq5xv8Dc/MrKwTR48e2rZt29oVK1Ysmjdv3qz/JHqW9Boy\nZMirvXnd2vGaP+5O/m5sXFZZUfaWW26pUc2T55jxzqv+z1L81RqDNH2YxH8SZXmv5cs98pe4\nJlury1KO5hoTXb5qx9SLkj1j8+2Ymi3ZM7bQb2+X+4ZOG5X1np8UbH3oJX9WVtwkb1Lq4d1b\n1678dv7MGZMnvzlq2MCB3bt3bNXqL54hVapWrWj2O7+oMtVv9nR7I2/3N3P3kGfJg62Ke667\nX42VFdXOq/6vdMEkpG/6yhebEwu8194F6jckCfeR3GMj+D7SI8q31YzQ6yd5vje/9Tt68uC2\nVUu/nj51wshh/XvUU1aU00ZyfZ26DRs2fLhp06cTE1t17dq1T//+/Ye+/vrr4yZMeKO8e0Xj\nr7/+eulKnvzzatuhXYfdKY/mfawM+Urlq2vuXlFf5T7S3uruJa+p/69k9j7ScVeGJH08vOga\nkIyMDV9I6T9+dzjgkB+k03HK92b30COWKStiP/10whsDe3VyNb339puqxgZTUqvO3fc0bZbY\nsutz/fsPl218kDydo/hp27ZfDy/md7TU7+F/VkVe0WCf37GAR+1eSWCs/Fi1IYeekIfcvVFt\nyaKbGYvvFfhAe0CmIH2xXZKGvX5yQ5tNxdeAZGBs2EKad5N8j2KI/zE/SKc8IroFfuaRLcvm\nfPTW4J5tH7+/UiCYqjfe1rDp0226vjj49fEfJM/7ZuXGpFvkHzbTVU9l70tPPPOV+tnuf3/4\n5wHf34G/R9o1bcav6kNo42crNX6PlL7xhwMaQ8xBajlNki6+07kff0aD9xqQDIwNV0hbFQPv\n+R30v2l3v2IjiW8f37R4xtuv9GjZ5M7flfFxU6VyHEto1OvVN9+f/vWKn/cEfyA85ajag9+m\ni5RnNmimeU6AFMaQXlQk3OJ30B/SmnJ8Rd1//v3+2sWPf5W/qeFjbXu9Oi5p3vKthziQo5r/\nlMOfIgRIlowNV0gtFBdl/A56IB1c9dlbL7Z98OYED56YWnc+2q7vGx98uWJnqol/CpAASXxs\nuELqoxip7Xvs8JZpbzz31zsqePzUaPBkt1cnzFq++5TYPwVIgCQ+NlwhrXPfbmOj3Tv7F78/\nIPFu5fkxrMIdj3V97cNFW7WeAao7QAIk8bGlBOnYpIHvBT647V8y/51J99XJQ1s1rKzcyqv7\naO/Xk1cctO4sPQESIImPLR1IK6/lN82WqKxIXzOxZaPa7l8Vlbn1iT7vfLX9TFj/GUXJAEkw\nQNKEdKqu+2fMjcEfGtj9xfCWd7ofxK58b8fXZ24q/qsBQAIk0aIK0lLPwwXzAo5nbkl+qZn7\niablGnYcNW934OcBEiCJFlWQZnsgfexz7NdZLzdzP5zwu8cHTtsc/FkwgARIokUVpO0eSGuV\n3dM/jHHdxPdvaj58jtrTaAAJkESLKkj0nNtRe3nr2NxBTfhzp69p9vLnWk9FAyRAEi66IKUN\nqMQqvLBvzouN+ONy9TpN3qjr6W6ABEiiRRckopPJL90fz1j8vc/P3K97LCABkmhRBWnTuMfl\nm3Oxdz//xTFDYwEJkESLGkipn3evw5/W3XPmIcNjAQmQRIsOSL++92Q5xir+7Z1tpsYCEiCJ\nFgWQ1g9vHMtY3X4LVV9rUS1AAiTRIgnSikcqX9v+F/9jG169Q75XdP/rG0QGAxIgiRZBkJS/\nXK3j8zjCpmF3Mpbw+KR9oT9JV4AESKJFEKRHlactDPbs/jr6D4yVeex99b+a0BUgAZJoEQRJ\nefE19gTfPvnx4/Es7tGJ1vwVESABkmgRBMnzWtmJREu6VGHsrlF7rZoMSIAkWgRB6qlAem9c\nfcZq9f3JwsmABEiiRRCklAbc0c3lWfw/Zmu8lY/BAAmQRIsgSHTq7SY1GLtp2B6rBwMSIIkW\nQZB+eakai3tqtsYrVJsJkABJtIiB9EPrMqzaUGNPRtUbIAGSaJEBKWP6/YzVG58aNm99qSdA\nAiTxsVZCOj3ldhbz6JeZYfQesnoCJEASH2sdpPTxv2fxbZWXYQAkQBLOmZBOjL6OJTy71bMH\nSIAknBMhpb9Zk5XvXfyEb0ACJOGcB+n0pBtZhf6+L70ASIAknNMgZSbfysp09X8+HSABknAO\ngzT3Dyyu446Ag4AESMI5CtKGx1iMq+RfvQISIAnnIEgHe8azB1cE+QAgAZJwjoF0akINduOU\noC+RCkiAJJxTIM2txyoND/H+k4AESMJFOaQTA2+r1WwJ7Ulksc+EfNV7QAIk4aIbUsYj7j9/\n7VmZNfox9CpAAiThohvSJ563NqoyVu3vjQAJkISLbkjPeyCpv9AjIAGScNENaYDiKPaE6ipA\nAiThohvS1wqkh9VXARIgCRfVkJbc4nZUQ+PdJAAJkISLYkgn+8fGdP28a/OhWi8/DEiAJFz0\nQlpal9VdrGslIAGScNEKKWN4fEzPVH1rAQmQhItSSPuasZpz9C4GJEASLjohfXUta6r/xfAB\nCZCEi0ZI6f1j41/J0D8WkABJuCiEtPmP7JaVRsYCEiAJF32QPq/C2h83NBaQAEm4aIOU8XJM\nwkSDYwEJkISLMkjHnmY3LDM6FpAASbjogrS6DnvI+JuRAxIgCRdVkJLKx/Q6ZXwsIAGScFEE\nKaMfqzjNzFhAAiThogfSib+zuutNjQUkQBIuaiDtbcTuP2BuLCBFEaRsrS5LOZprTHT5qh1T\nL0r2jM0P+aGNN7E2mSbHSqHHilRoy9QCe8bmX7Bjap50UWvJBYshXdYqV7qqucZEuXl2TM2R\n8u0Ye6Ug1Ee+rRIz/JLZsVLIsUIV2jK1QPtbxdTYK3ZMzZe0x1oMSfOnJG7aqdy0GxuXMNX8\nWNy0i6KbdprnBEghIWX0YrUM/xbWJ0ACJOGiANKpNqxe4Du1GAqQAEm4yIeU+lfWcH+Q4/oD\nJEASLuIhHf4T+8sxsbGABEjCRTqkvXexv4d4kwndARIgCRfhkDbexLob+FvY4AESIAkX2ZBW\nVWeDxccCEiAJF9GQVlSNGWfBWEACJOEiGdLKajHvWDEWkABJuAiGtLRy7AeWjAUkQBIuciEt\nqRw36dOhkw6KjwUkQBIuYiF9VyluzK2MsarzhMcCEiAJF6mQFlWMm/qg+01bqpv8K6TiAAmQ\nhItQSN+UKzN9h+d9LYXvKAESIAkXmZCWVSozk370QBojOhaQAEm4iIS0tlpsEtHRMgqkuaJj\nAQmQhItESJuvjfk3vx7odtQETxGyIEASLAIh7arNRrp3T79akcW3E36sAZAIkISLPEj767Gh\n3gMZO9MtGAtIgCRcxEE6WJ/1tngsIAGScJEG6XhD9lymxWMBCZCEizBIp5qxdsIPLgQGSIAk\nXGRByu3IHrbiXpF/gARIwkUWpCHs7hTrxwISIAkXUZAmst//asNYQAIk4SIJ0hfxNX62YSwg\nESAJF0GQfqhQbrX1UwmQeIAkWORA2lIzbqaedzU3HiABknARA+nA79m7Ot7V3EyABEjCRQqk\n9AdYfx3vam4qQAIk4SIF0rPs7xmARIAESEKNZA1StN/V3GSABEjCRQakuXHVtxEg8QAJkEy3\nvkrCIvdYQAIkQDLdwZvZJGUsIAESIJntVBM2wDMWkAAJkMzWiT9gp4wFJEACJJO9637AThkL\nSIAESOb6PqHqtqKxgARIgGSq/TfEzikeC0iABEhmOt2EDfMZC0iABEhmep496fMSDYAESIBk\nphkxdY/4jgUkQAIk4/1cudwqv7GABEiAZLjjdwS8ZQsgARIgGS6zOesVMBaQAAmQjPYWu/9U\nwFhAAiRAMtiaslV3BI4FJEACJGOl1IuZWWIsIAESIBmrLetTciwgARIgGWoCa1jyRb4BCZAA\nyVDrylfcFGQsIAESIBnoxB0sKdhYQAIkQDJQJ9Yt6FhAAiRA0t9H7M6TQccCEiABku62V66w\nIfhYQAIkQNJbxp/ZxBBjAQmQAElvw9jTocYCEiABks5+TKi5L9RYQAIkQNJX6m0xX4QcC0jO\nhlS4bNzUFGVz0Wvutkuf86tkQAqsO+sReiwgORtSUpcv3m533L25eaZckmuv9NIoeWM5IAU0\nL+a2VGVr7V8rVX7S/9E7QHI2pKxWOyVp5MTiA8nvSVLHHX5rNM/JGZAOXhv/g7K14xomV/0X\nv7GA5GhIy7oUStKP7Yr2j3U9L11wpQNSyVxsuGerPXP3jN9YQHI0pNmvyBf7XZe8+68tlKQj\nruQXX/k6n++m79u37+BZrS5KlzXXmOhirh1TsyVzYyezB37zbN6lQGro++HzeYLnFTwp35ax\nBbZMzS+0Z+x5O6bmSppjzxuA9OEo+SLVlebZ3d0lR5LWuyZu+rbTFL4/vHHjxo9rDnFCKZUr\nH/NuP6hAeqQ0zwfZX0HRljakT0bIF0ddGZ7d4bPli3NH5Iu1idzjN2PHjp14RaurUp7mGhPl\n5tsxNUcyM/byE+yjop2xCqTxfmMLRE8saJI9YwvtmSrZMrYgx46p+ZL2WAOQvnpJvtjluqrs\nnUk87f1AluuAd1Pz5qYD7iONZ49kFu2cbsYdPXnGbyzuIzn6PtLOlvLdozn9PHtfDuKXSxbI\nF8dcRWM0zyn6Ie2oVHmnz27mtJ49Z2T6rQAkZ0Mq/Od0Keu5byTpF/5I3aCZ/NiKxI352a+P\nLFqjeU5RDymzWajnqhaPBSRHQ5KOde3aaqJ8r6rlNEm61GKzG9eclu0SRxRP0TynqIf0ru8N\nuxBjAcnZkKTClHP8aq987+jKnhzlWO7RLJ8VmucU7ZACbtgFHwtIDoekneY5RTkk+YbdJO2x\ngARIgKTau6yp1g07QOIBEiCppOeGHSDxAAmQQpfZVMcNO0DiARIghe4DPTfsAIkHSIAUsoM1\nym3VNRaQAAmQQteajdQ3FpAACZBCNpc1OKW9igCJB0iAFKLUOnErdY4FJEACpFD1Zv30jgUk\nQAKkEC2Pq52idywgARIgBe/0H9k83WMBCZAAKXjDWQf9YwEJkAApaJvLVT+gfywgARIgBe0R\nNtXAWEACJEAKVhJramQsIAESIAXp2HUJG42MBSRAAqQg9WSDDI0FJEACpJKtia+damgsIAES\nIJUo4z4W8q2Qgo8FJEACpBKNZ80NjgUkQAKkwA5WL7/D4FhAAiRACqwj+z+jYwEJkAApoMUx\n9fX9FZLPWEACJEDy7/RdMd8YHgtIgARI/v0f62h8LCABEiD5tbdyVf1PVi0aC0iABEh+tWfj\nTIwFJEACJN+WGn+kgQCJB0iAVFzGvWyhmbGABEiA5NME1sLUWEACJEAq7nCNcgaf0+AZC0iA\nBEjF9WTDzI0FJEACpKLWlamTZm4sIAESIBX1F/aZybGABEiA5O1TQ6/T4DcWkAAJkDydrF3G\nyOs0+I0FJEACJE+DdL/Ud8mxgARIgKS0s1zNY6bHAhIgAZJSazbB/FhAAiRAcvd9TP3T5scC\nEiABEi/zT+wrgbGABEiAxJvKnhYZC0iABEjkfuh7k8hYQAIkQJIbwp4XGgtIgARIRLsrVDsk\nNBaQAAmQiDqwt8XGAhIgARKtir3dxN+X+44FJEACJHqQzRUcC0iABEjT2F9FxwISIDkeUnrd\n+PWiYwEJkBwP6U32nM9e2oKkVcbHAhIgOR3SkeoVfy3eW1aHMdbsiNGxgARITofUz/cFT47c\nwHitjY4FJEByOKQdZa87Ubw3xe2IxR40OBaQAMnhkNqw9332RiiQ2FqDYwEJkJwNaVXsXWd8\ndj9WHMUbvJMESIDkcEiPsHm+u6m3uyH1MDoWkADJ0ZA+Zw/7L9jYWHb0zEmjYwEJkJwM6Uz9\n2MBfGmVu/naf8bGABEga/aaVDElzjYmyr9gx9ayU47s7nj1jzdirlowJTMqzZWyBLVPzC20Z\nm5dlx1QZktaSsxZDytMqXyrQXGOifHum+p1s1u8qpFgzttCSMYFJNo21ZWqhTWPtmap9slct\nhqT5UzKCb9oNZgMtGoubdrhp51xI+ytVN/2SkAFjAQmQnAvpn+wtq8YCEiA5FtK2hNrpVo0F\nJEByLKQWLMmysYAESE6FtCKmQYZlYwEJkJwKKeDJQWJjAQmQHAppLmti4VhAAiRnQsq4O2aZ\nhWMBCZCcCekj1tLKsYAESI6ElH5zmc1WjgUkQHIkpDcN/8mR+lhAAiQnQjpew/eVgywYC0iA\n5ERIr7DB1o4FJEByIKSDla16tqp3LCABkgMh9WGjLR4LSIDkPEg7E643+qIMWmMBCZCcB6kT\nm2T1WEACJMdB2hR/62mrxwISIDkOUnOWbPlYQAIkp0Gy8s8nisYCEiA5DVJTC/98omgsIAGS\nwyCtYA/aMBaQAMlhkB5gS20YC0iA5CxIC9nf7BgLSIDkKEjUINbgWx/pCpAAyVmQprJOdowF\nJEByFKRTdctY+ucT3gAJkBwFaTzrY89dL0ACJAdBSruhbCogARIgCTaaPV/izZgtCZAAyUGQ\nUmpWOAhIgESAJNZQ9nKJN2O2JkACJOdAOlL1msOARIBEgCTUS2x4iXc1tyhAciKkPGdCOlCp\nxnFA4gGSCKRe65TrU21WOxNSTzY28D1kLQuQnAOpQ+yAy5JU+EnV8nscCUl5xRNAIkAiIUj5\n4yvctu7Qo+xvxww5ihpIz7CJBEjuAEnsPtLRx2PLXjfHGKOogbQ5vi5/xRNAIkAiQUi7H4ip\nVHuxQyG1ZVP5FSARIJEQpCuvlam3Ov0frMMZJ0JaF1vf/YongESAREKQEuMGXZavZlStdsCB\nkJqz6e5rQCJAIiFI/TYr1+n/cODD32ti7850bwASARJZ9AvZHOdBeoJ9qWwAEgES4SlCJlse\nc4/yAwmQeIAESOZ6hM33bAESARIBkrm+K35NSEAiQCJAMtdDbJF3E5AIkAiQTDWX/bVoG5AI\nkAiQTNWY/VC0DUgESARIZprJ/lG8A0gESARIJsq4y/dFigGJAIkAyUQfs1Y+e4BEgESAZLwz\nt8dt8NkFJAIkAiTjfcA6+O4CEgESAZLhTt0ct8l3H5AIkAiQDPce6+K7u7PXo60tf0tzAiQe\nIEUxpPTfJ2z32V1dgcn1sGBwQIAESFEN6W1/NXczd4tCLTcdIAFSNENKu6HsLz67hxRH7GXx\nyQEBktMhpS5ef8WzuW623Cr/YxEOaTTr67t7wANpgPjkgADJ4ZBWth7ao99ZZfuVXoMHD/6P\n/7HIhnSiVoV9fgfqKZDmCk8ODJCcDSm382Ipp99nyk7XTSWPRTakEay//4Hv3I5aCg8uESA5\nG9LatrmS9E1XBZArpcSxyIaUUqPigYBDP7Ws/+dxp0UHlwyQnA3pq5fliz2uXL6d6tqYNHO/\n/7HIhjSEDS5xDL+QJUAiyyEljZQvUlwZfHur64VP3kxc6nvszWbNmrUq1EzSXlIaZVWr+t+S\nRyV7ztamqTjZ0jvZfAOQpipo0vn2sYU5kjS3XY7PsQmJiYnP5mtVIBVorjFRQaHggGFsdJCp\nkujYoAmfbPDsOdl8yZaphTaNteW7S8fJFr+NmDakLwfJF3tdxQ92Z7gOBx7T/CkZpjftDlaq\nfrzkUdy0I9y0I8tv2q1pL//8Wvyse3vXXvnijOu077FIhtSPjQxyFJAIkMhySFfar5LyX54m\nSQWF0px2v0lSct/ComORDenX8temBjkMSARIZP0vZFe1e/elF89JUstp0uU32ozs/+yB4mOR\nDakXeyvYYUAiQCIbniKUsnj9Jflqzk5JKvzlmzXZPsciGtKectedDHYckAiQCE9a1V03Nj7o\ncUAiQCJA0tuOhNrpQT8ASARIBEh668QmB/8AIBEgESDpbIvyzstBAiQCJAIknXneeTlIgESA\nRICkr/Wxt54J8SFAIkAiQNJXIgv5SkGARIBEgKSr1bENMkJ9DJAIkAiQdPU3NivkxwCJAIkA\nSU8rYhplhvwgIBEgESDpqRmbF/qDgESARICkoyXsTyofBSQCJAIkHTVh36h8FJAIkAiQtPuK\nNVP7MCARIBEgaXcf+17tw4BEgESApNks9pTqxwGJAIkASavMP8asUl0ASARIBEhaTWOJ6gsA\niQCJAEmjjLvi1quvACQCJAIkjT5i7TVWABIBEgGSeqfqxm/RWAJIBEgESOpN8H/n5WABEgES\nAZJqaTf6vfNy0ACJAIkASbU3WW/NNYBEgESApFaJN7oMFiARIBEgqTVcz5ssAxIBEgGSSker\nVTmkvQqQCJAIkFQayIbpWAVIBEgESKEL/r5iJQIkAiQCpND1YWP0LAMkAiQCpJCFehuXwACJ\nAIkAKWRd2b91rQMkAiQCpFBtD/U2LoEBEgESAVKo2rEP9S0EJAIkAqQQrY2tH+pV8wMCJAIk\nAqQQPcU+07kSkAiQCJCC90PMPaFfpNg/QCJAIkAK3l/YAr1LAYkAiQApaPPYo7qnAhIBEgFS\nsDLvjVmmeyogESARIAUrmTXXPxWQCJAIkIJ05natl+DyDZAIkAiQgjSJdTYwFZAIkAiQSpZ2\nk/YrnvgESARIBEglG836GJkKSARIBEglSqlZUfsVT3wCJAIkAqQSvcoGGZoKSARIBEiB7a9U\n/aihqYBEgESAFFhP9qbP3vHF8w9qfAIgESARIAXk//d8n1RnrNzr6p8BSARIBEgBtfb9e74f\nyzLex6qfAUgESARI/q2OvTOjeK+T2xG7R/VTAIkAiQDJv0fZHJ+9pgqkmqqfAkgESARIfn3D\nHvLd7ahAaqT6OYBEgESA5FtmQ7bUd3+5ch/pI9VPAiQCJAIk3z5mLv8DH1VlLGGo+icBEgES\nAZJPp+rGbQg4dHT+HK3nCwESARIBkk/jWDcTUwGJAIkAqbiUa8v9YmIqIBEgESAVN0TP+/OV\nDJAIkAiQivq1osFnq3oCJAIkKhVIeVrlSwWaa0yUrzq1J5tkbqpNJ1tox9Q8yaaxtkwttGms\nPVO1T/aqxZB+00r+iaS5xkTZV1Q++HN8nVOmpp6VcsydjsbYq3ZM/U3Ks2VsgS1T8wttGZuX\nZcdU+SeS1pKzFkPS/ClZGjftnmDTzU3FTTvCTTvCfSRP37L79L7Yd0CARIBEgKSU2YgtMTkV\nkAiQCJCUpgY+OUh/gESARIDkLr1Omc1mpwISARIBkruRrJfpqYBEgESAxDtcrdJ+01MBiQCJ\nAInXm2m8wIlagESARIBE/JWDbkg1PxWQCJAIkORcbIrAVEAiQCJAIloa88eMoB/QFyARIBEg\nUUYjtkhkKiARIBEg0fssUWgqIBEgESClXJewRWgqIBEgESC9zF4SmwpIBEjkeEi/lK9p6u9i\niwMkAiRyPKTWbKLgVEAiQCKnQ/o+psEZwamARIBEDoeU+Sf2lehUQCJAIodDSmJPC08FJAIk\ncjakk7UTNglPBSQCJHI2pCGsr/hUQCJAIkdD2lm+xhHxqYBEgESOhpQo/NA3D5AIkMjJkBbH\n3C3yrG9vgESARA6GdOaumO+smApIBEjkYEjjWAdLpgISARI5F9Kh6pX2WDIVkAiQyLmQurKR\n1kwFJAIkciykVXF1062ZCkgESORUSJkPsS8tmgpIBEjkVEgfsr9ZNRWQCJDIoZBSrrPgSXae\nAIkAiRwKqZ+5910OGiARIJEzIa0rc2OKZVMBiQCJHAkpswmbad1UQCJAIkdC+oA9ZuFUQCJA\nIidCOlyz3DYLpwISARI5EdJz7DUrpwISARI5ENKKuFssek6DEiARIJHzIGXcy+ZaOhWQCJDI\neZDGsdbe/VU9nupr+j2YiwIkAiRyHKR911Ta7dn9gMmVFf7xBEgESOQ4SG3ZWM/evgocEquV\nJjgVkAiQyGmQFsT84bRnL5kpLRacCkgESOQwSCfrxn7v3ZvqgbRAcCogESCRwyD1Zz2L9jYp\njhIOC04FJAIkchakrWWu83kzpH+5IY0NvV5fgESARI6CdPYB9pnPbsbEe2rdP014KiARIJGj\nII0TfNvl4AESARI5CdLOiv1DZt0AAA0lSURBVNfs1l5lOEAiQCInQXqKfWTDVEDiAZJzIP2H\nNbls/VRAcgdIjoF0qFbCFnu+4wEJkMg5kDqxYcHe1Vw8QCJAIsdA+jrmznRAAiQeIJnveO24\nZUHejNmKAIkAiZwCqRvrH+xdza0IkAiQyCGQFsTcdhKQeIAESOY7cXPsEgIkHiBFDKTClHPF\n26lp+fJV6h6546UH6Tn2L34FSIDEiwxIh7u3dk3MU7a3d2zTsvdhSRrukhtRapAWxNQ7ya8B\nCZB4EQGpsO/k3BOdl7i3M1rOvHppVJ9CqedPV69ezS8tSPINO+WvYAEJkHgRAWlnYrYkzern\n3l76TKEkpbhO5rc45LdG85yshfQce1HZACRA4kUEpIX/ki92JLp//OxYJF/sd5044zq5Ycel\n0oK00HPDDpB4gBQhkJL5XaGjrv9693NHvlC4y/XMS12e3ct3vxwyZMjoHK2uSnmaa/RGN8et\n8Y7Nt2yqT7mSPWML7JiaI9kzttCeqZI9Y3PtmFogaY81AOnDUfJFquukZ/foy11PSLvGnJTy\n3u7Bf0oNb9y48eOaQ6ysOxv8P/33EApVQdGWNqRZr0r85ly28omzW77nfSz8hCtVvvxvWlra\nqSytLkqXNdfo7HN2x2nv9oUcq6b6dl6yZey5PDumZkn2jC2wZWp+oT1jz9oxNVc6p7Wk+PdC\n2pC+7ypf/NS20L3zXnf37bmznFWW67h3jebNTevuI+2vmbC6aAf3kXAfiRcR95GyWuyXpHET\n3NvbE9Pc10k9ciVpRce8/z2kzCfYqOI9QAIkXkRAkpJ6Lk3qcEyShi2WPuw2m3fudLdBX33c\nZkXREs1zsgzS2+zBM8V7gARIvMiAVLh0XNJR+XrEUmn6a+4ypHOzxkzZW7xE85ysgrSxfJUd\nPruABEi8yICkI81zsgjS6XvZf3z3AQmQeIBksAGsud8+IAESD5CMtTTu+kN+BwAJkHiAZKhj\ndWLn+x8BJEDiAZKhWrEXAo4AEiDxAMlIE1mjwDcvByRA4gGSgdaVr7gp8BggARIPkPSX1oAl\nlTgISIDEAyT9dWFdSx4EJEDiAZLukln91JJHAQmQeICkt61Vyq4JchiQAIkHSDo7dR+bHOw4\nIAESD5B01oe1CXockACJB0j6mhZzy3H5au+8ZSf9PwBIgMQDJF39XLnsj0SZ/cowduM8v48A\nEiDxAElPKXewD+SrMYxXebvvhwAJkHiApKPMRNaTX9/ghsQG+H4MkACJB0g6eoPdx59ilxGj\nQGrr+zFAAiQeIGm3ML7WbvfG9fiJFDRAAiQd7f1d/LfK1ki3o4pbfT8KSIDEAyStTj3ARns2\nM3rJjn432+/DgARIPEDSqjtLzCza+eWLbwOebwdIgMQDJI3Gs/rH1T4OSIDEAyT1vkuotll1\nASABEg+QVNteI36++gpAAiQeIKl1/E72rsYSQAIkHiCplPEk66G1BpAAiQdIKvVlDwW+aFCJ\nAAmQeIAUug9YnYOaiwAJkHiAFLIlCZXXaa8CJEDiAVKoNlWPn6e9CpAIkHiAFKJ9dTQfsHMH\nSIDEA6Tgpd7n/yTvkAESIPEAKWhnnmatM7WXESDxAAmQQtSdNdF84FsJkACJB0jBeoXdeUTn\nUkACJB4gBWkSu36X3rWABEg8QCrZnPhrdPwCyRMgARIPkEq0uHzCt/pW8gAJkHiAFNjqqnHJ\nuhYqARIg8QApoE21YibpWecNkACJB0j+7byJjdHzFRYFSIDEAyS/9tRhI3R9hUUBEiDxAMm3\ng3fofGJQcYAESDxA8ulYI/ZPnV9hUYAESDxAKu7ofayjvifY+QRIgMQDpKKO3MtanNH9JXoD\nJEDiAZK3I/ewGx7+1wH9X6QSIAESD5A8Hb7H/Qr51XcY+DJ5gARIPEBSOtxIecsW9oSRr5MA\niQdIgOTpwF2sngIpweDDDYAESDxA4u2rz7p0VCCVyTD2pQISIPEASW7Hraxb5mQFUhODXyog\nARIPkIjWXc/6ZNKZP3NHFfT/JZISIAESD5BoeXXWn1+ffK3x7e3V38MlSIAESDxAml8p7t/G\nv8KiAAmQeI6H9J8yCZ+a+AqLAiRA4jkd0juxFfS8MHHoAAmQeM6GlDmQVV9u7kv0BkiAxHM0\npLSWrPZGk1+iN0ACJJ6TIR18kDXeZ/ZL9AZIgMSLGkhXtLoq5fkf2FOPtfyv5qdplZsvPCJI\nOZI9YwvsmHpFsmdsoT1Ttb9VzFSQY8fUfEl7rMWQsrW6LOX47S+6hvU9p/lZ2mOvis8o2UXJ\nnrH5dkzNluwZW2jL1AJ7xuZfsGNqnnRRa8kFiyFp/pQMuGk3uUz8ePM/covDTTvctONFzU07\nzXPyg3S6D6u6QOhL9AZIgMRzJqT9Tditog/XeQIkQOI5EtKPtdkTet+2RStAAiSeEyFNKRfT\n3+BfHYUOkACJ5zxI6b1Y5ZniX6I3QAIknuMg7WzM7txiwZfoDZAAiec0SDOqspYnrPgSvQES\nIPEcA+nMt0kL09O6s7LjrPkSvQESIPGcAml7A8ZY7Xqs3k8WfYneAAmQeA6BlHmf8tImbVKs\n+hK9ARIg8RwCabXn5R/F/ogvWIAESDyHQPrKA+kDq77CogAJkHgOgbTVA2mJVV9hUYAESDyH\nQKLWbkePWPaEhqIACZB4ToF0rGMMY82F/x62ZIAESDynQCJK3WL5I3Y8QAIknnMg6XpXc+MB\nEiDxAEkwQAIkHiAJBkiAxAMkwQAJkHiAJBggARIPkAQDJEDiAZJggARIPEASDJAAiQdIggES\nIPEASTBAAiQeIAkGSIDEAyTBAAmQeIAkGCABEg+QBAMkQOIBkmCABEg8QBIMkACJB0iCARIg\n8QBJMEACJB4gCQZIgMSLGkjntDq58ZjmGhNduGzH1MyN++wYm51jx9RzG3fZMvaqLVO3/2zL\n2NzzdkzduzFTa0m2xZA0+6HxF/+bf8iKjjYeXdqnYKD7upX2GRio/cOlfQYGGtY4Xf9iQCoZ\nINkWIAkGSLYFSHYFSIIBkm0BkmBnVpz43/xDVnRxxa+lfQoGWrG5tM/AQD+vKu0zMNDuFZf1\nL/4fQUIougMkhCwIkBCyILshpQzr8MrBINthWeaYLl0nnVe2R7nkRpTu+ag2n59gS2W7cFaP\n7sn5pXs+am1wuZvk3vE98XBs5WF+ub5/p7fOeo74bofKZkiXO07c91GbcyW2w7JL3Ubs3drn\nDWXnhenbt28/XLonpNpHI+UT3KFsz+m8cUvX5NI9H7XOyqe6fUvnn9w7vicehmU/u0G+3Ju4\nYM/QgcoR3+2Q2Qzpu+4FUmHfOSW2w7JVbS9L0m7Xb3y7sM3e0j4djUZ+WbSZ33WZJK3pcKUU\nz0ZHCzz/ifI58fDr9OQuLg7prfGS9Fui8uit73bIbIb01mT54uMRJbbDsuX/li9OuI7w7SwX\nZeeU8vmo1/enK95neqW4SJIuuML7QfvMzmeUDZ8TD78yFy5sxSF15o/T/0v5r77vdshshvQK\n/z3s/H4ltsO3T5/J5Vf7XINciSN/K+2zCV1hq0GJrhf2ubd3JBbIl23Xl+4ZaTQpSbn2PfGw\nrJ0M6aprt7w1yn3KvtuhsxlS3wXyxdIuJbbDtStJicq345r2Ky4d7T+0lE9Hpd/afHqW3u3s\nvsu5ph2/fHZp6Z6ReuntPP9V8j3xsIxDynLxu8fvvsP3fbdDZzOkQbPli/l9S2yHab/8s/ce\nn929roxSOxVd5bT7kV9tSyyUL9uuKeWzUe3Dd333PCcelnFIucpPoSl833c7dDZDGsP/+eRh\nJbbDsyWt5uT57p9zhfWj9XIvzOeXR11ZknTZ/f93uJbbcavfvnLiYRmHJHXg/1UaMNt9wHc7\nZDZDWtRL/m9l/zkltsOylBabiraTx8gXO1tcKr2z0Wj9C+cl6VJb97dnfhf5v+8b24fzo3Yb\nOnh/zeV74mGZGxJ/ZCw7UXnk1nc7ZDZDyu48/fTs9iRlzkor2g7bpnffvUcuR/piu7S7xccH\nfu75SWmfUuiyu4zcuXfEwAJppXzXaHb3A4d6fVrap6TWFOWxb/lki048XHND2tVqddqYAZK0\nZ1Ze0bZqdj+z4fhQ97MZ9rm2Fm2HbaOV37+fkFpOk/8XHNqh79wwfrKAlPlWl26TsyVpxEBJ\nKvysR/dPw/dbU663csuIn6z3xMM1NyRpbf9OY89J0lzXlaJt1fBcO4QsCJAQsiBAQsiCAAkh\nCwIkhCwIkBCyIEBCyIIACSELAiSELAiQIr273ijtM0ASIEV+lV4s7TNAEiBFeudWl2+1Opyf\n9u2UACmyW83kjpT2WSBAivhw0y4sAqRID5DCIkCK9AApLAKkSA+QwiJAivQAKSwCpEgPkMIi\nQIr0KvUu7TNAEiBFfg1qdj9T2ueAACniW/9sizB/PVhHBEgIWRAgIWRBgISQBQESQhYESAhZ\nECAhZEGAhJAFARJCFgRICFkQICFkQYCEkAUBEkIWBEgIWRAgIWRB/w+b2wgWgSTB0wAAAABJ\nRU5ErkJggg==" }, "metadata": { "image/png": { "height": 300, "width": 420 } }, "output_type": "display_data" } ], "source": [ "# intervalos de tempo: uma sequência de zero a dez em passos de 0,5\n", "time <- seq(0, 10, by = 0.5)\n", "# condição inicial\n", "x0 <- 0.1\n", "## A função a ser integrada (expressão à direita da derivada acima)\n", "f <- function(x){x * (1.-x)}\n", "## Um vetor vazio para armazenar os resultados\n", "x <- c()\n", "## Armazena a condição inicial na primeira posição do vetor\n", "x[1] <- x0\n", "\n", "\n", "# loop ao longo do tempo: aproxima a função a cada passo de tempo\n", "for (i in 1:(length(time)-1)){\n", " x[i+1] = x[i] + 0.5 * f(x[i])\n", "}\n", "\n", "## plotando com ggplot2\n", "library(ggplot2)#carrega cada biblioteca uma vez por sessão R\n", "p <- ggplot(data = data.frame(x = x, t = time), aes(t, x)) + geom_point()\n", "analytic <- stat_function(fun=function(t){0.1 * exp(t)/(1+0.1*(exp(t)-1.))})\n", "print(p+analytic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Por que usar bibliotecas científicas?\n", "\n", "O método que acabamos de usar acima é chamado de *método de Euler* e é o mais simples disponível. O problema é que, embora funcione razoavelmente bem para a equação diferencial acima, em muitos casos não funciona muito bem. Há muitas maneiras de melhorá-lo: de fato, existem muitos livros inteiramente dedicados a isso. Embora muitos estudantes de matemática ou física aprendam a implementar métodos mais sofisticados, o tópico é realmente profundo. Felizmente, podemos contar com a experiência de muitas pessoas que já criaram bons algoritmos que funcionam bem na maioria das situações. Estes algoritmos estão já disponíveis na maioria das linguagens de programação. Aqui vamos usar uma ótima implementação disponível em R." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Então, como... ?\n", "\n", "Vamos demonstrar como usar bibliotecas científicas para integrar equações diferenciais. Embora os comandos específicos dependam do software, o procedimento geral é geralmente o mesmo:\n", "\n", "1. Defina os valores dos parâmetros e a condição inicial \n", "2. Escolha um intervalo de tempo ou uma sequência de tempos em que deseja a solução calculada\n", "3. Definir a função derivada na linguagem computacional (o lado direito da equação diferencial)\n", "4. passar a função, seqüência de tempo, parâmetros e condições iniciais para uma rotina de computador que executa a integração.\n", "\n", "### Uma única equação\n", "\n", "Então, vamos começar com a mesma equação acima, a equação logística, agora expressa com os parâmetros para taxa de crescimento ($r$) e capacidade de carga ($K$):\n", "\n", "$$ \\frac{dx}{dt} = f(x) = r x \\left(1 - \\frac{x}{K} \\right) $$\n", "\n", "E vamos usar o caso em que $r=2$, $K=10$ e a condição inicial $x(0) = 0.1$. Mostramos como integrá-lo usando R abaixo:\n", "\n", "#### 1. Defina os valores dos parâmetros e condição inicial" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [], "source": [ "# paramentros: devem estar em um vetor\n", "parameters <- c(r = 1.5, K = 10)\n", "\n", "# condições iniciais: também devem estar em um vetor, mesmo que seja um só valor\n", "state <- c(x = 0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2. Escolha um intervalo de tempo ou uma sequência de tempos em que deseja a solução calculada\n", "Note que estes são momentos no tempo para os quais você quer a solução. **Não são os intervalos de integração**, que são definidos internamente pela função de integração." ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "## seqüência de tempo\n", "time <- seq(from=0, to=10, by = 0.01)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3. Defina uma função em R para a EDO ser integrada\n", "\n", "Vamos definir uma função em R para o lado direito da equação diferencial.\n", "Para esta função ser reconhecida pelas rotinas de integração da biblioteca que vamos usar (*deSolve*), esta função em R deve calculas os valores da derivada em um tempo $t$.\n", "Há muitas maneiras de fazer isso, mas o formato recomendado é:\n", "* Faça uma função com três argumentos: seqüência de tempo, variáveis de estado e parâmetros, nesta ordem.\n", "* A função deve retornar uma lista com resultados da função a ser integrada.\n", "Para fazer isso use `with(as.list(c(state, parameters){ ... }` dentro da função R.\n", "Inclua entre parênteses a(s) função(ões) a ser(em) integrada(s)\n", "e então feche retornando a lista dos valores calculados." ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [], "source": [ "## A ODE logística a ser integrada\n", "logistic <- function(t, state, parameters){\n", " with(\n", " as.list(c(state, parameters)),{\n", " dx <- r*x*(1-x/K)\n", " return(list(dx))\n", " }\n", " )\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 4. Integrar a função\n", "\n", "Agora chame a função da biblioteca \"deSolve\" `ode`. Para realizar a integração, os argumentos básicos da função 'ode' são\n", "* `y`: o vetor de condições iniciais\n", "* `times`: o vetor com a sequência de tempo\n", "* `func`: a função R como descrito acima\n", "* `parms`: vetor de valores de parâmetro (nomeado)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [], "source": [ "library(deSolve)# Carrega a biblioteca para integração, basta chamar uma vez por seção do R\n", "## Executa a integração\n", "out <- ode(y = state, times = time, func = logistic, parms = parameters)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "O objeto resultante tem os valores da integração\n", "em cada ponto de tempo no vetor de tempos" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
time | x |
---|---|
0.00 | 0.1000000 |
0.01 | 0.1014960 |
0.02 | 0.1030141 |
0.03 | 0.1045547 |
0.04 | 0.1061181 |
0.05 | 0.1077046 |