6.1.4 Equações Diferenciais Ordinárias — Métodos de Runge-Kuta — Métodos de Runge-Kuta de Segunda Ordem

1. Introdução

Os métodos de Runge-Kuta possuem a seguinte fórmula geral

$$ y_{i+1} = y_i + \phi (x_i, y_i, h) h $$

onde $\phi(x\_i, y\_i, h)$ é a função incremento, que assim como no método de Euler, representa a inclinação no i-ésimo intervalo. O objetivo da utilização dessa função é aproximar numericamente uma série de Taylor sem calcular as derivadas de ordem superior. Essa função incremento possui forma geral

$$ \phi = a_1 k1 + a_2 k_2 + \cdots + a_n k_n $$

onde os $a\_i$ são constantes a serem determinadas e os termos $k\_i$ são combinações lineares da função diferencial a ser resolvida

$$ f(x,y) = \dfrac{dy}{dx} $$

Assim, os $k$ ‘s são

$$k_1 = f(x_1, y_i) $$
$$k_2 = f(x_i + p_1 h, q_{11} k_1 h) $$
$$k_3 = f(x_i + p_2 h, y_i + q_{21} k_1 h + q_{22} k_2 h) $$
$$\vdots $$
$$k_n = f(x_i + p_{n-1} h, y_i + q_{n-1,1} k_1 h + q_{n-1,2} k_2 h + \cdots + q_{n-1,n-1} k_{n-1} h )$$

onde $q\_{i,j}$ e $p\_i$ são constantes a serem determinadas de forma que obedeçam as relações de recorrência dos $k_i$.

Podemos utilizar vários níveis de aproximações. Isto é, utilizar os $k_i$’s de $1$ à $n$. Além disso, existem infinitos valores para os $p$’s e $q$’s que podem ser utilizados nessas relações de recorrência. Sendo assim, vários métodos de Runge-Kuta podem ser deduzidos usando-se um número diferente de termos na função incremento.

Cada novo $n$ escolhido corresponde à ordem do método de Runge-Kuta. Ou seja, para $n=2$ temos o método de Runge-Kuta de segunda ordem, para $n=3$ temos o de terceira ordem e assim por diante. Das relações de recorrência, podemos observar que se utilizarmos apenas $n=1$, temos o método de Euler, que equivale ao método de Runge-Kuta de primeira ordem.

Uma vez que a ordem do método é escolhida, devemos nos preocupar em encontrar os vetores $a$, $p$ e $q$. Cada elemento desses vetores pode ser obtido ao igualarmos termo a termo a equação

$$ y_{i+1} = y_i + \phi (x_i, y_i, h) h $$

com a expansão em série de Taylor.

 

2. Métodos de Runge-Kuta de Segunda Ordem

Tomando $n=2$ , temos o método de RK de segunda ordem. Assim, com apenas $2$ elementos de $k$ escolhidos, temos que

$$ y_{i+1} = y_i + (a_i k_1 + a_2 k_2)h $$

onde

$$k_1 = f(x_i, y_i) $$

e

$$k_2 = f(x_i + p_1 h, y_i + q_{11} k_1 h)$$

O próximo passo é determinarmos os valores de $a\_1$, $a\_2$, $p\_1$ e $q\_{11}$. Para isso, decompomos $y\_{i+1}$ em uma série de Taylor de segundo grau em termos de $y\_i$ e $f(x\_i, y\_i)$:

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

em que $f'(x\_i, y\_i)$ deve ser determinado pela regra da cadeia

$$ f'(x_i, y_i) = \dfrac{\partial f(x,y)}{\partial x} + \dfrac{\partial f(x,y)}{\partial y} \dfrac{dy}{dx} $$

Utilizando as duas últimas equações acima, temos que

$$ y_{i+1} = y_i + f(x_i, y_i) h + \left(\dfrac{\partial f(x,y)}{\partial x} + \dfrac{\partial f(x,y)}{\partial y} \dfrac{dy}{dx} \right) \dfrac{h^2}{2!} $$

A série de Taylor para uma função de duas variáveis é definida por

$$ g(x+r, y+s) = g(x,y) + r \dfrac{\partial g}{\partial x} + s \dfrac{\partial g}{\partial y} + \cdots $$

Aplicando esse método para expandir $k\_2 = f(x\_i + p\_1 h, y\_i + q\_{11} k\_1 h)$, isto é,

$$ k_2 = f(x_i + p_1 h, y_i + q_{11} k_1 h) = f(x_i, y_i) + p_1 h \dfrac{\partial f}{\partial x} + q_{11} k_1 h \dfrac{\partial f}{\partial y} + O(h^2) $$

Substituindo agora $k\_1$ e $k\_2$ na equação de segunda ordem, temos que

$$ y_{i+1} = y_i + \left( a_1 f(x_i, y_i) + a_2 f(x_i, y_i) \right) h + \left( a_2 p_1 \dfrac{\partial f}{\partial x} + a_2 q_{11} f(x_i, y_i) \dfrac{\partial f}{\partial y}\right) h^2 + O(h^3) $$

onde os termos abaixo devem valer o seguinte

$$a_1 + a_2 = 1 $$, $$a_2 p_1 = \dfrac{1}{2} $$, e $$a_2 q_{11} = \dfrac{1}{2} $$

Como temos três equações com quatro incógnitas, devemos escolher um valor para uma das incógnitas para determinar as outras três. Suponha que especifiquemos um valor para $a\_2$ , com isso podemos determinar $a\_1$ , $p\_1$ e $q\_{11}$. Isto é,

$$a_1 = 1 – a_2 $$, $$p_1 = q_{11} = \dfrac{1}{2 a_2}$$

Note que é possível escolher infinitos valores para $a_2$ , o que determina infinitas fórmulas para resolver o Runge-Kuta de segunda ordem. Segue abaixo as escolhas mais frequentemente usadas.

 

2.1. Método de Heun

Escolhendo-se $a\_2 = \frac{1}{2}$, temos que $a\_1 = \frac{1}{2}$ e $p\_1 = q\_{11} = 1$. Esses parâmetros, quando substituídos na equação do método, nos fornece

$$ y_{i+1} = y_i + \left( \dfrac{1}{2} k_1 + \dfrac{1}{2} k_2 \right) h $$

onde

$$ k_1 = f(x_i, y_i) $$

e

$$ k_2 = f(x_i + h, y_i + k_1 h) $$

Observe que $k\_1$ é a inclinação no início do intervalo e $k\_2$ é a inclinação no final do intervalo. Assim, o Runge-Kuta de segunda ordem com o parâmetro $a_2$ escolhido como $0.5$ equivale ao método de Heun sem iteração.

 

2.2. Método do Ponto Médio

Quando escolhemos $a\_2 = 1$ , temos que $a\_1 = 0$ , $p\_1 = q\_{11} = \frac{1}{2}$ o que nos dá

$$ y_{i+1} = y_i + k_2 h $$

onde

$$ k_1 = f(x_i, y_i) $$

e

$$ k_2 = f\left( x_i + \frac{1}{2} h, y_i + \frac{1}{2} k_1 h \right) $$

que equivale ao método do ponto médio.

 

2.3. Método de Ralston

Ralston e Rabinowitz[1] demonstram que a escolha de $a\_2 = \frac{2}{3}$ fornece um limitante inferior para o erro de truncamento para o método de segunda ordem. Com esse valor, temos que $a\_1 = \frac{1}{3}$, $p\_1 = q\_{11} = \frac{3}{4}$ e temos

$$ y_{i+1} = y_i + \left( \frac{1}{3} k_1 + \frac{2}{3} k_2 \right) h $$

onde

$$ k_1 = f(x_i, y_i) $$

e

$$ k_2 = f\left( x_i + \frac{3}{4} h, y_i + \frac{3}{4} k_1 \right) h $$

 

3. Implementação

 

def integral_runge_kuta_order_2(f, xi, xe, y0, n=1e6, method=ralston):
    """
    Numerical integration for solve ODE's using Second Order Runge-Kuta

    integral = integral_runge_kuta_order_2(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
      * method: method used to solve a unique step (ralston, heun or midpoint)

    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 integrator(x, y, h, n):
        for i in xrange(n):
            (x, y) = method(f, x, y, h)
        return y

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

Onde temos as variantes locais ralston, heun e midpoint, conforme demonstrado na seção anterior. Essas funções são passadas no parâmetro method da função acima implementada. Essas funções são implementadas abaixo.

def ralston(f, x, y, h):
    k1 = f(x, y)
    k2 = f(x + 0.75 * h, y + 0.75 * h * k1)
    y1 = y + (0.333 * k1 + 0.6667 * k2) * h
    x1 = x + h
    return (x1, y1)

def midpoint(f, x, y, h):
    k1 = f(x, y)
    k2 = f(x + 0.5 * h, y + 0.5 * k1 * h)
    y1 = y + k2 * h
    x1 = x + h
    return (x1, y1)

def heun(f, x, y, h):
    k1 = f(x, y)
    k2 = f(x + h, y + k1 * h)
    y1 = y + (0.5 * k1 + 0.5 * k2) * h
    x1 = x + h
    return (x1, y1)

Observe que outras variantes podem ser passadas como parâmetro na função integral_runge_kuta_order_2. Assim, se desejarmos utilizar outras fórmulas decorrentes de outras escolhas de $a_2$, podemos implementá-las de forma que sempre tenham a assinatura de entrada (f, x, y, h) e retorno (x1, y1).

Um exemplo de utilização dessas funções pode ser obtido em http://www.sawp.com.br/code/ode/runge_kuta_order_2.py.

 

4. 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).