Métodos Runge-Kutta

Os métodos de Taylor, são construídos para aproximar a solução do PVI padrão $$\begin{array}{}\text{(1)}& y^{\prime }(t)=f(t,y),t>t_0,y(t_0)=y_0,\end{array}$$ onde $y:\mathbb{R}\to {\mathbb{R}}^m,$ $f:\mathbb{R}\times {\mathbb{R}}^m\to {\mathbb{R}}^m$ e $y_0\in {\mathbb{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 $\text{(1)}$, entre $t_n$ e $t_{n+1},$ temos que $$y(t_{n+1})-y(t_n)={\int }_{t_n}^{t_{n+1}}{y}^{\prime }(t_n)dt={\int }_{t_n}^{t_{n+1}}f(t,y(t))dt,$$ ou seja $$\begin{array}{}\text{(2)}& y(t_{n+1})=y(t_n)+{\int }_{t_n}^{t_{n+1}}f(t_n,y(t))dt.\end{array}$$ 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 $$\{\begin{array}{rcl}\hat{y}& =& y_n+h_{n}f(t_n,y_n),\\ y_{n+1}& =& y_n+h_{n}f(t_{n+1},\hat{y}).\end{array}$$ 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 $\text{(2)}$ 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+\frac{h_n}{2}[f(t_n,y_n)+f(t_{n+1},y_{n_{+}1})].$$ 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 $$\{\begin{array}{rcl}\hat{y}& =& y_n+h_{n}f(t_n,y_n),\\ y_{n+1}& =& y_n+\frac{h_n}{2}[f(t_n,y_n)+f(t_{n+1},\hat{y})].\end{array}$$

Um último exemplo de método de Runge-Kutta de segunda ordem é o método do ponto médio, obtido quando a integral em $\text{(2)}$ é 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 $$\{\begin{array}{rcl}\hat{y}& =& y_n+\frac{h_n}{2}f(t_n,y_n),\\ y_{n+1}& =& y_n+h_{n}f(t_n+\frac{h_n}{2},\hat{y}).\end{array}$$ 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+\frac{h_n}{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{array}{rl}{Y}_1& =y_n,\\ Y_2& =y_n+\frac{h_n}{2}f(t_n,Y_1),\\ Y_3& =y_n+\frac{h_n}{2}f(t_{n+1/2},Y_2),\\ Y_4& =y_n+h_{n}f(t_{n+1/2},Y_3),\end{array}$$ $$y_{n+1}=y_n+\frac{h_n}{6}[f(t_n,Y_1)+2f(t_{n+1/2},Y_2)+2f(t_{n+1/2},Y_3)+f(t_n,Y_4)].$$

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^{\prime }=t^2-y,$ $y(0)=1.$ [ $y(t)=-e^{-t}+t^2-2t+2$ ]
2. $y^{\prime }=3y+3t,$ $y(0)=1.$ [ $y(t)=\frac{4}{3}{e}^{3t}-t-\frac{1}{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^{\prime }=y/t^2,t>1,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^{\prime }=1+(t-y{)}^2,$ $2\le t\le 3$, $y(2)=1$, com $h={10}^{-1}.$ [ $y(t)=t+(1-t{)}^{-1}$ ]
2. $y^{\prime }=1+y/t,$ $1\le t\le 2$, $y(1)=2$, com $h={10}^{-2}.$ [ $y(t)=t\ln (t)+2t$ ]