Derivadas e EDO Renato Assuno DCC UFMG Derivada

Derivadas e EDO Renato Assunção DCC, UFMG

Derivada numerica �Lembre da definição de derivada �Como no caso da integral, a definição e’ uma operação de limite quando h 0 �Uma estimativa simples e’ então tomar h ≈0 e usar a aproximação �Este e’ chamado o método das diferença sucessiva (forward difference method)

Precisão da diferença sucessiva �Se f e’ diferenciavel duas vezes, podemos escrever sua expansao de Taylor de 2ª ordem: �Onde ε (x, x+h) �Entao �Portanto, o erro de Dhf e’

Exemplo �Considere a função e calcule �Pelo método da diferença sucessiva: �Por outro lado, sabemos do calculo que f ’(x)=1/(1+x 2) e portanto

Um método mais preciso �Considere as seguintes DUAS expansões de Taylor: �Isolando f ’ (x) encontramos: �Isto produz o método da diferença simétrica

Diferença simétrica �Este método e’ uma media dos métodos de diferença sucessiva e diferença retroativa. �Qual a precisão desta media? �Como �e �O método da diferença simétrica tem um erro de aproximação igual a

Extrapolação de Richardson �Extrapolação de Richardson pode ser usada para melhorar qualquer método numérico que tenha a ordem de grandeza de seu erro conhecida. �Podemos usa-la para melhorar: �Diferença sucessiva (O(h)) �Diferença simétrica (O(h 2)) �Nos podemos ser bastante específicos sobre o impacto da extrapolação de Richardson nestes dois casos.

�Se nos tivéssemos tomado a expansão completa de Taylor quando derivamos o método das diferenças simétricas teríamos: �ou �onde k 2, k 4, . . são constantes independentes de h.

�Agora podemos olhar a precisao da extrapolacao de Richardson para diferencas simetricas (p=2): � 4/3 Dh/2 – 1/3 Dh �Do slide anterior: �Multiplicando pelos fatores apropriados temos �Substituindo (2) em (1) temos �Assim, diferenças simétricas com extrapolação de Richardson tem erro O(h 4)

Estimando a derivada segunda �Vamos considerar de novo as duas expansões de Taylor: �Isolando f’’(x) encontramos �E assim temos a aproximação

Equações diferenciais ordinárias de 1ª ordem

Calculo diferencial e integral �Newton e Leibniz inventaram o calculo diferencial e integral por volta de 1670. �O objetivo era ter ferramentas matemáticas apropriadas para lidar com o movimento e mudanças no tempo. �Entender, modelar e predizer o movimento dos corpos celestes era um dos maiores objetivos da ciência naqueles tempos. �Para os homens daquele tempo, compreender o movimento dos céus era ouvir a voz de Deus.

Calculo diferencial e integral �Logo depois ocorre uma explosão cientifica revolucionaria: Os irmãos Bernoulli, Euler, Lagrange, Laplace, etc. �A ciência e engenharia modernas nascem, vicejam, crescem e criam o que temos hoje em dia. �Derivadas e integrais aparecem em todos os modelos científicos para descrever a natureza se um processo de mudança estiver envolvido no fenômeno estudado. �A mais famosa lei da fisica, a 3ª lei de Newton, envolve uma segunda derivada: F = m * d 2 y(t)/dt 2

Um passo alem: Equações diferenciais �Equação não linear usual: �achar os valores de t para os quais a igualdade f(t)=0 e’ valida �Equação diferencial ordinária: achar as funções y(t) para as quais a igualdade g(t, y(t), y’(t)) = 0 PARA TODO t �Por exemplo: achar y(t) tal que seja valida a equação �y’(t)– 3 y(t) = 0 OU SEJA y’(t) = 3 y(t). �Isto e’, queremos achar as todas as funções y(t) tais que a sua função derivada y’(t) seja igual a 3 vezes a própria função y(t).

EDO de 1ª ordem �Achar y(t) tal que y’(t) = 3 y(t) �Geometricamente: �Desenhe o gráfico da função y(t) �Calcule a inclinação da reta tangente y’(t) em cada ponto t �A inclinação deve ser igual a 3 vezes o valor da função y(t) �Existe alguma função que satisfaz esta condição? �Se existem, e’ possível encontra-las? �Técnicas de solução ANALITICA de EDO: solução exata

EDO de 1ª ordem �Achar y(t) tal que y’(t) = 3 y(t) �Que tal y(t) = t 2 ? ? �Neste caso, y’(t) = 2 t t 2 = y(t) �Que tal y(t) = cos(t) ? �Neste caso, y’(t) = -sen(t) cos(t) = y(t) �OU ainda y(t) = log(t) ? ? �y‘(t) = 1/t que não e’ a própria função y(t)=log(t)

EDO de 1ª ordem �Achar y(t) tal que y’(t) = 3 y(t) �Que tal y(t) = e 3 t ? ? �De fato, para esta função, temos y’(t) = 3 e 3 t = 3 y(t) � e’ uma solução da equação diferencial. �Existe alguma outra função que também seja solução? �Sim: todas as funções da forma y(t) = c e 3 t onde c R também e’ uma solução. �Existem outras? Não, estas são todas, não existem mais funções para as quais temos y’(t) = 3 y(t) �PARA TODO t

EDO 1ª ordem com valor inicial �Achar y(t) tal que �y’(t) = 3 y(t) �E ALEM DISSO, y(0) = 2 �Agora, colocamos uma restrição adicional, uma condição sobre o valor inicial da função y(t). �No tempo t=0, a função deve valer y(0)=2 �Como todas as soluções de y’(t) = 3 y(t) são da forma y(t) = c e 3 t temos de encontrar alguma que satisfaça a condição inicial. � 2 = y(0) = c e 3*0 = c. 1 y(t) = 2 e 3 t

Notação: ordinária? �Equações diferenciais: ok �Mas por que Equações diferenciais ORDINÁRIAS? �Existem Equações diferenciais EXTRAORDINÁRIAS? �Não. O palavra “ordinária” e’ usada para diferenciar das equações diferenciais PARCIAIS. �Parciais: equações que envolvem funções de mais de uma variável e suas derivadas parciais. �Exemplo: Equação de difusão do calor numa barra de densidade homogênea. �Seja u(t, x) a temperatura no ponto x no tempo t �Então onde c depende do material

EDO de 1ª ordem �Equações diferenciais ordinárias de 1ª ordem são equações envolvendo apenas a derivada y’(t) e a funcao y(t) e possivelmente outras funções FIXAS e conhecidas (tais como sin(t) ou exp(t)). �Por exemplo: �y’(t) = p(t) * y(t) + g(t) �Casos particulares: �y’(t) = 3 * y(t) + sin(t) �y’(t) = (3*t 2 + 2 t -1) * y(t) + sin(t)

EDO de ordem n �EDO ordem n são equações envolvendo : �As derivadas yn(t), yn-1(t), . . . , y’(t) �a função y(t) �e possivelmente outras funções FIXAS e conhecidas (tais como sin(t) ou exp(t)). �Por exemplo, uma EDO de 2ª ordem: �y”(t) = sin(t) * y’(t) + y(t) + 3 t �Qual a (ou as) FUNCAO y(t) tal que a sua FUNCAO derivada segunda y’’(t) obedece a equação acima? �Mas isto e’ so’ um exercício de matemáticos sem ter o que fazer, certo?

Exemplos de EDOs famosas �Decaimento radioativo: proporção carbono 14/carbono-12 presente na matéria orgânica viva é constante. �No entanto, na matéria orgânica morta a quantidade de 14 C diminui com o tempo, a uma taxa proporcional à quantidade existente. �Se designarmos essa quantidade por Q, teremos: �Q’(t) = -c Q(t) onde c > 0 e’ uma constante

Exemplos de EDOs famosas �Corpo em queda livre com atrito devido a resistência do ar: �Mv’(t) = mg – k v(t) ou v’(t) + k/m v(t) – g = 0 �Engenharia Química: balanço de massa ou volume ou energia num reator químico. �O volume de líquido num tanque e a concentração de uma solução A mudam com o tempo. �Entra e sai líquido a taxas constantes e diferentes. �Os líquidos possuem concentrações de A diferentes. �Descrever a concentração de A em cada instante : terminamos em uma EDOs

Exemplos de EDOs famosas �Oscilador harmônico amortecido: �y’’(t) + a y’(t) + b y(t) = 0 �Para descrever a física do átomo de hidrogênio: �Legendre: (1 -t 2)y’’(t) – 2 t y’(t) + k(k+1) y(t) = 0 �Para descrever o comprimento de onda no átomo de hidrogênio: �Laguerre: t y’’(t) + (1 -t) y’(t) + k y(t) = 0 �Membranas vibratórias: �Bessel : t 2 y’’(t) + t y’(t) + (t 2 – k 2) y(t) = 0 �Mecânica quântica: �Hermite: y’’(t) – 2 t y’(t) + 2 k y(t) = 0 �Arco-íris �y’’(t) + t y(t) = 0

EDO linear de 1ª ordem �Suponha que y’(t) = a(t) * y(t) + b(t) �Solução: �Com o fator integrante �Se y’(t) = - y 2(t) EDO NÃO-LINEAR. �Esta EDO particular pode ser resolvida pelo método de separação de variáveis.

EDO de 1ª ordem �Vamos considerar problemas do seguinte tipo: �Existe alguma solução? Quando ela e’ única? �TEOREMA:

Métodos numéricos: Euler �Nosso problema de EDO de 1ª ordem: �Euler e’ o método mais simples. �Acha uma aproximação para a solução y(t) num intervalo [t 0, t. N] �Divide o intervalo em N subintervalos de comprimentos iguais: t 0, t 1, . . . , t. N onde tk=t 0 + k * h com h=(t. N-t 0)/N �Se h e’ pequeno, temos

Método de Euler �Nosso problema de EDO de 1ª ordem: �t 0, t 1, . . . , t. N onde tk=t 0 + k * h com h=(t. N-t 0)/N �Se h e’ pequeno, temos �Algoritmo: �Temos uma aproximação y 0, y 1, y 2, . . . , y. N para y(t)

Método de Euler – 1ª iteracao y(t) Valor de y(t 0+ h) Slope (t 0, y 0) Valor aproximado y 1 Passo h t t 0 t 1 Interpretação gráfica: primeiro passo do método de Euler 29

Método de Euler – 2ª iteração Valor verdadeiro y(x 2) y(t) Note que y 2 e’ o Valor em x 2 de uma reta que passa por (t 1, y 1) e que tem inclinação f(t 1, y 1). y 2 Valor aproximado y 1 h Tamanho do passo t 1 t 2 t Segunda iteração do método de Euler 30

Método de Euler – 2ª iteração NÃO ESTAMOS USANDO Valor verdadeiro y(t 2) y(t) y*2 = y(t 1) + f(t 1, y(t 1))h uma reta passando por y(t 1) e com inclinação f(t 1, y(t 1)) pois NOS NÃO TEMOS y(t 1) Na 1ª iteração obtivemos uma aproximação y 1 para y(t 1) y 2 Valor aproximado y 1 h Tamanho do passo t 1 t 2 t 31

Método de Euler – iteração i y Valor verdadeiro y(ti+1) yi+1, Valor aproximado yi h Tamanho do passo ti ti+1 t Passo genérico do método de Euler 32

Erros em Euler �Assim, na n-ésima iteração, gostaríamos de aproximar y n+1 pelo valor em tn+1 = tn+h da reta tangente a y(t) no ponto (tn, y(tn)) �Entretanto, NÃO TEMOS y(tn) �mas somente uma aproximação yn �Assim, temos dois erros acumulando-se �em cada iteração do método de Euler. �Existe um erro em aproximar y(tn) �por yn , a n-ésima iteração �Além disso, gostaríamos de ter f(tn, y(tn)) �mas usamos f(tn, yn)

Exemplo Uma esfera possui temperatura de 1200 K (ou 920 o. C ) e comeca a resfriar ‘a temperatura ambiente de 300 K (ou 27 o. C). Assumindo que o calor e’ dissipado sem interferência, a equação diferencial que reflete a temperatura (t) da esfera no tempo t deve satisfazer a seguinte EDO Encontrar a temperatura em t=480 segundos usando o método de Euler. Assuma um passo de tamanho h=240 segundos 34

Solução Primeira iteração: e’ a temperatura aproximada em 35

Solução - continuação Iteração 2: Para e’ a temperatura aproximada em 36

Solução – continuação A solução exata da EDO e’ dada pela raiz da equação não-linear A solução (480) desta equação não-linear em t=480 segundos e’ Bem diferente da aproximação: 37

Comparação das soluções exata e numérica Euler 38

Efeito do tamanho do passo h Temperatura aos 480 segundos como uma função do passo h Step, h (480) Et |єt|% 480 240 120 60 30 − 987. 81 110. 32 546. 77 614. 97 632. 77 1635. 4 537. 26 100. 80 32. 607 14. 806 252. 54 82. 964 15. 566 5. 0352 2. 2864 (valor exato) 39

Comparação com resultado exato Apenas h=480, 240 e 120 40

Efeito do tamanho do passo h em Euler Valor exato 41

Mais um exemplo �Considere a EDO �Como x 0=0 então xn=nh �Iteração: �Usando h=0. 1 e 0. 001 �E comparando com a �solução exata �temos a tabela ao lado

Erros no método de Euler Vimos que o método de Euler PODE ter erros grandes. Para entender a ordem de grandeza desses erros, vamos fazer a expansão de Taylor em torno de x i Isto e’: Os dois primeiros termos da serie de Taylor e’ o método de Euler O erro na aproximação e’ dado por 43

Runge - Kutta �Euler fez a seguinte aproximação �Que tal usar uma aproximação melhor para a integral? �Por exemplo, podemos usar a regra do trapézio: �Neste caso, teremos então a aproximação �E o algoritmo

Runge-Kutta �Encontramos a equação de iteração: �Existe um problema no entanto: yn+1 aparece dos dois lados da equação acima. Não conseguimos isolar yn+1. �Uma possibilidade e’ substituir yn+1 NO LADO DIREITO por sua aproximação baseada em Euler: yn+1 = yn + f(tn, yn)h �Este e’ o metodo de Runge-Kutta de 2ª ordem

Runge Kutta de 2ª ordem �Equação de iteração: �ou simplesmente �onde �Assim, este e’ um método de Euler com inclinação (s 1+s 2)/2

Runge – Kutta de 2ª ordem �E’ possível uma interpretação gráfica-geométrica deste método de Runge-Kutta. Temos �com �Isto corresponde ao seguinte esquema em dois passos: �Tome um passo preliminar de Euler com inclinação s 1 em t n: �Com isto, obtenha uma segunda inclinação s 2 em tn+h �A atualização de Euler realmente dada usa a média das inclinações s 1 em tn e s 2 em tn+h

Um segundo método de Runge-Kutta �O método de Runge-Kutta que acabamos de estudar começou aproximando uma integral pela regra do trapézio: �Podemos usar alguma outra regra: Simpson ou midpoint �Vamos usar midpoint: �Neste caso �Note que y(t+h/2) no lado direito não e’ conhecido.

2º. Método de Runge - Kutta �Temos a aproximação �Usamos a aproximação de Euler para o termo y(tn+h/2): �y(tn+h/2) ≈y(tn)+h/2 * f(tn, yn) �Substituindo a iteração para yn+1 temos �Este método e’ conhecido como método de Euler modificado ou método do ponto médio

2º metodo de Runge-Kutta �Também podemos ver este novo método de Runge. Kutta como um processo em dois estágios. �Escrevemos �como �onde

2º metodo de Runge-Kutta �Também podemos ver este novo método de Runge. Kutta como um processo em dois estágios. �Tome um passo de Euler mas apenas com metade do comprimento do intervalo h/2 �Isto corresponde ao tempo tn+h/2 = tn+1/2 �A seguir, de mais um passo de Euler de comprimento h usando a inclinação no ponto médio (tn+1/2, yn+1/2)

Resumo dos 2 métodos de R-K �Primeiro: o método clássico de 2ª ordem de R-K (ou método de Euler melhorado) �yn+1 = yn + h (s 1+s 2)/2 �com �Segundo: Método de Euler modificado (método do ponto médio) �yn+1 = yn + h s 2 �com �O que eles tem em comum?

Comparando os dois R-K �Os dois métodos usam dois estágios intermediários s 1 e s 2 para obter uma iteração. �Os estágios correspondem a diferentes estimativas para a inclinação da solução. � No método clássico de RK (Euler melhorado) nós damos um passo completo yn+1 = yn + h (s 1+s 2)/2 tomando a media das inclinações s 1 em tn e s 2 em tn+h �No método de Euler modificado (ponto médio), nós usamos s 1 em tn para dar um meio-passo ate tn+h/2. A seguir, calculamos s 2, a estimativa da inclinação no ponto médio, e então tomamos o passo completo yn+1 = yn + h s 2

Exemplo �Considere a EDO �Euler modificado: yn+1 =yn+hs 2 �Temos s 1=x 2 n+ y 2 n �e s 2=(xn+h/2)2+(yn+s 1/2)2 �Exemplo numérico na tabela ao lado

Runge-Kutta 2ª ordem geral �Podemos imaginar varias outras maneiras alternativas de calcular s 1 e s 2. �O método geral de Runge-Kutta de 2ª ordem e’ da forma �onde �com (esta notação vem de uma teoria mais avançada ligada a métodos implícitos) �Clássico RK (Euler melhorado): � Euler modificado (ponto médio): γ 1=0, γ 2=1 e α 2= β 21=1/2

Tabela de Butcher �E’ costume arranjar os coeficientes αi, βij e γi em uma tabela chamada tabela de Butcher �Onde α 2 = β 21 �Para o método ser de segunda ordem e ter certas propriedades desejáveis impomos também as condições

Tabela de Butcher �RK Clássico (Euler melhorado) � α 2 = β 21 �RK : Euler modificado (ponto médio) Método de Ralston �RK: Método de Heun 0 0 α 2=3/4 β 21=3/4 0 Γ 1=1/3 0 Γ 2=2/3

Exemplo Uma esfera possui temperatura de 1200 K (ou 920 o. C ) e comeca a resfriar ‘a temperatura ambiente de 300 K (ou 27 o. C). Assumindo que o calor e’ dissipado sem interferência, a equação diferencial que reflete a temperatura (t) da esfera no tempo t deve satisfazer a seguinte EDO Encontrar a temperatura em t=480 segundos usando o método EULER MELHORADO (ou metodo classico de Runge-Kutta de segunda ordem) Assuma um passo de tamanho h=240 segundos 58

Solução Iteração 1: 59

Solução - continuação Iteração 2: 60

Solução - continuação A solução exata da EDO e’ dada pela solução de uma equação não -linear: A solução para esta equação não-linear em t=480 segundos e’ 61

Comparação com resultado exatos Euler melhorado (ponto médio) para diferentes valores de h 62

Efeito do tamanho do passo h Temperatura em t=480 segundos como uma funcao do tamanho do passo h Passo h (480) Erro = Et |єt|% 480 240 120 60 30 − 393. 87 584. 27 651. 35 649. 91 648. 21 1041. 4 63. 304 − 3. 7762 − 2. 3406 − 0. 63219 160. 82 9. 7756 0. 58313 0. 36145 0. 097625 (exact) 63

Efeito do tamanho do passo h 64

Comparação de Euler e RK de 2 a ordem Passo h 480 240 120 60 30 q(480) Euler Melhorado Ponto Medio Ralston − 987. 84 110. 32 546. 77 614. 97 632. 77 − 393. 87 584. 27 651. 35 649. 91 648. 21 1208. 4 976. 87 690. 20 654. 85 649. 02 449. 78 690. 01 667. 71 652. 25 648. 61 (exato) 65

Comparação de Euler e RK de 2 a ordem Passo h 480 240 120 60 30 Euler Modificado Ponto Médio Ralston 252. 54 82. 964 15. 566 5. 0352 2. 2864 160. 82 9. 7756 0. 58313 0. 36145 0. 097625 86. 612 50. 851 6. 5823 1. 1239 0. 22353 30. 544 6. 5537 3. 1092 0. 72299 0. 15940 (exato) 66

Comparação de Euler e RK de 2 a ordem Ponto Medio Modificado 67

Runge-Kutta de 4ª ordem �E’ o mais famoso método de Runge-Kutta �com �E tabela de Butcher

Para a prova �Memorizar apenas os dois métodos mais simples de Runge-Kutta: �Euler melhorado (RK clássico de 2ª ordem) �Euler modificado (ponto médio)