Laboratório de quadratura · Métodos Numéricos

#Objetivo

Implemente a regra composta dos trapézios e verifique a convergência em um caso com solução conhecida. Você precisa de Python 3 e de noções de integrais e funções; este exemplo usa apenas a biblioteca padrão.

O notebook é fornecido sem execução e sem saídas salvas. Execute todas as células, de cima para baixo, em Jupyter. Os valores esperados descritos no texto são resultados analíticos, não registros de uma execução.

Implemente o método.
Compare com uma integral exata.
Examine o efeito do refinamento.
Resolva um exercício de verificação.

Consulte a derivação nas notas do curso.

#Método

Com $h = (b - a) / n$ , usamos

T_{n} = h [\frac{f (a) + f (b)}{2} + \sum_{i = 1}^{n - 1} f (a + i h)] .

O parâmetro $n$ é o número de subintervalos, não o número de pontos.

In [ ]

Código

def trapezoid(f, a, b, n):
    if not isinstance(n, int) or isinstance(n, bool) or n < 1:
        raise ValueError("n must be a positive integer")
    h = (b - a) / n
    interior = sum(f(a + i * h) for i in range(1, n))
    return h * (0.5 * f(a) + interior + 0.5 * f(b))

#Comparação com uma solução exata

Para $f (x) = x^{2}$ em $[0, 1]$ , a integral é $1 / 3$ e o erro em aritmética exata é

T_{n} - \frac{1}{3} = \frac{1}{6 n^{2}} .

Por exemplo, com quatro subintervalos esperamos $T_{4} = 0.34375$ . A tolerância abaixo permite diferenças de arredondamento em ponto flutuante.

In [ ]

Código

from math import isclose

for n in [1, 2, 4, 8, 16]:
    approximation = trapezoid(lambda x: x * x, 0.0, 1.0, n)
    error = approximation - 1.0 / 3.0
    expected_error = 1.0 / (6 * n * n)
    assert isclose(error, expected_error, rel_tol=1e-10, abs_tol=1e-15)
    print(f"n={n:2d}  T_n={approximation:.10f}  error={error:.10f}")

#Refinamento

Calcule $E_{n} / E_{2 n}$ . A expressão analítica prevê uma razão igual a quatro; ela é uma propriedade deste exemplo suave, não uma garantia para qualquer integrando.

In [ ]

Código

errors = [
    abs(trapezoid(lambda x: x * x, 0.0, 1.0, n) - 1.0 / 3.0)
    for n in [1, 2, 4, 8, 16]
]
ratios = [coarse / fine for coarse, fine in zip(errors, errors[1:])]
print(ratios)

#Exercício

Troque o integrando por $f (x) = x$ . Explique por que a integral deve ser $1 / 2$ e o método deve ser exato em aritmética exata. Complete a célula abaixo e use isclose para conferir o resultado.

In [ ]

Código

# Complete este exercício antes de executar:
# result = trapezoid(...)
# assert isclose(result, 0.5)

#Limites e extensão

Para $n$ muito grande, o erro de arredondamento pode interferir na análise do erro de discretização. Não conclua a ordem de um método a partir de apenas uma malha.

Como extensão, integre $e^{x}$ em $[0, 1]$ e compare com $e - 1$ , usando exp e expm1 da biblioteca math. Depois volte ao índice de aulas.