 |
| Janela principal do script. |
Neste script, o usuário escreve a equação diferencial (EDO) em uma caixa de edição de texto. Após clicar no botão 'RUN', os parâmetros de entrada são analisados e se estiver tudo certo, a equação é resolvida numericamente usando o método de Runge-Kutta de 6a. ordem. O resultado é mostrado graficamente em uma nova janela. Depois disso, o usuário tem a opção de editar a EDO e calcular uma nova solução ou 'fechar' (sair) do programa. Esse script apresenta dois diferenciais:
- a equação é dinamicamente analisada;
- é usado o método de Runge-Kutta de 6a. ordem.
Em geral, para resolver uma EDO, é usado o método de RK de 4a. ordem. Exemplos de saída:
Exemplo de saída 1.
Exemplo 2.
 |
| Exemplo de erro, o valor de h deve ser positivo. |
Código Scilab:
// ===============================================
// rk4_gui.sce
//
// Resolve numericamente uma Equacao Diferencial Ordinaria de 1a ordem
// dy/dt = f(t,y)
// pelo metodo de Runge-Kutta de 4a ordem (RK-6), com uma interface
// grafica onde o usuario informa:
// - a equacao diferencial f(t,y) (ex: -2*y + t)
// - o instante inicial t0
// - o valor inicial y0
// - o tempo final de integracao tf
// - o passo de integracao h
//
// Ao clicar em "RUN", a solucao numerica e calculada e plotada dentro
// da propria janela.
//
// COMO USAR:
// 1) Abra este arquivo no SciNotes (editor do Scilab).
// 2) Execute (Executar > Executar sem eco no Scilab, ou tecla F5),
// ou no console digite: exec('rk6_gui.sce', -1)
// 3) A janela sera aberta. Edite os campos e clique em RUN.
//
// EXEMPLOS DE EQUACOES QUE PODEM SER DIGITADAS (use as variaveis t e y):
// -2*y + t -> equacao linear
// y -> crescimento exponencial (dy/dt = y)
// -0.5*y -> decaimento exponencial
// y*(1-y) -> equacao logistica
// -y + sin(t) -> oscilador amortecido forcado
// t^2 - y -> nao-linear simples
//
// OBSERVACAO: o programa resolve EDOs escalares de 1a ordem da forma
// dy/dt = f(t,y). Sistemas de equacoes e EDOs de ordem superior nao
// sao suportados nesta versao.
// =========================================================
// -------------------------------------------------------------------
// Funcao que monta a interface grafica (janela, campos, botoes, eixo)
// -------------------------------------------------------------------
function rk6_gui()
f = figure('figure_name', 'EDO - Método de RK-6a. Dr. Aquino.', ...
'position', [80, 60, 340, 630], ...
'background', [0.95 0.95 0.95], ...
'menubar', 'none', ...
'toolbar', 'none', ...
'resize', 'off');
// ---- Painel esquerdo: entradas do usuario ----
uicontrol(f, 'style', 'frame', 'position', [10, 10, 320, 610]);
uicontrol(f, 'style', 'text', ...
'string', 'Parâmetros do Metodo RK-6a.', ...
'position', [20, 590, 300, 25], ...
'fontsize', 11, 'fontweight', 'bold', ...
'horizontalalignment', 'center');
uicontrol(f, 'style', 'text', ...
'string', 'Prof. Dr. Francisco J. A. de Aquino.', ...
'position', [20, 570, 300, 25], ...
'fontsize', 11, 'fontweight', 'bold', ...
'horizontalalignment', 'center');
uicontrol(f, 'style', 'text', 'string', 'dy/dt = f(t, y) :', ...
'position', [20, 550, 300, 20], 'horizontalalignment', 'left');
edit_eq = uicontrol(f, 'style', 'edit', 'string', '-y/2 + 2*exp(-t/2)*cos(2*t)', ...
'position', [20, 525, 300, 25]);
uicontrol(f, 'style', 'text', ...
'string', '(use t e y; funções: sin, cos, exp, sqrt, log,...)', ...
'position', [20, 503, 300, 18], 'fontsize', 10, ...
'horizontalalignment', 'left');
uicontrol(f, 'style', 'text', 'string', 'Instante inicial t0 :', ...
'position', [20, 465, 170, 20], 'horizontalalignment', 'left');
edit_t0 = uicontrol(f, 'style', 'edit', 'string', '0', ...
'position', [200, 465, 120, 25]);
uicontrol(f, 'style', 'text', 'string', 'Valor inicial y0 :', ...
'position', [20, 425, 170, 20], 'horizontalalignment', 'left');
edit_y0 = uicontrol(f, 'style', 'edit', 'string', '0', ...
'position', [200, 425, 120, 25]);
uicontrol(f, 'style', 'text', 'string', 'Tempo final tf :', ...
'position', [20, 385, 170, 20], 'horizontalalignment', 'left');
edit_tf = uicontrol(f, 'style', 'edit', 'string', '10', ...
'position', [200, 385, 120, 25]);
uicontrol(f, 'style', 'text', 'string', 'Passo de integracao h :', ...
'position', [20, 345, 170, 20], 'horizontalalignment', 'left');
edit_h = uicontrol(f, 'style', 'edit', 'string', '0.01', ...
'position', [200, 345, 120, 25]);
btn_run = uicontrol(f, 'style', 'pushbutton', 'string', 'RUN', ...
'position', [20, 280, 135, 42], ...
'fontsize', 12, 'fontweight', 'bold', ...
'callback', 'rk6_run_callback()');
btn_clear = uicontrol(f, 'style', 'pushbutton', ...
'string', 'Limpar grafico/Sair', ...
'position', [165, 280, 155, 42], ...
'fontsize', 12, 'fontweight', 'bold', ...
'callback', 'close');
txt_status1 = uicontrol(f, 'style', 'text', 'string', '', ...
'position', [20, 45, 300, 50], ...
'horizontalalignment', 'left', 'fontsize', 12, ...
'background', [1 1 1]);
txt_status2 = uicontrol(f, 'style', 'text', 'string', '', ...
'position', [20, 105, 300, 50], ...
'horizontalalignment', 'left', 'fontsize', 12, ...
'background', [1 1 1]);
txt_status3 = uicontrol(f, 'style', 'text', 'string', '', ...
'position', [20, 165, 300, 50], ...
'horizontalalignment', 'left', 'fontsize', 12, ...
'background', [1 1 1]);
txt_status4 = uicontrol(f, 'style', 'text', 'string', '', ...
'position', [20, 225, 300, 50], ...
'horizontalalignment', 'left', 'fontsize', 12, ...
'background', [1 1 1]);
btn_sobre = uicontrol(f, 'style', 'pushbutton', 'string', 'Sobre', ...
'position', [120, 15, 100, 22], ...
'fontweight', 'bold',...
'callback', 'b_sobre');
// Guarda os handles dos controles no userdata da figura,
// para que o callback consiga acessa-los depois
handles = struct('eq', edit_eq, 't0', edit_t0, 'y0', edit_y0, ...
'tf', edit_tf, 'h', edit_h, ...
'status1', txt_status1,...
'status2', txt_status2,...
'status3', txt_status3,...
'status4', txt_status4);
f.userdata = handles;
f.tag = 'RK4_GUI';
endfunction
// ----- Sobre
function b_sobre()
// Buttons labels + "modal" replaces title
msg = ["Este Script em Scilab resolve uma EDO de 1a. ordem usando RK-6."
"A EDO precisa estar em um formato legível."
" "
"Script desenvolvido pelo prof. Dr. Francisco Aquino,"
"Departamento de Telemática - IFCE."];
messagebox(msg, "Sobre este script", "modal", "info", ["OK"])
endfunction
// -------------------------------------------------------
// Callback do botao RUN: le os parametros,
// resolve a EDO com RK-6
// e plota o resultado
// --------------------------------------------------------
function rk6_run_callback()
f = gcf();
handles = f.userdata;
eq_str = handles.eq.string;
t0 = evstr(handles.t0.string);
y0 = evstr(handles.y0.string);
tf = evstr(handles.tf.string);
h = evstr(handles.h.string);
// --- validacao basica das entradas ---
if isempty(t0) | isempty(y0) | isempty(tf) | isempty(h) then
messagebox('Preencha todos os campos numericos.', 'Erro de entrada', 'error');
return;
end
if h <= 0 then
messagebox('O passo h deve ser um numero positivo.', 'Erro de entrada', 'error');
return;
end
if h > 0.5 then
messagebox('O passo h deve ser um numero menor que 0.5.', 'Erro de entrada', 'error');
return;
end
if tf <= t0 then
messagebox('O tempo final tf deve ser maior que t0.', 'Erro de entrada', 'error');
return;
end
// --- monta dinamicamente a funcao f_user(t,y) a partir do texto digitado ---
erro_eq = %f;
try
deff('dydt = f_user(t,y)', 'dydt = ' + eq_str);
teste = f_user(t0, y0); // chamada de teste para forcar avaliacao
if ~isreal(teste) | isnan(teste) then
erro_eq = %t;
end
catch
erro_eq = %t;
end
if erro_eq then
messagebox('Equacao invalida. Use somente as variaveis t e y ' + ...
'(ex: -2*y + t , sin(t)-y^2 , y*(1-y) ).', ...
'Erro na equacao', 'error');
return;
end
// --- metodo de Runge-Kutta de 6a ordem ---
n = floor((tf - t0) / h); // número de pontos
tt = zeros(1, n+1); // vetor de tempo
yy = zeros(1, n+1); // resposta/solução da EDO
tt(1) = t0; // valores iniciais
yy(1) = y0; // valores iniciais
t = t0;
y = y0;
for i = 1:n
K1 = f_user(t, y);
K2 = f_user(t + h/9, y + K1*h/9);
K3 = f_user(t + h/6, y + (K1 + 3*K2)*h/24);
K4 = f_user(t + h/3, y + (K1 - 3*K2 + 4*K3)*h/6);
K5 = f_user(t + h/2, y + (-5*K1 + 27*K2 - 24*K3 + 6*K4)*h/8);
K6 = f_user(t + 2*h/3, y + (221*K1 - 981*K2 + 867*K3 - 102*K4 + K5)*h/9);
K7 = f_user(t + 5*h/6, y + (-183*K1 + 678*K2 - 472*K3 - 66*K4 + 80*K5 + 3*K6)*h/48);
K8 = f_user(t + h, y + (716*K1 - 2079*K2 + 1002*K3 + 834*K4 - 454*K5 - 9*K6 + 72*K7)*h/82);
y = y + h*(41*K1 + 216*K3 + 27*K4 + 272*K5 + 27*K6 + 216*K7 + 41*K8)/840;
t = t + h;
tt(i+1) = t;
yy(i+1) = y;
end
// --- atualiza a caixa de status ---
linha1 = ' Passos calculados: ' + string(n)+ '. ';
linha2 = ' Instante final atingido: t = ' + string(tt($))+ '. ';
linha3 = ' y(tf) = ' + string(yy($));
[ymax, tmax] = max(yy);
tmax = tt(tmax);
linha4 = ' y(max) = '+string(ymax)+', tmax = '+string(tmax);
handles.status2.string = linha3;
handles.status3.string = linha2;
handles.status4.string = linha1;
handles.status1.string = linha4;
// --- plota o resultado em uma nova janela ---
figure;
plot(tt, yy, 'b-','thickness',2);
xlabel('tempo','fontsize',4);
ylabel('y(t)','fontsize',4);
title('dy/dt = ' + eq_str + ' -- (RK-6, h = ' + string(h) + ')',...
'fontsize',4,'color','blue');
xgrid;
// Obtém o handle da figura atual
fig = gcf();
set(fig, 'background', 8);
endfunction
// --------------------------------------------
// Lanca a interface ao executar o script
// --------------------------------------------
clc; close; // limpa o console e fecha alguma janela
rk6_gui();**********************************
****************
******
***
*
Nenhum comentário:
Postar um comentário