\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