Integração numérica 1: A regra dos trapézios

Um exemplo simples da integração numérica usando a regra dos trapézios:

A área de um único trapézio entre os pontos $x_j$ and $x_{j+1}$ é dada por: $$ \frac{h}{2} (f(x_j) + f(x_{j+1}). $$

Somando isso sobre $n$ trapezóides em um dado intervalo [a,b] temos: Summing this up over all the trapezoids gives: $$ h\left(\frac 1 2 f(x_0) + f(x_1) + f(x_2) + \cdots + f(x_{n-2}) + \frac 1 2 f(x_{n-1})\right) = h\sum_{j=0}^{n-1} f(x_j) - \frac h 2 \left(f(x_0) + f(x_{n-1})\right) = h\sum_{j=0}^{n-1} f(x_j) - \frac h 2 \left(f(a) + f(b))\right). $$


In [36]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

Primeiro definimos a função para qual desejamos obter o valor da integral:


In [37]:
def f(x):
    return (x-3)*(x-5)*(x-7)+85

Para podermos fazer o gráfico mostrando a função e os trapézio utilizados na integração, amostramos a função entre 0 e 10 utilizando 200 pontos:


In [38]:
x = np.linspace(0, 10, 200)
y = f(x)

defina a região onde se quer integrar e o numero de trapézios a usar


In [54]:
a, b = 1, 8                     # define os limnites de integração

ndiv = 5                        # numero de intervalos a serem utilizados
xint = np.linspace(a,b,ndiv)    # cria um vetor com ndiv elementos entre os limites de integração
yint = f(xint)                  # calcula o valor da função nos pontos do vetor acima

print(xint)
print(yint)


[ 1.    2.75  4.5   6.25  8.  ]
[  37.         82.609375   86.875      81.953125  100.      ]

plotamos a curva e os trapézios usados para a integração


In [55]:
plt.plot(x, y, lw=2,color='r')
plt.axis([a-1, b+1, 0, max([f(a),f(b)])])
plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)
plt.text(0.5 * (a + b), 30,r"$\int_a^b f(x)dx$", horizontalalignment='center', fontsize=20);

plt.plot(xint,yint,'bo-')


Out[55]:
[<matplotlib.lines.Line2D at 0x7f9e047def50>]

Usando as funções nativas do SciPy fazemos a conta usando quadratura e os trapézios


In [56]:
from scipy.integrate import quad, trapz
integral, error = quad(f, a, b)
print("Integral por quadratura:", integral, "+/-", error)
print("Integral usando ", len(xint), "trapezóides:", trapz(yint, xint))


Integral por quadratura: 565.25 +/- 6.27553564669e-12
Integral usando  5 trapezóides: 559.890625

Exercícios

Exercício 1: Escreva uma função em python que calcule a integral pelo método dos trapézios, utilizando um looping "FOR" para calcular a somatória. Com essa função calcule as integrais das funções abaixo:

$$a)~ (x-3)*(x-5)*(x-7)+85 ~~~com~a=1~e~b=9\\ b)~f_2(x) = 1 + x^3 + \sin(kx)~~~com~a=0~e~b=2$$

onde $k$ é um parametro.

  • Faça uma análise dos erros cometidos a medida que se aumenta o numero de intervalos. Para o numero de intervalos utilize a seginte sequencia: [5, 10, 20, 40, 80, 160, 320]. Utilize gráficos para mostrar seus resultados.

  • Tente usar a função SUM para escrever uma outra versão da função de calculo via trapézios.

Sua função deve ser algo como:


In [57]:
def trapezoidal(f, a, b, N):
    """Essa função integra uma função f de a até b usando N intervalos."""
    # seu código aqui.

In [ ]:
# voce deve poder chamar sua função com algo como
trapezoidal(func, 0, 10., 5)