Curso MATLAB 6

1. Curso MATLAB 6 Instrutor: Marcelo Escobar Métodos Númericos para Engenheiros

2. Métodos Numéricos Álgebra Matricial Sistemas Lineares Sistemas não lineares Equações Integrais Equações Diferenciais Otimização Manipulação Simbólica

3. Álgebra Matricial: Tópicos de Ajuda: >>help matfun >>help elmat >>help sparfun Multiplicação Matricial: [Produto Interno] Dadas as Matrizes A e B: A * B = C [n x m ] [ m x p] [n x p] >> A*B Divisão Matricial: [Produto Externo] B= C/A >>C\A

4. 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.

5. 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)

6. 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: >>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 L y = b U x

7. Exemplo Método de Gauss: 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 x1=-1 x2=2 x3=0

8. Exemplo Método de Crammer: X1=det(Ax)/det(A) X2=det(Ay)/det(A) X3=det(Az)/det(A) Linha1=b Ax Linha2=bAy Linha3=bAz x1=-1 x2=2 x3=0

9. Métodos Indiretos: Métodos Indiretos: ( Principais) Iterações de Jacobi 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:

10. 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:

11. 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)

12. 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: >>help norm A norma é utilizada como critério de parada em loops multivariaveis.

13. 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)

14. 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])

15. Métodos para Sistemas Não Lineares Método de Substituição Sucessiva: O processo iterativo é aplicado à equação algébrica na forma modificadada 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.

16. Métodos para Sistemas Não Lineares: Método da Bisseção: Os pontos iniciais devem satisfazer a condição: onde a função sign(f(x)) fornece o sinal da função f(x).

17. 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:

18. Newton para Sistemas não Lineares: Para sistemas não lineares: Onde:

19. 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:

20. 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)

21. Equações Integrais : >>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]

22. 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)

23. Equações Diferenciais: Aproximação com Δx pequeno: >>help diff 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]

24. Exemplo Equações Diferenciais: Simples: y x t 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

25. Exemplo Equações Diferenciais: Sistema: 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);

26. 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

27. f(x) f(x2) f(x1) f(xo) x2 x1 xo Método de Euler Implicito: Condicoes Iniciais: t=0 f(0)=fo Passo de Integracao:

28. Método de Euler Implícito : Para sistemas : Sistema de Equações Algébricas que devem ser resolvidas simultaneamente:

29. Método de Runge Kutta Implicito :

30. 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 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.

31. 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.

32. 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

33. 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: 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

34. 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.

35. 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]

36. Otimização sem Restrição: Exemplo:min 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,....)

37. 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]

38. Otimização: Programação Linear: S,g e h devem ser linear. Max Devemos Escrever na Forma: s.a: min f'*x s.a.: A.x <= b Aeq.x=Beq Lb-limite inferior Ub-Limite superior >>[x, S]=linprog(f,A,b,Aeq,Beq,lb,ub) x=[-121.8936 ; -1.8723; 25.9787 ] S=258.0213

39. 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.

40. 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.a. A função objetivo deve ser escrita na forma: S=0.5*x'*H*x + f'*x [x,S]=quadprog(H,f,A,b,Aeq,beq) x=[0.4812 2.4962 0.5263 -0.6023 0.7514] S=17.7989

41. Otimização: Programação Não Linear: Exemplo: min s. a. 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

42. 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

43. 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.

44. 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]

45. 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]

46. 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.

47. 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

48. 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)