\[ \frac{dy}{dt} = f(t,y) \]
A solução de Euler (mais simples) é dada por:
\begin{align*} y_{n+1} & = y_n + \Delta t f(t_n,y_n) \\
t_{n+1} & = t_n + \Delta t
\end{align*}
Sendo que $y_0$ deve ser conhecido para um instante $t_0$ e passo $\Delta t$ não pode ser muito grande. Se $\Delta t$ for um valor muito grande o algoritmo terá problemas de estabilidade. A solução pelo método RG4 é dada pelo seguinte algoritmo:
\begin{align*} k_1 &= f(t_n, y_n)\\ k_2 &= f(t_n + \Delta t/2, y_n + k_1 \Delta t/2)\\
k_3 &= f(t_n + \Delta t/2, y_n + k_2 \Delta t /2)\\
k_4 &= f(t_n + \Delta t, y_n + k_3 \Delta t)\\
y_{n+1} &= y_n + \Delta t(k_1 + 2k_2 + 2k_3 + k_4)/6\\
t_{n+1} &= t_n + \Delta t \\
\end{align*}
Esse algoritmo leva a resultados muito melhores que o método de Euler. Já uma equação de 2a. ordem é da forma:
\[ \frac{d^2y}{dt^2} + f(t,y) \frac{dy}{dt} = g(t,y)\]
E sua solução requer o uso de uma variável auxiliar:
\[ z(t,y) = \frac{dy}{dt} \]
Assim, teremos que resolver simultaneamente
\[ \frac{dy}{dt} = z(t,y) \] e
\begin{align*} \frac{dz}{dt} = & g(t,y) - f(t,y) z(t,y) \\
= & h(t,y,z)
\end{align*}
Já uma equação de 3a. ordem é da forma:
\[ \frac{d^3y}{dt^3} + f(t,y) \frac{d^2y}{dt^2} + g(t,y) \frac{dy}{dt} = h(t,y)\]
E sua solução requer o uso de duas variáveis auxiliares:
\[ z(t,y) = \frac{dy}{dt}, \text{ } w(t,y,z) = \frac{dz}{dt}\]
Assim, teremos que resolver simultaneamente
\begin{align*} \frac{dy}{dt} = & z(t,y) \\
\frac{dz}{dt} = & w(t,y,z) \\
\text{e } &\\
\frac{dw}{dt} = & h(t,y) - g(t,y)z - f(t,y)w \\
= & s(t,y,z,w)\\
\end{align*}
Um exemplo. Seja a EDO descrita por:
\[ \frac{d^3y}{dt^3} + \frac{1}{3}\frac{d^2y}{dt^2} + \frac{dy}{dt} = \frac{1}{3}y + \frac{1}{9}exp(-t/3)\cos(t) \]
Com as seguintes condições inicias: $y(0) = 0$, $z = dy/dt = 1$, $w = dz/dt = -2/3$, $t_0 = 0$, $\Delta t = 0.05$. Então
\begin{align*}
\frac{dy}{dt} = & z\\ \frac{dz}{dt} = & w\\
\frac{dw}{dt} = & -w/3 - z + y/3 + (1/9)exp(-t/3)cos(t)
\end{align*}
Logo, o algoritmo necessário é dado por:
\begin{align*}
ky1 = &\Delta t*z; kz1 = w*\Delta t; kw1 = \Delta t*ftyzw(t,y,z,w);\\
ky2 = &\Delta t*(z + kz1/2); kz2 = \Delta t*(w + kw1/2); \\
kw2 = &\Delta t*(ftyzw(t+dt/2,y+ky1/2,z+kz1/2,w+kw1/2));\\
ky3 = &\Delta t*(z + kz2/2); kz3 = \Delta t*(w + kw2/2); \\
kw3 = &\Delta t*(ftyzw(t+dt/2,y+ky2/2,z+kz2/2,w+kw2/2));\\
ky4 = &\Delta t*(z + kz3); kz4 = \Delta t*(w + kw3);
\\ kw4 = &\Delta t*(ftyzw(t+dt,y+ky3,z+kz3,w+kw3));\\
y = &y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6; \\
z = &z + (kz1 + 2*kz2 + 2*kz3 + kz4)/6; \\
w = &w + (kw1 + 2*kw2 + 2*kw3 + kw4)/6;\\
t = &t + \Delta t;
\end{align*}
com $ftyzw(t,y,z,w) = -w/3 - z + y/3 + (1/9)exp(-t/3)cos(t)$. Então, sua solução numérica usando RG4 é dada pelo código Scilab:
ky1 = dt*z; kz1 = w*dt; kw1 = dt*ftyzw(t,y,z,w);
ky2 = dt*(z + kz1/2); kz2 = dt*(w + kw1/2);
kw2 = dt*(ftyzw(t+dt/2,y+ky1/2,z+kz1/2,w+kw1/2));
ky3 = dt*(z + kz2/2); kz3 = dt*(w + kw2/2);
kw3 = dt*(ftyzw(t+dt/2,y+ky2/2,z+kz2/2,w+kw2/2));
ky4 = dt*(z + kz3); kz4 = dt*(w + kw3);
kw4 = dt*(ftyzw(t+dt,y+ky3,z+kz3,w+kw3));
y = y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6;
z = z + (kz1 + 2*kz2 + 2*kz3 + kz4)/6;
w = w + (kw1 + 2*kw2 + 2*kw3 + kw4)/6;
t = t + dt;
onde ftywz é uma função descrita no código.
Erro de Aproximação
Pro sua própria natureza, os métodos numéricos resultam em soluções aproximadas. Um certo grau de erro ocorre em casa iteração e esses tendem a se acumular. Em alguns casos, após várias iterações, a distância da solução real poderá ser tão grande que novos cálculos se tornam inúteis. Os erros cometidos são de duas naturezas: erros de arredondamentos e erros de truncamento.
Erros de Arredondamentos:
Quando usamos um computador ou uma calculadora para aproximar a solução de uma equação diferencial, por apresentarem casas decimais finitas trabalha-se com, erro de arredondamento.
Erros de truncamento:
Erros de truncamentos são causados pelo tipo de técnica empregada para a atualização do valor de $y$. Os erros de truncamento de Euler são maiores que os de Runge-Kutta.
Sobre Runge e Kutta
"Carl David Runge (1856-1927), matemático e físico alemão, trabalhou muitos anos em espectroscopia. A análise de dados o levou a considerar problemas em computação numérica e o método de Runge-Kutta tem origem em seu artigo sobre soluções em 1901 por M. Wilhelm Kutta (1867-1944). Kutta era um matemático alemão que trabalhava com aerodinâmica e é, também, muito conhecido por suas contribuições importantes à teoria clássica de aerofólio." (Valle, 2012)
Bibliografia
BOYCE, E. William, DIPRIMA, C. Richard. Equações diferenciais elementares e problemas de valores de contorno. Sétima edição, editora LTC.
BARROSO, Leonidas Conceição; et al. Calculo numérico (com aplicações). 2. Ed. Minas Gerais, editora HARBRA.
VALLE, Karine Nayara F. Métodos Numéricos de Euler e Runge-Kutta. Monografia (Especialização) - Programa de Pós-graduação em Matemática para Professores com Ênfase em Cálculo da UFMG. 2012.
Código Scilab:
function ydot=ftyzw(t, y, z, w)
ydot = -w/3 - z + y/3 + (1/9)*exp(-t/3).*cos(t);
//ydot = -6*y + cos(4*t)
endfunction
clc;
close;
/////////////////////
// y = exp(-t/3)sin(t)
// y''' + (1/3)y'' + y' = (1/3)y + (1/9)exp(-t/3)cos(t)
// y' = z
// z' = w
// w' = f(t,y,z,w)
////////////////////
y = 0; z = 1; w = -2/3;
t = 0; dt = 0.05;
vy = zeros(1,200);
for k=1:200
vy(k)=y;
ky1 = dt*z; kz1 = w*dt;
kw1 = dt*ftyzw(t,y,z,w);
ky2 = dt*(z + kz1/2);
kz2 = dt*(w + kw1/2);
kw2 = dt*(ftyzw(t+dt/2,y+ky1/2,z+kz1/2,w+kw1/2));
ky3 = dt*(z + kz2/2);
kz3 = dt*(w + kw2/2);
kw3 = dt*(ftyzw(t+dt/2,y+ky2/2,z+kz2/2,w+kw2/2));
ky4 = dt*(z + kz3);
kz4 = dt*(w + kw3);
kw4 = dt*(ftyzw(t+dt,y+ky3,z+kz3,w+kw3));
y = y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6;
z = z + (kz1 + 2*kz2 + 2*kz3 + kz4)/6;
w = w + (kw1 + 2*kw2 + 2*kw3 + kw4)/6;
t = t + dt;
end
ty=0:0.05:(k*dt-0.05);
yy = exp(-ty/3).*sin(ty);
erro3 = abs(yy-vy);
figure; plot(ty,vy,ty,yy,ty,1e6*erro3);
legend('Solução analítica','Solução numérica','|Erro x 10e6|');
title('Método de Runge-Kutta 4 - EDO de 3a. ordem');
xlabel('tempo (s)'); ylabel('Amplitude');
Nenhum comentário:
Postar um comentário