Parte III: Uma Introdução ao Cálculo Numérico

Integração Numérica

Nesse capítulo veremos um método iterativo que, a partir de uma função conhecida \(f\) e de valores \(a,b \in \mathcal{D}(f)\), encontra um valor aproximado para a área entre a curva de \(f\) e o eixo horizontal limitada no intervalo \([a,b]\), representada por cinza na Figura abaixo.

Dependendo da função \(f\), essa conta pode ser feita de forma precisa a partir de integrais definidas, como vocês aprenderão em cálculo. Mas algumas vezes isso não é possível, ou seja, para algumas funções não é possível encontrar o valor dessa área de forma exata, mesmo usando os recursos de cálculo. Nesses casos a alternativa é encontrar uma aproximação a partir de métodos numéricos, como veremos nesse capítulo.

Antes de seguirmos com a ideia por trás do métodos vale destacar que nesse contexto a área entre a curva de \(f\) e o eixo horizontal, a qual queremos encontrar, considera área acima do eixo como positiva e área abaixo como negativa. Dessa forma, considerando a Figura acima, a área total entre a curva de \(f\) e o eixo, denominada \(A\), é \[ A = A1 - A2, \qquad A_1>0 \ \text{ e } \ A_2>0 \] onde \(A1\) e \(A2\) são os valores de área, ilustrados na mesma figura, ambos positivos.

Aproximação por retângulos e trapézios

Para fazer a aproximação por retângulos primeiro vamos dividirmos o intervalo \([a,b]\) em \(n\) subintervalos de mesmo tamanho. O tamanho de cada subintervalo será \[\delta = \frac{b-a}{n}\] e os subintervalos são definidos pelos segmentos \([x_{i-1}, x_i]\), com \(i=1,\ldots,n\) e \(x_0=a\), \(x_1=a+\delta\), \(x_2=x_1+\delta\), \(\ldots\), \(x_n=b\).

Para cada subintervalo \([x_{i-1},x_i]\) vamos definir um retângulo cuja base é o próprio subintervalo e altura pode ser definida de diferentes formas, serão apresentada 4 maneiras para fazer isso.

Primeiro considere que a altura é o valor de \(f\) avaliado no ponto mais à esquerda desse subintervalo, nesse caso a altura do retângulo de base \([x_{i-1},x_i]\) é \(f(x_{i-1})\) (figura (a) abaixo). Também podemos definir a altura do retângulo como o valor de \(f\) avaliado no ponto mais à direita desse subintervalo, nesse caso a altura do retângulo de base \([x_{i-1},x_i]\) é \(f(x_i)\) (figura (b) abaixo).

Podemos definir também a altura do retângulo de base \([x_{i-1},x_i]\) pelo maior ou pelo menor valor entre \(f(x_{i-1})\) e \(f(x_i)\). Nesse caso a altura do retângulo de base \([x_{i-1},x_i]\) pode ser o \(\min(f(x_{i-1}),f(x_i))\) (figura (c) abaixo) ou pode ser o \(\max(f(x_{i-1}),f(x_i))\) (figura (d) abaixo). A figura abaixo exemplifica esses quatro diferentes possibilidades.

Veja que podemos aproximar a área entre a curva de \(f\) e o eixo horizontal pela soma das áreas dos \(n\) retângulos. Veja também que quanto maior o número retângulos, ou seja, quanto maior \(n\), mais próximo está essa aproximação do real valor da área.

Uma observação interessante é que se definirmos a altura do retângulo de base \([x_{i-1},x_i]\) como o menor valor entre \(f(x_{i-1})\) e \(f(x_i)\) (figura (c) acima), então a soma das áreas dos retângulos será sempre menor que o real valor da área. Dessa forma, dizemos que a soma das áreas desses retângulos, denominada \(\hat{A}_I\), é um limite inferior para o real valor da área \(A\), pois \(\hat{A}_I \le A.\)

Por outro lado, se definirmos a altura do retângulo de base \([x_{i-1},x_i]\) como o maior valor entre \(f(x_{i-1})\) e \(f(x_i)\) (figura (d) acima), então a soma das áreas dos retângulos será sempre maior que o real valor da área. Dessa forma, dizemos que a soma das áreas desses retângulos, denominada \(\hat{A}_S\), é um limite superior para o real valor da área \(A\), pois \(A \le \hat{A}_S.\)

Juntando as duas ideias podemos afirmar que \[ \hat{A}_L \le A \le \hat{A}_U \] Dessa forma, escolher como aproximação para o valor de \(A\) a média entre \(\hat{A}_I\) e \(\hat{A}_S\) nos possibilita controlar o erro da aproximação. Como \(\hat{A}_I \le A \le \hat{A}_S\) isso significa que \(A \in [\hat{A}_I,\hat{A}_S]\) e por isso a distância do ponto médio desse intervalo, que é \(\frac{\hat{A}_S + \hat{A}_I}{2}\), dista de \(A\) no máximo o valor do raio desse intervalo, que é \(\dfrac{\hat{A}_U - \hat{A}_L}{2}\). Ou seja, \[ \left| A - \dfrac{\hat{A}_S + \hat{A}_I}{2} \right| < \dfrac{\hat{A}_S - \hat{A}_I}{2} \]

Também é possível mostrar que a média desses dois valores coincide com a aproximação caso escolhêssemos usar a área dada pelo trapézio de base \([x_{i-1},x_i]\) e alturas \(f(x_{i-1})\) e \(f(x_i))\). Veja a demonstração de que para um subintervalo qualquer de base \([x_{i-1},x_i]\), a média das áreas dos retângulos com aulturas \(\max\{f(x_{i-1}),f(x_i)\}\) e \(\min\{f(x_{i-1}),f(x_i)\}\) é igual a área do trapézio de base \([x_{i-1},x_i]\) e alturas \(f(x_{i-1})\) e \(f(x_i)\).

Sejam \(\hat{A}_{I_i} = \delta \min\left(f(x_{i-1}),f(x_{i})\right)\) e \(\hat{A}_{S_i} = \delta \max\left(f(x_{i-1}),f(x_{i})\right)\) as áreas dos \(i\)-ésimos retângulos abaixo e acima da curva de \(f\), respectivamente. Veja que o valor médio entre essas duas áreas pode ser calculado da seguinte maneira: \[\begin{align*} \dfrac{\hat{A}_{I_i} + \hat{A}_{S_i}}{2} & = \frac{1}{2} \left( \delta \min\left(f(x_{i-1}),f(x_{i})\right) + \delta \max\left(f(x_{i-1}),f(x_{i})\right) \right)\\ & = \frac{\delta}{2} \left( \min\left(f(x_{i-1}),f(x_{i})\right) + \max\left(f(x_{i-1}),f(x_{i})\right) \right)\\ & = \frac{\delta}{2} \left( f(x_{i-1}) + f(x_{i}) \right) \end{align*}\]

O valor \(\frac{\delta}{2} \left( f(x_{i-1}) + f(x_{i}) \right)\) é exatamente a área do trapézio de base \([x_{i-1},x_i]\) e alturas \(f(x_{i-1})\) e \(f(x_i)\).

Veja na figura a seguir que para o mesmo valor de \(n\), no caso do exemplo \(n=8\), a aproximação para \(A\) dada pela área dos trapézios parece ser bem melhor que a aproximação pela somas das áreas dos retângulos, apresentados na figura anterior.

Pseudocódigo

Como já foi discutido vamos usar a área dos trapézios para aproximar o valor de \(A\) Essa área será calculada a partir do valor médio das áreas dos retângulos com altura máxima e mínima. Precisamos então escolher o valor de \(n\), números de retângulos, e calcular a soma das áreas dos retângulos. Mas como determinar o valor de \(n\) que irá fornecer uma aproximação relativamente boa?

Veja que quanto mais subdivisões tem o intervalo \([a,b]\), mais precisa será a aproximação. E para esse caso, temos o controle do erro. Então, passado como entrada um valor de erro \(\varepsilon\) vamos aumentando o número de retângulos, começando com \(n=100\), por exemplo, até que o erro da aproximação seja menor que \(\varepsilon\), isto é, até que \[ \dfrac{\hat{A}_S - \hat{A}_I}{2} < \varepsilon \quad \text{ ou, o que é equivalente, } \quad (\hat{A}_S - \hat{A}_I) < 2\varepsilon \] Quando isso acontecer podemos afirmar que o ponto médio entre \(\hat{A}_S\) e \(\hat{A}_I\), definido por \[ \frac{\hat{A}_S + \hat{A}_I}{2}, \] é uma aproximação para a área \(A\) com erro menor que \(\varepsilon\).

Entrada: a, b, f e e.

Saída: uma aproximação para a área entre a curva de \(f\) e o eixo horizontal, limitada pelo intervalor \([a,b]\), com erro menor que e.

Nome: IntegralNumérica

1. Se o intervalo [a,b] não estiver contido no domínio de f, pare e retorne erro. 
2. Defina n = 100.
3. Defina d = (b-a)/n.
4. Defina A_S = 0 e A_I = 0.
5. Inicie x = a.
6. Faça A_S = A_S + d*max(f(x),f(x + d)).
7. Faça A_I = A_I + d*min(f(x),f(x + d)).
8. Incremente x = x + d.
9. Se x < b, volte para a linha 5.
10. Se (A_S - A_I) < 2*e, retorne (A_S + A_I)/2.
11. Caso contrário faça n = n+10 e volte para a linha 3.

Podemos também fazer de forma recursiva, mas nesse caso o valor de \(n\) tem que ser um parâmetro de entrada.

Entrada: a, b, f, e e n=100.

Saída: uma aproximação para a área entre a curva de \(f\) e o eixo horizontal, limitada pelo intervalor \([a,b]\), com erro menor que e.

Nome: IntegralNuméricaRec

1. Se o intervalo [a,b] não estiver contido no domínio de f, pare e retorne erro. 
2. Defina d = (b-a)/n.
3. Defina A_S = 0 e A_I = 0.
4. Inicie x = a.
5. Faça A_S = A_S + d*max(f(x),f(x + d)).
6. Faça A_I = A_I + d*min(f(x),f(x + d)).
7. Incremente x = x + d.
8. Se x < b, volte para a linha 5.
9. Se (A_S - A_I) < 2*e, retorne (A_S + A_I)/2.
10. Caso contrário retorne IntegralNuméricaRec(a,b,f,e,n+10).