Métodos Runge-Kutta

$\newcommand{\R}{\mathbb{R}}$

Os métodos de Taylor, são construídos para aproximar a solução do PVI padrão \begin{equation}\label{standard} y'(t) = f(t,y), \quad t>t_0, \qquad y(t_0)=y_0, \end{equation} onde $y:\R\rightarrow \R^m,$ $f:\R\times\R^m\rightarrow \R^m$ e $y_0\in\R^m.$ Isso é feito construindo o polinômio de Taylor que aproxima $y(t)$ para $t$ próximo de $t_n$ e usando-o para estimar $y(t_{n+1}).$ O problema desta estratégia é que derivadas de $y(t_n)$ precisam ser computadas em termos da função $f$ e isso leva ao cálculo de diversas derivadas de $f$ apenas para montar a fórmula de iteração do método.

Para superar essa dificuldade, os métodos de Runge-Kutta trocam as avaliações de múltiplas derivadas de $f$ por diversas avaliações de $f.$ Essa é a forma de incorporar mais informação sobre a função $f,$ o que permite atingir alta acurácia.

Existem várias maneiras de construir métodos de Runge-Kutta. Neste curso, vou exibir uma forma simples, baseada em integração numérica. Integrando a equação diferencial \eqref{standard}, entre $t_n$ e $t_{n+1},$ temos que $$ y(t_{n+1})- y(t_n) = \int_{t_n}^{t_{n+1}} y'(t_n)\: dt = \int_{t_n}^{t_{n+1}} f(t,y(t))\: dt, $$ ou seja \begin{equation}\label{int} y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t_n,y(t))\: dt. \end{equation} Esta equação é exata, porém, sem conhecer $y(t)$ não há como avaliar a integral acima. Diferentes métodos de Runge-Kutta surgem de diferentes abordagens para aproximar esta integral.

O primeiro método de Runge-Kutta sai de aproximar a integral pela valor da área do retângulo de base $h_n$ e altura $f(t_n,y(t_n)$ (figura abaixo).

Com essa aproximação, temos que $$ \int_{t_n}^{t_{n+1}} f(t_n,y(t))\: dt \approx (t_{n+1}-t_n)f(t_n,y(t_n)) \approx h_n f(t_n,y_n), $$ e portanto o método numérico obtido tem iteração $$ y_{n+1} = y_n + h_n f(t_n,y_n). $$ Esse nada mais é que o método de Euler. Portanto, o método de Euler também pode ser entendido como um método tipo Runge-Kutta.

Claro que essa não é a única forma de aproximar a integral. Poderíamos aproximar essa área por um retângulo cuja altura é dada pelo valor de $f$ em $t_{n+1},$ como na figura a seguir.

Neste caso, a fórmula de iteração do método seria $$ y_{n+1} = y_n + h_n f(t_{n+1},y_{n+1}). $$ Assim como o método de Euler, este também é um método de primeira ordem, porém implícito, uma vez que a equação para $y_{n+1}$ depende do próprio valor de $y_{n+1}.$ Esse método é conhecido como método de Euler reverso.

Utilizar métodos implícitos é mais complicado que utilizar métodos explícitos. Perceba que a determinação de $y_{n+1}$ no método de Euler reverso equivale a resolver a equação não linear $F(z) = 0,$ $$ F(z) \equiv z - y_n - h_n f(t_{n+1},z) = 0. $$ Poderíamos considerar utilizar o método de Newton para isso, mas lembre-se que o objetivo inicial era evitar o cálculo de derivadas da função $f.$ Se o método de Newton for empregado, então precisaremos da derivada de $F$ e, portanto, da derivada de $f.$

Uma maneira de converter o método de implícito para explícito é computar aproximações intermediárias e empregá-las na fórmula de iteração. Por exemplo, o método de Euler reverso poderia ser adaptado da seguinte forma $$ \left\{ \begin{array}{rcl} \hat{y} &=& y_n + h_nf(t_n,y_n),\\ y_{n+1} &=& y_n + h_n f(t_{n+1},\hat{y}). \end{array} \right. $$ Neste caso, uma aproximação $\hat{y}$ para $y_{n+1}$ é computada usando o método de Euler e depois empregada para substituir $y_{n+1}$ na fórmula de iteração do método de Euler reverso, tornando-o assim um método explícito.

Todos os métodos vistos até aqui são de primeira ordem. Para conseguir um método melhor, a integral em \eqref{int} deve ser melhor aproximada. Por exemplo, considere a aproximação ilustrada a seguir.

Isto dá origem ao método do trapézio, cuja fórmula de iteração é $$ y_{n+1} = y_n + {h_n\over 2} \left[ f(t_n,y_n) + f(t_{n+1},y_{n_+1})\right]. $$ Novamente, este é um método implícito. Sua versão explícita, conhecida como método de Heun, pode ser obtida da mesma forma e ficaria $$ \left\{ \begin{array}{rcl} \hat{y} &=& y_n + h_nf(t_n,y_n),\\ y_{n+1} &=& y_n + {h_n\over 2} \left[ f(t_n,y_n) + f(t_{n+1},\hat{y})\right]. \end{array} \right. $$

Um último exemplo de método de Runge-Kutta de segunda ordem é o método do ponto médio, obtido quando a integral em \eqref{int} é aproximada pela área de um retângulo cuja altura é o valor $f$ no ponto médio do intervalo $[t_n,t_{n+1}].$ A figura abaixo ilustra esta aproximação.

O método, em sua versão explícita, seria $$ \left\{ \begin{array}{rcl} \hat{y} &=& y_n + {h_n\over 2}f(t_n,y_n),\\ y_{n+1} &=& y_n + h_n f(t_n+{h_n\over 2},\hat{y}). \end{array} \right. $$ Perceba $\hat{y}$ foi obtido com o método de Euler, usando um passo de metade de $h_n.$ Desta forma, $\hat{y}$ é uma aproximação para $y$ em $t_{n+1/2}\equiv t_n+{h_n\over 2}.$

É possível construir métodos de Runge-Kutta de ordens superiores. Como exemplo, exibo aqui um método de Runge-Kutta de quarta ordem, dado por \begin{equation*} \begin{split} Y_1 & = y_n,\\ Y_2 & = y_n +{h_n\over 2}f(t_n,Y_1),\\ Y_3 & = y_n +{h_n\over 2}f(t_{n+1/2},Y_2),\\ Y_4 & = y_n + h_n f(t_{n+1/2},Y_3), \end{split} \end{equation*} $$y_{n+1} = y_n +{h_n\over 6}\Big[f(t_n,Y_1) + 2f(t_{n+1/2},Y_2) +2f(t_{n+1/2},Y_3)+ f(t_n,Y_4)\Big].$$

Para concluir, é importante entender que a qualidade da aproximação está sujeita à ordem do método escolhido e ao tamanho do passo de discretização $h_n.$ Fixado o método, a redução no valor do passo de discretização melhora a qualidade da estimativa $y_n,$ porém não de forma ilimitada. Quanto menor for $h,$ maior a quantidade de passos necessários para cobrir um intervalo de interesse e, portanto, mais erros se acumularão. Além disso, se $h$ for pequeno demais, os erros de computação em precisão finita podem tornar-se significativos. Usar métodos de ordem superior, reduz a necessidade de empregar passos muito pequenos.

Os melhores métodos para PVI realizam controle automático do passo, de modo a garantir que a qualidade desejada para as estimativas é sempre atingida. Esse recurso é indispensável. Por isso, devemos sempre que possível utilizar implementações que tenham controle de passo embutido. No Octave, as funções lsode e ode45, por exemplo, têm essa funcionalidade.

1. Para os problemas de valor inicial abaixo, estime $y(1)$ por um método de RK de ordem 2 usando (i) $h=0.1$ e (ii) $h=0.001.$ Compute o erro absoluto nas duas aproximações, usando a solução exata apresentada.

1. $y'=t^2-y,$ $y(0)=1.$ [ $y(t) = -e^{-t}+t^2-2t+2$ ]
2. $y'=3y+3t,$ $y(0)=1.$ [ $y(t) = {4\over 3}e^{3t}-t-{1\over 3}$ ]

2. Aplique um método de Runge-Kutta de segunda ordem com $h=0.1$ para aproximar o valor solução do PVI abaixo em $t=2.$ $$y'= y/t^2, \quad t>1,\qquad y(1) =1.$$ Usando o comando lsode do Octave, refaça este exercício com $h=10^{-2}.$ Use o resultado do Octave como se fosse o valor real para estimar o erro da aproximação que você obteve com o RK de segunda ordem.

3. Aplique um método de Runge-Kutta de segunda ordem para resolver os problemas de valor inicial abaixo. Compare a solução obtida com a solução real. (Procure variar o método de Runge-Kutta empregado.)

1. $y'=1+(t-y)^2,$ $2\le t\le 3$, $y(2) = 1$, com $h=10^{-1}.$ [ $y(t)=t+(1-t)^{-1}$ ]
2. $y'=1+y/t,$ $1\le t\le 2$, $y(1) = 2$, com $h=10^{-2}.$ [ $y(t)= t\ln(t)+2t$ ]