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