Projeto Controlduino: Erros na integração

No artigo introdutório sobre Cálculo falei que a integração numérica é passível de erros. Isso deveria ser óbvio: Temos um \delta x pequeno, mas significativo e, portanto, o cálculo das áreas está longe de tomar, como base, um valor infinitesimal. O resultado são espaços vazios embaixo da curva, quando ela sobe, e áreas extras, quando ela desce:

Mesmo fazendo δx cada vez menos, temos erros...
Mesmo fazendo δx cada vez menos, temos erros…

Observe o triângulo que “sobra” na descida da curva e o triângulo que “falta”, na subida… Mesmo que \delta x seja pequenininho, ainda teremos esses triângulos e, portanto, erros. Esse é o motivo de não usarmos o método retangular para integreação, mas o método trapezoidal… É como se estivessemos somando o triângulo que falta ou sobra na área do retângulo.

Mas, mesmo o método trapezoidal introduz erros na integração numérica:

Mesmo fazendo δx cada vez menor, o erro ainda está lá...
Mesmo fazendo δx cada vez menor, o erro ainda está lá…

Alguém poderia argumentar que o método trapezoidal é suficientemente preciso, dado o \delta x de, digamos, 20 ms… Mas, muita coisa pode acontecer em 20 ms! Temos que compensar possíveis erros causados pela inclinação entre os valores de duas amostras. Como fazer isso?

Voltando ao caso do método retangular,  basta determinar se, nas duas amostras, a curva está subindo ou descendo. Isso nos dará o sinal do erro (se ele deve ser acrescentado ou retirado da área retangular). Depois calculamos a área do triângulo entre as amostras e a “acrescentamos” a área do retângulo:

diferenca = y1 - y0;
area_triangulo = diferença * dx / 2.0;
area_total = (diferenca * dx) - (sinal(diferenca) * area_triangulo)

Isso ai dá, quase exatamente, o método do trapézio. Mas um erro ainda menor continua (como mostrado na animação do método trapezoidal, acima). Só que existe um jeito de medir o “resíduo” deixado no cálculo da área… A dedução da fórmula do erro do método do trapézio, abaixo, é complicada de explicar…

\displaystyle error=-\frac{(b-a)^3}{12N^2}f''(\xi)

Os valores a e b correspondem ao intervalo \delta x = (b - a) e o valor ξ é um valor entre duas amostras (usaremos a média entre os valores y_n e y_{n-1}). O valor N é a quantidade de amostras. Esse ƒ” significa “derivada segunda”, ou seja, é uma derivada de uma derivada da função ƒ.

A ideia geral é que, a derivada te dará a taxa de variação da curva e a derivada segunda te dará a variação da variação. Esse valor será mais ou menos proporcional ao erro no ponto intermediário entre duas amostras (por isso usarei a média). A constante, antes de ƒ” ajusta a taxa de erro ao seu valor mais real… Como eu disse: deduzir essa fórmula é complicado pacas, exige a aplicação de vários outros conceitos de cálculo que não há espaço aqui para explicar.

Há alguns fatos interessantes que podemos deduzir da fórmula acima: A derivada segunda terá o mesmo sinal da variação do sinal de entrada… Se o sinal sobe temos ƒ” ≥ 0, se desce, o contrário. Mas, como vimos, quando a curva sobe, “falta” um pedaço da área e vice-versa, por isso o sinal negativo na frente da fórmula.

Como esse erro será acumulado, tudo o que temos que fazer é usar \delta x, ao invés de (b-a) e N=2 (duas amostras). Mas, temos um problema: Como calcular a derivada segunda? Não temos uma “função”, apenas um conjunto de pontos! Felizmente, podemos guardar a informação da taxa de variação do sinal anterior, assim, de posse de duas taxas de variação temos dois valores de onde calcular a nova taxa, ou seja, a derivada segunda. A fórmula do erro fica:

\displaystyle \begin{array}{rl}  erro &\approx \frac{\delta x^3}{12\cdot2^2}\cdot\frac{\frac{y_n - y_{n-1}}{\delta x} - \frac{y_{n-1} - y_{n-2}}{\delta x}}{\delta x}\\ \\ &\approx \frac{\delta x^3}{48}\cdot\frac{\frac{y_n - 2y_{n-1} + y_{n-2}}{\delta x}}{\delta x} \\ \\ &\approx \frac{\delta x^3}{48}\cdot\frac{y_n - 2y_{n-1} + y_{n-2}}{\delta x^2}\\ \\ &\approx \frac{y_n - 2y_{n-1} + y_{n-2}}{48}\delta x \end{array}

Mesmo que estejamos usando 3 amostras aqui o valor de N continua sendo 2. A terceira amostra só é usada para obtermos a segunda derivada. Assim, nosso processo de integração pelo método trapezoidal acumulará as pequenas áreas sob a curva usando o cálculo do trapézio, ainda exigindo duas amostras, adicionada ao erro acima (que exige 3 amostras). Então, cada passo da integração usará 3 amostras e teremos:

\displaystyle area(y_n \dots y_{n-2}, dx)=\frac{y_n+y_{n-1}}{2}\delta x + \frac{y_n - 2y_{n-1}+y_{n-2}}{48}\delta x

O erro tem que ser subtraído, ao invés de somado (por isso a inversão do sinal), porque quando a curva é ascendente, temos uma perda da área. Quando é descendente, um ganho…

Por que não usar a “regra de Simpson”?

Ao invés de usar o método do trapézio, poderíamos usar a regra de Simpson:

\displaystyle\int_{a}^{b}f(x)dx \approx \frac{b-a}{6} \left[ f(a) + 4f(\frac{a+b}{2}) + f(b) \right]

Acontece que essa regra também adiciona um pequeno erro… A inviabilidade do uso da regra de Simpson, pelo menos nos microcontroladores AVR de 8 bits, é que o erro usa uma derivada quarta:

\displaystyle \frac{1}{90}\left(\frac{b-a}{2}\right)^5\left|f^{(4)}(\xi)\right|

Para calcularmos a derivada quarta precisamos manter as 5 últimas amostras e você já viu o qual complicado por ficar o calculo envolvendo a derivada segunda! Sem contar com os problemas que deríamos com a precisão limitada do tipo float e a enorme quantidade de ciclos gastos com tanta conta pra fazer…

Mesmo assim, se você preferir a interpolação do método de Simpson, fique à vontade. Sugiro que ignore o erro residual e conviva com algum problema na acumulação das pequenas áreas…

Nossas rotinas de derivação e integração:

Assim, posso agora mostrar pra vocês as duas rotinas de derivação e integração… Elas ficam assim:

float deriv(float y1, float y0, float dx)
{
  return (y1 - y0)/dx;
}

float integr(float y2, float y1, float y0, float dx)
{
  // Simplificação da função area(), acima.
  return ((25.0 * y2 - 26.0 * y1 + y0) / 48.0) * dx;
}

Bem direto, huh? Não há muito o que fazer com relação à derivação. O máximo que obteremos é uma secante da curva, entre dois pontos. Sem os pontos intermediários fica difícil (mas não impossível) estimar um erro… No entanto, a secante é uma aproximação muito boa da tangente. E agora temos nossa integração, com os pequeninos erros corrigidos… A correção é importante porque precisaremos acumular esses valores no decorrer do tempo…

A única preocupação que temos que ter é quanto aos valores de y_2, y_1, y_0, e dx. Eles não devem ser diferentes o suficiente para que existam erros de truncamento. Como o tipo float nos dá 7 algarismos de precisão e provavelmente dx será da ordem de 10^{-3}, se usarmos valores percentuais para y_n e, considerando que o degrau entre duas amostras é bem próximo de 0.1% (ou próximo da ordem de 10^{-3}) não deveríamos ter problemas, especialmente com a rotina de integração.

Porém (e sempre tem um porém, não é mesmo?), como veremos mais tarde, não estamos interessados em integrar e derivar um sinal cujos valores sejam os pontos do valor de entrada do ADC… Estmos interessados em fazer isso com o desvio desses valores em relação a um ponto preestabelecido, ou seja, um valor que é a diferença entre o valor real, obtido do ADC e um valor de ajuste (set point). Essa diferença pode fazem com que os valores de y sejam bem pequenos, além da precisão de floats… Teremos que transformar esses valores, escalonando-os antes da operação e invertendo o escalonamento depois dela. Só assim não ocorrerão truncamentos… Lembre-se que não podemos usar tipos como doublelong double no Arduino… O tipo double é mapeado diretamente ao tipo float e estou assumindo que as rotinas da libm sigam à risca, a especificação IEEE-754. Ou seja, o tipo tem precisão de 7 algarismos e pode apresentar erro de ±FLT_EPSILON.

Anúncios

Deixe um comentário

Faça o login usando um destes métodos para comentar:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s