{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fórmulas do Trapézio e de Simpson\n", "# Nelson Kuhl (nmkuhl@usp.br)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fórmula dos $n$-Trapézios\n", "Dada uma função $f:[a,b]\\to\\mathbb{R}$, a aproximação da sua integral em $[a,b]$ pela fórmula dos $n$-Trapézios, $T_n$, é dada por\n", "\n", "$$\\int_a^b f(x)\\,dx \\approx T_n = h\\left[\\frac{1}{2}f(x_0) + \\sum_{j=1}^{n-1}f(x_j) + \\frac{1}{2}f(x_n)\\right],\\quad \\text{com } h=\\frac{b-a}{n}, \\quad x_j = a + jh,\\ 0 \\le j \\le n.$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def n_trapezios(a, b, n, func):\n", " # Aproximação da integral de a até b da função func usando n-Trapézios\n", " h = (b - a)/n # espaçamento\n", " x = np.linspace(a, b, n+1) # malha\n", " y = func(x) # valores da função na malha\n", " return h * (0.5*y[0] + np.sum(y[1:n]) + 0.5*y[n]) # y[1:n] são os valores y[1], ... , y[n-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exemplo\n", "Vamos aproximar\n", "\n", "$$\n", " \\ln2 = \\int_1^2 \\frac{1}{x}\\,dx\n", "$$\n", "\n", "pela fórmula dos $n$-Trapézios usando vários valores de $n$ e\n", "\n", "$$\n", " a = 1.0,\\ b = 2.0\\ \\text{ e }\\ f(x) = \\frac{1}{x}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ln2 = np.log(2)\n", "print('Valor \"exato\" de ln(2):',ln2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = 1.0; b = 2.0 # Extremos do intervalo\n", "def f(x): return 1.0/x # Integrando" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cálculo das aproximações com vários valores de $n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "print(\" n T_n Erro\")\n", "print()\n", "for n in range(1, 365):\n", " ntrap = n_trapezios(a, b, n, f)\n", " erro = abs(ln2 - ntrap)\n", " print(\"{:10d} {:20.16f} {:10.2e}\".format(n, ntrap, erro))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos verificar o comportamento do erro com $h$ da seguinte forma: calcule o erro $E_k$ para $n=2^k$ e veja o comportamento da razão $\\vert E_{k-1}/E_k\\vert$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\" n T_n Erro Razão\")\n", "print()\n", "n = 1; ekm1 = np.nan # ekm1 aramazena E_{k-1}; inicialmente, não está definido\n", "for k in range(1, 22):\n", " ntrap = n_trapezios(a, b, n, f)\n", " ek = abs(ln2 - ntrap) # E_k\n", " razao = abs(ekm1/ek) # E_{k-1}/E_k\n", " print(\"{:10d} {:20.16f} {:10.2e} {:7.2f}\".format(n, ntrap, ek, razao))\n", " n = 2*n\n", " ekm1 = ek" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quando dividimos $h$ por $2$ (ou duplicamos $n$), o erro cai aproximadamente por um fator $4$. Isto é coerente com a teoria e será demosndtrado em aula." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fórmula de $n$-Simpsons\n", "A aproximação da integral pela fórmula de $n$-Simpsons, $S_n$, é dada por\n", "\n", "$$\n", " \\int_a^b f(x)\\,dx \\approx S_n = \\frac{h}{3}\\left[f(x_0) + 4f(x_1) + 2f(x_2) + \\cdots + 2f(x_{2n-2})\n", " + 4f(x_{2n-1}) + f(x_{2n})\\right],\\ \\text{com} \\ h=\\frac{b-a}{2n},\\ x_j = a + jh,\\ 0 \\le j \\le 2n.\n", "$$\n", "\n", "Note que $n$ é o **número de repetições**. Para não termos de lidar com os índices onde aparecem multiplicação por $2$ ou por $4$, podemos usar o seguinte resultado (verifique como exercício):\n", "\n", "$$\n", " S_n = \\frac{4T_{2n} - T_n}{3},\n", "$$\n", "\n", "onde $T_{2n}$ é a fórmula dos $2n$-Trapézios (mesmo $h$ que $S_n$) e $T_n$ é a fórmula dos $n$-Trapézios (que usa o dobro do valor de $h$ usado por $S_n$). Além disso, podemos explorar a identidade (exercício!)\n", "\n", "$$\n", " T_{2n} = \\frac{1}{2}T_n + h\\sum_{j=0}^{n-1} f\\left[a + (2j+1)h\\right],\\quad \\text{com}\\quad h = \\frac{b-a}{2n},\n", "$$\n", "\n", "para evitar a duplicação de cáculos de valores da função $f$ (pense)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def n_simpsons(a, b, n, func): #ATENÇÃO: n é o número de REPETIÇÕES; h = (b-a)/2n\n", " # Aproximação da integral de a até b da função func usando n-Simpsons\n", " h = (b - a)/(2*n) # espaçamento\n", " Tn = n_trapezios(a, b, n, func) # n-Trapézios (usa 2h)\n", " x_new = np.linspace(a+h, b-h, n) # pontos onde é preciso calcular func para T_2n\n", " T2n = 0.5*Tn + h*np.sum(func(x_new)) # T_2n\n", " return (4*T2n - Tn)/3 # n-Simpsons" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vejamos resultados para alguns valorea de $n$ (número de repetições):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\" n S_n Erro\")\n", "print()\n", "for n in range(1, 9):\n", " nsimp = n_simpsons(a, b, n, f)\n", " erro = abs(ln2 - nsimp)\n", " print(\"{:10d} {:20.16f} {:10.2e}\".format(n, nsimp, erro))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Se quisermos comparar os métodos dos Trapézios e de Simpson, é adequado observar os resultados com o mesmo espaçamento, que envolve também a mesma quantidae de avaliações de função. Olhando os resultados tabelados acima, vemos que o erro para $S_8$ é $4.72\\times10^{-7}$ enquanto que o erro para $T_{16}$ é $2.44\\times10^{-4}$. Além disso, o mesmo erro por $n$-Trapézios só é atingido após 364 repetições." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para termos uma ideia do comportamento do erro de $n$-Simpsons, vamos fazer um procedimento parecido com o que foi feito para trapézios: calcule o erro $E_k$ para $n=2^k$ repetições e veja o comportamento da razão $\\vert E_{k-1}/E_k\\vert$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\" n S_n Erro Razão\")\n", "print()\n", "n = 1; ekm1 = np.nan # ekm1 aramazena E_{k-1}; inicialmente, não está definido\n", "for k in range(1, 11):\n", " nsimp = n_simpsons(a, b, n, f)\n", " ek = abs(ln2 - nsimp) # E_k\n", " razao = abs(ekm1/ek) # E_{k-1}/E_k\n", " print(\"{:10d} {:20.16f} {:10.2e} {:7.2f}\".format(n, nsimp, ek, razao))\n", " n = 2*n\n", " ekm1 = ek" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Aparentemente, o erro agora cai por um fator 16, o que é coerente com a teoria (decaimento do erro proporcionalmente a $h^4$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparando fórmulas de erro com erro exato\n", "Fórmulas de erro são extremamente úteis para demonstrações de convergência de métodos numéricos. Na prática, são menos usadas pois são pessimistas e, em geral, difíceis de calcular. Mas podem dar informações úteis quando é possível obter valores numéricos." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fórmula do erro para $n$-Trapézios\n", "Se $f\\in C^2([a,b])$, então o erro $E_{T_n} = I - T_n$ admite a expressão\n", "\n", "$$\n", " E_{T_n} = -\\frac{(b-a)^3}{12n^2}f^{''}(t) = -\\frac{b-a}{12}f^{''}(t)h^2\\ \\left(\\text{pois}\\ h=\\frac{b-a}{n}\\right),\n", "$$\n", "\n", "onde $t$ é um ponto do intervalo $[a,b]$. Como $f^{''}$ é limitada em $[a,b]$, pois é contínua, a expressão acima nos dá as seguintes informações: para *qualquer* $f\\in C^2([a,b])$, a aproximação por $n$-Trapézios converge para a integral e o erro tende a zero proporcionalmente a (pelo menos) $h^2$. Note que podemos afirmar isso sem especificar a função e o valor da integral.\n", "\n", "Como o ponto $t$ é em geral desconhecido, se quisermos obter um valor para a estimativa de erro, podemos usar uma das seguintes desigualdaes equivalentes:\n", "\n", "$$\n", " \\vert E_{T_n}\\vert \\le \\frac{(b-a)^3}{12n^2}M_2\\quad \\text{ou}\\quad \\vert E_{T_n}\\vert \\le \\frac{b-a}{12}M_2h^2,\\quad \\text{onde}\\quad M_2 = \\max_{a\\le x\\le b}\\vert f^{''}(x)\\vert.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exemplo\n", "Para a aproximação de $\\ln 2$ usando $n$-Trapézios na integral $\\int_1^2 \\frac{1}{x}\\,dx$ estime, com a fórmula do erro: **(i)** o erro para $6$-Trapézios; **(ii)** o número mínimo de repetições para se garantir um erro menor do que $10^{-6}$. Compare as respostas com os cálculos feitos acima.\n", "\n", "**Solução**\n", "\n", "Temos $a = 1$, $b = 2$ e $f(x) = \\frac{1}{x}$. Portanto, $f^{''}(x) = \\frac{2}{x^3}$ e $\\vert f^{''}(x)\\vert = \\frac{2}{x^3} \\le 2 = M_2$, $\\forall x\\in[a,b]$. Logo,\n", "\n", "$$\n", " \\vert E_{T_n}\\vert \\le \\frac{1}{6n^2}.\n", "$$\n", "\n", "**(i)** Para $n=6$ obtemos\n", "\n", "$$\n", " \\vert E_{T_6}\\vert \\le \\frac{1}{216} \\approx 4.63 \\times 10^{-3},\n", "$$\n", "\n", "que é cerca de 2.7 vezes maior que o erro exato.\n", "\n", "**(ii)** Basta determinar o menor valor de $n$ tal que\n", "\n", "$$\n", " \\frac{1}{6n^2} < 10^{-6} \\Longrightarrow n > \\frac{1000}{\\sqrt{6}} = 408.248\\dots\n", "$$\n", "\n", "Logo, segundo a estimativa, precisamos de no mínimo $409$ repetições. Na verdade, $251$ repetições bastam." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fórmula do erro para $n$-Simpsons\n", "Agora precisamos de mais regularidade para a função: $f\\in C^4([a,b])$. Nesse caso, o erro $E_{S_n} = I - S_n$ tem a forma\n", "\n", "$$\n", " E_{S_n} = -\\frac{(b-a)^5}{2880\\,n^4}f^{(4)}(t) = -\\frac{b-a}{180}f^{(4)}(t)\\,h^4\\ \\left(\\text{pois}\\ h = \\frac{b-a}{\\boldsymbol{2}\\boldsymbol{n}}\\right),\n", "$$\n", "\n", "onde $t$ é um ponto do inervalo $[a,b]$. Note agora que o erro tende a zero proporcionalmente a (pelo menos) $h^4$. Analogamente ao que foi feito para trapézios, podemos deduzir as estimativas equivalentes\n", "\n", "$$\n", " \\vert E_{S_n}\\vert \\le \\frac{(b-a)^5}{2880\\,n^4}\\quad \\text{ou}\\quad \\vert E_{S_n}\\vert \\le \\frac{b-a}{180}M_4h^4,\\quad \\text{onde} \\quad M_4 = \\max_{a\\le x\\le b}\\vert f^{(4)}(x)\\vert.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exemplo\n", "Para a mesma integral do exemplo anterior, use a fórmula do erro de $n$-Simpsons para estimar: **(i)** o erro para $3$-Simpsons; **(ii)** o número mínimo de *repetições* para se garantir um erro menor do que $10^{-6}$.\n", "\n", "**Solução**\n", "\n", "Novamente temos $a=1$, $b=2$ e $f(x)=\\frac{1}{x}$. Portanto, $f^{(4)}(x) = \\frac{24}{x^5}$ e $\\vert f^{(4)}(x)\\vert = \\frac{24}{x^5} \\le 24 = M_4$, $\\forall x\\in[a,b]$. Logo,\n", "\n", "$$\n", " \\vert E_{S_n}\\vert \\le \\frac{1}{120\\,n^4}.\n", "$$\n", "\n", "**(i)** Para $n=3$ (**repetições**) obtemos\n", "\n", "$$\n", " \\vert E_{S_6}\\vert \\le \\frac{1}{9720} \\approx 1.03 \\times 10^{-4},\n", "$$\n", "\n", "cerca de $70$ vezes maior do que o erro exato.\n", "\n", "**(ii)** Precisamos calcular o menor valor de $n$ tal que\n", "\n", "$$\n", " \\frac{1}{120\\,n^4} < 10^{-6} \\Longrightarrow n > \\sqrt[4]{\\frac{10^6}{120}} \\approx 9.55\\dots\n", "$$\n", "\n", "Segundo esta estimativa, precisamos de $10$ repetições ($h = 0.05$) para garantir um erro menor do que $10^{-6}$. Das simulações, observamos que $8$ repetições ($h=0.125$) bastam." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.12" } }, "nbformat": 4, "nbformat_minor": 2 }