Derivadas e EDO

Renato AssunçãoDCC, UFMG

Derivada numéricaLembre 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 sucessivaSe 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’

ExemploConsidere a função e

calcule

Pelo método da diferença sucessiva:

Por outro lado, sabemos do calculo que f ’(x)=1/(1+x2) e portanto

)2('f

Um método mais precisoConsidere as seguintes DUAS expansões de

Taylor:

Isolando f ’ (x) encontramos:

Isto produz o método da diferença simétrica

Diferença simétricaEste 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 RichardsonExtrapolaçã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(h2))

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 k2, k4, .. 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(h4)

Estimando a derivada segundaVamos considerar de novo as duas expansões

de Taylor:

Isolando f’’(x) encontramos

E assim temos a aproximação

Calculo diferencial e integralNewton 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 integralLogo 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 física, a 3ª lei de Newton, envolve uma segunda derivada: F = m * d2 y(t)/dt2

Um passo alem: Equações diferenciaisEquação não linear usual:

achar os valores de t para os quais a igualdade f(t)=0 e’ válida

Equação diferencial ordinária: achar as funções y(t) para as quais a igualdade g(t, y(t), y’(t)) = 0 e’ válida PARA TODO t

Por exemplo: achar y(t) tal que seja válida a equação

y’(t)–3y(t) = 0 OU SEJA y’(t) = 3y(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ª ordemAchar y(t) tal que y’(t) = 3y(t)Geometricamente:

Desenhe o gráfico da função y(t)Calcule a inclinação da reta tangente y’(t) em

cada ponto tA 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ª ordemAchar y(t) tal que y’(t) = 3y(t)Que tal y(t) = t2 ??Neste caso, y’(t) = 2 t t2 = 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ª ordemAchar y(t) tal que y’(t) = 3y(t)Que tal y(t) = e3t ??De fato, para esta função, temos y’(t) = 3 e3t = 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 e3t 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 inicialAchar y(t) tal que

y’(t) = 3y(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)=2Como todas as soluções de y’(t) = 3y(t) são da

forma y(t) = c e3t temos de encontrar alguma que satisfaça a condição inicial.

2 = y(0) = c e3*0 = c.1 y(t) = 2 e3t

Notação: ordinária?Equações diferenciais: okMas 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 2

2

x

uc

t

u

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*t2 + 2t -1) * y(t) + sin(t)

EDO de ordem nEDO 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) + 3t 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 famosasDecaimento 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 14C 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 famosasCorpo 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 famosasOscilador harmônico amortecido:

y’’(t) + a y’(t) + b y(t) = 0Para descrever a física do átomo de hidrogênio:

Legendre: (1-t2)y’’(t) – 2 t y’(t) + k(k+1) y(t) = 0Para 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 : t2 y’’(t) + t y’(t) + (t2 – k2) y(t) = 0

Mecânica quântica:Hermite: y’’(t) – 2 t y’(t) + 2 k y(t) = 0

Arco-írisy’’(t) + t y(t) = 0

EDO linear de 1ª ordemSuponha que y’(t) = a(t) * y(t) + b(t) Solução:

Com o fator integrante

Se y’(t) = - y2(t) EDO NÃO-LINEAR. Esta EDO particular pode ser resolvida pelo

método de separação de variáveis .

EDO de 1ª ordemVamos considerar problemas do seguinte

tipo:

Existe alguma solução? Quando ela e’ única?TEOREMA:

Métodos numéricos: EulerNosso 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 [t0, tN] Divide o intervalo em N subintervalos de

comprimentos iguais: t0, t1, ..., tN onde tk=t0 + k * h com h=(tN-t0)/N

Se h e’ pequeno, temos

Método de EulerNosso problema de EDO de 1ª ordem: t0, t1, ..., tN onde tk=t0 + k * h com h=(tN-t0)/NSe h e’ pequeno, temos

Algoritmo:

Temos uma aproximação y0, y1, y2, ..., yN para y(t)

29

Método de Euler – 1ª iteracao

Passo h

t

y(t)

(t0, y0)

Valor de y(t0+ h) 00,)(,)(' yytytfty

Slope 00 , ytf

010001 , ttytfyy

)()(, 10000 tyhtyhytfy

Interpretação gráfica: primeiro passo do método de Euler

Valor aproximado y1

t0 t1

30

Método de Euler – 2ª iteração

Tamanho do passo

h

Valor verdadeiro y(x2)

y2 Valor aproximado

y1

t

y(t)

t1 t2

Segunda iteração do método de Euler

hytfyy 1112 ,

Note que y2 e’ o valor em t2 de uma reta que passa por (t1, y1) e que teminclinação f(t1,y1).

31

Método de Euler – 2ª iteração

Tamanho do passo

h

Valor verdadeiro y(t2)

y2 Valor aproximado

y1

t

y(t)

t1 t2

hytfyy 1112 ,

NÃO ESTAMOS USANDO

y*2 = y(t1) + f(t1,

y(t1))h

uma reta passando por y(t1) e com inclinação f(t1,y(t1)) poisNOS NÃO TEMOSy(t1)

Na 1ª iteração obtivemos uma aproximação y1 para y(t1)

32

Método de Euler – iteração i

Tamanho do passo

h

Valor verdadeiro y(ti+1)

yi+1, Valor aproximado

yi

t

y

ti ti+1

Passo genérico do método de Euler

hytfyy iiii ,1

ii tth 1

Erros em EulerAssim, na n-ésima iteração, gostaríamos de

aproximar yn+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-seem cada iteração do método de Euler.Existe um erro em aproximar y(tn)por yn , a n-ésima iteraçãoAlém disso, gostaríamos de ter f(tn,y(tn)) mas usamos f(tn,yn)

34

Uma esfera possui temperatura de 1200K (ou 920oC ) e comeca a resfriar ‘a temperatura ambiente de 300K (ou 27oC). 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

Kdt

d12000,1081102067.2 8412

Encontrar a temperatura em t=480 segundos usando o método de Euler.

Assuma um passo de tamanho h=240 segundos

35

K

f

htf

htf iiii

09.106

2405579.41200

24010811200102067.21200

2401200,01200

,

,

8412

0001

1

Primeira iteração:

1 e’ a temperatura aproximada em 240240001 httt

K09.106240 1

8412 1081102067.2

dt

d

8412 1081102067.2, tf

36

Solução - continuaçãoPara 09.106,240,1 11 ti

K

f

htf

32.110

240017595.009.106

240108109.106102067.209.106

24009.106,24009.106

,

8412

1112

Iteração 2:

2 e’ a temperatura aproximada em 48024024012 httt

K32.110480 2

37

Solução – continuação

A solução exata da EDO e’ dada pela raiz da equação não-linear

9282.21022067.0)(00333.0tan8519.1300)(

300)(ln92593.0 31

tt

t

t

A solução (480) desta equação não-linear em t=480 segundos e’

K57.647)480(

Bem diferente da aproximação: K32.110480 2

38

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

0

200

400

600

800

1000

1200

1400

0 100 200 300 400 500

Te

mp

era

ture

,

Time, t(sec)

h=240

Exact Solution

θ(K

)

Euler

Step, h (480) Et |єt|%

4802401206030

−987.81110.32546.77614.97632.77

1635.4537.26100.8032.60714.806

252.5482.96415.5665.03522.2864

39

Efeito do tamanho do passo hTemperatura aos 480 segundos como uma função do passo h

K57.647)480( (valor exato)

40

Comparação com resultado exato

-1500

-1000

-500

0

500

1000

1500

0 100 200 300 400 500

Time, t (sec)Tem

per

atu

re,

Exact solution

h=120h=240

h=480

θ(K

)

Apenas h=480, 240 e 120

41

Efeito do tamanho do passo h em Euler

-1200

-800

-400

0

400

800

0 100 200 300 400 500

Step size, h (s) Te

mp

era

ture

, θ(K

)Valor exato

Mais um exemploConsidere a EDO

Como x0=0 então xn=nh Iteração:

Usando h=0.1 e 0.001E comparando com a solução exata temos a tabela ao lado

43

Erros no método de EulerVimos 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 xi

...!3

1

!2

1 31

,

3

32

1

,

2

2

1,

1 ii

yx

ii

yx

iiyx

ii xxdx

ydxx

dx

ydxx

dx

dyyy

iiiiii

...),(''!3

1),('

!2

1),( 3

12

111 iiiiiiiiiiiiii xxyxfxxyxfxxyxfyy

Os dois primeiros termos da serie de Taylor e’ o método de Euler

hyxfyy iiii ,1

O erro na aproximação e’ dado por

...

!3

,

!2

, 32

hyxf

hyxf

E iiiit

2hEt

Isto e’:

Runge - KuttaEuler 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-KuttaEncontramos 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ª ordemEquação de iteração:

ou simplesmente

onde

Assim, este e’ um método de Euler com inclinação (s1+s2)/2

Runge – Kutta de 2ª ordemE’ 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 s1 em tn:

Com isto, obtenha uma segunda inclinação s2 em tn+hA atualização de Euler realmente dada usa a média das

inclinações s1 em tn e s2 em tn+h

48


Kdt

d12000,1081102067.2 8412

Encontrar a temperatura em t=480 segundos usando o método EULER MELHORADO (ou metodo classico de Runge-Kutta de segunda ordem)


8412 1081102067.2

dt

d

8412 1081102067.2, tf

hssii

211 2

1

2

1

49

Iteração 1: Kti 1200)0(,0,0 00

5579.4

10811200102067.2

1200,0

,

8412

01

f

tfs o

017595.0

108109.106102067.2

09.106,240

2405579.41200,2400

,

8412

1002

f

f

hshtfs

K

hss

16.655

2402702.21200

240017595.02

15579.4

2

11200

2

1

2

12101

50

Solução - continuaçãoIteração 2: Khtti 16.655,2402400,1 101

38869.0

108116.655102067.2

16.655,240

,

8412

111

f

tfs

20206.0

108187.561102067.2

87.561,480

24038869.016.655,240240

,

8412

1112

f

f

hkhtfs

K

hss

27.584

24029538.016.655

24020206.02

138869.0

2

116.655

2

1

2

12112

51

Solução - continuação

A solução exata da EDO e’ dada pela solução de uma equação não -linear:

9282.21022067.0)(0033333.0tan8519.1300)(

300)(ln92593.0 31

tt

t

t

A solução para esta equação não-linear em t=480 segundos e’

K57.647)480(

52

Comparação com resultado exatos

Euler melhorado (ponto médio) para diferentes valores de h

-400

0

400

800

1200

0 100 200 300 400 500

Time, t(sec)

Tem

per

atu

re,θ

(K) Exact h=120

h=240

h=480

53

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|%

4802401206030

−393.87584.27651.35649.91648.21

1041.463.304

−3.7762−2.3406

−0.63219

160.829.7756

0.583130.36145

0.097625

K57.647)480( (exact)

54

Efeito do tamanho do passo h

-400

-200

0

200

400

600

800

0 100 200 300 400 500Step size, h

Tem

per

atu

re,

θ(48

0)

Um segundo método de Runge-KuttaO 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. Vamos usar Euler de novo para este valor.

2º. Método de Runge - KuttaTemos 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/(,2

))(,(1 htyh

thfydyfyy nnn

ht

t

nn

n

n

2º metodo de Runge-KuttaTambém podemos ver este novo método de

Runge-Kutta como um processo em dois estágios.

Escrevemos

como

onde

Resumo dos 2 métodos de R-KPrimeiro: o método clássico de 2ª ordem de

R-K (ou método de Euler melhorado) yn+1 = yn + h (s1+s2)/2com

Segundo: Método de Euler modificado (método do ponto médio)yn+1 = yn + h s2

com

O que eles tem em comum?

Comparando os dois R-KOs dois métodos usam dois estágios intermediários s1 e s2

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 (s1+s2)/2 tomando a media das inclinações s1 em tn e s2 em tn+h

No método de Euler modificado (ponto médio), nós usamos s1 em tn para dar um meio-passo ate tn+h/2. A seguir, calculamos s2, a estimativa da inclinação no ponto médio, e então tomamos o passo completo yn+1 = yn + h s2

ExemploConsidere a EDO

Euler modificado: yn+1 =yn+hs2

Temos s1=x2n+ y2

n

e s2=(xn+h/2)2+(yn+s1/2)2

Exemplo numérico na tabela ao ladoy(xi) e’ o valor exato e yi e’ a aproximação numérica

61


Kdt

d12000,1081102067.2 8412

Encontrar a temperatura em t=480 segundos usando o método EULER MELHORADO (ou método clássico de Runge-Kutta de segunda ordem) e EULER MODIFICADO (ou ponto médio)


Passoh

(480)

EulerEuler

MelhoradoPonto Medio

Ralston(ignore)

4802401206030

−987.84110.32546.77614.97632.77

−393.87584.27651.35649.91648.21

1208.4976.87690.20654.85649.02

449.78690.01667.71652.25648.61

62

Comparação de Euler e RK de 2a ordem

K57.647)480( (exato)

63

Passo hEuler

Euler Melhorado

Ponto Médio

Ralston (ignore)

480

240

120

60

30

252.54

82.964

15.566

5.0352

2.2864

160.82

9.7756

0.58313

0.36145

0.097625

86.61250.8516.58231.1239

0.22353

30.5446.55373.1092

0.722990.15940

K57.647)480( (exato)

%t


64

500

600

700

800

900

1000

1100

1200

0 100 200 300 400 500 600

Tem

per

atu

re,

Time, t (sec)

Analytical

Ralston

Midpoint

Euler

Heun

θ(K

)


Melhorado

Ponto Médio

Para a provaMemorizar apenas:

Método de Euler e os dois métodos mais simples de Runge-Kutta: Euler melhorado (RK clássico de 2ª ordem) Euler modificado (ponto médio)

Pode ignorar o resto dos slides

Runge-Kutta 2ª ordem geralPodemos imaginar varias outras maneiras alternativas

de calcular s1 e s2.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 ButcherE’ 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 ButcherRK Clássico (Euler melhorado)

RK : Euler modificado (ponto médio)

RK: Método de Heun

α2 = β21

Método de Ralston0 0 0

α2=3/4

β21=3/4

0

Γ1=1/3

Γ2=2/3

Runge-Kutta de 4ª ordemE’ o mais famoso método de Runge-Kutta

com

E tabela de Butcher

Documents

Derivadas e EDO