2-Zero de funções.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Zeros de funções\n",
    "\n",
    "Neste *notebook* vamos aprender sobre alguns métodos para encontrar zeros de funções. Este é um problema clássico que aparece em diversas áreas: temos uma função f(x) e queremos encontrar o valor de x que faz com que a função f possua um determinado valor y. Em alguns casos, é possível simplesmente isolar x e obter a resposta. Em outros, f(x) pode ser muito complicado e precisamos encontrar a resposta numericamente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Definimos a nossa função de teste e o valor de y\n",
    "def f(x):\n",
    "    return x*x + 2*x + 1\n",
    "\n",
    "# Queremos encontrar o valor de x para o qual f(x) = y0\n",
    "y0 = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Para entendermos um pouco melhor o nosso problema, vamos fazer um gráfico. Em Python, podemos utilizar as \n",
    "bibliotecas `matplotlib` e `numpy`. A primeira é uma biblioteca de gráficos propriamente dita enquanto a \n",
    "segunda é uma biblioteca de computação científica com vários recursos como vetores, matrizes, funções \n",
    "especiais, números aleatórios, etc.\n",
    "\n",
    "O primeiro passo é importar as duas bibliotecas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "O numpy (agora chamado de `np` no nosso código) possui várias funções para gerar sequências (ou vetores) de valores numéricos. Queremos gerar uma sequência de valores para x para calcularmos os respectivos y gerando\n",
    "os pontos que aparecerão no gráfico. A função que faz isto se chama `linspace`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "X = np.linspace(0, 10)   # cria uma sequência de números começando em zero e terminando em 10\n",
    "Y = f(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Aqui utilizamos `np.linspace` pois importamos o módulo numpy com o \"apelido\" `np`. Poderíamos ter importando a \n",
    "função linspace diretamente (`from numpy import linspace`) ou utilizado o nome completo do módulo (`import numpy`). De todo modo, os apelidos que escolhemos, `np` e `plt`, são considerados padrão na comunidade Python.\n",
    "\n",
    "Podemos imprimir X e Y para ver os resultados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "print('X:', X)\n",
    "print('Y:', Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Agora vamos aos gráficos! O comando mais simples possível é a função `plt.plot(X, Y)`, que faz um gráfico simples da sequência de valores X pela sequência Y. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.plot(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Podemos acrescentar outros elementos chamando vários comandos de gráficos na mesma célula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.plot(X, Y)\n",
    "plt.title('Busca de zero de funções')  # define um título\n",
    "plt.xlabel('argumento')                # nome para eixo x\n",
    "plt.ylabel('valor de f(x)')            # nome para eixo y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Podemos facilmente desenhar várias linhas. Basta chamar a função `plt.plot` várias vezes. Queremos desenhar \n",
    "$f(x)$ com uma linha horizontal no valor $y_0$.\n",
    "\n",
    "Para criar os pontos para a linha horizontal, escolhemos como valores para X os pontos 0 e 10 e para Y o valor constante $y_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.plot(X, Y)                # f(x)\n",
    "plt.plot([0, 10], [y0, y0])   # função horizontal\n",
    "plt.title('Ponto de intersecção')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "É possível verificar visualmente que o valor de $x$ que faz com que $f(x) = y_0$ está entre 5 e 6. Precisamos refinar um pouquinho mais a nossa busca para encontrar qual é o valor correto.\n",
    "\n",
    "Um método pouco elegante é fazer uma busca por força bruta: percorremos vários valores no intervalo e tentamos encontrar aquele que a diferença entre $f(x)$ e $y_0$ seja igual ou muito próxima de zero. O código abaixo faz exatamente isto. Dividimos o intervalo em 25 valores e buscamos qual x produz o melhor resultado."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "for i in range(25):\n",
    "    x = 5 + i / 25\n",
    "    print('%2d) x: %.2f;   f(x)-y0 = %6.3f' %(i, x, f(x) - y0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Podemos ver que o melhor resultado desta lista foi $x = 5.48$, onde $f(x)$ difere de $y_0$ por um fator de apenas $0.010$. É lógico que é possível melhorar este resultado. Uma abordagem seria repetir um passo similar ao anterior, mas agora em um intervalo próximo de $x=5.48$ (quem sabe entre $x=5.44$ e $5.52$). A partir daí encontraríamos uma estimativa mais refinada e poderíamos continuar até obter um resultado satisfatório.\n",
    "\n",
    "Este método de busca por força bruta, por mais que funcione, não parece muito elegante. De fato, este método é ruim sob vários aspectos. Primeiramente (fora Temer!), exige um grande número de avaliações de $f(x)$. Isto pode ser irrelevante no nosso caso, mas existem situações em que a função $f(x)$ pode envolver cálculos extremamente complexos e um grande número de requisições para $f(x)$ significaria um programa lento. Além deste problema, o método exige uma escolha manual dos intervalos de início e fim em cada etapa. Isto é inaceitável e queremos métodos em que o computador consiga realizar todos os cálculos automaticamente.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Método da bisseção\n",
    "\n",
    "Iniciamos a nossa discussão com o método da bisseção pois este se trata de uma versão mais sofisticada do raciocínio que fizemos anteriormente.\n",
    "\n",
    "Começamos mudando ligeiramente o nosso problema. Ao invés de encontrar o valor de x para o qual $f(x) = y_0$, passamos $y_0$ para o lado esquerdo da igualdade e buscamos o valor de $x$ para o qual $g(x) = f(x) - y_0 = 0$. Em matemática este problema é conhecido como encontrar o zero de uma função. É lógico que é trivial transformar um problema em que buscamos $x$ que produz um valor arbitrário (ex.: 42) em um problema de zero de funções. Utilizar o zero, no entanto, simplifica alguns cálculos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Definimos g(x), a função que queremos encontrar os zeros.\n",
    "def g(x):\n",
    "    return f(x) - y0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "O método da bisseção, assim como o método de força bruta mostrado anteriormente, exige um intervalo de busca inicial. Este intervalo deve ser escolhido para conter o zero da função e deve conter um valor para o qual $g(x)$ seja positivo e outro valor para o qual $g(x)$ seja negativo. Chamamos estes valores de $x_a$ e $x_b$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Nosso intervalo de busca inicial é entre 0 e 10\n",
    "# (No método de força bruta, usamos entre 5 e 6. Vamos utilizar\n",
    "# um intervalo maior só por segurança.)\n",
    "x0, x1 = 0, 10\n",
    "\n",
    "# Verificamos os valores de g(x) para descobrir quem é xa e xb.\n",
    "# No nosso caso sabemos que xa = x0 e xb = x1, mas caso a função\n",
    "# seja decrescescente no intervalo, esta relação se inverte\n",
    "g0, g1 = g(x0), g(x1)\n",
    "if g0 < g1:\n",
    "    xa, xb = x0, x1\n",
    "    ga, gb = g0, g1\n",
    "else:\n",
    "    xa, xb = x1, x0\n",
    "    ga, gb = g1, g0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Fazemos um gráfico dos resultados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.plot(X, Y - y0)\n",
    "plt.plot([xa, xb], [0, 0], 'ko-')     # 'k-' é uma linha (-) preta (k, de black) com círculos (o)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "O próximo passo consiste em avaliar o valor de g(x) no ponto central do intervalo e investigar qual será o novo intervalo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xmed = (xa + xb) / 2\n",
    "gmed = g(xmed)\n",
    "print(gmed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vimos que o valor no ponto médio é igual à -6.0, e portanto negativo. Deste modo, sabemos que o zero de g(x) deve estar entre este valor e xb. Atualizamos nossas variáveis e mostramos o resultado em um gráfico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xa = xmed\n",
    "ga = gmed\n",
    "\n",
    "plt.plot(X, Y - y0)\n",
    "plt.plot([xa, xb], [0, 0], 'ko-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que o intervalo reduziu pela metade. Agora repetimos outra vez o mesmo raciocínio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xmed = (xa + xb) / 2\n",
    "gmed = g(xmed)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Observe que desta vez o valor de g(xmed) ficou positivo. Isto significa que devemos substituir $x_b$ e não $x_a$. Novamente, o intervalo de valores aceitáveis para o zero da função reduziu pela metade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xb = xmed\n",
    "gb = gmed\n",
    "\n",
    "plt.plot(X, Y - y0)\n",
    "plt.plot([xa, xb], [0, 0], 'ko-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Agora juntamos o raciocínio realizado nas duas etapas anteriores em um único passo em que o computador decide\n",
    "automaticamente qual dos dois valores ($x_a$ ou $x_b$) deve ser atualizado:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xmed = (xa + xb) / 2\n",
    "gmed = g(xmed)\n",
    "if gmed <= 0:\n",
    "    xa, ga = xmed, gmed\n",
    "else:\n",
    "    xb, gb = xmed, gmed\n",
    "\n",
    "print('intervalo:', xa, xb)\n",
    "plt.plot(X, Y - y0)\n",
    "plt.plot([xa, xb], [0, 0], 'ko-')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Observe que o intervalo aceitável diminui, mas sempre contêm o valor do zero da função. Podemos verificar isto \n",
    "executando várias vezes a célula acima.\n",
    "\n",
    "É lógico que quanto mais repetições realizarmos, melhor será a estimativa do intervalo da função. Cada repetição reduz o intervalo pela metade. Deste modo, após 10 repetições teríamos um intervalo $2^{10} = 1024$ vezes menor que o intervalo inicial. Nada mal!\n",
    "\n",
    "Não queremos executar a célula acima várias vezes manualmente. Vamos então programar o computador para fazer isto automaticamente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Reseta o intervalo inicial (vamos começar do zero)\n",
    "x0, x1 = 0, 10\n",
    "\n",
    "# Calcula xa/ga e xb/gb\n",
    "g0, g1 = g(x0), g(x1)\n",
    "if g0 < g1:\n",
    "    xa, xb = x0, x1\n",
    "    ga, gb = g0, g1\n",
    "else:\n",
    "    xa, xb = x1, x0\n",
    "    ga, gb = g1, g0\n",
    "\n",
    "# Atualiza 20 vezes o intervalo\n",
    "for i in range(20):\n",
    "    xmed = (xa + xb) / 2\n",
    "    gmed = g(xmed)\n",
    "    if gmed <= 0:\n",
    "        xa, ga = xmed, gmed\n",
    "    else:\n",
    "        xb, gb = xmed, gmed\n",
    "        \n",
    "    print('%2d) centro: %.6f,  g(x): %9.2e' % (i, xmed, gmed))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que o valor de x converge para 5.487... e o valor de g(x) se aproxima de zero. Podemos refazer o cálculo acima acompanhando o valor do erro, para depois mostrarmos em um gráfico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Reseta o intervalo inicial (vamos começar do zero)\n",
    "x0, x1 = 0, 10\n",
    "\n",
    "# Calcula xa/ga e xb/gb\n",
    "g0, g1 = g(x0), g(x1)\n",
    "if g0 < g1:\n",
    "    xa, xb = x0, x1\n",
    "    ga, gb = g0, g1\n",
    "else:\n",
    "    xa, xb = x1, x0\n",
    "    ga, gb = g1, g0\n",
    "\n",
    "# Atualiza 20 vezes o intervalo\n",
    "g_lista = []\n",
    "for i in range(20):\n",
    "    xmed = (xa + xb) / 2\n",
    "    gmed = g(xmed)\n",
    "    if gmed <= 0:\n",
    "        xa, ga = xmed, gmed\n",
    "    else:\n",
    "        xb, gb = xmed, gmed\n",
    "        \n",
    "    g_lista.append(abs(gmed))\n",
    "    \n",
    "plt.plot(g_lista)\n",
    "plt.xlabel('número de iterações')\n",
    "plt.ylabel('|g(x)|')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "\n",
    "### Critério de parada\n",
    "\n",
    "Vimos no gráfico anterior que o erro claramente reduz com o número de iterações. Mas quantas iterações devemos \n",
    "realizar? A resposta é sempre \"depende\".\n",
    "\n",
    "O método da bisseção, assim como vários outros métodos numéricos atinge a resposta correta apeanas após um número infinito de iterações. Se por um lado é impossível esperar estas infinitas iterações, por outro, raramente precisamos do valor \"completo\" da solução com todas suas infinitas casas decimais. Na prática precisamos apenas de um valor \"próximo o suficiente\" do correto e é lógico que o que é \"suficiente\" depende muito da aplicação.\n",
    "\n",
    "No método da bisseção escolhemos tipicamente dois critérios para definir o que é \"bom o suficiente\". Eles \n",
    "se traduzem em uma margem de tolerância para y ou para x.\n",
    "\n",
    "1. $|g(x)| < y_{tol}$\n",
    "2. $|x_b - x_a| < x_{tol}$\n",
    "\n",
    "Tipicamente, paramos quando um dos dois critérios for atingido. É lógico que podemos também adotar apenas um dos dois critérios ou, se quisermos ser mais rigorosos, paramos apenas quando os dois critérios forem atingidos. O código abaixo implementa a bisseção com o critério de parada."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Define tolerâncias\n",
    "x_tol = 1e-6\n",
    "y_tol = 1e-6\n",
    "\n",
    "# Reseta o intervalo inicial (vamos começar do zero)\n",
    "x0, x1 = 0, 10\n",
    "\n",
    "# Calcula xa/ga e xb/gb\n",
    "g0, g1 = g(x0), g(x1)\n",
    "if g0 < g1:\n",
    "    xa, xb = x0, x1\n",
    "    ga, gb = g0, g1\n",
    "else:\n",
    "    xa, xb = x1, x0\n",
    "    ga, gb = g1, g0\n",
    "\n",
    "# Atualiza o intervalo até atingir o critério de parada\n",
    "n_iter = 0\n",
    "while True:\n",
    "    n_iter += 1\n",
    "    xmed = (xa + xb) / 2\n",
    "    gmed = g(xmed)\n",
    "    if gmed <= 0:\n",
    "        xa, ga = xmed, gmed\n",
    "    else:\n",
    "        xb, gb = xmed, gmed\n",
    "        \n",
    "    if abs(xb - xa) < x_tol or abs(gmed) < y_tol:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Verificamos os resultados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xmed = (xa + xb) / 2\n",
    "print('número de iterações: %s' % n_iter)\n",
    "print('xa, xb: %.7f, %.7f' % (xa, xb))\n",
    "print('meio do intervalo: %.7f' % xmed)\n",
    "print('g(x) no meio do intervalo: %.3e' % g(xmed))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Convertendo em uma função\n",
    "\n",
    "Seria bom reutilizarmos a lógica do método da bisseção com qualquer função arbitrária, não? Isto é fácil fazer em Python. Vamos criar uma função que recebe uma funçao g e um intervalo x0, x1 e a partir disto calcula o zero de g(x) contido neste intervalo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def bissect(g, x0, x1, x_tol=1e-6, y_tol=1e-6):\n",
    "    \"\"\"\n",
    "    Calcula o zero de g(x) dentro do intervalo (x0, x1).\n",
    "    \n",
    "    Args:\n",
    "        g: uma função de uma única variável\n",
    "        x0, x1: intervalo inicial para a busca do zero de g(x)\n",
    "        x_tol: tolerância em x (retorna quando intervalo for menor que x_tol)\n",
    "        y_tol: tolerância em y (retorna quando |g(x)| < y_tol)\n",
    "    \n",
    "    Retuns:\n",
    "        Retorna um zero de g(x) (valor de x em que g(x) = 0).\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calcula xa/ga e xb/gb\n",
    "    g0, g1 = g(x0), g(x1)\n",
    "    if g0 < g1:\n",
    "        xa, xb = x0, x1\n",
    "        ga, gb = g0, g1\n",
    "    else:\n",
    "        xa, xb = x1, x0\n",
    "        ga, gb = g1, g0\n",
    "\n",
    "    # Atualiza o intervalo até atingir o critério de parada\n",
    "    n_iter = 0\n",
    "    while True:\n",
    "        n_iter += 1\n",
    "        xmed = (xa + xb) / 2\n",
    "        gmed = g(xmed)\n",
    "        if gmed < 0:\n",
    "            xa, ga = xmed, gmed\n",
    "        elif gmed > 0:\n",
    "            xb, gb = xmed, gmed\n",
    "        else:\n",
    "            return xmed\n",
    "\n",
    "        if abs(xb - xa) < x_tol or abs(gmed) < y_tol:\n",
    "            break\n",
    "    if abs(ga) < abs(gb):\n",
    "        return xa\n",
    "    else:\n",
    "        return xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Agora usamos a função"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Encontra o zero de g(x)\n",
    "bissect(g, 0, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Acha o zero do cos(x) no intervalo de 0 a 4 (pi / 2)\n",
    "bissect(np.cos, 0, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Encontra o ponto onde cos(x) é igual a 0.5 no intervalo de 0 a 4\n",
    "def g2(x):\n",
    "    return np.cos(x) - 0.5\n",
    "\n",
    "bissect(g2, 0, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Encontra o ponto onde cos(x) é igual a sin(x) no intervalo de 0 a 4\n",
    "# (sabemos que a resposta é pi/4)\n",
    "def g3(x):\n",
    "    return np.cos(x) - np.sin(x)\n",
    "\n",
    "bissect(g3, 0, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Método da posição falsa\n",
    "\n",
    "O método da posição falsa é um outro método para encontrar o zero de uma função e foi desenvolvido no Egito antigo cerca de 1.650 a.C!\n",
    "\n",
    "Assim como o método da bisseção, a idéia é começar com um intervalo inicial (xa, xb) e reduzí-lo gradualmente até atingir uma tolerância pré-estabelecida. A única diferença com relação ao método da bisseção é que o valor\n",
    "escolhido para a divisão do intervalo não é mais o centro. Na realidade o valor é calculado aproximando a função\n",
    "por uma reta e escolhendo o ponto onde esta reta atinge o zero.\n",
    "\n",
    "Esta estratégia de escolha pode ser mais eficaz e muitas vezes faz com que o método da posição falsa seja bem mais eficiente que o método da bisseção. Se a função for muito próxima de uma reta, a convergência se dá em apenas uns poucos passos. Mostramos o resultado em um gráfico:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x0, x1 = 0, 10\n",
    "g0, g1 = g(x0), g(x1)\n",
    "\n",
    "# A expressão abaixo fornece o zero da reta\n",
    "x_falsa = (x0 * g1 - x1 * g0) / (g1 - g0)\n",
    "X = np.linspace(0, 10)\n",
    "Y = g(X)\n",
    "\n",
    "plt.plot(X, Y, 'b-', label='função')\n",
    "plt.plot([x0, x1], [g0, g1], 'kx-', label='aproximação linear')\n",
    "plt.plot([x_falsa, x_falsa], [0, g(x_falsa)], 'ko--')\n",
    "plt.legend(loc='best')   # mostra legenda com os 'labels' dados\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Faremos algumas iterações manuais do algoritmo para ver como se dá a convergência. Note que como o novo valor \n",
    "x_falsa é menor que zero, substituímos xa no intervalo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xa, xb = x0, x1\n",
    "ga, gb = g0, g1\n",
    "\n",
    "for _ in range(5):\n",
    "    x_falsa = (gb * xa - ga * xb) / (gb - ga)\n",
    "    g_falsa = g(x_falsa)\n",
    "    if g_falsa <= 0:\n",
    "        xa, ga = x_falsa, g_falsa\n",
    "    else:\n",
    "        xb, gb, = x_falsa, g_falsa\n",
    "\n",
    "    # Fazemos o gráfico\n",
    "    x_falsa = (gb * xa - ga * xb) / (gb - ga)\n",
    "    plt.plot(X, Y, 'b-', label='função')\n",
    "    plt.plot([xa, xb], [ga, gb], 'kx-', label='aproximação linear')\n",
    "    plt.plot([x_falsa, x_falsa], [0, g(x_falsa)], 'ko--')\n",
    "    plt.legend(loc='best')   # mostra legenda com os 'labels' dados\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "    print('xa: %.7f, ga: %.2e' % (xa, ga))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que apesar da escolha inicial de novo intervalo ser boa, o método da posição falsa gera uma convergência muito lenta para certas classes de funções (infelizmente isto inclui a nossa $g(x)$ considerada). Isto ocorre porque um dos pontos extremos se mantêm fixo em cada iteração do método, enquanto o outro converge lentamente para a posição correta.\n",
    "\n",
    "Este comportamento torna o método da posição falsa muito mais ineficiente que a bisseção para algumas funções. Como nem sempre sabemos se nossas funções vão se comportar bem ou mal com este método, é necessário criar uma estratégia para contornar este problema. Uma estratégia comum é utilizar o algoritmo de Illinois, que atualiza a posição da nova divisão com uma regra um pouco diferente da posição falsa:\n",
    "\n",
    "$$x' = \\frac{\\frac12 g_b x_a - g_a x_b}{\\frac12 g_b - g_a}$$\n",
    "\n",
    "(note a presença dos fatores de 1/2)\n",
    "\n",
    "Esta regra deve ser utilizada apenas se houver duas atualizações seguidas com o mesmo ponto.\n",
    "\n",
    "Quando fazemos isto, o resultado melhora significativamente:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xa, xb = 0, 10\n",
    "ga, gb = g(xa), g(xb)\n",
    "ultimo_ponto = None\n",
    "illinois = False\n",
    "\n",
    "for _ in range(5):\n",
    "    if illinois:\n",
    "        x_falsa = (0.5 * gb * xa - ga * xb) / (0.5 * gb - ga)\n",
    "        illinois = False\n",
    "    else:\n",
    "        x_falsa = (gb * xa - ga * xb) / (gb - ga)\n",
    "    g_falsa = g(x_falsa)\n",
    "    \n",
    "    if g_falsa <= 0:\n",
    "        xa, ga = x_falsa, g_falsa\n",
    "        illinois = ultimo_ponto == 0\n",
    "        ultimo_ponto = 0\n",
    "    else:\n",
    "        xb, gb, = x_falsa, g_falsa\n",
    "        illinois = ultimo_ponto == 1\n",
    "        ultimo_ponto = 1\n",
    "        \n",
    "    # Fazemos o gráfico\n",
    "    x_falsa = (gb * xa - ga * xb) / (gb - ga)\n",
    "    plt.plot(X, Y, 'b-', label='função')\n",
    "    plt.plot([xa, xb], [ga, gb], 'kx-', label='aproximação linear')\n",
    "    plt.plot([x_falsa, x_falsa], [0, g(x_falsa)], 'ko--')\n",
    "    plt.legend(loc='best')   # mostra legenda com os 'labels' dados\n",
    "    plt.show()\n",
    "    print('xa: %.7f, ga: %.2e' % (xa, ga))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Finalmente, criamos a nossa função para encontrar o zero pela posição falsa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def pos_falsa(g, x0, x1, x_tol=1e-6, y_tol=1e-6, illinois=True):\n",
    "    \"\"\"\n",
    "    Calcula o zero de g(x) dentro do intervalo (x0, x1) utilizando o método da posição falsa.\n",
    "    \n",
    "    Args:\n",
    "        g: uma função de uma única variável\n",
    "        x0, x1: intervalo inicial para a busca do zero de g(x)\n",
    "        x_tol: tolerância em x (retorna quando intervalo for menor que x_tol)\n",
    "        y_tol: tolerância em y (retorna quando |g(x)| < y_tol)\n",
    "        illinois: se verdadeiro, usa correção de Illinois\n",
    "    \n",
    "    Retuns:\n",
    "        Retorna um zero de g(x) (valor de x em que g(x) = 0).\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calcula xa/ga e xb/gb\n",
    "    g0, g1 = g(x0), g(x1)\n",
    "    if g0 < g1:\n",
    "        xa, xb = x0, x1\n",
    "        ga, gb = g0, g1\n",
    "    else:\n",
    "        xa, xb = x1, x0\n",
    "        ga, gb = g1, g0\n",
    "\n",
    "    # Atualiza o intervalo até atingir o critério de parada\n",
    "    n_iter = 0\n",
    "    use_illinois = False\n",
    "    ultimo_ponto = -1\n",
    "    \n",
    "    while True:\n",
    "        n_iter += 1\n",
    "        \n",
    "        if use_illinois:\n",
    "            x_falsa = (0.5 * gb * xa - ga * xb) / (0.5 * gb - ga)\n",
    "            use_illinois = False\n",
    "        else:\n",
    "            x_falsa = (gb * xa - ga * xb) / (gb - ga)\n",
    "        g_falsa = g(x_falsa)\n",
    "\n",
    "        if g_falsa < 0:\n",
    "            xa, ga = x_falsa, g_falsa\n",
    "            use_illinois = illinois and ultimo_ponto == 0\n",
    "            ultimo_ponto = 0\n",
    "        elif g_falsa > 0:\n",
    "            xb, gb, = x_falsa, g_falsa\n",
    "            use_illinois = illinois and ultimo_ponto == 1\n",
    "            ultimo_ponto = 1\n",
    "        else:\n",
    "            return x_falsa\n",
    "\n",
    "        if abs(xb - xa) < x_tol or abs(g_falsa) < y_tol:\n",
    "            break\n",
    "    \n",
    "    # Resultado\n",
    "    if abs(ga) < abs(gb):\n",
    "        return xa\n",
    "    else:\n",
    "        return xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Testamos a função"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x1 = pos_falsa(g, 0, 10)\n",
    "x2 = bissect(g, 0, 10)\n",
    "print('posição falsa: x=%.7f, g=%.2e' % (x1, g(x1)))\n",
    "print('bisseção:      x=%.7f, g=%.2e' % (x2, g(x2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "pos_falsa(np.cos, 0, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Vemos que ambos os métodos retornam valores parecidos dentro da faixa de tolerância aceitável. Queremos descobrir qual dos dois métodos é o mais rápido:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 pos_falsa(g, 0, 10)\n",
    "%timeit -n1000 pos_falsa(g, 0, 10, illinois=False)\n",
    "%timeit -n1000 bissect(g, 0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Estes valores mudam com diferentes tolerâncias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 pos_falsa(g, 0, 10, x_tol=1e-10, y_tol=1e-10)\n",
    "%timeit -n1000 pos_falsa(g, 0, 10, illinois=False, x_tol=1e-10, y_tol=1e-10)\n",
    "%timeit -n1000 bissect(g, 0, 10, x_tol=1e-10, y_tol=1e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "E com diferentes funções ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 pos_falsa(np.cos, 0, 4)\n",
    "%timeit -n1000 pos_falsa(np.cos, 0, 4, illinois=False)\n",
    "%timeit -n1000 bissect(np.cos, 0, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Método da secante\n",
    "\n",
    "O método da secante utiliza uma filosofia um pouco diferente dos métodos anteriores. Ao invés de analisar a convergência de um intervalo que contêm a raiz, o método da secante monitora um único ponto que (sob condições favoráveis) deve convergir para a raiz. \n",
    "\n",
    "Assim como os métodos anteriores, começamos com dois pontos que definem um intervalo (mas estes não precisam conter a raiz). A regra de iteração é a mesma do método do ponto falso, mas ao invés de escolher qual dos dois\n",
    "valore xa e xb iremos escolher, fazemos sempre a substituição:\n",
    "\n",
    "$$\\begin{align} x_b' & = \\frac{x_a g_b - x_b g_a}{g_b - g_a} \\\\ x_a' & = x_b \\end{align}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "X = np.linspace(0, 10)\n",
    "Y = g(X)\n",
    "xa, xb = 0, 10\n",
    "ga, gb = g(xa), g(xb)\n",
    "\n",
    "for i in range(6):\n",
    "    # Propriedades da reta aproximadora para desenhar no gráfico\n",
    "    m = (gb - ga) / (xb - xa)\n",
    "    Y_reta = ga + m * (X - xa)\n",
    "    x_falsa = xa - ga / m\n",
    "    \n",
    "    plt.plot(X, Y, 'b-', label='função')\n",
    "    plt.plot(X, Y_reta, 'k-', label='aproximação linear')\n",
    "    plt.plot([xa, xb], [ga, gb], 'ko')\n",
    "    plt.plot([x_falsa, x_falsa], [0, g(x_falsa)], 'ro--')\n",
    "    plt.legend(loc='best')\n",
    "    plt.show()\n",
    "    print('%s) g(x): %.2e' % (i, gb))\n",
    "    \n",
    "    # Regra de atualização\n",
    "    xa, xb = xb, (xa * gb - xb * ga) / (gb - ga)\n",
    "    ga, gb = gb, g(xb)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def secante(g, x0, x1, x_tol=1e-6, y_tol=1e-6):\n",
    "    \"\"\"\n",
    "    Calcula o zero de g(x) dentro do intervalo (x0, x1) utilizando o método da \n",
    "    secante. O método não exige que a raiz de g(x) esteja neste intervalo e\n",
    "    nem garante que o resultado estará entre x0 e x1.\n",
    "    \n",
    "    Args:\n",
    "        g: uma função de uma única variável\n",
    "        x0, x1: intervalo inicial para a busca do zero de g(x)\n",
    "        x_tol: tolerância em x (retorna quando intervalo for menor que x_tol)\n",
    "        y_tol: tolerância em y (retorna quando |g(x)| < y_tol)\n",
    "    \n",
    "    Returns:\n",
    "        Retorna um zero de g(x) (valor de x em que g(x) = 0).\n",
    "    \"\"\"\n",
    "    \n",
    "    xa, xb = x1, x2\n",
    "    ga, gb = g(x1), g(x2)\n",
    "    \n",
    "    while True:\n",
    "        xa, xb = xb, (xa * gb - xb * ga) / (gb - ga)\n",
    "        ga, gb = gb, g(xb)\n",
    "        if abs(xb - xa) < x_tol or abs(gb) < y_tol:\n",
    "            break\n",
    "    return xb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Funciona como as outras funções"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "secante(g, 0, 10), secante(np.cos, 0, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Fazemos uma comparação de velocidade e vemos que o método da secante é bem mais rápido!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit pos_falsa(g, 0, 10)\n",
    "%timeit bissect(g, 0, 10)\n",
    "%timeit secante(g, 0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Método de Newton\n",
    "\n",
    "O método de Newton é semelhante ao método da secante, mas ao calcular a inclinação da reta a partir de dois pontos consecutivos, utilizamos uma derivada fornecida explicitamente. Desta forma, a regra de atualização fica\n",
    "\n",
    "$$x' = x - \\frac{f(x)}{f'(x)}$$\n",
    "\n",
    "O método de Newton tende a ser extremamente rápido, mas o usuário precisa fornecer a derivada explicitamente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def g_linha(x):\n",
    "    return 2 * x + 2\n",
    "\n",
    "X = np.linspace(0, 10)\n",
    "Y = g(X)\n",
    "x = 5\n",
    "\n",
    "for i in range(6):\n",
    "    # Propriedades da reta aproximadora para desenhar no gráfico\n",
    "    m = g_linha(x)\n",
    "    Y_reta = g(x) + m * (X - x)\n",
    "    x_zero = x - g(x) / m\n",
    "    \n",
    "    plt.plot(X, Y, 'b-', label='função')\n",
    "    plt.plot(X, Y_reta, 'k-', label='aproximação linear')\n",
    "    plt.plot([x], [g(x)], 'ko')\n",
    "    plt.plot([x_zero, x_zero], [0, g(x_zero)], 'ro--')\n",
    "    plt.legend(loc='best')\n",
    "    plt.show()\n",
    "    print('%s) g(x): %.2e' % (i, g(x)))\n",
    "    \n",
    "    # Regra de atualização\n",
    "    x = x - g(x)/g_linha(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def newton(g, g_linha, x0=0, x_tol=1e-6, y_tol=1e-6):\n",
    "    \"\"\"\n",
    "    Calcula o zero da função g(x) com derivada g_linha(x) partindo do valor inicial\n",
    "    x0.\n",
    "    \n",
    "    Args:\n",
    "        g: uma função de uma única variável\n",
    "        g_linha: derivada de g(x)\n",
    "        x0: valor inicial\n",
    "        x_tol: tolerância em x (retorna quando intervalo for menor que x_tol)\n",
    "        y_tol: tolerância em y (retorna quando |g(x)| < y_tol)\n",
    "    \n",
    "    Retuns:\n",
    "        Retorna um zero de g(x) (valor de x em que g(x) = 0).\n",
    "    \"\"\"\n",
    "    \n",
    "    x = x0\n",
    "    \n",
    "    while True:\n",
    "        value = g(x)\n",
    "        diff = g_linha(x)\n",
    "        x = x - value / diff\n",
    "        if abs(value / diff) < x_tol or abs(value) < y_tol:\n",
    "            break\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Funciona como as outras funções"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def menos_seno(x):\n",
    "    return -np.sin(x)\n",
    "\n",
    "newton(g, g_linha, 5), newton(np.cos, menos_seno, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "O método de Newton é normalmente o mais rápido de todos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 pos_falsa(np.cos, 0, 10)\n",
    "%timeit -n1000 bissect(np.cos, 0, 10)\n",
    "%timeit -n1000 secante(np.cos, 0, 10)\n",
    "%timeit -n1000 newton(np.cos, menos_seno, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Mas isto depende um pouco da função..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 pos_falsa(g, 0, 10)\n",
    "%timeit -n1000 bissect(g, 0, 10)\n",
    "%timeit -n1000 secante(g, 0, 10)\n",
    "%timeit -n1000 newton(g, g_linha, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Ponto fixo\n",
    "\n",
    "O problema do ponto fixo consiste em encontrar o valor de x para o qual $f(x) = x$. Este problema geralmente ocorre em sistemas dinâmicos em que $f(x)$ calcula o novo estado do sistema a partir do estado anterior $x$. Por exemplo, se $x$ representa a posição de um sistema em uma animação $f(x)$ representaria a posição no frame posterior. O ponto fixo portanto representaria uma posição estacionária neste sistema: se $f(x) = x$, então o próximo estado é idêntico ao primeiro e assim sucessivamente. É lógico que se reorganizarmos um pouco os termos temos $f(x) - x = 0$, que é um problema de encontrar as raízes (zeros) de uma função.\n",
    "\n",
    "Deste modo, se quisermos encontrar o ponto fixo de uma função (digamos $f(x) = cos(x)$), então podemos utilizar qualquer método de encontrar raízes para a função $g(x) = f(x) - x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Ponto fixo de cos(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def g_fixo(x):\n",
    "    return np.cos(x) - x\n",
    "\n",
    "x_fixo = secante(g_fixo, 0, 2)\n",
    "print('x: %.8f,  f(x): %.8f' % (x_fixo, np.cos(x_fixo)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Ponto fixo como mapa\n",
    "\n",
    "Também podemos encontrar um ponto fixo explicitamente aplicando $f(x)$ sucessivamente nos valores de x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "valores = [1]\n",
    "x = 1\n",
    "for i in range(20):\n",
    "    x = np.cos(x)\n",
    "    valores.append(x)\n",
    "    print('%2d) x: %.7f' % (i, x))\n",
    "    \n",
    "plt.plot(valores)\n",
    "plt.xlabel('número de iterações')\n",
    "plt.ylabel('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Podemos visualizar os dados anteriores em um gráfico com a função"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "points = [(1, 0)]\n",
    "x = 1\n",
    "for i in range(20):\n",
    "    x, x_old = np.cos(x), x\n",
    "    # Ponto na funçao cos(x)\n",
    "    points.append((x_old, x))\n",
    "\n",
    "    # Ponto projetado na reta y = x\n",
    "    points.append((x, x))\n",
    "   \n",
    "# Desenha y = cos(x) e y = x\n",
    "X = np.linspace(0, 2)\n",
    "Y = np.cos(X)\n",
    "plt.plot(X, Y, 'b-', label='cos(x)')\n",
    "plt.plot(X, X, 'k-')\n",
    "\n",
    "# Desenha mapa de iteração\n",
    "X, Y = np.array(points).T\n",
    "plt.plot(X, Y, 'r-', label='x_n')\n",
    "\n",
    "# Define os limites do gráfico [xmin, xmax, ymin, ymax]\n",
    "plt.axis([0.5, 1.1, 0.4, 1.0]) \n",
    "\n",
    "# Elementos do gráfico\n",
    "plt.legend(loc='best')\n",
    "plt.ylabel('cos(x)')\n",
    "plt.xlabel('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Método do ponto fixo\n",
    "\n",
    "Vimos que é possível resolver um problema de ponto fixo como um problema de encontrar raízes. Também é possível fazer o contrário e resolver um problema de raízes como se fosse um ponto fixo. Para isto, considere a que queremos encontrar as raízes tais quais $g(x) = 0$. Com um pouquinho de álgebra, temos $x + k g(x) = x$, o que pode ser reescrito como $h(x) = x$, onde definimos $h(x) = x + k g(x)$ e $k$ é uma constante arbitrária.\n",
    "\n",
    "Deste modo, encontrar a raiz de $g(x)$ é equivalente a encontrar o ponto fixo de $h(x) = x + k g(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def h(x):\n",
    "    return x - 0.1 * g(x)  # definimos k = -0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Desenhamos um gráfico como no caso anterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "x = 5\n",
    "points = []\n",
    "for i in range(20):\n",
    "    x, x_old = h(x), x\n",
    "    \n",
    "    # Ponto na funçao h(x)\n",
    "    points.append((x_old, x))\n",
    "\n",
    "    # Ponto projetado na reta y = x\n",
    "    points.append((x, x))\n",
    "   \n",
    "# Desenha y = g(x) + x e y = x\n",
    "X = np.linspace(0, 10)\n",
    "Y = h(X)\n",
    "plt.plot(X, Y, 'b-', label='x - g(x)')\n",
    "plt.plot(X, X, 'k-')\n",
    "\n",
    "# Desenha mapa de iteração\n",
    "X, Y = np.array(points).T\n",
    "plt.plot(X, Y, 'r-', label='x_n')\n",
    "\n",
    "# Define os limites do gráfico [xmin, xmax, ymin, ymax]\n",
    "plt.axis([5, 6, 5, 6]) \n",
    "\n",
    "# Elementos do gráfico\n",
    "plt.legend(loc='best')\n",
    "plt.ylabel('h(x)')\n",
    "plt.xlabel('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Podemos definir uma função chamada ponto fixo que calcula as raízes como as outras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def ponto_fixo(g, x0, k=1, x_tol=1e-6, y_tol=1e-6):\n",
    "    \"\"\"\n",
    "    Calcula o zero da função g(x) usando o método do ponto fixo e iterando \n",
    "    sobre a função h(x) = x + k * g(x)\n",
    "    \n",
    "    Args:\n",
    "        g: uma função de uma única variável\n",
    "        x0: valor inicial\n",
    "        k: valor usado para criar h(x) = x + k * g(x)\n",
    "        x_tol: tolerância em x (retorna quando intervalo for menor que x_tol)\n",
    "        y_tol: tolerância em y (retorna quando |g(x)| < y_tol)\n",
    "    \n",
    "    Retuns:\n",
    "        Retorna um zero de g(x) (valor de x em que g(x) = 0).\n",
    "    \"\"\"\n",
    "    x = x0\n",
    "    \n",
    "    while True:\n",
    "        value = g(x)\n",
    "        x_new = x + k * value\n",
    "        x, dx = x_new, x_new - x\n",
    "        if abs(dx) < x_tol or abs(value) < y_tol:\n",
    "            break\n",
    "        x = x_new\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Calculamos alguns exemplos como antes..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ponto_fixo(g, 5, k=-0.08), ponto_fixo(np.cos, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Testamos a performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 bissect(g, 0, 10)\n",
    "%timeit -n1000 pos_falsa(g, 0, 10)\n",
    "%timeit -n1000 secante(g, 0, 10)\n",
    "%timeit -n1000 newton(g, g_linha, 5)\n",
    "%timeit -n1000 ponto_fixo(g, 5, k=-0.08)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%timeit -n1000 bissect(np.cos, 0, 4)\n",
    "%timeit -n1000 pos_falsa(np.cos, 0, 4)\n",
    "%timeit -n1000 secante(np.cos, 0, 4)\n",
    "%timeit -n1000 newton(np.cos, menos_seno, 1)\n",
    "%timeit -n1000 ponto_fixo(np.cos, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "O ponto fixo é considerado um método bastante rápido. No entanto, é um dos métodos mais frágeis dos considerados. Em muitos casos não há convergência ou a velocidade de convergência pode ser fortemente dependente\n",
    "do parâmetro k. Como regra geral **nunca use o método do ponto fixo!**. A única exceção é quando é possível provar que a função que você deseja encontrar a raiz possui as propriedades necessárias para que o método funcione.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Cálculo da raiz quadrada\n",
    "\n",
    "Podemos definir a função y = raiz(x) como a solução da equação $x^2 - y = 0$. Deste modo, pensamos em encontrar as raizes para $g(x) = x^2 - y$, com $y$ forncecido. Esta função aceita o método do ponto fixo. Deste modo, criamos a função $h(x) = x + k (x^2 - y)$ e iteramos este mapa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Vamos calcular raiz de 2\n",
    "y = 2\n",
    "x = y   # começamos no final do caminho\n",
    "k = - 1 / (2 * y)\n",
    "for i in range(10):\n",
    "    x = x + k * (x * x - y)\n",
    "    print('%2d) x: %.6f' % (i, x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Método Babilônico\n",
    "\n",
    "Apesar de obter a raiz quadrada do número desejado, o método converge de forma muito lenta. Compare com o método Babilônico, baseado na regra de Newton.\n",
    "\n",
    "Se queremos calcular a raiz de y, definimos $g(x) = x^2 - y = 0$. A derivada de g é $g'(x) = 2x$. Juntando os dois resultados, temos a regra de iteração:\n",
    "\n",
    "$$x' = x - \\frac{x^2 - y}{2x} = \\frac12\\left(x + \\frac yx\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "y = 2\n",
    "x = 1\n",
    "for i in range(20):\n",
    "    x = 0.5 * (x + y / x)\n",
    "    print('%2d) x: %.8f' % (i, x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "A convergência é extremamente rápida e o resultado é muito preciso. Compare com o valor calculado usando o numpy."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}