Curso MATLAB 6
Instrutor: Marcelo Escobar
Métodos Númericos para Engenheiros
Métodos Numéricos
Álgebra Matricial
Sistemas Lineares
Sistemas não lineares
Equações Integrais
Equações Diferenciais
Otimização
Manipulação Simbólica
Álgebra Matricial:
Tópicos de Ajuda:
>>help matfun
>>help elmat
>>help sparfun
Multiplicação Matricial: [Produto Interno]
Dadas as Matrizes A e B:
A
*
[n x m ]
B =
[ m x p]
C
[n x p]
>> A*B
Divisão Matricial: [Produto Externo] B= C/A
>>C\A
Conceitos Importantes:
Conceitos Importantes:
Matriz Transposta: B=AT se b(j,i)=A(i,j)
Matriz Identidade: I(i,j)=1 se i==j e I(i,j)=0 se i~=j
Matriz Inversa: se B*A=I, B é a inversa da matriz A
Matriz Singular: se det(A)=0, A é singular
Matriz Simétrica: se A= AT
Diagonal Principal da Matriz : A(i,i) para i=1:n
Matriz Triangular Superior: A(i,j)=0 se i>j
Matriz Triangular Inferior: A(i,j)=0 se i<j
Ortogonalidade de Vetores: se a*b’=0
a[ 1 x n] e b[ 1 x n] a e b são ditos ortogonais.
Sistemas Lineares:
Sistemas Lineares: Forma Geral [ Ax=b ]
Classificação:
Possível e Determinado : se det(A)~=0
Possível e Indeterminado: se det(A)=0 e todos det(A(:,i)=b)=0 i=1:n
Impossível: se det(A)=0 e pelo menos um det(A(:,i)=b)~=0 i=1:n
Posto de uma Matriz: Número de Equações Independentes
>> rank(A)
Valores Característicos: A-λI=A para λ~=0
>>eig(A)
Vetores Característicos: A* (λ*I) = (λ*I) *V
>>[lambda V]=eig(A)
Métodos Diretos:
Métodos de Resolução de Sistemas Lineares:
Forma mais simples no Matlab: x=A\b
Mínimos Quadrados: x=lsqlin(A,b)
Métodos Diretos: ( Principais)
Eliminação Gaussiana:
Fatorização:
a kj

a

 kj
k  1,..., N
a kk

j  k ,..., N  1 
i  1,..., N (  k ) aij  aij  aik a kj


>>help lu
[ Decomposição LU]
>>help qr
[ Decomposição Ortogonal Triangular]
>>help svd [ Decomposição em Valores Singulares]
>>help schur [ Decomposição Schur]
Ex: A = L U
Ly=b
Ux
Exemplo Método de Gauss:
1 2 3   x1  3 
4 5 6. x   6

  2  
7 8 0   x3  9 
A.x  y
Exemplo:
Linha1=linha1/A(1,1)
Linha2=linha2-A(2,1)*linha1
Linha3=linha3-A(3,1)*linha1
Linha2=linha2/A(2,2)
Linha1=linha1-A(1,2)*linha1
Linha3=linha3-A(3,2)*linha3
Linha3=linha3/A(3,3)
Linha1=linha1-A(1,3)*linha1
Linha2=linha2-A(2,3)*linha2
3  3 
1 2
0  3  6    6 



0  6  21  12
1 0  1  1
0 1 2   2 

 
0 0  9  0 
1 0 0  1
0 1 0   2 

 
0 0 1  0 
x1=-1 x2=2 x3=0
Exemplo Método de Crammer:
1 2 3   x1  3 
4 5 6. x   6

  2  
7 8 0   x3  9 
A.x  y
Linha1=b Ax
Linha2=bAy
Linha3=bAz
 3 2 3
Ax  6 5 6
9 8 0
X1=det(Ax)/det(A)
1 3 3 
Ay  4 6 6
7 9 0
X2=det(Ay)/det(A)
1 2 3 
Az   4 5 6
7 8 9 
X3=det(Az)/det(A)
x1=-1 x2=2 x3=0
Métodos Indiretos:
Métodos Indiretos: ( Principais)
Iterações de Jacobi
x k 1  M x k  c , k  0,1,2,...
onde M = D-1 B, c = D-1 b, B = D - A. Sendo D a diagonal da matriz
A. O método escrito para cada elemento do vetor x apresenta a
seguinte forma:
bi 
xik 1

N
 aij x kj
j 1 ( i )
a ii
, i  1,..., N
e k  0,1,2,...
Métodos Indiretos:
Iterações de Gauss-Seidel :
Este método é uma modificação do método de Jacobi, cujo princípio
é de usar os novos valores de x tão logo eles estejam disponíveis.
Neste caso a matriz M = (D - L)-1 U e o vetor c = (D - L)-1 b, onde
D, L e U são as matrizes diagonal, triangular inferior e triangular
superior, respectivamente, extraídas da matriz A = D - L - U. O
método escrito para cada elemento do vetor x apresenta a seguinte
forma:
i 1
xik 1 
bi  
j 1
a ij x kj 1
a ii

N
 aij x kj
j i 1
, i  1,..., N
e k  0,1,2,...
Sistemas Esparsos:
Sistemas Esparsos: vários elementos nulos
>>help issparse [ teste de esparsidade]
>>help sparse [ conversão de matriz cheia para matriz esparsa]
>>help full
[ conversão de matriz esparsa para matriz cheia]
Geração de Matrizes Esparsas:
>>help sprand
[geração de matriz esparsa aleatória]
>>help sparndsym [geração de matriz esparsa simétrica aleatória]
Métodos para Sistemas Esparsos:
>> help pcg Conjugate Gradiente
>> help cgs Conjugate Gradient Squared (CGS)
>> help bicg BiConjugate Gradient (BiCG)
>>help bicgstab BiConjugate Gradient Stabilized (BiCGSTAB)
>>help gmres Generalized Minimum Residual (GMRES)
>>help qmr Quasi-Minimal Residual without lookahead (QMR)
Dicas –Sistemas Lineares:
Sistemas Sub-Determinados:
Numero de Equações (ne) menor que o numero de incógnitas(ni)
>>A\b assume (ni-ne) variáveis nulas
Sistemas Sobre-Determinados:
Numero de Equações (ne) Maior que o numero de incógnitas(ni)
>>A\b utiliza mínimos quadrados para minimizar os resíduos
Residuo=(A*x-b)
Conceito de Norma:
N
M
1
 max  mij
j i 1
N
M

 max  mij
i j 1
M
F

N
N
 m
i 1 j 1
2
ij
>>help norm
A norma é utilizada como critério de parada em loops multivariaveis.
Equações Transcendentais:
Equações sem solução analítica:
Ex: f(x)= x*exp(x/2) qual x / f(x)=0?
Resolução no matlab:
>>help optim
>>help fzero
>>help fsolve
Para utilizar as funções deve-se criar uma função com a equação:
function f=funcao_teste(x)
f=x*exp(x/2)
chute inicial
>>fsolve(’funcao_teste’, 0)
ou
>> fzero(’funcao_teste’, 0)
Sistemas Não Lineares:
Para resolver sistemas de equações:
function f=funcao_teste2(x)
x1=x(1);x2=x(2);
f1=x1*x2-6;
f2=x2+x1-5;
F=[f1;f2];
chute inicial
>>fsolve(’funcao_teste2’, [ 3; 4])
Métodos para Sistemas Não Lineares
Método de Substituição Sucessiva:
x  g( x)
O processo iterativo é aplicado à equação algébrica na forma
modificada da equação , que pode ser obtida por um rearranjo interno
desta equação ou pela simples adição de x em ambos os lados da
igualdade.
x k 1  g( x k )
k  0,1,2,
g(x)
45°
x*
x1
x0
x
Métodos para Sistemas Não Lineares:
Método da Bisseção:
x Rk
k

x

k 1

x R

x k 1


x Rk  x Lk

2
se sign f ( x k )  sign f ( x Rk1 )
k  0,1,2,
f(x)

caso contrá rio
x 0L

x
x Lk  
k 1

x L
k



se sign f ( x )  sign f
k
( x Lk 1 )

caso contrá rio
Os pontos iniciais devem satisfazer a condição:



sign f ( x 0L )   sign f ( x R0 )

onde a função sign(f(x)) fornece o sinal da função f(x).
x * x1
R
x 0R
x
Métodos para Sistemas Não Lineares:
Método de Newton Raphson:
O processo iterativo é aplicado diretamente sobre a equação algébrica
na forma:
x k 1  x k 
f (x k )
k  1,2,3,
f ( x )
k
f(x)
x*
x1
x0
x
Newton para Sistemas não Lineares:
Para sistemas não lineares:
f1 ( x1 , x2 , x3 ,..., xn )  0
f 2 ( x1 , x2 , x3 ,..., xn )  0
f3 ( x1 , x2 , x3 ,..., xn )  0
...
f n ( x1 , x2 , x3 ,..., xn )  0
x  x1
x2
x3
xn 
k
f
(
x
)
K 1
K
x x 
Ja( x k )
Onde:
 df1

 dx1
 df2
Ja   dx
 1
 ...

 ...

df1
dx2
df2
dx2
...
df1
dx3
df2
dx3
...
...
...
df1 

dxn 
df2 

dxn 
... 
dfn 
dxn 
 f1 
f 
f ( x)   2 
 f3 
 
 fn 
Métodos para Sistemas Não Lineares
Método de Newton Secante:
O método de Newton-secante baseia-se na aproximação da derivada da
função f(x), que aparece no método clássico de Newton, pela equação de
diferenças à esquerda:
df
f ( x ) 
dx
k
x
k 1
xk
 x  f (x )
k
k
f
f ( x k )  f ( x k 1 )


x
x k  x k 1
x k  x k 1
f (x )  f (x
k
k 1
k  1,2,3,
)
Dicas-Sistemas não lineares:
Uma vez definida a função e criado o arquivo contendo a mesma,
Podemos executar a subrotina criada, lembrando que a solução numérica
é sujeita a uma tolerância:
Se f(x)>tol, x é solução da equação
Para sistemas multivariaveis iterativos devemos usar a norma:
norma(xk+1 –xk)>tol
Ao usar os métodos do matlab podemos criar um vetor de opções:
op=optimset(‘metodo’)
Op=(‘Propriedade1’, valor, ‘Propriedade2’,valor)
Equações Integrais :
x max
 f ( x).dx
x min
>>help quad [ Método da quadratura de Simpson]
>>help quadl [Método da quadratura de Lobato]
>>help quad8 [Método de Alta ordem]
>>help trapz [ Método Trapezoidal]
>>help bdlquad [Método para Integrais Duplas]

y max
y min

x max
x min
f ( x, y).dx.dy
Métodos Integrais :
Regra dos Trapézios:
Regra de Simpson:
Ex: integral do sin(x)/(x+1) de 0 a 3.14
function f=funcao_01(x)
f=sin(x)./(x+1);
>> quad(‘funcao_01’,0,3.14)
Equações Diferenciais:
Aproximação com Δx pequeno:
>>help diff
dy
f ( x  x)  f ( x)
 lim x0
dx
x
Utilizando os métodos integrais:
>>help ode 23 [baixa ordem]
>>help ode23s [baixa ordem rígido]
>>help ode15s [ordem moderada rígido]
>>help ode 45 [ alta ordem]
>>help ode45s [alta ordem rígido]
>>help odeset [ set de propriedades dos métodos]
Exemplo Equações Diferenciais:
Simples:
dy
 0.1( y  10)
dt
yo  100
Function dy= df(t,y)
dy=-0.1*(y-10)
>>[t,y]=ode23(‘df’, [0 60] , 100)
Ou
>>ode23(‘df’, [0 60] , 100)
Mostra a evolução da Integração
yxt
Exemplo Equações Diferenciais:
Sistema:
dx
u
dt
du
 (1  u 2 )u  x
dt
function dydt = vdp1(t,y)
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
>>ode23(‘vdp1’,[0 20],[2 0]);
Ou
>>[t,y]= ode23(‘vdp1’,[0 20],[2 0]);
>>x=y(:,1);
>>u=y(:,2);
xo  2
u0  0
Equações Diferenciais:
Problemas de Valor inicial no Matlab:
[t,y]=odexx(‘func’,tspan,yo)
Yo-valor inicial
Tspan- valor inicial e final de tempo
Func-função a ser integrada
Método de Euler Implicito:
df ( x)
 g ( x, t )
dt
Condicoes Iniciais:
f(x)
t=0 f(0)=fo
df ( x0 , t0 ) f ( x) f ( x1 , t1 )  f ( x0 , t0 )


dt
t
x1  x0
f(x2)
f(x1)
Passo de Integracao:
f(xo)
xo
x1
x2
h  x1  x0
f ( x1, t1 )  f ( x0 , t0 )  h.g ( x0 , t0 )
Método de Euler Implícito :
Para sistemas :
df
f1 ( x1 , x2 , x3 ,..., xn ) 
dx1
f 2 ( x1 , x2 , x3 ,..., xn ) 
df
dx2
f 3 ( x1 , x2 , x3 ,..., xn ) 
df
dx3
...
df
f n ( x1 , x2 , x3 ,..., xn ) 
dxn
x  x1
x2
x3
xn 
Sistema de Equações
Algébricas que devem ser
resolvidas simultaneamente:
f1 ( x1 , t )  f1 ( x0 , t0 )  h. f1 ( x0 , t0 )
f 2 ( x1 , t )  f 2 ( x0 , t0 )  h. f 2 ( x0 , t0 )
f 3 ( x1 , t )  f 3 ( x0 , t0 )  h. f 3 ( x0 , t0 )
f n ( x1 , t )  f n ( x0 , t0 )  h. f n ( x0 , t0 )
Método de Runge Kutta Implicito :
Métodos Explícitos:
São os métodos implicitos só que a função g usa o valor de t e
y no novo ponto.
Exemplo: Euler explicito
f ( x1, t1 )  f ( x0 , t0 )  h.g ( x1, t1 )
Os métodos explícitos são mais estavéis, no entanto se a a
função g é não linear o calculo requer a solução de um sistema
não linear a cada iteração.
Para facilitar podemos usar o método preditor corretor:
1-usamos o método explicito para calcular o novo ponto
2-usamos o método implícito para corrigir.
Dicas Equações Diferenciais:
Estabilidade:
Os métodos explícitos requerem passos pequenos para manter a
estabilidade, o menor passo que pode ser dado pode ser
calculado pela expressão:
Onde p depende do método ( p=2 para Euler) e lambda é o valor
caracteristico do sistema.
Os valores característicos , são os valores que multiplicam a
variável t.
Dicas Equações Diferenciais:
O passo é limitado pela dinâmica mais rápida do sistema, uma
forma de medir essa limitação é dada pelo conceito de rigidez:
Rigidez
Problemas de Contorno:
Os métodos vistos para problemas de valor inicial, podem ser aplicados a
problemas de contorno substituindo por exemplo t por x.O único problema
é que para equações de ordem superior esses métodos requerem como
entrada a derivada no primeiro ponto de x o que em alguns casos não é
conhecido.
Para contornar esse tipo de problema podemos chutar valores para a
derivada e verificar se a solução satisfaz o outro ponto.[Shooting Method]
Ex:
d2y
x
 6y  0
dx2
2
d y
z
dx
d z 6y
 2
dx
x
C.C. Y=1 para x=1 e x=2
function F=teste(x,y)
F(1)=y(2);
F(2)=6*y(1)/x^2;
>>chute=-1.5
>> [x,y]=ode45(‘teste’,[1 2],[1 chute])
Interpolando: para x=2
y=1.1
Solução Precisa: chute=-1.516
Equações Diferenciais Parciais:
>>help pde [ toolbox de Equações Dif. Parciais]
>>pdetool [Ferramenta para Simulação]
Usando o pdetool:
1)Devemos desenhar os contornos do problema
2)Em PDE, devemos editar a equação a ser resolvida
3)Solve para resolver o problema.
Uma forma alternativa é usar Diferenças Finitas e transformar
o nosso problema em um sistema de equações algébricas.Esse
procedimento pode ser usado também para problemas de
contorno.
Otimização:
>>help optim [ toolbox de otimização]
>>help optimset [ set de propriedades dos métodos]
Otimização sem restrição:
Problema a ser resolvido:
1)podemos aplicar o conceito de derivada nula, gerando um
sistema de equações que podem ser resolvidos como visto
anteriormente usando fzero e fsolve.
2)Usando as funções:
>>help fmin
[monovariavel]
>>help fminbnd [monovariavel com limites]
>>help fminsearch [multivariavel]
Otimização sem Restrição:
Exemplo:min

S ( x)  100. x2  x12
  1  x 
2
function S=test(x)
S=100*(x(2)-x(1).^2).^2+(1-x(1)).^2;
>>x0 = [-1.2, 1]
>>[X ,S]= FMINSEARCH(‘test’,X0)
xo = [1, 1] [ ótimo encontrado]
S=0
[valor da função objetivo]
Podemos criar um vetor de Opções:
Op=(optimset,’Propriedade1’,valor,....)
2
1
Otimização:
Otimização com restrição:
Problema a ser resolvido:
>>help fmincon [ Restrições lineares e não lineares]
>>help constr [ Restrições lineares]
>>help linprog [ Programação Linear]
>>help quadprog [Programação Quadrática]
>>help lsqlin [Mínimos Quadrados]
Otimização:
Programação Linear:
S,g e h devem ser linear.
Max
s.a:
Devemos Escrever na Forma:
S  3x1  2 x2  4 x3
min f'*x
2x1  8x2  9x3  5
s.a.: A.x <= b
1.6x1  4.7 x2  8x3  4
Aeq.x=Beq
x1  5x3  8
x2  0
f   3 2  4
'
1.6 4.7 8 
A

  1 0  5
Aeq  2  8 9
Lb-limite inferior
4
b 
 8
beq  5
  Inf 
ub   0 
  Inf 
Ub-Limite superior
>>[x, S]=linprog(f,A,b,Aeq,Beq,lb,ub)
x=[-121.8936 ; -1.8723; 25.9787 ]
  Inf 
lb    Inf 
  Inf 
S=258.0213
Otimização:
Programação Não Linear:
[Programação Linear Sucessiva-SLP]
Quando um problema de otimização é não linear, seja na função
objetivo ou nas restrições, uma possibilidade para encontrar o
ótimo é através da linearização em torno do ponto ótimo.
Além disso podemos utilizar métodos que transformam um
problema com restrição em um problema sem restrição.
Exemplo:
[Multiplicadores de Lagrange]
[Função Penalidade]
Para maiores detalhes sobre os métodos, vide na referência
o material sobre otimização.
Otimização:
Programação Quadrática:
Se função objetivo é quadrática e as restrições são
lineares, podemos utilizar quadprog do Matlab.
Exemplo: min
S  x12  2x22  3x32  4x42  5x52
s.a.
2
0

H  0
0

0
2x1  x2  4x3  x4  x5  0
5x1  2x3  x4  x5  0
x1  2 x2  x3  6
0 0
4 0
0 6
0 0
0 0
0 0
0 0 
0 0
8 0

0 10
4 x3  x4  2 x5  0
A função objetivo deve ser escrita na forma: S=0.5*x'*H*x + f'*x
f  0 0 0 0 0
'
2 1  4 1  1
Aeq  

5 0  2 1  1
 1  2  1 0 0 
A
0
4 1  2
0
0
beq   
0
  6
b 
0
[x,S]=quadprog(H,f,A,b,Aeq,beq)
x=[0.4812 2.4962 0.5263 -0.6023 0.7514]
S=17.7989
Otimização:
Programação Não Linear:
Exemplo: min
s. a.


exp(x1 ). 4x12  2x22  4x1 x2  2x2  1
1.5  x1 x2  x1  x2  0
x1 x2  10
Podemos usar fmincon do Matlab, devemos criar um arquivo com a
função objetivo e se as restrições são não lineares, precisamos
criar outro arquivo com as restrições.
Function S=fob(x)
S=exp(x(1))*( 4*x(1)^2 +2*x(2)^2 +4*x(1)*x(2)+2*x(2)+1);
Function [G,H]=rest(x)
G(1)=1.5+x(1)*x(2)-x(1)-x(2);
G(2)=-x(1)*x(2)-10;
H(1)=0;
>>[x S]=fmincon('fob',[-1;2],[],[],[],[],[],[],'rest')
x=[ -9.5474 1.0474]
S=0.0236
Otimização:
Programação Inteira Mista:
Muitos problemas em operação, projeto, localização e
escalonamento de plantas envolvem variáveis que não são contínuas
e sim discretas, ou seja, variáveis que são inteiras.
Exemplo:
Um dos algoritmos numéricos mais empregados para PIM e
denominada Branch and Bound Technique.
O Matlab não possui uma rotina pronta para esse tipo de
problema. As rotinas podem ser encontradas no diretório rotinas
prontas/otimização/MILP e MINLP
Rotinas Prontas:
No diretório do CD-ROM
Rotinas\MetodosNumericos
Temos uma série de rotinas prontas separadas por tópicos.
As rotinas seguem um padrão bem similar às funções
embutidas do Matlab e um help nome da função explica o
seu funcionamento.
Variáveis Simbólicas:
>>help symbolic [ toolbox ]
Criando variaveis simbolicas:
>>help sym
>> help syms ex: >>syms x y [ cria x e y como var. simbólicas]
Manipulando variáveis simbólicas:
Uma vez criada as variáveis simbólicas podemos usar todas as
operações matemáticas do matlab.
>>f=x+2*y >>g=x*y >>f+g; >>f*g; >>f/g;
Além de algumas operações especificas para var. simbólicas:
>>finverse(f); [inversa da função f]
>>compose(f,g); [ função composta f(g(x))]
>>ezplot(f,2,3); [ plotagem de f entre os limites 2-3]
Variáveis Simbólicas:
O produto das operações pode resultar em expressões
matemáticas complicadas:
>>simple; [ coloca a expressão na forma mais simples]
>>simplify; [ simplifica a expressão]
>> pretty; [ exibe a expressão de uma forma mais visual]
Após a manipulação e simplificação pode-se desejar substituir
valores para as variáveis simbólicas:
>>subs( f,2) [ substitui em f x=2]
>>subs(f, x,2) [ substitui em f x=2 se f é função multivariavel]
>>subs(f,x,y) [ substitui em f x=y]
ex:
>> f=x+y
>>subs(subs(f,x,2),y,3) [ ans=5] [x=2 e y=3]
Resolução Simbólica:
Equações Algébricas:
>>help solve
>>[x1,x2,..xn]=solve( ‘eq1’,’eq2’,...’eqn’)
As equações podem ser escritas na forma:
‘x*y=2’ ou ‘x*y-2’
Exemplo:
>> syms x
>>f=x+4;
>>g=‘x+4’;
>>solve(f) [ ans=-4]
>>solve(g) [ ans=-4]
A vantagem é que o solve retorna todas as soluções do sistema,
no entanto, o solve não é muito robusto. Não resolvendo
sistemas muito complexos.
Resolução Simbólica:
Derivadas Simbólicas:
>> help diff
>>syms x y
>>f= 2*x + x*y + 2*y;
>>diff(f, x) [ derivada parcial em relação a x]
>>help jacobian:
>>jacobian( [f; g], [ x ; y])
Integrais Indefinidas:
>>help int
>>int(g) [ g(x), integra g em relação a x com constante de int=0]
>>int(g,x) [ g(x), integra g em relação a x com constante de int=0]
>>int(g,a,b,c) [integral definida entre a e b]
Se a constante de integração é diferente de zero, devemos somar
essa constante à solução obtida
Resolução Simbólica:
Equações Diferenciais:
>>help dsolve
>>dsolve(‘Dy=4*y’)
>>dsolve(‘Dy=4*y’, ‘y(0)=1’)
>>dsolve( eqdif 1, eqdif 2, ...., cond inicial 1,....)
Exemplo:
>>S = dsolve('Dx = y', 'Dy = -x', 'x(0)=0', 'y(0)=1')
As vezes, o matlab retorna a resposta em uma estrutura:
S.x= sin(t)
S.y=cos(t)