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 $$ 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 $$ \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].$$

Depois de ver este vários métodos, é 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.

Por fim, deixo duas opções de projetos para você se aprofundar no tema de métodos para problemas de valor inicial. Na primeira delas, você verá como podemos utilizar esses métodos em um problema de interesse real, na área de Farmacocinética. No segundo projeto, mais teórico, você verá que os métodos implícitos, como o método de Euler reverso ou método do trapézio, têm seu valor, apesar de serem mais trabalhosos. Escolha um, ou os dois, e siga em frente! Ou, se preferir, pode ir direto para o próximo tópico, e começar a estudar problemas de valor de contorno.

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$ ]