6.1.2 Equações Diferenciais Ordinárias — Métodos de Runge-Kuta — Método de Heun

O método de Euler possui uma velocidade de convergência baixa porque utiliza a derivada do início do intervalo como média para o intervalo todo. Para melhorar a estimativa da inclinação, o método de Heun determina duas derivadas para o intervalo. A primeira derivada é tomada no início do intervalo e equivale àquela usada no método de Euler. A segunda é tomada no ponto final do intervalo. A partir dessas duas derivadas, tomamos a média delas para obter uma estimativa melhorada da inclinação no intervalo.

Do método de Euler, tínhamos que

$$ y’_i = \dfrac{dy}{dx} = f(x_i, y_i) $$

era a inclinação do início do intervalo. Essa mesma equação pode ser usada para extrapolar a solução

$$ y_{i+1}^0 = y_i + f(x_i, y_i) h $$

No método de Euler, essa equação é a usada para solucionar o problema das EDOs. Contudo, no método de Heun, esta estimativa é usada como aproximação inicial do valor, chamado de “predita”. Assim, essa equação é chamada de equação preditora, pois ela fornece uma estimativa de $y\_{i+1}$ que permite o cálculo da inclinação na extremidade final do intervalo. A figura abaixo mostra graficamente o valor predito ($y\_{i+1}^0$) e o corrigido ($y_{i+1}$).

As duas inclinações podem ser combinadas para obter uma inclinação média do intervalo:

$$ \bar{y}’ = \dfrac{y_i’ + y_{i+1}’}{2} = \dfrac{f(x_i, y_i) + f(x_{i+1}, y_{i+1}^0)}{2} $$

Essa inclinação média é então usada para extrapolar a fórmula de Euler:

$$ y_{i+1} = y_i + \dfrac{f(x_i, y_i) + f(x_{i+1}, y_{i+1}^0)}{2} h $$

Esta abordagem de Heun é chamada do tipo preditor-corretor. Todos os métodos de passo múltiplos são dessa categoria. Uma característica interessante dessa abordagem é a capacidade de se acelerar a convergência das aproximações sucessivas. A figura abaixo ilustra a redução do erro de acordo com o número de divisões do intervalo de integração. Note que para o mesmo $n$, o método de Heun decai o erro mais velozmente.

 

1. Implementação

 

def integral_heun(f, xi, xe, y0, n=1e6):
    """
    Numerical integration for solve ODE's using Heun's Method

    integral = integral_heun(f, xi, xe, n=1e6)

    INPUT:
      * f: derivative function f(x, y) = dy/dx
      * xi: beginning of integration interval
      * xe: end of integration interval
      * y0: initial estimative for y
      * n: number of points used

    return: \int_{xi}^{xe} f(x)

    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br

    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/

    Sep 2012
    """
    def predictor(x, y, h):
        return y + f(x, y) * h

    def corrector(x, y, h):
        pred = predictor(x, y, h)
        corr = y + (f(x, y) + f(x + h, pred)) / 2.0 * h
        return corr

    def heun(x, y, h):
        return (x + h, corrector(x, y, h))

    def integrator(x, y, h, n):
        for i in xrange(n):
            (x, y) = heun(x, y, h)
        return y

    (y, x) = (y0, xi)
    h = abs(xe - xi) / n
    return integrator(x, y, h, int(n))

Um exemplo de utilização dessa função pode ser obtido em http://www.sawp.com.br/code/ode/heun.py.

 

2. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001).
[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall (2006).