aula-06-newton.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementação do método de Newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neste capítulo, utilizamos uma implementação própria do método de Newton para resolver equações não-lineares unidimensionais. O algoritmo é limitado para a entrada de funções matemáticas.\n",
    "\n",
    "Para ser executado, o método `newton` requer 5 parâmetros: \n",
    "\n",
    "- a estimativa inicial, representada pela variável `x0`;\n",
    "- a função $f(x)$ propriamente dita, representada por `f`;\n",
    "- a derivada $f'(x)$, representada por `df`;\n",
    "- o erro relativo assumido, representado por `tol`;\n",
    "- o número máximo de iterações $N$ para tentativa de solução, representado por `nmax`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import inspect, re, numpy as np\n",
    "from matplotlib.pyplot import plot\n",
    "from prettytable import PrettyTable as pt\n",
    "\n",
    "def newton(x0,f,df,tol,N):\n",
    "    \"\"\"Algoritmo para determinação de raízes pelo método de Newton.\n",
    "\n",
    "    Parâmetros: \n",
    "        x0: estimativa inicial\n",
    "        f: string dependendo de uma variável, i.e., a função-alvo\n",
    "            (e.g., 'x**2 - 1', 'x**2*cos(x)', etc.) \n",
    "        df: string dependendo de uma variável, i.e., a derivada da função-alvo\n",
    "        tol: erro desejado (tolerância)\n",
    "        N: número máximo de iterações a repetir\n",
    "\n",
    "    Retorno: \n",
    "        x: aproximação para a raiz da função    \n",
    "    \"\"\"\n",
    "\n",
    "    # construtor de tabela\n",
    "    table = pt()\n",
    "    \n",
    "    # substitui expressões da string pelas chamadas das funções do numpy\n",
    "    f = re.sub('(sin|sinh|cos|cosh|tan|tanh|exp|log|sqrt|log10|arcsin|arccos|arctan|arcsinh|arccosh|arctanh)', r'np.\\1', f)\n",
    "    df = re.sub('(sin|sinh|cos|cosh|tan|tanh|exp|log|sqrt|log10|arcsin|arccos|arctan|arcsinh|arccosh|arctanh)', r'np.\\1', df)\n",
    "    \n",
    "    # identifica a variável independente em f\n",
    "    var = re.search(r'([a-zA-Z][\\w]*) ?([\\+\\-\\/*]|$|\\))', f).group(1)\n",
    "    \n",
    "    # cria função\n",
    "    f = eval('lambda ' + var + ' :' + f)\n",
    "    \n",
    "    # checa se a função é de uma variável, senão lança erro        \n",
    "    try: \n",
    "        len(inspect.getfullargspec(f).args) - 1 > 0\n",
    "    except:\n",
    "        raise ValueError('O código é válido apenas para uma variável.')\n",
    "    finally:\n",
    "        # cria função derivada\n",
    "        df = eval('lambda ' + var + ' :' + df)\n",
    "    \n",
    "    it = 0 # contador de iteracoes\n",
    "    \n",
    "    # cabeçalho de tabela\n",
    "    table.field_names = ['i','x','f(x)','f\\'(x)','ER']\n",
    "\n",
    "    # imprime estimativa inicial\n",
    "    print(f'Estimativa inicial: x0 = {x0:.6f}\\n')  \n",
    "\n",
    "    # Loop \n",
    "    for i in range(0,N):\n",
    "        \n",
    "        x = x0 - f(x0)/df(x0) # funcao de iteracao \n",
    "        \n",
    "        e = abs(x-x0)/abs(x) # erro\n",
    "        \n",
    "        # tabela\n",
    "        # impressão de tabela\n",
    "        table.add_row([i,np.round(x,8),np.round(f(x),8),np.round(df(x),4),f'{e:e}'])\n",
    "        table.align = 'c';      \n",
    "        \n",
    "        if e < tol:\n",
    "            break\n",
    "        x0 = x                \n",
    "        \n",
    "    print(table)\n",
    "       \n",
    "    if i == N:\n",
    "        print(f'Solução não obtida em {N:d} iterações')\n",
    "    else:\n",
    "        print(f'Solução obtida: x = {x:.6f}')\n",
    "\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exemplo\n",
    "\n",
    "Resolva o problema $f(x) = 0$, para $f(x) = -\\text{arccos}(x) + 4\\text{sen}(x) + 1.7$, no intervalo $-0.2 \\le x \\le 1.0$ e $\\epsilon = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimativa inicial: x0 = -0.100000\n",
      "\n",
      "+----+-------------+-------------+--------+--------------+\n",
      "| i  |      x      |     f(x)    | f'(x)  |      ER      |\n",
      "+----+-------------+-------------+--------+--------------+\n",
      "| 0  |  0.00656144 |  0.16201076 | 3.4999 | 1.624055e+01 |\n",
      "| 1  | -0.03972877 | -0.06940881 | 3.4961 | 1.165156e+00 |\n",
      "| 2  | -0.01987529 |  0.02983116 | 3.499  | 9.989026e-01 |\n",
      "| 3  | -0.02840088 | -0.01278928 | 3.498  | 3.001876e-01 |\n",
      "| 4  | -0.02474469 |  0.00548778 | 3.4985 | 1.477564e-01 |\n",
      "| 5  | -0.02631332 |  -0.0023538 | 3.4983 | 5.961326e-02 |\n",
      "| 6  | -0.02564047 |  0.00100976 | 3.4984 | 2.624162e-02 |\n",
      "| 7  | -0.02592911 | -0.00043314 | 3.4983 | 1.113178e-02 |\n",
      "| 8  | -0.02580529 |  0.00018581 | 3.4983 | 4.798023e-03 |\n",
      "| 9  |  -0.0258584 |  -7.97e-05  | 3.4983 | 2.053976e-03 |\n",
      "| 10 | -0.02583562 |  3.419e-05  | 3.4983 | 8.818630e-04 |\n",
      "+----+-------------+-------------+--------+--------------+\n",
      "Solução obtida: x = -0.025836\n"
     ]
    }
   ],
   "source": [
    "# Chamada da função\n",
    "xm = newton(-0.1,\n",
    "            '-arccos(x) + 4*sin(x) + 1.7',\n",
    "            '4*cos(x) - 1/(1 - x**2)**1/2',\n",
    "            1e-3,\n",
    "            30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exemplo\n",
    "\n",
    "Resolva o problema $h(z) = 0$, para $h(z) = -z(1-2z)^{-1} - \\text{tan}(z+1)$, no intervalo $[1,8]$ e $\\epsilon = 10^{-5}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como no exemplo anterior, para utilizarmos o método de Newton é preciso saber a derivada da função $h(z)$. Vamos encontrá-la utilizando o módulo de computação simbólica `sympy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2*z/(1 - 2*z)**2 - tan(z + 1)**2 - 1 - 1/(1 - 2*z)\n"
     ]
    }
   ],
   "source": [
    "# Importa variável z como símbolo\n",
    "from sympy.abc import z \n",
    "import sympy as sym\n",
    "\n",
    "# Derivada\n",
    "dh = sym.diff(-z/(1 - 2*z) - sym.tan(z+1))\n",
    "\n",
    "# Impressão\n",
    "print(dh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A partir daí, utilizamos a expressão normalmente na função."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimativa inicial: x0 = 5.000000\n",
      "\n",
      "+---+------------+------------+---------+--------------+\n",
      "| i |     x      |    f(x)    |  f'(x)  |      ER      |\n",
      "+---+------------+------------+---------+--------------+\n",
      "| 0 | 4.75884953 | 0.01963205 | -1.3483 | 5.067411e-02 |\n",
      "| 1 | 4.77341064 | 0.00056164 | -1.3262 | 3.050462e-03 |\n",
      "| 2 | 4.77383412 | 1.173e-05  | -1.3256 | 8.870850e-05 |\n",
      "| 3 | 4.77384296 |  2.4e-07   | -1.3256 | 1.852869e-06 |\n",
      "+---+------------+------------+---------+--------------+\n",
      "Solução obtida: x = 4.773843\n"
     ]
    }
   ],
   "source": [
    "zm = newton(5,\n",
    "            'z/(1 - 2*z) - tan(z+1)',\n",
    "            '-2*z/(1 - 2*z)**2 - tan(z + 1)**2 - 1 - 1/(1 - 2*z)',\n",
    "            1e-5,30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare a quantidade de iterações obtidas com o mesmo exemplo resolvido com o algoritmo da bisseção."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tarefa\n",
    "\n",
    "- Generalize o código acima para que a expressão da derivada seja calculada diretamenta por computação simbólica sem intervenção manual."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}