aula-22-metodo-euler.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ipub": {
     "ignore": true
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Método de Euler \n",
    "\n",
    "As considerações anteriores sobre malhas numéricas são fundamentais para sugerir a notação que utilizaremos para expressar os métodos numéricos que estudaremos. O primeiro deles, chamado _método de Euler_, é considerado o método numérico mais simples para resolver PVIs. Embora não seja muito eficiente, é o ponto de partida para a compreensão de uma enorme família de métodos. \n",
    "\n",
    "Ao longo do texto, denotaremos por $y(t_n) = y_h(t_n) = y_n, \\ \\ n = 0,1,2,\\ldots,N$ a solução aproximada de um PVI.\n",
    "\n",
    "## Definição do método de Euler \n",
    "\n",
    "A derivação do método de Euler inicia-se com a seguinte aproximação para a derivada: \n",
    "\n",
    "$$y'(t) \\approx \\dfrac{y(t+h) - y(t)}{h},$$\n",
    "\n",
    "conhecida como _aproximação por diferença finita avançada (ou progressiva)_. Se aplicarmos esta definição ao nosso PVI padrão, em $t=t_n$ teremos $y'(t_n) = f(t_n,y(t_n))$, donde segue que \n",
    "\n",
    "$$\\dfrac{y(t_{n+1}) - y(t_n)}{h} \\approx f(t_n,y(t_n))$$\n",
    "$$y(t_{n+1}) \\approx y(t_n) + h \\, f(t_n,y(t_n)).$$\n",
    "\n",
    "O método de Euler toma esta aproximação como exata, de modo que o esquema numérico resultante é \n",
    "\n",
    "$$y_{n+1} = y_n + h \\, f(t_n,y_n), \\ \\ 0 \\le n \\le N-1.$$\n",
    "\n",
    "A estimativa inicial é $y_0 = Y_0$ ou alguma aproximação de $Y_0$. Quando $Y_0$ é obtido empiricamente, seu valor é conhecido apenas aproximadamente. A fórmula anterior permite o cálculo sequencial das iteradas do método de Euler $y_1,y_2,\\ldots,y_n$, aproximações para os valores exatos de $y$ nesses instantes.\n",
    "\n",
    "## Interpretação geométrica \n",
    "\n",
    "A figura a seguir ajuda-nos a interpretar o método de Euler geometricamente. A aproximação numérica da curva exata (azul) é feita por meio de retas tangentes. O valor $y(t_{n+1})$ é erroneamente computado (comprimento $BD$) e excede o valor exato (comprimento $BC$) por uma quantidade (comprimento $CD$). Isto é, a equação da reta tangente a $(t_n,y(t_n))$ é $r(t) = y(t_n) + f(t_n,y(t_n))(t-t_n)$. Na verdade, $r(t_{n+1})$ coincide com o ponto $D$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(2.8, 2.2, 'D')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD8CAYAAADUv3dIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAe1UlEQVR4nO3de5xVZb3H8c+Xu6ioiaIWinhDySS8JGZIysnMyLyQJ9Li5C1vaerRNCsvlXkNr5nl0cwsTT15QzMNykTUVEQCJPOChnfBC2iC/M4fz5ozwzAzzDCz11p77+/79VqvDbPWzP75cvv1mWf91vMoIjAzs3x0K7oAM7N64tA1M8uRQ9fMLEcOXTOzHDl0zcxy5NA1M8uRQ9fMLEcOXTOzHDl0zcxy5NA1M8uRQ9fMLEcOXTOzHDl0zcxy5NA1M8uRQ9fMLEcOXTOzBlLPSr+FQ9fMTOqOdCTwFNJ6lXwrh66Z1TdpW2AqcAkwB+hdybdz6JpZfZK6IV0EPAR8BPgy8Bkinqvo23qPNDOrW9I1wJvAqUS8mctbOnTNrG5ImwEXAt8mYjpSNyKW5lmCpxfMrPZJfZC+DzwBfBLYBCDvwAXokfcbmpnlShoNXAZsBvwWOI6IF4sqx6FrZrVuFCDSTbI/FlyL53TNrMZI3YHDgTlE3I3UB4CI9wqtK+M5XTOrHdJ2wIPAxcBYIIVtSQIXHLpmVgukNZAuJvXcfpjUc3tosUW1zKFrZrVgb+BI4FJgCBG/paRzp57TNbPqlHpuNyfiDqRuwFAinii6rBVx6JpZdUk3xk4CTgZeAjYjYnGxRbWfpxfMrHqkntvpwGnA/wIjqilwwX26ZlYtpI8CfwSeoiQ9tyvD0wtmVl6p53YHIh7I/j4WuK1MLWAd5dA1s3JK69xeDgwDtiTiqWIL6hqe0zWzckk9tw3r3A4EDgT+WWxRXccjXTMrj9SZMBvYkLRIzalELCi0pi7m0DWz4knrEfFS9udDgceIeLjYoirD0wtmVhypN9L3gGeR9gAg4opaDVxwy5iZFWX5dW6nFVpPTjzSNbN8SUK6hNRz27DO7ZeLXFg8T57TNbN8pJ7bIGIp0nhgEPDjau65XRkOXTOrvMae2yuJuLzocork6QUzq5zle25fLbiiwvlGmplVhvQ54BfAetRoz+3KcOiaWaUsBuYBe9VyC1hHeU7XzLqG1Ju0zm03Ik7LvtaNiKVFllU2ntM1K5CkUZKiyfGBpPmSZkj6paTPSlLRda5Q6rl9Ajgd2ISGmh24y/H0glk5/AaYSOpbXR3YAvgi8FXgHkljo4zzodIA4AJgHGlRmt2JuLvYosrNoWtWDo9GxLVNvyDpOOAc4DhSKO9RRGErsDbwBdIIt+56bleGpxfMSioiPoiI44G/Ap+VtHPRNQEgDUc6DYCImcBAIk5z4LaPQ9es/K7MXvcstAqpH9KFwMPAN5DWAXAbWMc4dM3Kb3r2unkh757WSvgSaZ3bo4GfAkOIqPsHHVaG53TNyu+t7LVfQe+/JukR3mdwz22neaRrVn4NYftWm1d1pbTO7aFZn+18YCRpg0gHbic5dM3K72PZ65O5vJu0G2lK42fAaAAiZhDxQS7vX+Mcumbld1D2ekdF30VaD+nXwD1Ad9xzWxGe0zUrKaX1Z88GdgYmRsT9lXwz4FZgG+AMUs/tuxV7vzrm0DUrh+GSDsj+3PSJtI2Au0lPfHU96ePAHCIWIh0NzCdiTkXeywAveGNWKEmjgElNvrQUeAd4Afgb8JuIuKsCb7wGcCZwJHAGEad3+XtYixy6ZvUkTSOMBSbgdW4L4RtpZvXlh8D1wIvAJ4g4yoGbL490zWpdWue2LxHzkbYitYFd6hawYjh0zWpZ6rm9DJhGxP5Fl2OeXjCrTdIApGtp7Ln9n4IrsoxHuma1Jo1ubwJWIfX5nuWe2/Jwn65ZrZB6ELGEtG3On4CTicjn0WFrN490zaqRNA74EbAhqaf378BqwC7el6zcPKdrVm1S4P6c9LSagIHAZ4F3gd4FVmbt4JGuWbWRniUFbnPPETEo32Ksoxy6ZtVGWkoa4TYXRHRrdu3kdCZGVbosax9PL5hVC2lXpLuA51u5Ym6e5djKceialV1jz+29wKbAxcCiZlctAk7JuzTrOLeMmZWV1A04DDiL1HN7Jg09t9I8GrsX5gKnEHFdCz/l3rzKtfbxnK5ZWaXQnQwsBo5wz21tcOialYnUDzgV+AkRLyKtCbyJ/0OtGZ7TNSsDSUhjgVnACcDuAEQs6FTgSnci3dklNVqXcOiaFU3aBJgI3AC8DOxIxNVd9NNXyQ4rCYeuWfG+A3wSOBbYgYiHii3HKslzumZFkHYFXiFiBtI6QE8i5lXgfSYDfjiiRDzSNcvTsj233wEg4tWKBK6Vkvt0zfIgdQcOpbHn9gzgxzm88+05vId1gKcXzPIgHQlcQlrn1j23dcyha1Ypqed2IyKeQFoF2BO4yT239c2ha9bVJAH7AROA94Atsh0diqhlMuAbaSXiG2lmXWn5ntsvFxa4Vkq+kWbWVaRtgKmktRKOAS5z4FpzDl2zzpIGEPEyaUPIc4AriPhXwVVZSXl6wWxlpZ7bXwFPIq1HxFIivu/AtbZ4pGvWUWnJxYae21WBs4E3C62pdTcUXYAty90LZh2RWr8mAZ/APbe2Ejy9YNYeUvqtMOJd0s2yA4DRpQ9cqS9S36LLsEYe6Zq1JfXc7gucB4wh4omCK+oY9+mWjke6Zq1p7Ln9HTAf//diXcAfIrOWSCcCM2hc53Z7Ih4vtCarCQ5ds5b1A24DtiTiwqIfcpDUV9Kxku6T9IakxZJeljRR0ng1zDlb6flflBmknls4H7iOiInA94hYWnBVAEjaFLgD2By4h9Sq9hqwLjAauArYCjixqBqt/Ry6Vt8ae25/DPQFHgQoUeCuQloTdzCwb0Tc3OySsyVtD2zfyo+4uoLl2Upw94LVL2kYcDkl7rmVdDRwEXB2RHy76Hqs8zzStXq2LbAxqef2upKuc7tf9nrFSn231B+AiNe6qB7rJI90rX409tz2IeLabGqhHxELii2sdZJeB3pGRL+V/AGTAffploi7F6w+LNtzexCSsgVqFhRb2Ar1A94quoh6JFUmHz29YLVN6g38N2nn3SWknttLSzqV0JK3gNWLLqJWZcG6ITAkO7Zo8ufLgDO7+j0dulbrdiT9h/M74FtVuOziDGCkpMER8XTRxVQriT6klrshwJbZ0RCyfZpcOh+YBdwJPFaJWhy6VntSz+0uRNxAxJ+RPk7EtKLLWkk3ASOBg4FTCq6l9CTWoDFUmx6DAWWXBfAMMBu4lxSyT2Z/fy2Civ4W5BtpVjvSjbFDSD23vYANiXi92KI6R2mFsEdJoTE2Im5p4ZptgU9ExGUt/ID9AYi4vrKV5ktibVKYbtXs+HCTy94nhemsJsdsYE4E7+ZacBMOXasNy/bcTiL13M4utKYu0uyJtLuBPwKvA+sAnwZ2B86pxT7eLFy3AoY2ObYCBjS5bCEpUGc2e30mgtLtUefQteqXelHnAm8Dx1HentuVlo14DyO1vA0FVgPeAP4G/Ba4LiI+aOEbBwIQ8Xxeta4MiX6kMP1ok2MosF6Ty94hBerfm70+H0EpniBsD4euVafUczuCiCnZ378A3EfE/ELrKpuS9elK9CbdwNqaZQN2oyaXLSIFavPj+UrPt+bBoWvVRxoMXALsAXyaiMnFFlRiBYWuhIBBpHD9WPa6NWmKpHt22WLSHOsM0k7KM0jh+mw1jVw7yt0LVj1Sz+0JwKmknttvAX8ttCZDYnVSsDY9tmbZ/uJnSMF6M40BOyeCxflWWzyHrlWHNJ1wL2lR8Wrtua1qTUavw4BtmhwbN7lsASlUrwGmZ3+eEcHbOZZaap5esHKT1gHeIOIDpLHAO0TcWXRZVWMlpxeyhwmGkgJ2GI0B27AGxFLgH8Dj2TE9e32hFuZdK8mha+W0bM/tyURcXnBF1UkaA0DEba1fwprAx0nh2vC6JY2/Cb9DY7hOy15nRLCoIjXXOIeulU/quf0p6RHeyaSe21lFllQrJNYjBevwJq9NpwfmkYL1sex1GvB0Ld/YyptD18pFOh44h9T8fxzw61rrue0KEuOAH5EWa5kLnBLBdU3OC9hgV+4d8wwbD32GwYNI6wev3+THPEUK10ez12kRvJzTP0Ldcuha8dJNsu5ELEEaTVq4+2T33LYsC9yfk7YXavAucCmpx3VbYDuyp7a68QFL6T6TFK6PAo8Aj0fwZp51W+LQtWI19tw+QcRJRZdTDSTmAgNbOb2U9AjsI8AjE9lj/E5MeWeNeHNkbgVam9wyZsVYvud2YrEFlZPEaqR51x1Im09uT+uBG8AaEbzT+APu2qfSNVrHOHQtf9J2wK9Ij4PeCBzrnluQ6E5q09qBtHDPDqRHZBt2MHgOeBhYE1irhR8xd5nAtVJy6FoRFmavexJRtyPcrJNgx+z4BGkUu2p2ej7wEPB7UtA+FMEr2fe1NKe7CK+3WxU8p2uV19hzO5yIw/7/axF104Yk0YvUojWCxqBtWORlMal74MHseAh4qq2HDLLgvQBYlxa6F5pcOBqAiHu66B/FOsmha5W1/Dq3exJR2ALSeclGsSOAnbLX7YDe2enngQeAqdnxWATvFVGn5c+ha5UhrQ6cDhxDjffcZpsbbkVaF6LhGJydfp/USfAAMAV4IIJ5XfS+wwAimNbGRcNIF7V+jeXKoWuVIa1LWmD6Rmqs5zZbl2B74FPAzqSR7JrZ6VeA+7NjCvBoBP+uUB2TASIY1cZFk0kXtX6N5co30qzrpJ7bw4GTiHgFabNaCFuJtUij151JQbsdaQ82SD2xvyOF7F9Jj8x6JGOtcuha5y3fc3sN6WGHqgxciXVJ4boLaSfej5F2kl1M2h7nQlLA3h9BVW98aflz6FrnSKNIi9MMIW0XfiwRLxRZUkdJrE8K2F2AUaR/FkhtWFOA04D7gAe9spZ1lkPXVp7UA7iC9Dmqmp7brLNgVJNji+zUW6QR7FXAX4BH6nFnA6ss30izjkk9twcCNxKxEGkL4HkiSjsCzLbx3gXYNTu2zE69RRrBTs6OaWXcsrs1EjsBRDCljYt2Il3U+jWWK4eutd+y69weXtaFxSVWJc3JjiaF7DDSnOxCUsj+idQzXFUha7XBoWsrtnzP7fHAtWXpuZXoQWrhGp0dI4CepB7ZKaSQ/RPwcATvF1VnV/NItzo5dG3FpOuBsaT521L03EoMBnYHPkMazfYjrbL1GHBPdtxfyze+3KdbnXwjzVombQy8R8SLwPeBnxAxtbhy6EcK189kxybZqbnADcDdwKQIXiumQrP2cejasqRepJ7b75LC7GtEzM6/DETqj90D+Czp4YQepE0SJwETSEH7Dz+MYNXEoWuNUs/tZaS7+zcD38n37VmDNIr9HGnqoGE/r2nAecAfgCm1NC9r9ceha4l0MGmN1mfIqec2G81uCexJCtqdSZ/J+aRR7F3AHyJ4sdK1mOXFN9LqWeq5XZuIV7MFag4Hzq1kz61Eb1LP7Bjg88Cg7NR04A7Stj1T3cq1Yl5lrDo5dOuVtA1pnVsBO1VyQXGJ/qSR7BjStMHqpN1r7yEL2gier9T7m5WJpxfqTWPP7TeBN0g9t13+f96speuL2fFJ0j5fLwK/AW4D7o2g5hczrySJ0QARtL4rhHeOKB2PdOuJNJR0M2oDurjnNpufHQbsTQrarbNTjwO3ZsejEdTNFj2V5j7d6uSRbj2QehCxBPgnaQ+uc7ui5zbbMWEEsC+wD2nPr6WkRWOOA26J4OnOvo9ZLXHo1rLUc3s8cCDSdtkNsn079yPpQVrXYD/SqHZ90uO2d5OmLW6P4NVO1W1Wwxy6tUrahbQ4TUPPbV9YuUdiJbqTOg6+RBrRrkO6ETaRtIbuHRG81QVVm9U8h26tkfqSwvarwLPA54m4o+M/hu6kG2D7k0a165JW6bqNtO/ZXREs7KKqzeqGb6TVmtR7ezdp7vaHHem5zW6GbQeMI4Xt+qQR7e2kR4In1vICMtVGSouvR/BkGxelBdojWr/GcuXQrQWp5/Zs4L+IeBGpW0f6biWGAF8mhe2mpDnaicBvSVMH71SgarO65OmFapZ6bk8jrXP7BrA58GJ7AjfbsmYc8BVgOKnr4E/AWcDNESyoTNHWVSTGAERwWxsXjSFd1Po1liuPdKuVtDdwEfARGntu32j7W+gL7EXabmd30gMLDwO/Bm7wGgfVxX261ckj3eq1D2l0+yUiHmjtomyediTpxtpY0iO4c0kj2l+1OR9oZl3OoVstGntuJxLxOHAksCh76KGFy/kI8DXgv0gLfr9D6jq4BviznwwzK4ZDtxos23PbB3iciOX6YrMVvL4AfJ20Lm030i63p5Pmad3iZVYwh26ZpeUWz2UFPbdZ69ChpJHt2sALwI+AqyP4Z271mtkK+UZamUnfJW2bcy7Nem6zUe0+wGGkp8WWALeQFiK/J4IP8i/Y8iQxEKDNZTGlgaSLvHRmSdRk6Epai7SMYG/gwIi4tuCS2i/13PYj4j6kVYCNmu5RJrEp8A1gPGlU+zQpaK+K4OUCKrZOUNoiaVKzL/8bmAf8GTgnImblXJZVUK1OL3wF6EXaeuYgoPyhu2zP7SNIOxLxLjA7W81rD9LNsz1oHNX+jLQurW+KVb/fkB5IAViFtCnnwcC+kraOiOeaf4PE/gARXN/qT5X2J13U+jWWq1oN3YNIo4dbgAmSNomIcs5tSiKt1nUhy/bchsSHSN0HRwCDSaP304Ar3FNbcx5t/huZpH+QPhf7AD9p4XsOz17bCtT2XGM56lZ0AV1N0nDSYtq/JDX9LyYFV1ntTlqp6w3StjmHiVhP4grSDbHzstf9gY0iON2BWzfmZa/e/biG1Fzokka5C4GbIuJ10h5cX1NaCKYcpF5I22Z/uxs44BGGbyuin8RdwN9JT41dC2wTwS4R3BDB4qJKtorrK6l/dgyUtAfwQ+A10v+UrUbU1PSCpD6khVtujIiGntRfkn593x24s6ja/l9jz+2HkQaJeI+01u3jwFbAS8CpwM8ieK24Qi1np2dHUzOBT0XESwXUYxVSntFf19gHWIsUtA3uAF4hPTBQHGkdpKtJDyus8hSbHCriSOA50jzu+6Q+20ER/NCBW3euAP4jO8YAJwH9gYmSNiqyMOtaNTXSJU0tvAq8IGnTJl//IzBWUv+IyD/M0kMOs4DVX2LAxVsxM+bzoSuBVUmj73OByRFdvyuvVY1/xLI79t4u6c/AVNKynf/Zwvfs146f255rLEc1E7qSNgY+DQiY08plBwAT8qoJaV0iXiHilVna8qrD+Nng+xj5DdJvGL8Fzolgem71WFWJiAclvQns2vL5dvw2VMQgw9pUM6FL6lAQcAi0uBbsD0gj4QkVr0RajdTaddQFOm7c8VywP8waC7xH6q09P4JnK16H1YIepId8liMxHiCCq1v9bmk86aLWr7Fc1UToZp0J44EnIuIXrVwzFDhN0vYR8XClCgG+SLbO7S184ZkfcOpNwNukpRQneKdcay9J/0Gagrq/lUvGZ69Xt/Fj2nON5agmQpe0otZA4Mo2rrmJNPo8iLRwd9dKgXsjsM8/GfzWAVzLVEZ8CDgDuDCCNhcYt7o3XNIB2Z97A0NJv7UtJnWzWI2oldA9KHu9ubULImKGpDnAf0r6VqRHbDtP6kHEEhFDJ3DMkLlsyEV884Ml9PwecLG3vbF2+nJ2QNo66XXSDeCzKvabmRWiJhe8yY20y3v0/sURXPb8VXx9FGka4TzSyHa59W7NupK366lOtTLSzZe0zgLWuHRNGPsi6/M8AweS2nrO9TSCmbXFI90OmqcNDlmTBRf1YEmf8zl+6Tmc+LMFrHVGBH5qyHKVbTRKBIvauKgv6aLWr7Fc1doTaRUzVH/vKXH0GXzvgofYoc9u3HvLKZy1yfxY6wgHrhUhgkVtBm66aJEDt1w80l2B17X2ak8z+NrrGDdiAt9aVyz9Uw+WnPB+9Hqs6NqsvkkcARDBZW1cdATpotavsVx5pNsaSbdpzPFL6PHa9vxtr82Z0x34fNBttAPXSuJL2dHZayxHDt0W/Fpf+fjDbPfsGG4/7zX69ziSS849gp+uF8EdXh/BzDqj7kNXYpzEsxJLJZ6TuPIavjp5CLM3vJzD/nIUlwy4NI48MYIlRddq1kDSOFh/R+i2i6Rn09+tGtT1nK7EONKmjn2bnZq5F78/+PfxxQcKKMusTVnANv/cLgIOiYjrml08GXCfbolULHTV8C97WTdExGVKbSwTWzh/dURcLak/6ZHa5n4aEdcrbSv9qxbOnx8Rt0nagrSwTHM/iIh7JA0DJsCIHaFXk8VEfgTsBEx+CT79ZAvff2xETJM0mpYfzTwsIp6UNAY4voXzB0bE80qbBR7ewvn9IuI1pUVKxrdw/nMRsUjp5shy83SR/Ycl6QTg881OvxsRe2Tnvwvs1uz86xGxb3b+LGBEs/MvRMQB2fkJpC2RmpoTEYdm568ANm92flpEHJudv5a0H1xTD0TEydn5m0g7HTd1b0ScmZ2/k7R5Y1O3R8R52fnJLK9kn73lnBIRUyTtRPogNtf0s3c7LS+C8xxwNE0+e9tm/54ega392evYZy8q9D+qOn84oleLqzdBzwFAS6FrVgatfG7ZsPkXHoFplS3FOqoupxckNvwYj/9yAWuOmkuLi/I/F8GgnMsyaxdJz0KLH9znImJQvtVYR9XVjbRHNbzfJTrqGmDmdLbZYXXeuh2WaxxfBJxSQHlm7XUKLPdQhD+3VaI+QlfSZI06cgPmvXoUlx44gimPAUNnxNZjQIeQ5sIiez0kguva+nFmRcpulrXwuQ1/bqtAzU8vvKV+G89jg1uH8ORHZzD0g6sZ/4PzOeF099uaWRFqOnSHaPagvzBydl8W9b6Eo2b9goM/81Rs+kLRdZlZ/arJ0F2iHjusztubv8cql3yGP/R+l1VOv4+RZ3t0a2ZFq63QldZZSN8Jq7Jo3Ne5kqv4+hTgaxE8VXRpZmZQKzfSpG5IBy+mx1O9eH/cjzlp6W2M+R4w0oFrZmVSEw9HBFwl+OpUdoyjufhfjzNsnwgeKrouM7Pmqjd0pdWApSL67MONm6zO21zDV28Jun09gvlFl2dm1pLqm9NNW51/EbhoFkMmbcWskcAGwH8DF/lmmZmVWXXN6UqDgFuBm19hnW6H8PNxpObwnSO40IFrZmVXPaGbVkeaGfDpizlq6of51wb3s/OdwHDP35pZtSj/nK7Ug4glwPSF9P3jMKat9xSb7Uha/u67ESwtuEIzs3Yrb+hK6wDnkBZq3l9EL2AbYAAwLoLfFFmemdnKKN/0QtZzC8wGDgCe7q1/7wdMIf1P4lMOXDOrVuUKXWkwcB9pK5IZi+kxTMTC9+n9O+AJYPsI/lZojWZmnVCu0IW3gf7A+FFM2rUXi48AzgSuBUZF8GKh1ZmZdVLxfbrSXsCBwP5EfIDUTURPUtDuB5wLnOR2MDOrBcWNdKWNkG4Ffk/aSG4AgIjVSBsH7gecEMGJDlwzqxX5h67UE+kkYCawK3ACsC0R8yQGAJOBkcCBEZyfe31mZhWU//SC1Id0U2wGcAwRc9OX2QT4A7A+sF8Ed+ZbmJlZ5eUTulJ/4NvA94lYiPQhIt5oPM0WpBFuT2DPCB6sfFFmZvmr7PRC6rk9CHgSOIY0bUCzwN0MmJTVMtKBa2a1rHKhK21N6rn9BfB3YBgRdy57CZuQArcnsFsEMytWj5lZCVRuekG6BxhGulH2S5q9kcTGwJ9Jj/nuGsH0yhRiZlYelQzdQcDbRLy+/Ck2IgVuP1LgTqtMEWZm5ZJ794LEQFLgrkWaUng01wLMzAqU6ypjWR/uJGBtYLQD18zqTW6hK7EqcBupD3e3CB7O673NzMoil9CV6A78GtgW2DuCqXm8r5lZ2eQ10r0A2As4OoJbc3pPM7PSqfjaCxLHAN8EfhLBJZV+PzOzMqto94LE3sBNwP8CY72fmZnVu4qFrsQnSJ0K00m9uIsq8kZmZlWkIqErMRiYStoJYkQEr3T5m5iZVaFKzeluAbwPfM6Ba2bWqJLTC30ieK8iP9zMrEoVv0eamVkdKdtuwGZmNc2ha2aWI4eumVmOHLpmZjly6JqZ5ciha2aWI4eumVmOHLpmZjly6JqZ5ciha2aWI4eumVmOHLpmZjly6JqZ5ciha2aWI4eumVmOHLpmZjmq69CV9BFJ+xddh1lz/mzWrroOXWA3YHjRRZi1wJ/NGlW32/VI2hm4BVhA2rV474h4ptCizPBns9bVbegCSLoLOCEiZhRdi1lT/mzWrnqfXtgCeLLoIsxa0K7PpqStJW3Y7GuDJV0p6caKVWcrrW5DV9LawJsRsTi7aTFJ0rckXV90bVbfOvjZ3BYY3PQLEfF0RByUS7HWYXUbusDGwLzsz9sAv4+InwBLiivJDGjHZ1PSVpIuB74GHCfpckkD8i/VOqpH0QUUaDbQX9IM4Drg5uzr9TvJbWWxws9mRMwEviFpPPBsREzOu0hbOXU70o2IdyJih4j4KLApMEdSf+ClgkuzOtfZz6aktbNR8MclnVzJWq3j6rp7wcwsb3U70jUzK4JD18wsRw5dM7McOXTNzHLk0DUzy5FD18wsRw5dM7McOXTNzHLk0DUzy9H/Ab0ekLIuehCFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0.5,4,50)\n",
    "\n",
    "plt.plot(x,np.log(x),'b')\n",
    "x0 = 1.0\n",
    "x1 = 3.0\n",
    "plt.plot(x,np.log(x0) + x/x0 -1,'r--')\n",
    "plt.plot(x0,np.log(x0),'ob')\n",
    "plt.plot(x1,np.log(x1),'ob')\n",
    "plt.plot(x1,0,'ok')\n",
    "plt.plot(x,0*x,'k--')\n",
    "plt.plot(x1,np.log(x0) + x1/x0 -1,'or')\n",
    "plt.vlines(x1+0.1,ymin=0,ymax=np.log(x0) + x1/x0 -1,colors='r',linestyles='dashed')\n",
    "plt.vlines(x1,ymin=0,ymax=np.log(x1),colors='b',linestyles='dashed')\n",
    "\n",
    "\n",
    "xt = [x0,x1]\n",
    "plt.box(False)\n",
    "locs, labels = plt.xticks()\n",
    "plt.xticks(xt, ('$t_n$','$t_{n+1}$'))\n",
    "plt.tick_params(axis='both',width=0.0,labelleft=False)\n",
    "\n",
    "fs = 18\n",
    "plt.annotate('A',xy=(x0-0.2,np.log(x0)+0.2),fontsize=fs)\n",
    "plt.annotate('B',xy=(x1-0.2,0+0.2),fontsize=fs)\n",
    "plt.annotate('C',xy=(x1-0.1,np.log(x1)+0.2),fontsize=fs)\n",
    "plt.annotate('D',xy=(x1-0.2,np.log(x0) + x1/x0 - 1 +0.2),fontsize=fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exemplo:** A solução exata do PVI \n",
    "\n",
    "$$\\begin{cases}\n",
    "y'(t) = -y(t)\\\\\n",
    "y(0) = 1\n",
    "\\end{cases}$$\n",
    "\n",
    "é $y(t) = e^{-t}$. O método de Euler é dado por \n",
    "\n",
    "$$y_{n+1} = y_n - h y_n = (1-h)y_n, \\ \\ n \\geq 0,$$\n",
    "\n",
    "com $Y_0 = 1$ e $t_n = nh$.\n",
    "\n",
    "Para $h=0.1$, temos, por exemplo\n",
    "\n",
    "$$y_1 = (1-h)y_0 = 0.9(1) = 0.9$$\n",
    "$$y_2 = (1-h)y_1 = 0.9(0.9) = 0.81$$\n",
    "$$y_2 = (1-h)y_1 = 0.9(0.9) = 0.81,$$\n",
    "\n",
    "cujos erros são \n",
    "\n",
    "$$y_{h,1} - y_1 = e^{-0.1} - y_1 = 0.004837$$\n",
    "$$y_{h,2} - y_2 = e^{-0.2} - y_2 = 0.008731$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementação computacional \n",
    "\n",
    "O seguinte código implementa o método de Euler explícito. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "\n",
    "def euler_expl(t0,tf,y0,h,fun):\n",
    "    \"\"\"\n",
    "    Resolve o PVI y' = f(t,y), t0 <= t <= tf, y(t0) = y0\n",
    "    com passo h usando o metodo de Euler explicito. \n",
    "    \n",
    "    Entrada: \n",
    "        t0  - tempo inicial\n",
    "        tf  - tempo final \n",
    "        y0  - condicao inicial \n",
    "        h   - passo \n",
    "        fun - funcao f(t,y) (anonima)\n",
    "        \n",
    "    Saida:\n",
    "        t   - nos da malha numerica \n",
    "        y   - solucao aproximada\n",
    "    \"\"\"\n",
    "    \n",
    "    n = round((tf - t0)/h) + 1\n",
    "    t = linspace(t0,t0+(n-1)*h,n)\n",
    "    y = linspace(t0,t0+(n-1)*h,n)\n",
    "    y = zeros((n,))\n",
    "    \n",
    "    y[0] = y0\n",
    "\n",
    "    for i in range(1,n):\n",
    "        y[i] = y[i-1] + h*f(t[i-1],y[i-1])\n",
    "\n",
    "    return (t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exemplo:** Resolva \n",
    "\n",
    "$$\\begin{cases}\n",
    "y'(t) = \\frac{y(t) + t^2 - 2}{t+1}\\\\\n",
    "y(0) = 2 \\\\\n",
    "0 \\le t \\le 6 \\\\ \n",
    "h = 0.1\n",
    "\\end{cases}$$\n",
    "\n",
    "Defina $y_h(t)$ como a solução numérica, calcule o erro relativo e plote o gráfico de $y_h(t)$ juntamente com o da solução exata $y(t) = t^2 + 2t + 2 - 2(t+1)\\ln(t+1)$\n",
    "\n",
    "**Solução:** O processo iterativo de Euler será dado por \n",
    "\n",
    "$$y_{n+1} = y_n + \\dfrac{h(y_n + t_n^2 - 2)}{t_n + 1}, \\ \\ n \\geq 0, \\ \\ y_0 = 2, t_n = nh.$$\n",
    "\n",
    "Entretanto, vamos usar o nosso programa `euler_expl`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd023af8b20>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqGElEQVR4nO3de3zO5ePH8de12YbJYeSc0IFEOQx9E4b4IjJETm3OqUTpW0rpyxeRtM4ncirTYUKi5DTHNKdIOWSMGjmfdzC7d/3+uJdfirDTZ/e99/Px2MN27959v+9v9v5ervv6XJex1iIiIp7Hx+kAIiKSMSpwEREPpQIXEfFQKnAREQ+lAhcR8VD5cvLJSpQoYStWrJiTTyki4vE2btx41Fp7/V9vz9ECr1ixIhs2bMjJpxQR8XjGmH2Xul1TKCIiHkoFLiLioVTgIiIeKkfnwC/l/PnzxMfHk5yc7HSUbJc/f37Kly+Pn5+f01FExAs4XuDx8fFcd911VKxYEWOM03GyjbWWY8eOER8fT6VKlZyOIyJewPEplOTkZIoXL+7V5Q1gjKF48eJ54l8aIpIzHC9wwOvL+w955XWKSM7IFQUuIuKtkpKSGDx4MIcOHcryx1aBi4hkoxdffJE333yTrVu3Zvljq8BFRLLJ2rVriYiIoH///tx7771Z/vgq8GuUlJRE48aNcblcgHsVzWeffQZASkoKjRo1IjU11cmIIpILJCUl0atXL8qVK8crr7ySLc+hAr9GU6ZMoUOHDvj6+gKwdOlSNm3aBIC/vz/NmjW7UOgikneNGDGCnTt38uGHH1K4cOFseQ4VOLB161YaNGhw4etNmzbRtGnTS943MjKSdu3aAbB69WqGDBnCrFmzqFmzJnFxcYSGhhIZGZkjuUUkd4qJiWHChAn07duXFi1aZNvzOH4hz1+FhPz9ts6d4dFHITERWrf++/d79nR/HD0KDzxw8feWL7/yc95+++3s3r0bl8uFr68vTz31FK+++urf7peSksKePXv4Y0vce+65h7p16zJhwgSqV68OgMvlYv369Vd+UhHxSsnJyfTq1YuyZcsyYcKEbH0ujcABHx8fbr/9dn7++We++OILKlSoQO3atQHo27fvhfsdPXqUokWLXvSzO3fupEqVKhe+9vX1xd/fnzNnzuRIdhHJXUaOHMn27duZNGkSRYoUweUCa7PnuXLdCPyfRswFC/7z90uUuLoR96XcddddrFmzhnfffZeFCxcCkJiYSGxsLM8//zzbtm1jypQpF11JeezYMYoUKfK3vU3OnTtH/vz5MxZERDzW+vXrGT9+PL1796Zly5YAjB4NGzdCVBQEBGTt82kEnu6uu+7ihRdeoH379pQrVw5wz4W3atWKMWPGEBgYSLFixXC5XBdKPC4ujrJly170OMeOHeP666/XhlUieUxSUhLh4eGULVv2whTsxo3uAi9cOOvLG1TgF1StWpWAgACGDh164bb169fTtm1bgAurTlq0aMHq1asv/MzRo0epXr063333HQDR0dG0vtREvYh4teHDh7N9+3YmT55M0aJFSU6G8HAoWRLeeit7nlMFnu6NN95g7NixBAYGXrht27ZtVKlShaNHj3L99e7j6AYOHMj06dMBKFSoEOvWreOnn37i7rvvBmDmzJn0798/51+AiDhm1apVREREMGDAgAurTn7+Gfbvhw8/hGLFsud583yB7969m6pVq17458+fTZo0CV9fX0qUKHHh3eRatWrRpEmTCxfy/FlKSgqhoaEXvakpIt7t7Nmz9OzZk4oVK150wU6dOhAXB61aZd9z57o3MXPaTTfdxI4dO67pZ3r37n3J2/39/QkLC8uKWCLiIYYOHUpcXBzLly+nUKFCJCTA55+7p0/+smgty+X5EbiISEYtWbKEd999lyeeeIJGjRoB8Oyz0Ls3/PBD9j+/ClxEJANOnTpF7969qVKlCmPGjAFgyRJ4+20YPNg9hZLd8vwUiohIRgwePJgDBw6wZs0aChQowIkT7ivCq1SBl17KmQwqcBGRazR79mymT5/O8OHDqV+/PgCPPw6HDsHcue6LDnOCClxE5BocPHiQ/v37U6dOHYYPH37h9l694F//guDgnMuiAhcRuUrWWvr27UtCQgIzZszAz88Plwt8faFZM/dHTtKbmCIiV2nSpEksWLCA8ePHU7VqVdLSoE0b9+XyTlCBi4hchdjYWJ588kmaN2/OY489BsA778DChe7L5Z2gAr+Cyx2hpuPTRPKO1NRUwsLC8Pf3Z+rUqfj4+LB9OzzzDNx3H/Tr50wuFfgVXO4INR2fJpJ3jBs3jrVr1/Lee+9Rrlw5zp2Dbt2gUCH3XifGOJPrigVujLnBGBNtjNlujPnZGDM4/fYgY8xiY8yu9D+zabuW7PdPR6r90xFqOj5NxPutW7eOESNG0LVrV7p06QLA99/Dtm0wZQqULu1cNmOvcFSEMaYMUMZau8kYcx2wEQgFegLHrbXjjDHPAsWstUMv/0gQHBxsN2zYcNFt27dv57bbbgPgiSeeYPPmzRl7JZdRs2ZNXn/99X+8T1paGmXLlmX//v34+vrSpEkTXn31VapXr06FChU4ePDghfu2bNnywhFqLpeL0qVLc+TIkavO8+fXKyK529mzZ6lVqxYpKSls2bLlohO54uOhfPmcyWGM2Wit/dsCxSuOwK21v1trN6V/fgbYDpQD2gHT0+82HXepe6TLHal2pSPUdHyaiHd74okn2LNnDzNmzKBo0aIcOwZffun+Xk6V9z+5pnXgxpiKQC0gBihlrf0d3CVvjMn0+7BXGilnp0sdqVagQIErHqGm49NEvNMXX3zB5MmTGTZsGA0bNsRa95uV8+dDbCxUqOB0wmt4E9MYUwj4AnjCWnv6Gn6uvzFmgzFmw7VMNeS0Sx2pdqUj1HR8moh32r9/P/369SM4OJgRI0YA7vnuOXPc+5zkhvKGqyxwY4wf7vKOtNbOTr/5UPr8+B/z5Icv9bPW2onW2mBrbfAfp9rkRpc6Ug3++Qg1HZ8m4n3S0tIIDw/n3LlzREZG4ufnxy+/wKBB0LQpDBnidML/dzWrUAwwGdhurY3407fmAX8cYRMOfJn18XLOpY5Ug38+Qk3Hp4l4n4iICJYuXcrrr7/OrbfeSmoqdO3qPpR4+nTwyUWLr68mSgPgIaCpMWZz+kdrYBzQ3BizC2ie/rXH+acj1eDyR6jp+DQR77NhwwaGDRtG+/bt6du3LwD58sGTT8K0abnjjcs/u+Iywqx0pWWEeUFee70inuLMmTPUrl2b5ORktmzZQlBQECkp4O/vdLJMLCMUEckLHn/8cfbs2UNkZCRBQUEcPAi33gpRUU4nuzwVuIjkeZGRkUyfPp0XXniBRo0akZYGYWHuAxqqVXM63eXliv3ArbUYpzYTyEE5OV0lIldn9+7dPPLIIzRo0ODCAQ2vvgqLF8P778Pttzsc8B84PgLPnz8/x44d8/pys9Zy7NgxXfQjkoukpKTQtWtXfH19iYyMJF++fKxfD8OGQceOkNsXmTk+Ai9fvjzx8fHXtJ+Ip8qfPz/lc9vb2CJ52PDhw1m/fj1RUVHceOONACxfDuXKwaRJzu0yeLUcX4UiIuKEhQsX0qpVK/r3788HH3xw0fdOn4bChR0KdglahSIiku7AgQM89NBD1KhR48IeTDNmwNq17u/npvL+JypwEclTXC4X3bp1IzExkc8//5wCBQqweTP06QNjxzqd7to4PgcuIpKT/ve//7FixQqmT59O1apVOXMGOneGEiVg8mSn010bFbiI5BnLli1j1KhRhIeHExYWhrXw8MOwezdER0Mu3m/vklTgIpInHDp0iO7du1OlShXeeecdABYsgE8+gdGjoVEjhwNmgApcRLyey+XioYce4uTJkyxatOjCrqOtWrl3GOzRw+GAGaQCFxGvN3r0aBYvXsykSZOoUaMGp07B2bPu9d5hYU6nyzitQhERr7Z48WJGjhxJWFgYffr0wVro3Rvq1oWEBKfTZY4KXES81v79++nevTvVqlXj3XffxRjDG2/A7Nnuk3X+cn6Lx1GBi4hXOn/+PF26dCExMZFZs2YRGBjI2rXw9NMQGgpPPeV0wszTHLiIeKXnn3+e1atXM3PmzPTzbN3rvStUgKlTc/8+J1dDI3AR8Trz5s3jlVde4ZFHHqFr164A+PlBw4buAxqKFnU2X1bRZlYi4lViY2MJDg7m5ptvZs2aNQQEBOByga+v08kyTptZiYjXS0xMpGPHjvj6+jJr1iwCAgKYPx/q1YMDB5xOl/VU4CLiFay1DBgwgK1btxIZGUnFihWJjf3/i3SKFXM2X3bQm5gi4hXef/99Pv74Y0aOHEnLli1JTHSfquPrC198AQUKOJ0w66nARcTjff/99wwePJjWrVvzwgsvXNikautW+PprqFjR6YTZQ1MoIuLRDh8+TKdOnShfvjwff/wxPj4+nDoFmzfDiBHQsqXTCbOPRuAi4rFSU1Pp0qULR48e5bvvviMoKAhwLxOMiQFvP0NcI3AR8VhDhw4lOjqa9957j1q1arF/PzzyiHujqoIFwcfLG87LX56IeKuZM2cSERHBwIED6dmzJ8nJ0KGD+2zL+Hin0+UMTaGIiMf54Ycf6Nu3Lw0bNiQiIgJr4bHHYN0690ZVVas6nTBnaAQuIh7l6NGjtG/fnuLFixMVFYWfnx/vvQdTpsALL0D79k4nzDkagYuIx0hNTeXBBx/k4MGDrFq1ilKlSpGQ4D4S7b77YORIpxPmLBW4iHiMoUOHsmzZMqZMmULdunUB957e333nXnni7W9a/lUee7ki4qk++ugjIiIieOyxx+jVqxeJifDBB2Ct+0Idb9lh8FqowEUk14uJiaF///6EhITw2muvYS306uVeMrhxo9PpnKMCF5Fcbf/+/bRv356yZcteeNNyzBj4/HMYPx6C/7bJat6hOXARybWSkpJo3749Z86cYdGiRZQoUYI5c2D4cPcug95wLFpmqMBFJFey1tKvXz/Wr1/P3LlzqV69OqdPu0+Ur1cPJk3yjmPRMkMFLiK50iuvvEJkZCSjRo2iXbt2ABQuDHPnwi23eP8+J1dDc+AikuvMmzePZ599lk6dOvH888+TnAyLF7u/17gxlC3rbL7cQgUuIrnKli1b6NatG3Xq1GHatGmAoW9f97awO3c6nS53uWKBG2OmGGMOG2N++tNtI4wx+40xm9M/WmdvTBHJCw4ePEjbtm0pWrQo8+bNo2DBgowZA5GR8L//QZUqTifMXa5mBD4NuNSW6K9Za2umf3ydtbFEJK9JSkqiXbt2HDt2jK+++ooyZcoQFfX/K06GDXM6Ye5zxTcxrbUrjTEVcyCLiORR1lp69erF+vXrmT17NrVq1WLfPggLg7vv1oqTy8nMHPhAY8yP6VMslz3v2RjT3xizwRiz4ciRI5l4OhHxViNHjuSzzz5j7NixhIaGAnDjjfDOOzBnjlacXI6x1l75Tu4R+HxrbfX0r0sBRwELjALKWGt7X+lxgoOD7YYNGzIVWES8S2RkJD169CA8PJypU6dy+rRh3z644w6nk+UexpiN1tq/XXOaoRG4tfaQtdZlrU0DJgH1MhtQRPKelStX0rt3b0JCQpg4cSKpqYYHHoCQEDh1yul0uV+GLuQxxpSx1v6e/mV74Kd/ur+IyF/t3LmT0NBQKleuzOzZs/Hz86dfP1iyxH04Q5EiTifM/a5Y4MaYT4AQoIQxJh74LxBijKmJewplL/Bw9kUUEW9z5MgRWrduTb58+ViwYAHFihXjpZdg8mT3qpNevZxO6BmuZhVK10vcPDkbsohIHpCUlMT999/PgQMHWL58OZUrV2bJEnj+eejePe+dqpMZ2gtFRHJMWloa4eHhxMTEEBUVRf369QH35fGvvAKPP67lgtdCl9KLSI55+umniYqKYvz48XTs2JEdO+DQIfDzg//8BwICnE7oWVTgIpIjXnvtNSIiInj88cd56qmniI+H5s3hgQfcx6LJtVOBi0i2+/zzzxkyZAgdOnTgtdde49QpQ6tW7qWCb7+taZOM0hy4iGSrFStW8NBDD9GgQQNmzJhBaqovoaHunQW/+QbuvNPphJ5LBS4i2ebnn3++sNZ73rx5FChQgKefhhUrYMYMaNbM6YSeTQUuItkiPj6eVq1aUaBAARYuXEhQUBAAQ4dCzZruJYOSOZoDF5Esd/z4cf79739z8uRJvv76a2688UbmzIGUFChRQuWdVVTgIpKlEhISaNOmDbGxscybN4+aNWsycSJ06ABvveV0Ou+iKRQRyTLnz5+nU6dOFy7UCQkJ4Ysv4JFHoHVrGDTI6YTeRQUuIlkiLS2N3r1788033zBx4kQ6dOhAdDR06wb160NUlPuCHck6mkIRkUyz1vL0008zY8YMRo8eTb9+/UhJgZ494eabYf58KFjQ6ZTeRyNwEcm0cePGERERwaBBgxiWfnilvz/MmwfFi0P6AhTJYhqBi0imvPPOOwwbNowePXrw2muvsXevufBm5Z13QvnyzubzZhqBi0iGzZgxg4EDB9KuXTumTp3KwYM+3HsvnDwJnTpB6dJOJ/RuKnARyZC5c+fSs2dPmjVrxqeffsqpU/lo0QIOH4alS1XeOUEFLiLXbOnSpTz44IMEBwczd+5cUlLy06oVxMbCwoVQT6fk5gjNgYvINVm7di3t2rWjSpUqfP311xQqVIgVK2DLFvdSwZAQpxPmHRqBi8hV27BhAy1btqRs2bIsWrTowv4mbdvC7t16wzKnaQQuIldly5YttGjRgqCgIJYuXUqxYqXp2NE9ZQIqbyeowEXkirZt20bz5s0JDAxk2bJllC59Aw8+CLNnw/79TqfLuzSFIiL/aNeuXTRr1gxfX1+WLVvGDTdUont3+PJL92k6ffo4nTDvUoGLyGXFxcXRrFkzUlNTWbFiBZUr30JYGHz+OUyYAI895nTCvE0FLiKXtHfvXpo0acLZs2eJjo6mWrVqpKVBYCCMHQtPPeV0QlGBi8jf7N27l5CQEE6fPs2SJUuoXv1Ofv8dypSBDz7QIcS5hd7EFJGL7Nu3jyZNmnDq1CmWLFnCnXfWJjwc7rrLfYq8yjv3UIGLyAX79u0jJCSEkydPXlTekZHw8MNQpIjTCeXPNIUiIgD8+uuvNGnShJMnT7J48WLuvLMOYWEwcyaMHg3pu8RKLqIRuIiwZ88eGjVqxPHjx1m8eDHBwcGMGuUu77Fj4fnnnU4ol6IRuEge98svv9C0aVOSkpJYtmwZtWvXBuDJJ92n6Tz0kMMB5bI0AhfJw7Zt20bjxo1JSUkhOjqa226rzYsvQlISFC2q8s7tVOAiedSPP/5ISPrWgcuXL+emm+6gbVv3fHd0tLPZ5OqowEXyoI0bN9KkSRP8/f1ZsWIF5cpVo2VLd3FPmwatWzudUK6GClwkj1m9ejVNmzZN38d7BUFBt9K0KcTEwCefQFiY0wnlaqnARfKQb7/9lhYtWlC6dGlWr17NTTfdxNGjcOSIe3Oqzp2dTijXQqtQRPKI2bNn06VLF26//Xa+/fZbjCmJtVC1KuzaBQEBTieUa6URuEgeMH36dDp16kRwcDDR0dEcOVKSO++EcePc31d5eyYVuIiXe+utt+jZsydNmzZl0aJFbNtWlIYNwVr3UWjiuVTgIl7KWsuLL77IoEGDCA0N5auvvmLlykLcey8ULw7ffQfVqzudUjLjigVujJlijDlsjPnpT7cFGWMWG2N2pf9ZLHtjisi1cLlcDBgwgFGjRtG3b1+ioqI4eTI/HTvCbbfBmjVQqZLTKSWzrmYEPg1o+ZfbngWWWmtvAZamfy0iuUBycjKdO3dm4sSJDBs2jIkTJ5IvXz5Kl3afYRkdDSVLOp1SssIVC9xauxI4/peb2wHT0z+fDoRmbSwRyYjTp0/TqlUrZs+ezeuvv86oUWN47jnDnDnu77dqBYULO5tRsk5G58BLWWt/B0j/87L/f26M6W+M2WCM2XDkyJEMPp2IXMmBAwdo3Lgxq1evJjIykgEDBtOjB7z8Mqxc6XQ6yQ7Zvg7cWjsRmAgQHBxss/v5RPKibdu20apVK44fP878+fOpX//ftGrlni4ZOxaGDnU6oWSHjBb4IWNMGWvt78aYMsDhrAwlIldvxYoVhIaGkj9/flauXMnNN9fi7rth5074+GPo0cPphJJdMjqFMg8IT/88HPgya+KIyLX47LPPaNGiBWXKlOH777+nVq1aFCoEbdrAN9+ovL3d1Swj/ARYC1QxxsQbY/oA44DmxphdQPP0r0Ukh1hrmTBhAl26dKF+/fqsXr2abdtu5Mcf3YcOjx0LzZo5nVKy2xWnUKy1XS/zLf31EHHA+fPnGThwIBMnTqRTp0589NFHTJ6cn0GD4L77YN48pxNKTtGVmCIe5OTJk7Ru3ZqJEyfy3HPPERn5KUOH5mfgQHd5z5zpdELJSdqNUMRDxMXF0aZNG3bt2sWUKVPo3LkXHTrA/PnwxBMwYQL4+jqdUnKSClzEA6xdu5Z27dqRmprKokWLCAkJISUFzp2Dd96BRx91OqE4QVMoIrncRx99RJMmTShcuDBr167F3z+Ew4fB3x8WLlR552UqcJFcyuVy8fTTTxMeHs7dd99NTEwMq1dXISQEnnnGfR8f/QbnafrPL5ILnTp1irZt2zJhwgQeffRRFiz4llGjitO3L4SEwGuvOZ1QcgPNgYvkMrt27eL+++8nNjaW9957j86dB9CuHSxe7H6z8pVXIJ9+cwUVuEiu8s0339CtWzd8fX1ZvHgxISEhHDsGv/0GkydD795OJ5TcRFMoIrlAWloao0eP5r777uPGG29k3bp1JCS4V5oULw4//qjylr9TgYs47PTp03To0IHhw4fTrVs3Vqz4jjffrEybNu4lggB+fs5mlNxJUygiDtq+fTvt27cnNjaW119/nU6dBnH//YaVK93z3QMHOp1QcjMVuIhDPv/8c/r06UPBggVZunQpfn6NqV0bzpyBGTOge3enE0pupykUkRyWkpLCoEGDePDBB6lRowYbN26kcePGFCkC5cpBTIzKW66OClwkB/366680atSIt956iyeffJK5c1cwZ055rIXbb4cNG6B6dadTiqfQFIpIDvnmm2/o0aMH58+fZ9asWZQr15H69WH/fvfe3dWquffyFrlaGoGLZLPz588zdOhQWrduTfny5Vm/fiO7d3ekYUP391etcpe3yLXSCFwkG8XFxdG1a1diYmJ4+OGHiYiIoE+fgnz6KTzwAEyaBEWLOp1SPJUKXCSbREVF0a9fP8C94qRTp04AdOni3s+kf39NmUjmaApFJIslJiYyYMAAOnfuTNWqVfn++x9Yv74TERHu77drBw8/rPKWzNMIXCQLbdq0ie7du7Njxw6eeeYZuncfTdeufmzeDI8/7nQ68TYagYtkAZfLxbhx46hfvz5nzpxhyZKlVKjwMvXr+xEfD19+CW++6XRK8TYagYtk0r59+wgLC2PlypV06tSJ999/n/j4IFq0gBYtYOpUKF3a6ZTijVTgIhlkreWjjz5i8ODBuFwupk2bRo0aYQQFGYKCYM0aqF9fc92SfTSFIpIBBw8epF27dvTs2ZMaNWoQHb2FBQvCqVPHEBPjvs9dd6m8JXtpBC5yjT777DMeffRREhISePXVV6lYcTD33efLiRMwdizUqeN0QskrNAIXuUpHjhzhwQcfpEuXLtx8881s3ryZHTuG0LGjL2XKuPcxefZZHXcmOUcFLnIF1loiIyO57bbbmDNnDqNHj2H16jVUrVqVO+6A//4X1q2DO+5wOqnkNSpwkX/w22+/0aZNG3r06MEtt9zCwoU/EBMzjKgo9zB74EAYMQL8/Z3NKXmT/rEncglpaWl88MEHDB06FJfLRUTE6/j5DaR9e1/On4c2bZxOKKIRuMjfbN26lYYNG/Loo49Sv359oqJ+4rPPBvP4477Ures+YLh/f6dTiqjARS5ISEhg6NCh1K5dm19++YVp06axaNEizp2rRFyc+5izxYvh5pudTiripikUEWD+/PkMHDiQffv20bt3Hxo3fpnExOIYA6GhcO+9cN11TqcUuZhG4JKnxcXFERoaStu2bQkMDGTGjJXEx39IeHhxpkwBl8t9MY7KW3IjjcAlT0pMTOTll1/m5ZdfJl++fIwcOY7Tp5+kVy9/ChSAN96ARx8FX1+nk4pcngpc8hRrLXPmzGHIkCHs27ePrl278sorr3DoUDnq1YOwMPfVlKVKOZ1U5MpU4JJnbNmyhSFDhrBs2TJq1KjB228v59SpxpQrB+XKQWwsVKzodEqRq6c5cPF6v//+O3379qVWrVps2bKFkSPfomrVTQwc2Jh334XTp933U3mLp9EIXLxWYmIiERERjBs3jpSUFB57bAi+vs/z0kvF8PGBF1+EZ56BwECnk4pkTKYK3BizFzgDuIBUa21wVoQSyYzU1FSmTp3KiBEjOHDgAB06dGD8+PEEBt7ErbdC584wZgzccIPTSUUyJytG4E2stUez4HFEMuWPNyiHDRvGzp07ueuufxEe/il79jSkcmX3csA9e6BECaeTimQNzYGLV1i2bBn/+te/6NixIz4+PvznP3M4cWINY8c2ZM8eOH7cfT+Vt3iTzBa4BRYZYzYaYy65O4Qxpr8xZoMxZsORI0cy+XQiF1u1ahVNmjShWbNmxMfH89JLH5Iv349MmBCKr68hKgpiYqB4caeTimS9zBZ4A2ttbaAV8JgxptFf72CtnWitDbbWBl9//fWZfDoRt7Vr19K8eXMaNWrEjh07GDnyDWJjY3nyyT4UL56PyEj3plMPPKBjzcR7ZarArbUH0v88DMwB6mVFKJHLWbVqFS1btuTuu+9my5YthIe/SsmSu5k6dRC+vvnJnx+io6FbN11FKd4vwwVujAk0xlz3x+dAC+CnrAom8gdrLYsWLaJx48Y0atSITZs20aXLyxQvHsf06UNITi7IiBEaaUvek5lVKKWAOcb9W5MPmGmtXZglqUQAl8vFvHnzeOmll9iwYQPlypXjjTfeoHz5vnTsWJDq1eGTT6BTJ422JW/KcIFba/cAd2ZhFhEAkpKS+Oijj3j11VfZtWsXFSpUok2biTRoEMagQQGkpcGCBdCyJfhoHZXkYfrrL7nG0aNHGTlyJBUqVGDAgAEEBBTh3ns/49ChX5g/vx87dgQA7tJu3VrlLaJL6cVxmzdv5q233iIyMpJz587Rpk0bSpf+Dx9+2Ihduww9esBTT8FttzmdVCR3UYGLI1JTU5k7dy5vvvkmq1atomDBggQH92LYsEG0bn0b69a5N5fq3x+0+lTk0lTgkqPi4+OZPHkyH374IfHx8ZQtW4m6dV9l69ZerFlTjB9/dE+P1Kvn/hCRy1OBS7ZzuVwsXLiQDz74gAULFmCtpUWLFhQu/DbbtrXh1ClfwsLg4Yehdm2n04p4DhW4ZJvdu3czbdo0pk+fzm+//UaxYqWoV+9ZIiP7UrlyJcaOhaJFoXt3KFzY6bQinkcFLlnqzJkzREVFMW3aNFatWoWPjw+VK7fg+usjOHKkHdu3++Hv777vc885m1XE02khlmRaSkoK8+bNo2vXrpQqVYo+ffpw+PBhevYcS1rar+ze/Q133vkAM2f68fvvUL6804lFvING4JIhLpeLlStX8sknnzBr1ixOnDhBoULFKV48nNDQMN588y4SEw3Vq8ODD6q0RbKDClyu2vnz54mOjuaLL75g7ty5HD58mICAQIKC2uPj05WzZ5tTsqQf1aq59yUJDHSv3xaR7KECl3+UkJDA4sWL+fLLL/nyyy85ceIEAQGBhIa2oWPHjkyadB9xcQUZOtS9dWutWtpUSiSnqMDlb+Lj45k/fz5fffUVS5cu5dy5cwQEFCEgoC3wAOfOteDttwtQogS0aOFeQaLSFsl5KnAhOTmZVatW8e2337Jo0SK2bt0KwE033USDBo+wbFlbUlMbcs89frRrB23b/v/RZEWKOBhcJI9TgedBLpeLH374gejoaJYtW8aKFStISkrCx8efgICGwHhGj27DsGFVOXjQ8P330KyZ1mqL5DYq8DwgNTWVLVu2sGrVKqKjo1m+fAWnT58C4Oabq5CU1BdoScGCjQkJCaRpU2jf3j0tUqaM+3MRyX1U4F7oxIkTxMTEsGbNGtasWcPatTEkJycCkC/fzaSmdiI4uAlffhlC2bJleecdqFMHgoMhn/5GiHgM/bp6uISEBDZv3sz69euJiVnHd9+t59dfYwHw8fGhZs2aFCnSh+TkBhQv3oDGjctzzz3QpAmULet+jMcec/AFiEiGqcA9hLWW+Ph4tm7dml7Ym9m2bQu7du3CWguAMeWxti7Qm8KF6/Hbb/UpXLgQmzZBsWLu7Vm1WkTEe6jAc5mUlBTi4uLYtWsXO3bsYNu2bWzc+DO//LKd5OQzf7pnJe67rybdu3dn796abN9el7vuKkOtWu612FWrgp+f+57a4U/EO6nAc5i1liNHjrBv3z727t3L3r172bFjLz/8sJtff93F8eN7sTbtwv1Lly5N0aLVSE4Ox9+/GhUr3k6dOndSq1YR+vZ1j6xFJG9SgWcRl8vFsWPHOHToEL//fpi4uEOcOHGY48cPEBsbz48/7ufYsf2cObMflyv5op8tVKgYZ89WBurh59edkiVvoVKlWxgz5lYaNQri5ElISHDPWWsKRET+kGcL3OVyce7cOc6dO0dycjInT57jzJlkTpxI4MSJBE6dSsSYBPLnT+D48dOsW3eaU6dOc/r0Kc6cOc2ZM8fx8zuOy3Wco0ePcfLkScD+7XkCAgK4/vpyxMeXw5i6BAaGEhRUgXLlKtK/f0Xat78Ra4uwcydUquQ+PuyvJV20qPtDROTPPKLAR40axfjxn5CScvHtBQtaSpWyWGv57TdLaqoFLNa6ABcBAS6uuy4tfXR8HmvdH5AKpP39ia7IHygCFMaYIG66KYi6dW8iMDCIpUuDKFasJCVKlKR06VLccEMp/v3vktxzTzFcLsOpU+7pjsudpF6/fgbiiEie5hEFXqZMGcqUqU5S0sW3Fy0K1aoZjDGkpRlSUgxg8PHxxcfHl5IlfbnjDl98fHxYu9aPtDQ/8uVzf/j65qNSpfzUqxdA/vz5iY4OwN8/gEKFArnuukAKFw6kSpVAatcO5LrrChMfX5igoAAKFXLvspc//9VPZ+TLB8WLZ/X/KiKS15k/lqDlhODgYLthw4Ycez4REW9gjNlorQ3+6+06kUdExEOpwEVEPJQKXETEQ6nARUQ8lApcRMRDqcBFRDyUClxExEOpwEVEPFSOXshjjDkC7Mvgj5cAjmZhHCfpteQ+3vI6QK8lt8rMa7nRWnv9X2/M0QLPDGPMhktdieSJ9FpyH295HaDXkltlx2vRFIqIiIdSgYuIeChPKvCJTgfIQnotuY+3vA7Qa8mtsvy1eMwcuIiIXMyTRuAiIvInKnAREQ/lEQVujGlpjNlpjIk1xjzrdJ6MMsZMMcYcNsb85HSWzDDG3GCMiTbGbDfG/GyMGex0powyxuQ3xqwzxmxJfy0jnc6UGcYYX2PMD8aY+U5nyQxjzF5jzFZjzGZjjEefAmOMKWqMmWWM2ZH+O/OvLHvs3D4HbozxBX4BmgPxwHqgq7V2m6PBMsAY0wg4C3xkra3udJ6MMsaUAcpYazcZY64DNgKhHvrfxACB1tqzxhg/YDUw2Fr7vcPRMsQYMwQIBgpba9s4nSejjDF7gWBrrcdfxGOMmQ6sstZ+aIzxBwpaa09mxWN7wgi8HhBrrd1jrU0BPgXaOZwpQ6y1K4HjTufILGvt79baTemfnwG2A+WcTZUx1u1s+pd+6R+5e1RzGcaY8sB9wIdOZxE3Y0xhoBEwGcBam5JV5Q2eUeDlgN/+9HU8HloW3sgYUxGoBcQ4HCXD0qcdNgOHgcXWWk99La8DzwBpDufIChZYZIzZaIzp73SYTKgMHAGmpk9tfWiMCcyqB/eEAr/U2e8eOULyNsaYQsAXwBPW2tNO58koa63LWlsTKA/UM8Z43PSWMaYNcNhau9HpLFmkgbW2NtAKeCx9+tET5QNqA+9Za2sBCUCWvY/nCQUeD9zwp6/LAwccyiLp0ueLvwAirbWznc6TFdL/abscaOlskgxpANyfPnf8KdDUGDPD2UgZZ609kP7nYWAO7qlUTxQPxP/pX3WzcBd6lvCEAl8P3GKMqZT+BkAXYJ7DmfK09Df+JgPbrbURTufJDGPM9caYoumfFwDuBXY4GioDrLXPWWvLW2sr4v4dWWat7eFwrAwxxgSmvzlO+nRDC8AjV25Zaw8CvxljqqTf1AzIsjf782XVA2UXa22qMWYg8C3gC0yx1v7scKwMMcZ8AoQAJYwx8cB/rbWTnU2VIQ2Ah4Ct6XPHAMOstV87FynDygDT01c7+QCfW2s9egmeFygFzHGPE8gHzLTWLnQ2UqY8DkSmD0D3AL2y6oFz/TJCERG5NE+YQhERkUtQgYuIeCgVuIiIh1KBi4h4KBW4iIiHUoGLiHgoFbiIiIf6PwzOOP22cSekAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# define funcao\n",
    "f = lambda t,y: (y + t**2 - 2)/(t+1)\n",
    "\n",
    "# invoca metodo\n",
    "t0 = 0.0\n",
    "tf = 6.0\n",
    "y0 = 2.0\n",
    "h = 0.1\n",
    "t,y = euler_expl(t0,tf,y0,h,f)\n",
    "\n",
    "# plota funcoes \n",
    "yex = t**2 + 2*t + 2 - 2*(t+1)*log(t+1)\n",
    "plt.plot(t,y,'b--',label='$y_h(t)$')\n",
    "plt.plot(t,yex,'k',label='$y(t)$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agora, vamos computar a curva do erro relativo e plotá-la."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fd023bee880>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAia0lEQVR4nO3deXRUVbbH8e8mgAzODEIbFHQpNIoGjbBEm4eKzSAoTigOPNQWARUj0s7dYoNta4OigAgCyjwIKqNTgyOiEgQJiAwPFaIgCNrMEMh5f5ygCQmmklRyq279PmvVSlWdW7d2abFzsu8ZzDmHiIiEV7mgAxARkdKlRC8iEnJK9CIiIadELyISckr0IiIhVz7oAApSvXp1V7du3aDDEBGJG4sWLfrJOVejoLaYTPR169YlPT096DBEROKGmX13uDaVbkREQk6JXkQk5JToRURCToleRCTklOhFREJOiV5EJOSU6EVEQi4mx9GXRIsW+Z/r2BF69IBdu6Bt2/ztXbr4208/wTXX5G/v3h2uuw7Wr4ebb87fft990L49rFwJd9yRv/3RR6FlS1iyBNLS8rf/85/QrBl88gk8/HD+9oEDISUF/vMf6Ncvf/uwYVC/PsycCQMG5G8fOxbq1IHJk2Ho0PztU6dC9erwyiv+dqg5c6BKFXjhBZgyJX/7++/7n/37w6xZedsqV4Y33/T3+/aFuXPztlerBtOm+fsPPQQLFuRtT06GceP8/bQ0/98wt9NPh+HD/f2uXWHVqrztKSn+vx/ATTdBZmbe9vPPhyef9Pevvhq2bHGQtR/27IGsLC754w/8rdXnsH07bQZfxu6sJDDzN4x2DdfSu+1XcNxxtHiyFVQoD+XLQ4WKUL68vnv67gGRfvd++zzRFrpELxIZB7t2w7K1kDYCVq+Gj9JgZyXIPvDbYRlzYcrBDHc6UDnvab6ZBbMPZrj38rYlJcGqD2D2Aqh9Kqzr7bNW1SpQqXLOLwuR0mexuPFIamqq08xYiaoDB+Dzz2H2bN91W7QI/vtf31a5su+W1q372+3kk6FmTTj6aDjqqN9u5ctDdrY/38Hbzp3w88+wdetvPzdu9N3wdet+u23Z8ls8FSvCaadBw4bQuDGcc46/1ShwBrtIocxskXMutaA29eglvLZvh3fe8XWFOXNg82bfy05JgU6dIDUVzjvPJ9vyRfinkJTkbwcdeSSccELhr/vlF19jWbECvv7a/0xPh1df/e2Y5GQ491z/d32zZj7GypUPe0qRSCjRS7g45wvOw4b5BLpnDxx3HLRp44vZrVr5x0E49lho2tTfcvv5Z18A/uILWLzY/+UxfbpvK1/e9/SbNYP/+R9o3hyOP76sI5c4p9KNhMN//+uv/A0bBsuW+TLLjTf6nnuzZkXrsceCn37yJaYFC/wvrs8/h927fV3/7LPhoov8yIMWLXx5SRLe75VulOglvm3d6odcPP+8r5Wfe64fftKpky+phMXevbBwIbz3nr998ol/rnx5X+b585/97dxz85aVJGEo0Uv4bNsGzz3nk/z27X4MYu/ePtElgj17fG//3Xfh7bd92Qd8WefPf4bLLoPWrf3YRUkISvQSHnv3wqBB8K9/+VEsHTrAP/4BjRoFHVmwNm/2g93fessPHt+8GcqV89cDLrsMLr8czjxTQzpDTIlewmHBAvjLX+Crr/xF1b59/agZySs724/mmT3b3xYt8s/Xq+cT/hVXwJ/+FH/XLeR3KdFLfNu+HR55BAYP9sMPhw71vVSJzIYNftro9Om+1793rx951K4dXHWVL/VUqRJ0lFJCSvQSv958E7p185OP7rzTz9k/6qigo4pfO3b4uv706TBjhh/aWaWKr+dfdZVP/sccE3SUUgy/l+i1qJnEpn37oGdPv0BM1arw8ce+Nq8kXzJHHglXXukXlvnxR9/D79LFl8VuusnPBm7XDkaP9r8EJBTUo5fY88MPcO21fgjhPffAU0/BEUcEHVW4ZWfDp5/6Vb6mTvVLNlSo4FdEu/Zaf9E7qIlmEpES9+jNrLWZrTSzNWb2YAHtDcxsgZntNbPeBbQnmdliM5t1aJtIHh984GeCfvklTJrkl/9Tki995cr5iWUDBsC338Jnn/klG1esgFtv9Us8tGsHY8b8tkaQxI1CE72ZJQFDgDZAQ6CTmTU85LCtQE+g/2FOcw+wogRxStg555PMJZf4GvFnn/mx8VL2zKBJE3j6aVi71s/KveceyMiA//1fX9654gqYONHX/CXmRdKjbwKscc6tdc7tAyYBV+Q+wDm3yTm3EMg69MVmlgxcBoyIQrwSRgcO+AW9e/f2CWThQjjjjKCjEvBJ/7zz4N//9j39Tz/1F8UXLYIbbvBJv2NHeO01v0SDxKRIEv2JwPpcjzNznovUQOB+IPv3DjKzrmaWbmbpmzdvLsLpJa7t3QvXXw8jRvghlFOnau2WWGXmJ2A984yv4X/4oS/rfPCB3znjhBOgc2c/UiorX59PAhRJoi9oKl1EV3DNrB2wyTm3qLBjnXPDnXOpzrnUGlqTOzHs2OFXlJw61SePfv00czNelCvnJ10NHgzff++HbF57rR+y2bYt/OEPfnusDz/0F3olUJEk+kygTq7HycAPEZ7/AuByM/sWX/K52MzGFSlCCaetW/2Ijrlz4eWX4d57g45Iiqt8ef//cuRIP2Rz+nT/eMwYv7TyySfDX//ql2COwVF+iSCSRL8QOM3M6plZReB6YEYkJ3fOPeScS3bO1c153Tzn3E3FjlbCYcMGv6764sV+OF+XLkFHJNFyxBF+mYWJE2HTJpgwwe+gNXCgH031xz/6tYlWrw460oRSaKJ3zu0H7gLexo+cmeKcW25m3cysG4CZ1TKzTKAX8KiZZZqZCq2S39atfsr9t9/6Wm6HDkFHJKWlalW/XPSMGb6nP3w41K4Nffr4nbXPOw+efdbPm5BSpQlTUnZ27oRLL/UjNubM8UMpJfF8/z1Mnux7+4sW+esyLVr4UTxXX62JWcWkJRAkePv2wTXX+PHxEycqySeyE0+EXr38Cptffw1/+5tfy+j22/3InQ4d/DaQGq4ZNUr0Uvqys/1Em7fe8n++X3VV0BFJrKhfHx5/HFat8vMn7rrLT9Dq2NGP0e/c2X9vNFyzRJTopXQ5B3ff7ZczeOopuO22oCOSWGQGqal+mO369X401nXXwcyZfmP3E0/0E7Xmz9dwzWJQopfS9a9/wQsvwP33+5tIYZKS4OKL/SS6jRvhjTf8ZuijRsGFF/oNVB580K+HFIPXGGORLsZK6Zk920+I6tQJxo3TZCgpme3bfdKfOBHeeccvndGwof9+deoEp54adISB0sYjUvZWrfLD50491a8lrx2MJJo2b/YzqidM8N8v8Auxderk6/t/+EOw8QVAo26kbG3b5kdOVKwIr7+uJC/RV6OGX2Lho4/gu+/8SptZWX6GdXKyL/UMH+43kBcleomy7Gw/UmLVKpgyxU9/FylNJ53kl1j44gs/XPPvf/eTsO64A2rV8vsLjx3rOyAJSoleoqtfP7/WyYABvlclUpbq1/czb7/+2if+Xr1g2TLf+ahZ00/ImjIFdu0KOtIypRq9RM+sWf7ia+fOfk9SXXyVWOCcX0d/4kQ/EWvjRr88w+WX+yGcrVuHYhczXYyV0rdhAzRqBHXq+L1eK1cOOiKR/A4c8EsnT57sL+Zu2eL3P+jQwSf9li39taU4pIuxUrqys+GWW/xaNhMmKMlL7EpK8iXFF1/0nZO33vLlnBkzfC2/Vi0/qe+dd0I1G1eJXkpu8GB4+21fl//jH4OORiQyFSpAq1Z+ItaPP/rSY7t2vrzTqpVP+rffHoqkr9KNlMyyZX7qesuWfrq66vIS7/bs8cl9yhQ/sGDHDqhWDa680u+iddFF/pdEjFGNXkrH3r1+ksqGDZCR4VceFAmT3bv9X6tTpviOzI4dcPzxvqZ/7bV+FdYYSfqq0UvpeOQRWLrU/+mrJC9hVLmyT+oTJvgds954wy+y9uqr/ucJJ/gd0mbN8h2fGKUevRTP3Lm+XNO9u1+0TCSR7N3ryzvTpvnyzi+/+NE77dr5i7utW5f5jHCVbiS6du6EM87wY48XL9YSB5LY9u2DefP8cM033vBDNitX9j3+q6/2o3mOOabUw1DpRqKrTx+/vsiIEUryIhUr+h78wWWV586FW2+FBQvgxhv9jNy2bX37pk2BhKgevRTN4sV+Vcpbb/WLRolIwbKz/daZ06b5xf3WroVy5eCCC/wIniuvhLp1o/Z2Kt1IdBw4AE2bQmYmrFihTZxFIuWcH7jw+uv+tnSpfz4lxV/s7dABzjqrRMOTVbqR6Bg0CBYtgueeU5IXKQozOPtsX/b88ktYvRr694cjj/R75qak+L0bevUqlclZ6tFLZNat87v5NG/ud47SxCiR6PjxRz9G/403/P2FC4t1GpVupGSc8yv9zZsHy5dHta4oIrlkZ/s6fjH8XqIvX6KgJDFMm+YnhAwYoCQvUpqKmeQLPW2pnFXCY+dOSEuDxo2hZ8+goxGRYlCPXn5f//7w/fcwaRKU19dFJB6pRy+Hl5kJTz0FHTvChRcGHY2IFJMSvRzegw/6i0NPPRV0JCJSAkr0UrBPP4Xx4+G++3QBViTOKdFLfs75C7C1asFDDwUdjYiUkK6uSX4TJ/o1Ol5+2c/cE5G4ph695LVzJzzwAJxzDnTuHHQ0IhIF6tFLXv37+9E2EyaU2uQNESlb+pcsv9m0Cf79b79Zwp/+FHQ0IhIlSvTymyefhD174Ikngo5ERKIookRvZq3NbKWZrTGzBwtob2BmC8xsr5n1zvV8JTP73My+NLPlZvZ4NIOXKFq3zu/92qUL1K8fdDQiEkWF1ujNLAkYAlwKZAILzWyGc+6rXIdtBXoCHQ55+V7gYufcDjOrAHxsZm865z6NSvQSPX37+p9//3uwcYhI1EXSo28CrHHOrXXO7QMmAVfkPsA5t8k5txDIOuR555zbkfOwQs4t9tZFTnSrVvmhlN27w0knBR2NiERZJIn+RGB9rseZOc9FxMySzGwJsAl41zn32WGO62pm6WaWvnnz5khPL9Hw2GNQqZImR4mEVCSJvqCthCLulTvnDjjnUoBkoImZnXmY44Y751Kdc6k1atSI9PRSUl9+6VemTEuDE04IOhoRKQWRJPpMoE6ux8nAD0V9I+fcL8D7QOuivlZK0SOPwLHHQu/ehR4qIvEpkkS/EDjNzOqZWUXgemBGJCc3sxpmdmzO/cpAS+DrYsYq0TZ/vt//9YEHfLIXkVAqdNSNc26/md0FvA0kAaOcc8vNrFtO+4tmVgtIB44Gss0sDWgI1AZG54zcKQdMcc7NKp2PIkX26KO+XHP33UFHIiKlKKIlEJxzc4A5hzz3Yq77G/ElnUMtBRqXJEApJR9+CO+/DwMHQtWqQUcjIqVIM2MTVb9+ULMm3H570JGISClTok9En30G777rL8BWqRJ0NCJSypToE1HfvnD88X6ClIiEnhJ9olm82I+0ufdebSoikiCU6BNNv35wzDEaaSOSQJToE8myZfDaa9Czp0/2IpIQlOgTyRNP+HJNWlrQkYhIGVKiTxQrV8LkyXDnnf5CrIgkDCX6RPHkk36Fyl69go5ERMqYEn0iWLcOxo2Drl39JCkRSShK9Ing2WfBTL15kQSlRB92P/8ML70EnTpp9yiRBKVEH3ZDh8LOnVpvXiSBKdGH2Z498Pzz0Lo1nHVW0NGISECU6MNs7Fj48Ue4//6gIxGRACnRh9WBA9C/P6SmQosWQUcjIgGKaOMRiUMzZsCqVX6SlBW0v7uIJAr16MPIOXj6aTjlFLjqqqCjEZGAqUcfRvPnw6efwpAhUF7/i0USnXr0YfT001C9OnTpEnQkIhIDlOjDZuVKmDnTL16mbQJFBCX68HnuOTjiCOjRI+hIRCRGqIAbJlu3wujRcOONWrxMpABZWVlkZmayZ8+eoEMptkqVKpGcnEyFChUifo0SfZgMHw67dmljEZHDyMzM5KijjqJu3bpYHA47ds6xZcsWMjMzqVevXsSvU+kmLLKyYPBguOQSaNQo6GhEYtKePXuoVq1aXCZ5ADOjWrVqRf6LRD36sHj1Vfj+exg2LOhIRGJavCb5g4oTv3r0YeCcX3O+fn1o0yboaEQkxqhHHwaffALp6fDCC1BOv7tFJC9lhTB49lk47jjo3DnoSEQkBinRx7tvvoHXX/f7wVatGnQ0IhKBYcOGUbt2bVJSUn69ZWRkMGzYMGrVqsXZZ5/NqaeeypgxY6Lyfkr08W7QIL865V13BR2JiERo6dKl9OvXjyVLlvx6a9SoEUuXLqVPnz58+eWXTJw4kV5R2udZiT6ebd8OI0fCtddCcnLQ0YhIhDIyMkhJSSnw+fr16wNQr149KlasGJX308XYeDZ6NGzbpglSIsWRlgZLlkT3nCkpMHBgoYctX76cW265hXI5gyd69OhB165df030zjkGDx7ME088EZWwlOjjVXa2L9s0aQJNmwYdjYhEaP369dSsWZOlS5fme3779u20bduW77//nrPOOos+ffpE5T2V6OPVO+/4HaTGjQs6EpH4FEHPuzQsXbqUBg0aFPh88+bNmTdvHj///DNnnnkmCxYsoFmzZiV+T9Xo49WgQVCrlq/Pi0jcyMjIKDDRZ2Rk0LhxYwCOO+44brjhBmbPnh2V91Sij0erV8OcOdCtG0TpYo2IlI2MjAzGjh3767DKxo0bs2PHjjyJHqB9+/bMmTMnKu8ZUenGzFoDzwFJwAjn3L8OaW8AvAycAzzinOuf83wdYAxQC8gGhjvnnotK5Ils8GCoUAHuuCPoSESkiMaPHx/R882bN2fx4sVRec9CE72ZJQFDgEuBTGChmc1wzn2V67CtQE+gwyEv3w/c55z7wsyOAhaZ2buHvFaKYvt2ePlluO46X7oRESlEJKWbJsAa59xa59w+YBJwRe4DnHObnHMLgaxDnt/gnPsi5/52YAVwYlQiT1SjR/tkf/fdQUciInEikkR/IrA+1+NMipGszawu0Bj47DDtXc0s3czSN2/eXNTTJ4aDQyqbNvXDKkWkyJxzQYdQIsWJP5JEX9Dix0V6JzM7EpgGpDnnthV0jHNuuHMu1TmXWqNGjaKcPnEcHFLZs2fQkYjEpUqVKrFly5a4TfYHd5iqVKlSkV4XycXYTKBOrsfJwA+RvoGZVcAn+fHOudeKFJ3k9fzzvi5/zTVBRyISl5KTk8nMzCSeqwYH94wtikgS/ULgNDOrB3wPXA/cEMnJzW+FMhJY4Zx7pkiRSV6rV8Obb0KfPhpSKVJMFSpUKNJeq2FRaKJ3zu03s7uAt/HDK0c555abWbec9hfNrBaQDhwNZJtZGtAQOAu4GcgwsyU5p3zYORedwaGJZMgQDakUkWKJaBx9TmKec8hzL+a6vxFf0jnUxxRc45ei2LHDD6ns2FFDKkWkyDQzNh6MGeNXqdSa8yJSDEr0sc45PxM2NVWrVIpIsWj1ylg3bx6sWOEnSpmqYCJSdOrRx7pBg6B6dV+fFxEpBiX6WPbttzBzpt/4u4gTJEREDlKij2VDh/pyTbduQUciInFMiT5W7d4NI0ZAhw5Qp06hh4uIHI4SfayaMAG2btUqlSJSYkr0sejgkMpGjaB586CjEZE4p+GVsWj+fFiyBIYN05BKESkx9ehj0fPPw3HHwY03Bh2JiISAEn2sycyE116D226DqlWDjkZEQkCJPtYMHepr9HfeGXQkIhISSvSxZM8eGD4cLr8c6tYNOhoRCQkl+lgyaRL89JOGVIpIVCnRxwrn/EXYM86Aiy4KOhoRCRENr4wV8+fD4sUaUikiUacefazQkEoRKSVK9LFAQypFpBQp0ccCDakUkVKkRB80DakUkVKmRB+0CRM0pFJESpUSfZCcg4ED4ayzNKRSREqNhlcGae5cyMiAUaM0pFJESo169EF69lmoWRM6dQo6EhEJMSX6oHz9NcyZAz16aONvESlVSvRBee45OOII6N496EhEJOSU6IOwZQuMHu1nwdasGXQ0IhJySvRBGD4cdu+Ge+8NOhIRSQBK9GVt3z4YNAguvRTOPDPoaEQkAWh4ZVmbMgU2bICRI4OOREQShHr0Zck5P6SyQQNo1SroaEQkQahHX5Y++gi++AJefBHK6XesiJQNZZuy9PTTUL063Hxz0JGISAJRoi8ry5bB7Nl+8bIqVYKORkQSiBJ9WXn6aZ/gtea8iJSxiBK9mbU2s5VmtsbMHiygvYGZLTCzvWbW+5C2UWa2ycyWRSvouLNuHUycCF27QrVqQUcjIgmm0ERvZknAEKAN0BDoZGYNDzlsK9AT6F/AKV4BWpcszDj3zDP+Z69ewcYhIgkpkh59E2CNc26tc24fMAm4IvcBzrlNzrmFQNahL3bOfYj/RZCYtmyBl16CG26AOnWCjkZEElAkif5EYH2ux5k5z0WVmXU1s3QzS9+8eXO0Tx+cIUNg1y64//6gIxGRBBVJoi9oRwwX7UCcc8Odc6nOudQaNWpE+/TB2LkTnn8e2rWDM84IOhoRSVCRJPpMIHfNIRn4oXTCCZlRo3zp5oEHgo5ERBJYJIl+IXCamdUzs4rA9cCM0g0rBLKyYMAAuOACuPDCoKMRkQRW6BIIzrn9ZnYX8DaQBIxyzi03s2457S+aWS0gHTgayDazNKChc26bmU0EWgDVzSwTeMw5F/4VvSZPhu++8ytViogEyJyLerm9xFJTU116enrQYRTf/v2+Jn/EEbBkida1EZFSZ2aLnHOpBbVpUbPSMHEirFoF06YpyYtI4JSFom3/fvjHP+Dss6FDh6CjERFRjz7qxo+HNWvg9dfVmxeRmKBMFE3790PfvtC4MVxxReHHi4iUAfXoo2nsWPi//4Pp08EKmmcmIlL21KOPlqws35s/91xo3z7oaEREfqUefbSMGQPffOOXPFBvXkRiiHr00bBvH/TrB+edB5ddFnQ0IiJ5qEcfDa+8At9+61eqVG9eRGKMevQltX07PPYYnH8+tGkTdDQiIvmoR19STz0FGzfCG2+oNy8iMUk9+pJYt86vUHnDDdC0adDRiIgUSIm+JB5+2P988slg4xAR+R1K9MX1+ed+uYNeveCkk4KORkTksJToi8M5n+BPOAEefDDoaEREfpcuxhbH1Kkwfz689BIcdVTQ0YiI/C716Itqzx6/B2yjRnDLLUFHIyJSKPXoi2rgQL/UwbvvQlJS0NGIiBRKPfqiWLUKHn/cbyjSsmXQ0YiIRESJPlLZ2XDbbVCpErzwQtDRiIhETKWbSA0ZAh9/7Ne1qV076GhERCKmHn0kvvnGD6Ns0wY6dw46GhGRIlGiL4xz8Je/+Auvw4ZpPRsRiTsq3RTmpZdg3jyf5OvUCToaEZEiU4/+96xfD717w8UXw+23Bx2NiEixKNEfTlaWr8cfOOB79SrZiEicUunmcO67D95/H0aPhlNOCToaEZFiU4++ICNHwqBBfuEyjbIRkTinRH+oTz6B7t3h0kv97lEiInFOiT63zEy46iq/vvykSVBelS0RiX/KZAft3g1XXgk7d/rhlMcfH3REIiJRoUQPsHcv3HQTpKfD9OnQsGHQEYmIRI0S/bZtvlwzdy48+yxcfnnQEYmIRFViJ/qNG/36NcuWwZgxcPPNQUckIhJ1iZvoV6+GVq3gxx9h5kxo3TroiERESkViJvrPPoP27f2CZe+9B02aBB2RiEipSazhlevXw623QrNmULWq3+BbSV5EQi6iRG9mrc1spZmtMbMHC2hvYGYLzGyvmfUuymvLxNat8Ne/wmmnwfjxcO+9sGgRnH56IOGIiJSlQks3ZpYEDAEuBTKBhWY2wzn3Va7DtgI9gQ7FeG10OQebN/sa/Jo1/kLrSy/50TWdO/s9X08+udTeXkQk1kRSo28CrHHOrQUws0nAFcCvydo5twnYZGaXFfW1UXPgADRtSouUJf5+Lh2vPIMevSayq8GptB3fNt9Lu6R0oUtKF37a9RPXTLkmX3v31O5cd+Z1rP/vem5+Pf/InPvOv4/29duz8qeV3DHrjnztjzZ/lJantGTJxiWkvZWWr/2fl/yTZnWa8cn6T3h47sP52ge2HkhKrRT+s/Y/9PuwX772Ye2GUb96fWaunMmABQPytY+9cix1jqnD5GWTGZo+NF/71I5TqV6lOq8seYVXlrySr33OjXOoUqEKLyx8gSnLp+Rrf7/L+wD0/6Q/s1bNytNWuUJl3rzxTQD6ftCXud/MzdNerUo1pnWcBsBD/3mIBZkL8rQnH53MuKvGAZD2VhpLNi7J0356tdMZ3n44AF1ndmXVllV52lNqpTCw9UAAbnrtJjK3ZeZpPz/5fJ5s+SQAV0+5mi27tuRpv6TeJfztf/4GQJvxbdidtTtPe7vT29G7mf8jtsUrLThUxzM60uO8HuzK2qXvnr57edoL+u4d/DzRFkmiPxFYn+txJtA0wvNH/Foz6wp0BTjppJMiPH0uSUnQoAHU2gCVq0Dlyv5WqRKceR00agRZu4p+XhGROGfOud8/wOxaoJVz7i85j28Gmjjn7i7g2D7ADudc/6K+NrfU1FSXnp5ejI8jIpKYzGyRcy61oLZILsZmArn30EsGfojwvUvyWhERiYJIEv1C4DQzq2dmFYHrgRkRnr8krxURkSgotEbvnNtvZncBbwNJwCjn3HIz65bT/qKZ1QLSgaOBbDNLAxo657YV9NpS+iwiIlKAQmv0QVCNXkSkaEpaoxcRkTimRC8iEnJK9CIiIadELyIScjF5MdbMNgPfFfPl1YGfohhOkMLyWcLyOUCfJRaF5XNAyT7Lyc65GgU1xGSiLwkzSz/cled4E5bPEpbPAfossSgsnwNK77OodCMiEnJK9CIiIRfGRD886ACiKCyfJSyfA/RZYlFYPgeU0mcJXY1eRETyCmOPXkREclGiFxEJudAk+pjYhDwKzGyUmW0ys2VBx1JSZlbHzN4zsxVmttzM7gk6puIys0pm9rmZfZnzWR4POqaSMLMkM1tsZrMKPzp2mdm3ZpZhZkvMLK5XQjSzY81sqpl9nfNv5vyonTsMNfqcTchXkWsTcqBTqW5CXkrMrDmwAxjjnDsz6HhKwsxqA7Wdc1+Y2VHAIqBDnP5/MaCqc26HmVUAPgbucc59GnBoxWJmvYBU4GjnXLug4ykuM/sWSHXOxf2EKTMbDXzknBuRs39HFefcL9E4d1h69L9uQu6c2wcc3IQ87jjnPgS2Bh1HNDjnNjjnvsi5vx1Ygd9HOO44b0fOwwo5t7jsJZlZMnAZMCLoWMQzs6OB5sBIAOfcvmgleQhPoi9oE/K4TChhZWZ1gcbAZwGHUmw55Y4lwCbgXedcvH6WgcD9QHbAcUSDA94xs0Vm1jXoYErgFGAz8HJOSW2EmVWN1snDkuitgOfisrcVRmZ2JDANSHPObQs6nuJyzh1wzqXg9z5uYmZxV1ozs3bAJufcoqBjiZILnHPnAG2AO3NKn/GoPHAOMNQ51xjYCUTtWmNYEr02IY9ROfXsacB459xrQccTDTl/Ur8PtA42kmK5ALg8p7Y9CbjYzMYFG1LxOed+yPm5CXgdX8aNR5lAZq6/EqfiE39UhCXRaxPyGJRzAXMksMI590zQ8ZSEmdUws2Nz7lcGWgJfBxpUMTjnHnLOJTvn6uL/ncxzzt0UcFjFYmZVcy7yk1Pm+DMQl6PVnHMbgfVmVj/nqUuAqA1aKHRz8HhwuA3MAw6rWMxsItACqG5mmcBjzrmRwUZVbBcANwMZObVtgIedc3OCC6nYagOjc0Z4lQOmOOfiemhiCJwAvO77E5QHJjjn3go2pBK5Gxif01ldC9wSrROHYniliIgcXlhKNyIichhK9CIiIadELyISckr0IiIhp0QvIhJySvQiIiGnRC8iEnL/D9//GsTmuLh1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def erro_relativo(y,yh):\n",
    "    return abs(y-yh)/abs(yh)\n",
    "\n",
    "e = erro_relativo(yex,y)+0.1\n",
    "\n",
    "plt.plot(t,e,'r',label='$ER$')\n",
    "plt.legend()\n",
    "emax = np.max(e)\n",
    "emin = np.min(e)\n",
    "\n",
    "plt.plot(t,np.ones(t.shape)*emax,'b--')\n",
    "plt.plot(t,np.ones(t.shape)*emin,'g--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercício:** Resolva o PVI do exemplo anterior usando por computador para $h = 0.2, 0.1, 0.05$. Produza um `script` para imprimir, para cada $h$, os dados da saída na forma da tabela a seguir: \n",
    "\n",
    "|$t$|$y_{h}$$(t)$|$EA$|$ER$|\n",
    "|--|--|--|--|\n",
    "|$\\vdots$|$\\vdots$|$\\vdots$|$\\vdots$|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Análise de erro\n",
    "\n",
    "A análise de erro para o método de Euler tem o propósito de entender como o método funciona para que se possa estimar o erro ao usá-lo e, possivelmente, acelerar sua convergência. Procedimentos similares são aplicáveis a métodos numéricos mais eficientes. \n",
    "\n",
    "Para a análise, assumiremos que o PVI padrão possui solução única $y(t)$ em $[t_0,b]$ e que esta solução tem uma segunda derivada $y''(t)$ limitada neste intervalo. \n",
    "\n",
    "Consideremos a série de Taylor para aproximar $y(t_{n+1})$:\n",
    "\n",
    "$$y(t_{n+1}) = y(t_{n}) + hy'(t_{n}) + \\dfrac{h^2}{2}y''(\\xi_{n}), \\ \\ t_n \\leq \\xi_n \\leq t_{n+1}$$\n",
    "\n",
    "Uma vez que $y(t)$ satisfaz a EDO, temos que \n",
    "\n",
    "$$y(t_{n+1}) = y(t_{n}) + hf(t_{n},y(t_n)) + \\dfrac{h^2}{2}y''(\\xi_{n}), \\ \\ t_n \\leq \\xi_n \\leq t_{n+1}$$\n",
    "\n",
    "O termo $T_{n+1} = \\dfrac{h^2}{2}y''(\\xi_{n})$, _erro de truncamento_ para o método de Euler, é o erro na aproximação\n",
    "\n",
    "$$y(t_{n+1}) \\approx y(t_{n}) + hy'(t_{n}).$$\n",
    "\n",
    "Assim, tendo em vista que \n",
    "\n",
    "$$y_{n+1} \\approx y_n + hf(t_{n},y_n),$$\n",
    "\n",
    "subtraímos as equações para obter \n",
    "\n",
    "$$y(t_{n+1}) - y_{n+1} \\approx y(t_{n}) - y_n + h[ f(t_{n},y(t_n)) - f(t_{n},y_n)] + \\dfrac{h^2}{2}y''(\\xi_{n}),$$\n",
    "\n",
    "mostrando que o error em $y_{n+1}$ consiste de duas partes: \n",
    "\n",
    "- **erro propagado:** $T_{n+1} = y(t_{n}) - y_n + h[ f(t_{n},y(t_n)) - f(t_{n},y_n)]$.\n",
    "- **erro de truncamento:** $T_{n+1} = \\dfrac{h^2}{2}y''(\\xi_{n})$.\n",
    "\n",
    "O erro propagado pode ser simplificado pelo teorema do valor médio: \n",
    "\n",
    "$$f(t_{n},y(t_n)) - f(t_{n},y_n) = \\dfrac{\\partial f(t_n,\\zeta_n)}{\\partial y}[y(t_n) - y_n], \\ \\ y(t_n) \\leq \\zeta_n < y_n$$\n",
    "\n",
    "Se definirmos o erro $e_k = y(t_k) - y_k, k \\geq 0$, podemos reescrever a última equação como \n",
    "\n",
    "$$e_{n+1} = \\left[ 1 + h \\dfrac{\\partial f(t_n,\\zeta_n)}{\\partial y} \\right]e_n + T_{n+1},$$\n",
    "\n",
    "assim obtendo uma equação geral para análise de erro do método de Euler.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Teoremas para limite de erro \n",
    "\n",
    "**Teorema:** Suponha que $f(t,y)$ seja definida em um conjunto convexo $D \\subset \\mathbb{R}^2$. Se existe uma constante $L>0$ com \n",
    "\n",
    "$$\\left| \\dfrac{\\partial f}{\\partial y}(t,y)\\right| \\leq L, \\ \\ \\forall(t,y) \\in D$$\n",
    "\n",
    "**Teorema (limite de erro):** Suponha que f seja contínua e satisfaça a condição de Lipschitz com constante $L$ em \n",
    "$D = \\{ (t,y); \\ a \\leq t \\leq b \\}$ e que exista uma constante $M$ com $|y''(t)| \\leq M, \\ \\ \\forall t \\in [a,b]$. Seja $y(t)$ única solução do PVI \n",
    "\n",
    "$$y' = f(t,y), \\quad y(a) = y_0, \\quad t_0 \\leq t \\leq b$$\n",
    "\n",
    "e sejam $w_0,w_1,\\ldots,w_N$ aproximações geradas pelo método de Euler para algum inteiro positivo $N$. Então, para cada $i=0,1,2,\\ldots,N$, \n",
    "\n",
    "$$|y(t_i)-w_i| \\leq \\frac{hM}{2L}[e^{L(t_i-a)} -1].$$\n",
    "\n",
    "O último teorema fornece um limitante de erro para o método de Euler. Ele evidencia a dependência linear do tamanho de passo $h$. Logo, à medida que $h->0$, uma maior precisão nas aproximações deve ser obtida.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exemplo:** Consideremos o PVI\n",
    "\n",
    "$$\\begin{cases}\n",
    "y'(t) = y - t^2 + 1 \\\\\n",
    "y(0) = 0.5 \\\\\n",
    "0 \\le t \\le 2 \\\\ \n",
    "h = 0.2\n",
    "\\end{cases}$$\n",
    "\n",
    "Uma vez que $f(t,y) = y - t^2 + 1$, temos $\\dfrac{\\partial f}{\\partial y}(t,y) = 1, \\forall y$. Assim, $L=1$. A solução exata deste problema é $y(t) = (t+1)^2 - 0.5e^t$, de modo que $y''(t) = 2 - 0.5e^t$ e $|y''(t)| \\leq 0.5e^2 - 2 = M, \\ \\ \\forall t \\in [0,2]$. Pelo teorema do limite de erro, temos que \n",
    "\n",
    "$$|y_i - w_i| \\leq 0.1(0.5e^2-2)(e^{t_i}-1)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imprimindo comparação...\n",
      "\n",
      " ti |    ei   |   Ei\n",
      "\n",
      "0.0 | 0.00000 | 0.00000\n",
      "\n",
      "0.2 | 0.02930 | 0.03752\n",
      "\n",
      "0.4 | 0.06209 | 0.08334\n",
      "\n",
      "0.6 | 0.09854 | 0.13931\n",
      "\n",
      "0.8 | 0.13875 | 0.20767\n",
      "\n",
      "1.0 | 0.18268 | 0.29117\n",
      "\n",
      "1.2 | 0.23013 | 0.39315\n",
      "\n",
      "1.4 | 0.28063 | 0.51771\n",
      "\n",
      "1.6 | 0.33336 | 0.66985\n",
      "\n",
      "1.8 | 0.38702 | 0.85568\n",
      "\n",
      "2.0 | 0.43969 | 1.08264\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def f(t,y):\n",
    "    return y - t**2 + 1\n",
    "\n",
    "# solucao numerica \n",
    "t,y = euler_expl(0,2,0.5,0.2,f)\n",
    "\n",
    "# solucao exata\n",
    "yex = (t+1)**2 - 0.5*np.exp(t)\n",
    "\n",
    "# erro\n",
    "erro = np.abs(y - yex)\n",
    "\n",
    "# limite de erro \n",
    "lim_erro = 0.1*(0.5*np.exp(1)**2 - 2)*(np.exp(t) - 1.)\n",
    "\n",
    "# tabela\n",
    "\n",
    "print(\"Imprimindo comparação...\\n\\n ti |    ei   |   Ei\\n\")\n",
    "for i in range(len(y)): \n",
    "    print(\"{0:.1f} | {1:0.5f} | {2:0.5f}\\n\".format(t[i],erro[i],lim_erro[i]))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problemas\n",
    "\n",
    "1) Resolva os seguintes problemas usando o método de Euler com passos $h = 0.2,0.1,0.05$. Compute o erro relativo usando a solução exata $y(t)$. Para valores selecionados de $t$, observe a razão com que o erro diminui à medida que $h$ é reduzido pela metade. \n",
    "a. $y'(t) = [\\cos(y(t))]^2, 0 \\leq t \\leq 10, \\quad y(0) = 0; \\quad y(t) = \\tan^{-1}(t)$\n",
    "\n",
    "b. $y'(t) = \\frac{1}{1+t^2} − 2[y(t)]^2, 0 \\leq t \\leq 10, \\quad  y(0) = 0; \\quad y(t) = \\frac{t}{1+t^2}$\n",
    "\n",
    "c. $y'(t) = \\frac{1}{4}y(t)\\left[1 - \\frac{1}{20}y(t)􏰍\\right], 0 \\leq t \\leq 20, \\quad  y(0) = 1; \\quad y(t) = \\frac{20}{1 + 19e^{-t/4}}$\n",
    "\n",
    "d. $y'(t)= -[y(t)]^2, y(t) = \\frac{1}{t}, 1 \\leq t \\leq 10, \\quad  y(1) = 1 \\quad y(t) = \\frac{1}{t}$\n",
    "\n",
    "2) Considere o problema linear \n",
    "\n",
    "$$y'(t) = \\lambda y(t) + (1 - \\lambda) \\cos(t) - (1 + \\lambda) \\textrm{sen}(t), \\quad y(0) = 1.$$\n",
    "\n",
    "A solução exata é $y(t) = \\textrm{sen}(t) + \\cos(t)$. Resolva este problema usando o método de Euler com vários valores de $\\lambda$ e $h$, para $0 \\leq t \\leq 10$. Comente sobre os resultados.\n",
    "\n",
    "a. $\\lambda = −1; \\ \\ h = 0.5, 0.25, 0.125$.\n",
    "\n",
    "b. $\\lambda =1; \\ \\ h=0.5,0.25,0.125.$\n",
    "\n",
    "c. $\\lambda = −5; h = 0.5, 0.25, 0.125, 0.0625.$\n",
    "\n",
    "d. $\\lambda = 5; h = 0.125, 0.0625.$\n",
    "\n",
    "\n",
    "3) Faça uma análise do erro obtido pelo método de Euler ao ser resolver o caso a. do problema 2, $h = 0.25$."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "ipub": {
   "listcode": false,
   "listfigures": false,
   "listtables": false,
   "titlepage": {
    "author": "Prof. Dr. Gustavo Oliveira",
    "email": "gustavo.oliveira@ci.ufpb.br",
    "institution": [
     "Universidade Federal da Paraíba"
    ],
    "tagline": "Aula 3.",
    "title": "Soluções Numéricas de EDOs"
   },
   "toc": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}