CURSO de NIVELAMENTO
Métodos Computacionais
PARTE III – PROGRAMAÇÃO NO MATLAB
Rogério Pagano (MatLab)
MATLAB

Função
Arquivos de texto
 extensão “.m”
 Começam com a palavra chave “function”
 Têm parâmetros de entrada e saída
 As variáveis internas são locais
 Permitem incorporar “help”
 Chamadas a funções
MATLAB

Exemplo
MATLAB

Exemplo

Load e Save
MATLAB
Save
save X
save arq1 X Y Z
save arq2.sai X Y Z -ascii
save arq3.sai X Y Z -ascii -double
load
load x
load arq1
load arq2.sai
load arq3.sai
Gravando dados num arquivo : FID,FPRINTF,FCLOSE
fid=fopen(‘nome do arquivo.extensão’, ‘modalidade’)
fprintf(fid, formato,dados)
fclose(fid)
MATLAB
Exemplo:
MATLAB

Reserva de memória
Exemplo:
N= 1000;
x = rand(1,N);
% Versão sem reserva de memória
t = cputime;
y = 0;
for n = 1:N,
y = [y x(n)];
end
tsem = cputime-t
clear y
% Versão com reserva de memória
t = cputime;
y = zeros(1,N+1);
for n= 1:N,
y(n+1)= x(n);
end
tcom = cputime-t
% Velocidade relativa
rel= tsem/tcom
MATLAB

Vetorização
Exemplo:
N= 10000;
nciclos= 100;
a= rand(1,N);
b= rand(1,N);
% Cálculo utilizando um ciclo for
t = cputime;
for c = 1:nciclos,
y = 0;
for n = 1:N,
y = y+a(n)*b(n);
end
end
tfor = cputime-t
% Cálculo de forma vectorizada
t = cputime;
for c = 1:nciclos,
y = a*b’;
end
tmat = cputime-t
% Velociade de cálculo relativa
rel = tfor/tmat
MATLAB

Boas práticas em programação

Comentar o código
Preferir funções a scripts
Ler preferencialmente dados de ficheiros
Utilizar preferencialmente funções de alto-nível
Vetorizar operações
Evitar ciclos, utilizar pré-alocação em memória





+ importante do que a velocidade:
código correto e legível
MATLAB

Oscilador de Van der Pol
t0 = 0;
tf = 20;
x0 = [0 0.25];
[t,x] = ode23('volpol', t0, tf, x0);
plot(t,x)
function xdot=volpol(t,x)
xdot=[0 0]
xdot(1)=x(1).*(1- x(2).^2) - x(2);
xdot(2)=x(1);
xdot = xdot’

Runge-kutta
MATLAB
% metodo de Runge Kutta de 4 ordem (explicito)
function [t,x] = RK(t0,tf,x0,dxdt,n)
% atribui o valor de n = 10 como padrao
if nargin < 5, n = 10; end
% armazenado as condicoes iniciais
t(1) = t0;
xs = x0;
x = x0;
% calculo do passo
h = (tf-t0)/n;
function y = f2(t,x)
y = -x.**2./(1+3.*exp(-t).*x);
end
% loop de integracao
for i = 1:n
g1 = h*feval(dxdt,t(i),xs);
g2 = h*feval(dxdt,t(i)+.5*h,xs+.5*g1);
g3 = h*feval(dxdt,t(i)+.5*h,xs+.5*g2);
g4 = h*feval(dxdt,t(i)+h,xs+g3);
xs = xs+(g1+g4)/6+(g2+g3)/3;
t(i+1) = t(i)+h;
x = vertcat(x,xs);
end
MATLAB

Modelo do CSTR
% Definição das constantes do modelo
U =50;
% BTU/(h.ft2.R)
A = 120;
% ft2
DH = -30000;
% BTU/lbm
Ro = 50;
% lb/ft3
Cp = 0.75;
% BTU/(lbm.R)
E = 30000;
% BTU/lbm
R = 1.99;
% BTU/(lbm.R)
k0 = 7.08e10;
% 1/h
V =48;
% ft3
Te = 580;
%R
Th = 550;
%R
Fe = 18;
% ft3/h
Cre = 0.48;
% lbm/ft3
MATLAB

Modelo do CSTR
% Tempo de simulação
t = 0.0 : 0.01 : 10.0;
%h
% Condições iniciais
Cr0 = 0.16;
T0 = 603;
% lbm/ft3
%R
% Simulação do modelo
[t,y] = ode45('dcstr',t,[Cr0 T0],[],[U A DH Ro Cp E R k0 V Te Th Fe Cre]);
MATLAB

Modelo do CSTR
% Visualização da simulação
figure(1);
plot(t,y(:,1));
title('CSTR com Reação Exotérmica');
xlabel('Tempo (h)');
ylabel('Concentração de Reagente (lbm/ft3)');
figure(2);
plot(t,y(:,2));
title('CSTR com Reação Exotérmica');
xlabel('Tempo (h)');
ylabel('Temperatura (R)');
MATLAB

Modelo do CSTR
function dy = dcstr(t,y,flag,par);
U = par(1);
DH = par(3);
Cp = par(5);
R = par(7);
V = par(9);
Th = par(11);
Cre = par(13);
A = par(2);
Ro = par(4);
E = par(6);
k0 = par(8);
Te = par(10);
Fe = par(12);
dy(1) = (Fe/V)*(Cre-y(1)) - k0*exp(-E/(R*y(2)))*y(1);
dy(2) = (Fe/V)*(Te-y(2)) + ((DH*k0*exp(-E/(R*y(2)))*y(1))/(Ro*Cp)) - ...
(U*A*(y(2)-Th)/(V*Ro*Cp));