dicas_ep_num.md

Dicas EP1 numérico

Só uma mão para quem nunca usou MATLAB/Scilab/Octave e está perdido com o EP de numérico repentino. O objetivo é ser uma introdução rápida a linguagem.

Escolha de software:

A disciplina pede MATLAB ou Scilab. Pessoalmente tenho mais experiência com o Octave, mas nos três dá para resolver o EP. Rapidamente, o que penso deles:

• MATLAB: O programa comercial que vocês devem ter usado alguma vez, tem os toolboxes e o conhecido Simulink, que permite criação de modelos por interface gráfica e tira um pouco do trabalho de codá-los.
• Scilab: Nunca usei, mas dando uma olhada na internet ele é razoavelmente compatível com o MATLAB e possui também uma interface gráfica para criação de modelos, no estilo MATLAB, apesar de ter menos funções.
• Octave: É bastante compatível com MATLAB na linguagem e sintaxe e código que roda em MATLAB geralmente roda em Octave com poucas modificações. Implementa parte das funções de toolboxes pagas do MATLAB, como funções de otimização. A desvantagem dele é que não possui editor gráfico de modelos.

Tendo em mente o nosso curso, Octave provavelmente não vai servir para controle, mas pode ser bom se vocês quiserem aprender algo mais próximo de MATLAB sem investir grana ou baixar pirata. O código aqui foi testado em Octave e MATLAB.

O que torna esses pacotes atraentes

Eles são conhecidos por serem de fácil manupulação de vetores e matrizes e resolução de problemas matemáticos, tendo funções já implementadas para a maioria das coisas que necessitaremos no curso. Abra um prompt em qualquer um deles e teste um pouco:

a = 2
a + 3
a = [1 2] % vetor linha (aliás, % é usado para comentários)
a + 3
a = (a + [1 2]) * 3
a = [1 ; 2] % vetor coluna, que geralmente é mais usado
M = [0 1 ; 1 1]
(M ^ 20) * [0 ; 1] % 20º e 21º números da seq. de Fibonacci

Vetores e matrizes possuem operatores diretos de adição, subtração, multiplicação, divisão e exponenciação e muitos outros. A maioria das funções que funcionam em reais funcionam também em matrizes. Além disso, a linguagem não é tipada e é muito flexível, então a função:

% Declara soma(a, b) retornando n. Coloque-a em um arquivo soma.m
function n = soma(a, b)
	n = a + b;
end

Funciona para matrizes, vetores e reais.

O EP

Você vai precisar de pelo menos duas funções, uma para realizar a função f = dY/dt e outro para a simulação. Crie o arquivo f.m e coloque nele:

function valor_f = f(t, Y)
%F Coloque aqui descrição.
%   Coloque aqui documentação.
%   O que estiver aqui é considerado documentação e mostrado no comando help.
%   Tente "help f" e veja isso aparecer. O formato padrão é começar com o nome
%   da função em MAIÚSCULAS, uma descrição e depois a documentação tabulada.
	% Parâmetros do seu processo
	alfa = 1;
	beta = 2;
	gama = 3;
	del = 4;
	psi = 5;

	% Apenas um exemplo sem sentido, faça o seu.
	valor_f = [gama * Y(1) * (Y(1) * alfa - beta * Y(2)) ; Y(1) * (del * Y(2) + psi * t)];
end

O arquivo possui documentação e recebe t e Y e calcula f = dY/dt no momento desejado. Note que indexação em MATLAB inicia em 1 e não 0, como geralmente é de costume.

Para o arquivo da simulação, você vai precisar de:

1 - display('mensagem'): Mostra uma mensagem na tela.

2 - A = input('mensagem'): Mostra uma mensagem na tela e espera entrada. A entrada vai ser interpretada pelo programa e pode ser reais, vetores ou matrizes (nessa função o usuário pode entrar 3.14, [1 0 ; 0 1] ou [2 ; 3], por exemplo, e seu programa os recebe corretamente).

3 - Inicializar:

dt = (tf - t0) / iteracoes;
historico_Y = zeros(length(Y), iteracoes); % inicializa vetor vazio com
% length(Y) linhas e iteracoes colunas. Nele salvaremos as amostras do processo
% para plotar posteriormente. É recomendado e mais rápido alocar o vetor com um
% tamanho fixo antes e atribuir durante o loop que dar "append" no loop.
t = 0;

4 - Realizar o processo:

for i = 1:iteracoes % sintaxe para intervalo de 1 até iteracoes
	% você se vira pra fazer o processo aqui.

	% Depois faça:
	historico_Y(:, i) = Y; % salva a amostra. O operador : indica um intervalo
	% (como no for), sendo vazio indica toda a dimensão da matriz. Ou seja,
	% salvaremos nosso vetor (2 linhas, uma coluna, ou 2x1) nas linhas da matriz
	% historico_Y, coluna i. Por acaso, como alocamos corretamente antes,
	% historico_Y possui 2 linhas e iteracoes colunas. Se as dimensões estiverem
	% erradas você recebe um erro "dimension mismatch", que vai se tornar um
	% fenômeno frequente conforme você usa MATLAB.
end

Veja se entende a notação de intervalos usada no for. Para indexação de vetores ou matrizes, há algumas formas alternativas de usá-la:

• (:): Retorna toda a dimensão. Se aplicada em uma matriz, retorna seus elementos como um vetor coluna. Esse é bem importante, então aqui vão alguns exemplos:

A = magic(3) % quadrado mágico de dimensão 3
A(:, 1) %  mostra todas as linhas da primeira coluna, ou seja, a primeira coluna
% inteira da matriz.
A(2, :) % mostra a segunda linha.
A(:) % todos os elementos em uma coluna

Além disso, há a notação start:step:end, que gera sequências com incrementos diferentes de 1:

1:0.5:3 % 1, 1.5, 2, 2.5, 3

Note que a indexação é inclusiva no final.

5 - Plotar o gráfico:

Com o que vimos dá pra matar esse EP. A função plot mostra uma série de pontos, recebendo suas coordenadas Y (e plotando sequêncialmente em X = 0, 1, 2...) ou X e Y. O único truque aqui é chamar hold on se você escolher um modelo que necessita que você plote dois gráficos sobrepostos. Esse comando instrui o MATLAB a, quando receber comandos para plotar dois gráficos na mesma janela, manter ambos ao invés de apagar o antigo e plotar o novo por cima.

hold on;
plot(t0:dt:(tf - dt), historico_Y(1, :)); % plota os valores da 1ª linha de Y de
% t0 até tf - dt, em passos dt.
plot(t0:dt:(tf - dt), historico_Y(2, :)); % mesmo para a segunda linha.

Resultado:

EDIT: Para quem for usar Octave:

Octave suporta operadores de operação e atribuição (+=, *=, ...) e MATLAB não. É um detalhe, mas se for usar Octave rode uma vez em MATLAB o código antes de entregar ou tenha certeza que não está usando nada que MATLAB não suporte.