PARTE II
PROGRAMAÇÃO BÁSICA NO
MATLAB
Além das operações matemáticas “usuais”, o MATLAB possibilita o uso
de operações relacionais e lógicas. A função destes operadores é fornecer
respostas do tipo falso e verdadeiro e ajudar no controle de fluxo em
ambientes de programação. A saída de todas as expressões lógicas
podruz 1 para verdadeiro e 0 para falso.
Os operadores relacionais do MATLAB são:
Os caracteres de strings são colocados entre apóstrofos.
PROGRAMAÇÃO BÁSICA NO
MATLAB
Laços: São construções que possibilitam executar
uma sequência de declarações mais de uma vez.
Veremos duas formas básicas de construção de
laços: laços for e laços while. No laço while a
sequência é repetida até que uma condição seja
atingida e no laço for , a sequência é repetida um
número determinado de vezes.
Laço While: A forma geral do laço while é:
while expressão
...

...bloco de códigos
...
end
O exemplo que se segue é extraído do livro
Programação em MATLAB para Engenheiros:
% PROGRAMA PARA CALCULO DE DESVIO PADRAO
% DE UM CONJUNTO DE DADOS COM
ARBITRARIO DE VALORES DE ENTRADA
NUMERO
% DEFINICAO DAS VARIAVEIS
% N NUMERO DE DADOS DE ENTRADA
%STD_DEV DESVIO PADRAO
% SUMX SOMATORIO DOS VALORES DE ENTRADA
%SUMX2 SOMA DOS QUADRADOS DOS VALORES DE
ENTRADA
% X UM VALOR DE ENTRADA
% X BAR MEDIA DOS VALORES DE ENTRADA
% INICIALIZAÇAO DAS SOMAS
n=0; sumx=0; sumx2=0;
% LEITURA DO PRIMEIRO VALOR
x=input('entre primeiro valor');
% vamos assumir que todas as medidas sao
positivas ou zero
while x>=0
n=n+1;
sumx=sumx+x;
sumx2=sumx2+x^2;
% LEITURA DO PROXIMO VALOR
x=input('entre com o proximo valor');
end
%
VERIFICAÇAO DE EXISTEM NUMERO DE DADOS
SUFICIENTES
if n<2 % INFORMAÇOES INSUFICIENTES
disp('no minimo dois dados devem ser fornecidos');
else
xbar=sumx/n;
stddev=sqrt((n*sumx2-sumx^2)/(n*(n-1)));
fprintf('A
%f\n',xbar);
MEDIA
DESTES
DADOS
VALE
fprintf('O DESVIO PADRAO VALE: %f\n',stddev);
fprintf('O NUMERO DE PONTOS VALE: %f\n',n);
end
:
LAÇO FOR
for índice  expr
declaração1 

.

.


.
corpo
declaraçãon 



end
Onde índice representa a variável do laço e
expr é a expressão de controle do laço.
Para ilustrar vejamos o programa para o cálculo de
um fatorial:
% PROGRAMA PARA CALCULO DO FATORIAL DE UM
NUMERO LIDO VIA TECLADO
% N REPRESENTA O NUMERO
clc;
fat=1;
n=input('entre com o numero');
for i=1:n
fat=fat*i;
end
fprintf('o fatorial de %d vale %d\n',n,fat);
EXPRESSÕES break e continue
Duas declarações adicionais podem ser utilizadas
para contorlar a operação dos laços while e for: as
declarações break e continue.
A declaração break encerra a execução de
um laço e passa o controle para a próxima
declaração logo após o fim do laço e a declaração
continue termina a passagem corrente pelo laço e
retorna o controle para o início do laço.
Um exemplo da declaração break em um laço for
é:
for i=1:5
if i==3
break;
end
fprintf('i =%d\n',i);
end
disp('fim do laço');
A execução deste programa produz:
>> teste
i =1
i =2
fim do laço
A declaração continue executada dentro de um laço
induz que a execução da passagem corrente pelo laço
seja interrompida e que o controle retorne ao início do
laço, sendo que a variável do laço for assumirá o seu
próximo valor, e o laço será executado novamente:
for i=1:5
if i==3
continue;
end
fprintf('i =%d\n',i);
end
disp('fim do laço');
>> teste
i =1
i =2
i =4
i =5
fim do laço
LAÇOS ANINHADOS
Um laço pode estar completamente dentro de outro
laço e são denominados aninhados.
Exemplo:
% METODO NUMERICO
%FACULDADE DE ENGENHARIA QUIMICA DE LORENA
%PROFESSOR OSWALDO LUIZ COBRA GUIMARAES
% [email protected]
% [email protected]
% PROGRAMA PARA METODO DE NR
function newton
xo=input('entre com o valor do de xo');
% ENTRADA DA TOLERANCIA
eps=input('entre com o valor da precisao');
% ENTRADA DA FUNÇAO COMO STRING
y=input('entre com a funçao','s');
der=input('entre com a derivada','s');
fprintf('iteraçao
i=1;
x(n)
x(n+1)
erro\n');
erro=1;
x=xo;
n=input('entre com o numero de termos');
for i=1:n
while erro >eps
f=eval(y);
d=eval(der);
x=xo-f/d;
erro=abs((x-xo)/xo);
fprintf('%4d
\n',i,xo,x,erro);
i=i+1;
xo=x;
end
end
%13.5f
%9.5f
%9.2e
val1 = input('Entre um número : ');%permite a entrada de
valores
%via teclado em modo interativo
val2 = input('Entre outro número : ');
if val1<val2
disp('O primeiro valor é menor que o segundo');
elseif val1>val2
else
end
disp('O segundo valor é maior que o primeiro');
disp('Os valores são iguais');
% Conversão entre unidades centímetros, polegas e pés
fprintf('\n\n');%pula duas linhas
disp('Conversão entre unidades: centímetros, polegadas e pés');
fprintf('\n');
x=input('Entre valor numérico a converter : ');
fprintf('\n\n');
disp('Os seguintes são sistemas válidos de conversão:');
disp('cen-pol ; pol-cen; cen-pes; pes-cen; pol-pes; pes-pol');
fprintf('\n');
sistema = input('Entre sistema de conversão (ex:pol-pes):
','s');
switch sistema
case 'cen-pol'
y=0.393701*x;
disp([num2str(x),' centímetros = ',num2str(y),'
polegadas']);
case 'pol-cen'
y=2.54*x;
disp([num2str(x),' polegadas = ',num2str(y),'
centímetros']);
case 'cen-pes'
y=0.0328084*x;
disp([num2str(x),' centímetros = ',num2str(y),' pés']);
case 'pes-cen'
y=30.48*x;
disp([num2str(x),' pés = ',num2str(y),' centimetros']);
case 'pol-pes'
y=x/12;
disp([num2str(x),' polegadas = ',num2str(y),' pés']);
case 'pes-pol'
y=12*x;
disp([num2str(x),' pés = ',num2str(y),' polegadas']);
otherwise disp('Unidade desconhecida');
end
Encontrar raízes de um polinômio é um problema usual em
Engenharia. Seja o polinômio dado por x5-x4+x-3. No
MATLAB o polinômio é criado:
>> p=[1 -1 0 0 1 -3];
Observer que os coeficientes nulos devem ser incluídos no
vetor p.
A determinação das raízes é realizada utilizando-se a função
roots( ).
>> r=roots(p)
r=
-0.8751 + 0.7925i
-0.8751 - 0.7925i
0.6718 + 1.0386i
0.6718 - 1.0386i
1.4068
A multiplicação entre dois vetores é realizada pela
função conv.
Consideremos os polinômios x5-x4+x-3 e x-1.
>> p=[1 -1 0 0 1 -3];
>> w=[0 0 0 0 1 -1];
>> conv(p,w)
ans =
0
0
0
0
1
-2
1
0
O resultado é x6-2x5+x4+x2-4x+3.
1
-4
3
Em muitas situações, é necessário dividir um
polinômio por outro.
Isto pode ser feito utilizando-se a função
deconv.
>> a=[1 -1 +3];
>> b=[1 1];
>> [q,r]=deconv(a,b)
q=
1
-2
r=
0
0
5
Nos dá o polinômio quociente q e o resto r.
Para a determinação de derivadas de polinômios o
MATLAB apresenta a função polyder.
Seja o polinômio x5-x4+x-3
>> p=[1 -1 0 0 1 -3];
>> polyder(p)
ans =
5
-4
0
0
1
Isto representaria
polinômio 5x4-4x3+1.
os
coeficientes
do
Descrever dados experimentais é de suma importância na
Engenharia. Na verdade, após um experimento, em geral,
desejamos descrever por meio de funções matemáticas o
comportamento do fenômeno, de modo a predizer dados
que não constem de nossos experimentos.
A teoria de ajuste de curvas é vista na disciplina
Métodos Numéricos e neste curso focaremos nossa
atenção na aplicação do Método dentro do ambiente
MATLAB. No MATLAB a função polyfit resolve o
problema do ajuste polinimial de curvas. Vejamos um
exemplo:
>> x=0:.1:1;
>> y=[0 0.11 0.199 0.31 0.4 0.55 0.6 0.7 0.82 0.9 1];
>> n=1; %grau do polinomio interpolador
>> p=polyfit(x,y,n)
p=
1.03
0.0080
A saída da função polyfit é um vetor linha com os
coeficientes do polinômio, do maior grau para o menor grau,
e, desta forma o polinômio de aproximação é dado por
y=1,03x+0,0080.
Agora, façamos uma comparação gráfica entre os
dados reais e o polinômio de aproximação:
>> xi=linspace(0,1,100);
>> z=polyval(p,xi);
>> plot(x,y,xi,z,':')
>> xlabel('x'),ylabel('y=f(x)')
>> title('Ajuste de Curvas de Primeira Ordem')
Ajuste de Curvas de Primeira Ordem
1.4
real
apx
1.2
1
y=f(x)
0.8
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
A interpolação polinomial é utilizada frequentemente, quando
estamos interessados em um dado que não conste da tabela, e
não especificamente na função geradora dos dados tabelados.
A situação gráfica abaixo, ilustra o que acontece quando
dados são interpolados linearmente:
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
1
2
3
4
5
6
7
8
9
10
O MATLAB possui diversas funções de
interpolação.
interp1
interpola dados unidimensionais
interp2
interpola dados bidimensionais
interpn
interpola dados de dimensões maiores
Vejamos um exemplo: Os dados abaixo mostram o
limiar da audição humana, ou seja, o nível mínimo de
som perceptível pela audição humana e que varia com
a frequência do som emitida:
>> hz=[20:10:100 200:100:1000 1500 2000:1000:10000];
%frequencia em hertz
>> nps=[76 66 59 54 49 46 43 40 38 22 ...
14 9 6 3.5 2.5 1.4 0.7 0 -1 -3 ...
-8 -7 -2 2 7 9 11 12]; %em decibeis;>> semilogx(hz,nps,'o');>> xlabel('Frequencia, hz')
>> ylabel('Nivel de press~ao, dB')
>> grid on
>> s=interp1(hz,nps,2.5e3) %iterpolaçao linear
>> s=interp1(hz,nps,2.5e3,'cubic') %iterpolaçao linear
s = -6.0488
>> s=interp1(hz,nps,2.5e3,'spline')
s = -5.8690
>>
A interpolação bidmensional é utilizada quando
temos z=f(x,y). Vamos ilustrar co a seguinte
situação: uma empresa está mapeando o solo
oceânico e utiliza um radar para isso. Mapeia uma
área com espaçamento de 0,5 km, numa malha
retangular.
>> x=0:.5:4;
>> y=0:.5:6;
DIGITE.
>> z=[99 87 99 86 66 88 99 89 99
88 88 99 96 78 88 87 88 100
88 88 99 96 66 88 87 88 85
88 88 84 94 66 88 87 88 56
99 66 84 94 66 88 87 99 95
88 88 99 96 78 88 87 88 100
99 66 84 94 66 88 87 99 95
88 88 99 96 78 88 87 88 100
88 88 99 96 78 88 87 88 100
99 66 84 94 66 88 87 99 95
88 88 99 96 78 88 87 88 100
99 66 84 94 66 88 87 99 77
99 55 84 94 66 88 87 99 77]
>> zi=interp2(x,y,z,3.5,5.1)
zi = 90.2000
>> zi=interp2(x,y,z,3.5,5.1,'spline')
zi = 88.2251
>> zi=interp2(x,y,z,3.5,5.1,'cubic')
zi = 88.9680
Em muitas aplicações é de interesse
encontrar o mínimo ou o máximo de uma função em
um dado intervalo.
Como exemplo, vamos utilizar a função
f(x)=2e-xsin(x), no intervalo [0,8].
>> f='2*exp(-x).*sin(x)'; %entrada da funçao como uma
string;
>> fplot(f,[0,8])
>> title(f), xlabel('x')
2*exp(-x).*sin(x)
0.7
0.6
0.5
De acordo com a figura,existe
um máximo próximo a x=0,7
e um mínimo próximo a 4.
0.4
0.3
0.2
0.1
0
-0.1
0
1
2
3
4
x
5
6
7
8
Para determinar o ponto de mínimo da função
utilizamos:
>> f='2*exp(-x).*sin(x)'; %entrada da funçao como
uma string;
>> xmin=fminbnd(f,2,5)
xmin=fminbnd(f,2,5)
xmin =
3.9270
>> x=xmin;
>> ymin=eval(f)
ymin = -0.0279
>> fx='-2*exp(-x).*sin(x)'; %entrada da funçao como
uma string;>> xmax=fminbnd(f,0,3);
xmax = 0.7854
>> ymax=-eval(f)
ymax = 0.6448
Uma vez que o máximo de f(x) é o mínimo de –f(x),
efetuamos a troca de sinais na definição das funções.
Para a determinação dos zeros de uma função
pode ser utillizada a função fzero, que procura o zero
de uma função unidimensional.
O MATLAB possui a função trapz( ) que aproxima a
curva por trapézios e a função quad( ) que aproxima a
função a ser integrada pro segmentos de parábola,
equivalendo ao Método de Simpson.
a)
Utilizando o Método dos Trapézios
>>
x=0:.1:1;
>>
y=exp(x);
>>
area=trapz(x,y)
area =1.7197
a)
Utilizando o Método de Simpson
Dedemos primeiro criar um arquivo m com a função a
ser integrada. Desta forma, camos criar um arquivo
denominado, por exemplo, int.m.
function y=f(x)
y=exp(x)
Na área de comando do MATLAB, digitamos:
area=quad(‘int’,0,1), para uma integração de x=0 a x=1.
Além da integração unidimensional, o MATLAB possibilita
a integração bidimensional com a utilização da função
dblquad. Para usar dblquad é necessário primeiro criar
uma função que calcule f(x).
Vamos a um exemplo:
1 2
e
x y
dydx
0 0
Inicialmente vamos criar um arquivo m e nomeá-lo
como intdupla.m:
function z=fun(x,y)
z=exp(x+y);
A função pode ser representada graficamente:
>> x=linspace(0,1,30);
>> y=linspace(0,2,40);
>> [X,Y]=meshgrid(x,y);
>> Z=intdupla(X,Y);
>> mesh(X,Y,Z)
25
20
15
10
5
0
2
1.5
1
0.8
1
0.6
0.4
0.5
0
0.2
0
O valor aproximado da integral é:
integral=dblquad('intdupla',0,1,0,2)
integral =10.9782
A
sintaxe
da
função
dblquad
DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX)
é
Crie um
arquivo m
x = -pi/2:.1:pi/2; for c=1:20,
y = sin(2*x+c*pi/10);
axis([-pi/2 pi/2 -1 1]);
movie(M,20,10);
plot(x,y);
M(c) = getframe; end;
% A última linha de comando significa que os gráficos
armazenados em são mostrados 10 vezes a uma taxa de
repetição de 20 figuras por segundo.
Crie um
arquivo m
x = -pi/2:.1:pi/2;
y = -pi/2:.1:pi/2;
[X,Y] = meshgrid(x,y);
for c=1:20,
Z = sin(2*X+c*pi/10)+1.5*cos(2*Y+c*pi/10);
surf(X,Y,Z);
M(c) = getframe;
end;
movie(M,20);