Razes e otimizao

Derivadas e EDORenato AssunoDCC, UFMG1Derivada numericaLembre da definio de derivada

Como no caso da integral, a definio e uma operao de limite quando h 0Uma estimativa simples e ento tomar h 0 e usar a aproximao

Este e chamado o mtodo das diferena sucessiva (forward difference method)

2Preciso da diferena 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

3ExemploConsidere a funo e calcule

Pelo mtodo da diferena sucessiva:

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

4Um mtodo mais precisoConsidere as seguintes DUAS expanses de Taylor:

Isolando f (x) encontramos:

Isto produz o mtodo da diferena simtrica

5Diferena simtricaEste mtodo e uma media dos mtodos de diferena sucessiva e diferena retroativa.Qual a preciso desta media?Como

e

O mtodo da diferena simtrica tem um erro de aproximao igual a

6Extrapolao de RichardsonExtrapolao de Richardson pode ser usada para melhorar qualquer mtodo numrico que tenha a ordem de grandeza de seu erro conhecida. Podemos usa-la para melhorar:Diferena sucessiva (O(h))Diferena simtrica (O(h2))Nos podemos ser bastante especficos sobre o impacto da extrapolao de Richardson nestes dois casos.

7Se nos tivssemos tomado a expanso completa de Taylor quando derivamos o mtodo das diferenas simtricas teramos:

ou

onde k2, k4, .. so constantes independentes de h.

8Agora podemos olhar a precisao da extrapolacao de Richardson para diferencas simetricas (p=2):4/3 Dh/2 1/3 DhDo slide anterior:

Multiplicando pelos fatores apropriados temos

Substituindo (2) em (1) temos

Assim, diferenas simtricas com extrapolao de Richardson tem erro O(h4)

9Estimando a derivada segundaVamos considerar de novo as duas expanses de Taylor:

Isolando f(x) encontramos

E assim temos a aproximao

10Equaes diferenciais ordinrias de 1 ordem11Calculo diferencial e integralNewton e Leibniz inventaram o calculo diferencial e integral por volta de 1670.O objetivo era ter ferramentas matemticas apropriadas para lidar com o movimento e mudanas no tempo.Entender, modelar e predizer o movimento dos corpos celestes era um dos maiores objetivos da cincia naqueles tempos. Para os homens daquele tempo, compreender o movimento dos cus era ouvir a voz de Deus.12Calculo diferencial e integralLogo depois ocorre uma exploso cientifica revolucionaria: Os irmos Bernoulli, Euler, Lagrange, Laplace, etc. A cincia e engenharia modernas nascem, vicejam, crescem e criam o que temos hoje em dia. Derivadas e integrais aparecem em todos os modelos cientficos para descrever a natureza se um processo de mudana estiver envolvido no fenmeno estudado.A mais famosa lei da fisica, a 3 lei de Newton, envolve uma segunda derivada: F = m * d2 y(t)/dt2

13Um passo alem: Equaes diferenciaisEquao no linear usual: achar os valores de t para os quais a igualdade f(t)=0 e valida Equao diferencial ordinria: achar as funes y(t) para as quais a igualdade g(t, y(t), y(t)) = 0 PARA TODO tPor exemplo: achar y(t) tal que seja valida a equao y(t)3y(t) = 0 OU SEJA y(t) = 3y(t).Isto e, queremos achar as todas as funes y(t) tais que a sua funo derivada y(t) seja igual a 3 vezes a prpria funo y(t).

14EDO de 1 ordemAchar y(t) tal que y(t) = 3y(t)Geometricamente: Desenhe o grfico da funo y(t)Calcule a inclinao da reta tangente y(t) em cada ponto tA inclinao deve ser igual a 3 vezes o valor da funo y(t)Existe alguma funo que satisfaz esta condio?Se existem, e possvel encontra-las?Tcnicas de soluo ANALITICA de EDO: soluo exata

15EDO 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 no e a prpria funo y(t)=log(t)

16EDO de 1 ordemAchar y(t) tal que y(t) = 3y(t)Que tal y(t) = e3t ??De fato, para esta funo, temos y(t) = 3 e3t = 3 y(t) e uma soluo da equao diferencial. Existe alguma outra funo que tambm seja soluo? Sim: todas as funes da forma y(t) = c e3t onde c R tambm e uma soluo.Existem outras? No, estas so todas, no existem mais funes para as quais temos y(t) = 3 y(t) PARA TODO t

17EDO 1 ordem com valor inicialAchar y(t) tal que y(t) = 3y(t)E ALEM DISSO, y(0) = 2Agora, colocamos uma restrio adicional, uma condio sobre o valor inicial da funo y(t).No tempo t=0, a funo deve valer y(0)=2Como todas as solues de y(t) = 3y(t) so da forma y(t) = c e3t temos de encontrar alguma que satisfaa a condio inicial.2 = y(0) = c e3*0 = c.1 y(t) = 2 e3t

18Notao: ordinria?Equaes diferenciais: okMas por que Equaes diferenciais ORDINRIAS? Existem Equaes diferenciais EXTRAORDINRIAS?No. O palavra ordinria e usada para diferenciar das equaes diferenciais PARCIAIS.Parciais: equaes que envolvem funes de mais de uma varivel e suas derivadas parciais. Exemplo: Equao de difuso do calor numa barra de densidade homognea. Seja u(t,x) a temperatura no ponto x no tempo t Ento onde c depende do material

19EDO de 1 ordem Equaes diferenciais ordinrias de 1 ordem so equaes envolvendo apenas a derivada y(t) e a funcao y(t) e possivelmente outras funes 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)

20EDO de ordem nEDO ordem n so equaes envolvendo :As derivadas yn(t), yn-1(t),..., y(t) a funo y(t) e possivelmente outras funes 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 equao acima?Mas isto e so um exerccio de matemticos sem ter o que fazer, certo?21Exemplos de EDOs famosasDecaimento radioativo: proporo carbono-14/carbono-12 presente na matria orgnica viva constante.No entanto, na matria orgnica 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 constante22Exemplos de EDOs famosasCorpo em queda livre com atrito devido a resistncia do ar:Mv(t) = mg k v(t) ou v(t) + k/m v(t) g = 0Engenharia Qumica: balano de massa ou volume ou energia num reator qumico. O volume de lquido num tanque e a concentrao de uma soluo A mudam com o tempo. Entra e sai lquido a taxas constantes e diferentes. Os lquidos possuem concentraes de A diferentes. Descrever a concentrao de A em cada instante : terminamos em uma EDOs23Exemplos de EDOs famosasOscilador harmnico amortecido:y(t) + a y(t) + b y(t) = 0Para descrever a fsica do tomo de hidrognio:Legendre: (1-t2)y(t) 2 t y(t) + k(k+1) y(t) = 0Para descrever o comprimento de onda no tomo de hidrognio:Laguerre: t y(t) + (1-t) y(t) + k y(t) = 0Membranas vibratrias:Bessel : t2 y(t) + t y(t) + (t2 k2) y(t) = 0Mecnica quntica:Hermite: y(t) 2 t y(t) + 2 k y(t) = 0Arco-risy(t) + t y(t) = 024EDO linear de 1 ordemSuponha que y(t) = a(t) * y(t) + b(t) Soluo:

Com o fator integrante

Se y(t) = - y2(t) EDO NO-LINEAR. Esta EDO particular pode ser resolvida pelo mtodo de separao de variveis .

25EDO de 1 ordemVamos considerar problemas do seguinte tipo:

Existe alguma soluo? Quando ela e nica?TEOREMA:

26Mtodos numricos: EulerNosso problema de EDO de 1 ordem: Euler e o mtodo mais simples.Acha uma aproximao para a soluo 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)/NSe h e pequeno, temos

27Mtodo 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 aproximao y0, y1, y2, ..., yN para y(t)

2829Mtodo de Euler 1 iteracaoPasso ht y(t)(t0, y0)Valor de y(t0+ h)

Slope

Interpretao grfica: primeiro passo do mtodo de Euler

Valor aproximado y1t0t130Mtodo de Euler 2 iteraoTamanho do passohValor verdadeiro y(x2) y2 Valor aproximadoy1ty(t)t1t2Segunda iterao do mtodo de Euler

Note que y2 e o Valor em x2 de uma reta que passa por (t1, y1) e que teminclinao f(t1,y1).

31Mtodo de Euler 2 iteraoTamanho do passohValor verdadeiro y(t2) y2 Valor aproximadoy1ty(t)t1t2

NO ESTAMOS USANDO

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

uma reta passando por y(t1) e com inclinao f(t1,y(t1)) poisNOS NO TEMOSy(t1)

Na 1 iterao obtivemos uma aproximao y1 para y(t1)32Mtodo de Euler iterao iTamanho do passohValor verdadeiro y(ti+1) yi+1, Valor aproximadoyitytiti+1Passo genrico do mtodo de Euler

Erros em EulerAssim, na n-sima iterao, gostaramos de aproximar yn+1 pelo valor em tn+1 = tn+h da reta tangente a y(t) no ponto (tn, y(tn))Entretanto, NO TEMOS y(tn)mas somente uma aproximao ynAssim, temos dois erros acumulando-seem cada iterao do mtodo de Euler.Existe um erro em aproximar y(tn)por yn , a n-sima iteraoAlm disso, gostaramos de ter f(tn,y(tn)) mas usamos f(tn,yn)

3334ExemploUma 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 interferncia, a equao 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 mtodo de Euler. Assuma um passo de tamanho h=240 segundos

35Soluo

Primeira iterao:

e a temperatura aproximada em

36Soluo - continuaoPara

Iterao 2:

e a temperatura aproximada em

37Soluo continuao A soluo exata da EDO e dada pela raiz da equao no-linear

A soluo (480) desta equao no-linear em t=480 segundos e

Bem diferente da aproximao:

38Comparao das solues exata e numrica

EulerStep, hq(480)Et|t|%4802401206030987.81110.32546.77614.97632.771635.4537.26100.8032.60714.806252.5482.96415.5665.03522.286439Efeito do tamanho do passo hTemperatura aos 480 segundos como uma funo do passo h

(valor exato)40Comparao com resultado exato

Apenas h=480, 240 e 12041Efeito do tamanho do passo h em Euler

Valor exatoMais um exemploConsidere a EDO

Como x0=0 ento xn=nh Iterao:

Usando h=0.1 e 0.001E comparando com a soluo exata temos a tabela ao lado

4243Erros no mtodo de EulerVimos que o mtodo de Euler PODE ter erros grandes. Para entender a ordem de grandeza desses erros, vamos fazer a expanso de Taylor em torno de xi

Os dois primeiros termos da serie de Taylor e o mtodo de Euler

O erro na aproximao e dado por

Isto e:Runge - KuttaEuler fez a seguinte aproximao

Que tal usar uma aproximao melhor para a integral?Por exemplo, podemos usar a regra do trapzio:

Neste caso, teremos ento a aproximao

E o algoritmo

44Runge-KuttaEncontramos a equao de iterao:

Existe um problema no entanto: yn+1 aparece dos dois lados da equao acima. No conseguimos isolar yn+1.Uma possibilidade e substituir yn+1 NO LADO DIREITO por sua aproximao baseada em Euler: yn+1 = yn + f(tn,yn)hEste e o metodo de Runge-Kutta de 2 ordem

45Runge Kutta de 2 ordemEquao de iterao:

ou simplesmente

onde

Assim, este e um mtodo de Euler com inclinao (s1+s2)/2

46Runge Kutta de 2 ordemE possvel uma interpretao grfica-geomtrica deste mtodo de Runge-Kutta. Temos

com Isto corresponde ao seguinte esquema em dois passos:Tome um passo preliminar de Euler com inclinao s1 em tn:

Com isto, obtenha uma segunda inclinao s2 em tn+hA atualizao de Euler realmente dada usa a mdia das inclinaes s1 em tn e s2 em tn+h

47Um segundo mtodo de Runge-KuttaO mtodo de Runge-Kutta que acabamos de estudar comeou aproximando uma integral pela regra do trapzio:

Podemos usar alguma outra regra: Simpson ou midpointVamos usar midpoint:

Neste caso

Note que y(t+h/2) no lado direito no e conhecido. Vamos usar Euler de novo para este valor.

482. Mtodo de Runge - KuttaTemos a aproximao

Usamos a aproximao de Euler para o termo y(tn+h/2):y(tn+h/2) y(tn)+h/2 * f(tn, yn)Substituindo a iterao para yn+1 temos

Este mtodo e conhecido como mtodo de Euler modificado ou mtodo do ponto mdio

492 metodo de Runge-KuttaTambm podemos ver este novo mtodo de Runge-Kutta como um processo em dois estgios.

Escrevemos

como

onde

502 metodo de Runge-KuttaTambm podemos ver este novo mtodo de Runge-Kutta como um processo em dois estgios.Tome um passo de Euler mas apenas com metade do comprimento do intervalo h/2

Isto corresponde ao tempo tn+h/2 = tn+1/2A seguir, de mais um passo de Euler de comprimento h usando a inclinao no ponto mdio (tn+1/2, yn+1/2)

51Resumo dos 2 mtodos de R-KPrimeiro: o mtodo clssico de 2 ordem de R-K (ou mtodo de Euler melhorado) yn+1 = yn + h (s1+s2)/2com

Segundo: Mtodo de Euler modificado (mtodo do ponto mdio)yn+1 = yn + h s2com

O que eles tem em comum?

52Comparando os dois R-KOs dois mtodos usam dois estgios intermedirios s1 e s2 para obter uma iterao.Os estgios correspondem a diferentes estimativas para a inclinao da soluo.

No mtodo clssico de RK (Euler melhorado) ns damos um passo completo yn+1 = yn + h (s1+s2)/2 tomando a media das inclinaes s1 em tn e s2 em tn+h

No mtodo de Euler modificado (ponto mdio), ns usamos s1 em tn para dar um meio-passo ate tn+h/2. A seguir, calculamos s2, a estimativa da inclinao no ponto mdio, e ento tomamos o passo completo yn+1 = yn + h s253ExemploConsidere a EDO

Euler modificado: yn+1 =yn+hs2

Temos s1=x2n+ y2ne s2=(xn+h/2)2+(yn+s1/2)2Exemplo numrico na tabela ao lado

54Runge-Kutta 2 ordem geralPodemos imaginar varias outras maneiras alternativas de calcular s1 e s2.O mtodo geral de Runge-Kutta de 2 ordem e da forma

onde

com (esta notao vem de uma teoria mais avanada ligada a mtodos implcitos)

Clssico RK (Euler melhorado): Euler modificado (ponto mdio): 1=0, 2=1 e 2= 21=1/2

55Tabela de ButcherE costume arranjar os coeficientes i, ij e i em uma tabela chamada tabela de Butcher

Onde 2 = 21 Para o mtodo ser de segunda ordem e ter certas propriedades desejveis impomos tambm as condies

56Tabela de ButcherRK Clssico (Euler melhorado)

RK : Euler modificado (ponto mdio)

RK: Mtodo de Heun

2 = 21Mtodo de Ralston

0002=3/421=3/401=1/32=2/35758ExemploUma 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 interferncia, a equao 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 mtodo EULER MELHORADO (ou metodo classico de Runge-Kutta de segunda ordem) Assuma um passo de tamanho h=240 segundos

59SoluoIterao 1:

60Soluo - continuaoIterao 2:

61Soluo - continuaoA soluo exata da EDO e dada pela soluo de uma equao no -linear:

A soluo para esta equao no-linear em t=480 segundos e

62Comparao com resultado exatosEuler melhorado (ponto mdio) para diferentes valores de h

63Efeito do tamanho do passo h Temperatura em t=480 segundos como uma funcao do tamanho do passo hPasso hq(480)Erro = Et|t|%4802401206030393.87584.27651.35649.91648.211041.463.3043.77622.34060.63219160.829.77560.583130.361450.097625

(exact)64Efeito do tamanho do passo h

Passohq(480)EulerEuler MelhoradoPonto MedioRalston4802401206030987.84110.32546.77614.97632.77393.87584.27651.35649.91648.211208.4976.87690.20654.85649.02449.78690.01667.71652.25648.6165Comparao de Euler e RK de 2a ordem

(exato)66Passo hEulerEuler ModificadoPonto MdioRalston4802401206030252.5482.96415.5665.03522.2864160.829.77560.583130.361450.09762586.61250.8516.58231.12390.2235330.5446.55373.10920.722990.15940

(exato)

Comparao de Euler e RK de 2a ordem67

Comparao de Euler e RK de 2a ordemModificadoPonto MedioRunge-Kutta de 4 ordemE o mais famoso mtodo de Runge-Kutta

com

E tabela de Butcher

68Para a provaMemorizar apenas os dois mtodos mais simples de Runge-Kutta: Euler melhorado (RK clssico de 2 ordem)Euler modificado (ponto mdio)69Chart112001200901.08106.09775.09110.32699.66647.57

h=240Exact Solution(K)Time, t(sec)Temperature,

Sheet101200120901.08240775.09360699.66480647.5701200240106.09480110.32

Sheet1

h=240Exact Solution(K)Time, t(sec)Temperature,

Sheet2

Sheet3

Chart11200120012001200653.05106.09-987.81901.08607.03110.32775.09573.22699.66546.77647.57

Exact solutionh=120h=240h=480(K)Time, t (sec)Temperature,

Sheet101200120653.05240607.03360573.22480546.7701200240106.09480110.3201200480-987.8101200120901.08240775.09360699.66480647.57

Sheet1

Exact solutionh=120h=240h=480(K)Time, t (sec)Temperature,

Sheet2

Sheet3

Chart1-987.81647.57110.32546.78614.97632.77

(K)Step size, h (s)Temperature,

Sheet1480-987.81240110.32120546.7860614.9730632.770647.57

Sheet1

(K)Step size, h (s)Temperature,

Sheet2

Sheet3

Chart112001200120012001200653.04903.521007950.43901.08607.03779.99856.095811.79775.09573.22704.13757.15726.25699.66546.77651.35690.2667.71647.57

AnalyticalRalstonMidpointEulerHeun(K)Time, t (sec)Temperature,

Sheet1TimeEulerHeunMidpointRalstonExact012001200120012001200120653.04903.521007950.43901.08240607.03779.99856.095811.79775.09360573.22704.13757.15726.25699.66480546.77651.35690.2667.71647.57

Sheet1

AnalyticalRalstonMidpointEulerHeun(K)Time, t (sec)Temperature,

Sheet2

Sheet3