Terceira_unidade Numerico

Material sujeito a correções

1

Interpolação......................................................................................................................................... 2 1.1 Interpolação Linear ..................................................................................................................... 2 1.2 Diferenças divididas.................................................................................................................... 3 1.3 Diferenças Finitas ....................................................................................................................... 5 1.4 Interpolação Parabólica............................................................................................................. 11 2 Integração Numérica......................................................................................................................... 12 2.1 Métodos de Newton-Cotes........................................................................................................ 12 2.1.1 Método Dos Trapézios...................................................................................................... 15 2.1.2 Método de Simpson (1ª regra) .......................................................................................... 17 2.1.2.1 Outra dedução da 1ª regra de Simpson ......................................................................... 18 2.1.3 Segunda Regra de Simpson .............................................................................................. 20 3 Diferenciação Numérica ................................................................................................................... 21 3.1 Solução numérica de um PVI de 1ª ordem ............................................................................... 24 3.2 Método de Euler........................................................................................................................ 25 3.3 Método de Runge-Kutta............................................................................................................ 29

Prof. José Augusto Lucas Matos

Página 1 de 31

Material sujeito a correções

1 Interpolação 1.1 Interpolação Linear Seja a função y = f ( x ) dada pele tabela abaixo:

y0 = f ( x 0 ) y2 = f ( x1 ) y3 = f ( x2 ) .................. yn = f ( xn ) tabela esta obtida através da expressão analítica de experimentais. Considerando a figura abaixo, temos:

____

Temos

DB ____

EC

f ( x ) , ou de determinações

____

=

BA ____

CA

, ou

f ( x ) − f ( x1 ) x − x1 = ∴ f ( x2 ) − f ( x1 ) x 2 − x1

f ( x ) = f ( x1 ) +

f ( x2 ) − f ( x1 ) ( x − x1 ) x 2 − x1

Exemplo: Interpolar o seno de 22º ( f ( 22º ) na tabela abaixo:

0 − 0,00000 10 − 0,17365 20 − 0,34202 30 − 0,50000 40 − 0,64279  22 − 20  sen22 = 0,34202 +  (0,5 − 0,34202) = 0,37461  30 − 20  Prof. José Augusto Lucas Matos

Página 2 de 31

Material sujeito a correções

1.2 Diferenças divididas x x0 x1 x2 x3 x4

------

f ( x) f ( x0 ) f ( x1 ) f ( x2 ) f ( x3 ) f ( x4 )

1ª Dif . f [ x0 , x1 ] f [ x1 , x2 ] f [ x2 , x 3 ] f [ x3 , x4 ]

2ª Dif .div . f [ x0 , x1 , x2 ] f [ x1 , x 2 , x3 ] f [ x 2 , x 3 , x4 ]

3ª Dif . Dividida f [ x0 , x1 , x 2 , x3 ] f [ x1 , x 2 , x 3 , x4 ]

4ª Dif . Dividida f [ x0 , x1 , x2 , x 3 , x4 ]

A primeira diferença dividida é definida como:

f [ x0 , x1 ] =

f ( x0 ) − f ( x1 ) x0 − x1

analogamente, as diferenças de ordem superior são definida como:

f [ x0 , x1 , x2 ] =

f [ x0 , x1 , x2 , x 3 ] =

f [ x0 , x1 ] − f [ x1 , x 2 ] x0 − x1

f [ x0 , x1 , x2 ] − f [ x1 , x2 , x3 ] x0 − x 3

.......................................................................................

f [ x0 , x1 ,... xn − 1 , xn ] =

f [ x0 , x1 ,...., xn − 1 ] − f [ x1 , x2 ,...., xn ] x0 − x4

Desenvolvendo as expressões das diferenças divididas dos numeradores, resulta:

p

f [ x0 , x1 , x2 , x 3 , x4 ] =

∑ (x

k =0

k

f ( xk ) − x0 )( xk − x1 ).....( xk − x k − 1 )( xk − x k + 1 )....( xk − x p ) Equação 1

Prof. José Augusto Lucas Matos

Página 3 de 31

Material sujeito a correções

Sendo por definição:

f [ x0 , x1 ] =

f ( x0 ) − f ( x1 ) x0 − x1

Resulta:

f ( x ) = f ( x0 ) + ( x − x0 ) f [ x , x0 ] Como:

f [ x , x0 ] = f [ x0 , x1 ] + ( x − x1 ) f [ x , x0 , x1 ] Se substituirmos teremos:

f ( x ) = f ( x0 ) + ( x − x0 ) f [ x0 , x1 ]( x − x0 )( x − x1 ) f [ x , x0 , x1 ] Generalizando, temos:

f ( x ) = f ( x0 ) + ( x − x0 ) f [ x0 , x1 ]( x − x0 )( x − x1 ) f [ x , x0 , x1 ] + + ( x − x0 )( x − x1 )( x − x2 ) f [ x0 , x1 , x2 , x3 ] + .......................... + + ( x − x0 )( x − x1 )( x − x2 )....( x − xn − 1 ) f [ x0 , x1 , x 2 ,........., xn ] + + ( x − x0 )( x − x1 )( x − x2 )....( x − xn ) f [ x0 , x1 , x2 ,........., x n ] Equação 2

Que é fórmula de interpolação, por diferenças divididas, de Newton, onde a última parcela se refere ao erro de arredondamento. Das Equações 1 e 2, é obtida a Fórmula de Lagrange: n

f (x) =

∑

Lk ( f ( x k )

k =0

onde:

Lk ( x ) =

Prof. José Augusto Lucas Matos

( x − x0 )( x − x1 ).........( x − x k − 1 )( x − x k + 1 )..........( x − xn ) ( x k − x0 )( x k − x1 )......( x k − x k − 1 )( x − xk + 1 )......( x k − xn )

Página 4 de 31

Material sujeito a correções

Exemplo: Seja a tabela:

x 1 2 4 5

P( x) 10 15 20 30

Interpolar o valor de P ( 3) .

P ( 3) =

( 3 − 2)( 3 − 4)( 3 − 5) ( 3 − 1)( 3 − 4)( 3 − 5) .10 + .15 + (1 − 2)(1 − 4)(1 − 5) ( 2 − 1)( 2 − 4)( 2 − 5)

+

( 3 − 1)( 3 − 2)( 3 − 5) ( 3 − 1)( 3 − 2)( 3 − 4) .20 + .30 = (4 − 1)(4 − 2)(4 − 5) (5 − 1)(5 − 2)(5 − 4)

= −1,66666 + 10 + 13,33333 − 5 = 16,66666...

1.3 Diferenças Finitas Uma diferença finita de uma função f ( x ) é o valor da função no ponto x1 menos o valor da função no ponto x 2 . Algebricamente é representada como f ( x1 ) − f ( x 2 ) . Geometricamente é:

Onde:

∆f ( x i ) → é chamada de difernça ascendene(ou progressiva ); ∇f ( xi ) → é chamada de difernça descendente(ou regressiva) e

δf ( x i ) é chamada de diferença central. Prof. José Augusto Lucas Matos

Página 5 de 31

Material sujeito a correções

Considerando uma função f ( x ) , tabulada com espaçamentos iguais em seus valores de x, teremos:

x 0 1 2 3 4 5

f ( x) 0 10 20 50 60 100

∆ 10 10 30 10 40

∆2 0 20 − 20 30

∆3 20 − 40 50

∆4 − 60 90

∆5 150

As diferenças ascendentes (as quais utilizaremos no curso) na tabela anterior para o ponto x0 são:

∆f (0) = 10 ∆2 f (0) = 0 ∆3 f (0) = 20 ∆4 f (0) = −60 ∆5 f (0) = 150 Para o ponto x = 3 , teremos:

Difernça ascendente ∆f (3) = 10 e Difernça descendente ∆f (3) = 30.

Prof. José Augusto Lucas Matos

Página 6 de 31

Material sujeito a correções

A fórmula de Newton pode ser apresentada como:

f ( x ) = a1 + a 2 ( x − x1 ) + a3 ( x − x1 )( x − x 2 ) + a4 ( x − x1 )( x − x2 )( x − x 3 ) + + ............................................... + an − 1 ( x − x1 )( x − x 2 ).......( x − xn ) se substituirmos x por x1 , x 2 , x 3 ,......, xn teremos:

f ( x1 ) = a1 f ( x2 ) = a1 + a2 ( x2 − x1 ) f ( x3 ) = a1 + a2 ( x3 − x1 ) + a3 ( x3 − x1 )( x3 − x2 ) ........................................................................... f ( xn + 1 ) = a1 + a 2 ( x n + 1 − x1 ) + a 3 ( x n + 1 − x1 )( xn + 1 − x 2 ) + ..... + + a n + 1 ( xn + 1 − x1 )( x n + 1 − x2 )...................................( xn + 1 − xn ) Os coeficiente a1 , a 2 a 3 ,......., a n + 1 podem ser achados por substituições sucessivas. Quando consideramos os intervalos iguais, temos:

h = x 2 − x1 = x3 − x 2 = ......... = x n + 1 − x n Substituindo h no sistema anterior vem:

f ( x1 ) = a1 f ( x1 + h) = a1 + a 2 ( h) f ( x1 + 2h) = a1 + a 2 ( 2h) + a3 ( 2h)( h) ........................................................................... f ( x1 + nh) = a1 + a 2 ( nh) + a3 ( nh)[( n − 1)h)] + ..... + a n + 1 ( nh)[( n − 1)h]...........( h)

Prof. José Augusto Lucas Matos

Página 7 de 31

Material sujeito a correções

que, quando resolvido resulta em:

a1 = f ( x1 ) f ( x1 + h) − f ( x1 ) h f ( x1 + 2h) − 2 f ( x1 + h) + f ( x1 ) a3 = 2h 2 ...........................................................

a2 =

É então possível simplificas estas expressões, notando que os denominadores podem ser expressos por:

h i − 1 ( i − 1)! Donde obtemos:

a1 = f ( x1 ) a2 =

∆f ( x1 ) 1! h

∆2 f ( x1 ) 2! h2 ....................... a3 =

∆n f ( x1 ) an +1 = n! hn Que para o coeficiente a i , é:

∆i − 1 f ( x1 ) ai = ( i − 1)! hi − 1 Logo:

∆f ( x1 ) ∆2 f ( x1 ) f ( x ) = f ( x1 ) + ( x − x1 ) + ( x − x1 )( x − x 2 ) + 1! h 2! h2 ................. +

∆3 f ( x1 ) ( x − x1 )( x − x2 )( x − x 3 ) + ................ + 3! h3

∆n f ( x1 ) ....................... + ( x − x1 )( x − x2 )................( x − xn ) n! hn

Prof. José Augusto Lucas Matos

Página 8 de 31

Material sujeito a correções

Trocando a variável x por h , usando convenientemente a expressão: x = x1 + uh , Teremos:

( x − x1 ) = uh ( x − x2 ) = ( x1 + uh) − ( x1 + h) = h( u − 1) ( x − x3 ) = ( x1 + uh) − ( x1 + 2h) = h( u − 2) .................................................................... ( x − xn ) = ( x1 + uh) − [ x1 + ( n − 1)h] = h[u − ( n − 1)] Substituindo, teremos:

∆f ( x1 ) ∆2 f ( x1 ) f ( x ) = f ( x1 ) + u+ u( u − 1) + 1! 2! + ............ +

∆n f ( x1 ) u( u − 1)( u − 2).......[u − ( n − 1)] n!

Exemplo: Consideremos a tabela de diferenças finitas, abaixo:

x 1 3 5 7 9

f ( x) 30 50 60 90 140

∆ 20 10 30 50

∆2 − 10 20 20

∆3 30 0

∆4 − 30

Aplicando a fórmula de Newton para intervalos iguais, teremos para o ponto x = 4 :

x = x1 + uh , e como h = 2, x1 = 1 e x = 4 , vem: u=

4−1 ∴ u = 1,5 2

Então:

20 ( −10) 30 .1,5 + .1,5(1,5 − 1) + .1,5(1,5 − 1)(1,5 − 2) + 1! 2! 3! ( −30) + .1,5(1,5 − 1)(1,5 − 2)(1,5 − 3) = 4! ( −10) 30 ( −30) = 30 + 20(1,5) + .1,5.0,5 + .1,5.0,5( −0,5) + .1,5.0,5( −0,5)( −1,5) = 2 6 24 = 30 + 30 − 3,75 − 1,875 − 0,703125 = 53,671875

f (1,5) = 30 +

Prof. José Augusto Lucas Matos

Página 9 de 31

Material sujeito a correções

Resolvendo o mesmo problema por diferenças divididas teremos:

x 1 3 5 7 9

f ( x) 30 50 60 90 140

∆ 10 5 15 25

∆2 − 1,25 2,5 2,5

∆3 0,625 0

∆4 − 0,078125

140 − 90 9−7 Aplicando a fórmula de Newton para diferenças divididas, temos para o ponto x = 4 :

f ( x ) = ( 30) + (10)(4 − 1) + ( −1,25)(4 − 1)(4 − 3) + (0,625)(4 − 1)(4 − 3)(4 − 5) + + ( −0,078125)(4 − 1)(4 − 3)(4 − 5)(4 − 7) = = 30 + 30 − 3,75 − 1,875 − 0,73125 = = 53,671875

Prof. José Augusto Lucas Matos

Página 10 de 31

Material sujeito a correções

1.4 Interpolação Parabólica Dada a tabela:

x 20 25 30

P( x) 0,34202 0,42252 0,50000

calcular o seno de 22º. Se usarmos a expressão P ( x ) = a1 x 2 + a 2 x + a 3 e nela substituirmos os valores dos pontos da tabela, vem:

400a1 + 20a 2 + a 3 = 0,34202 625a1 + 25a 2 + a 3 = 0,42262 900a1 + 30a 2 + a 3 = 0,50000 E resolvendo o sistema, teremos:

a1 = −0,0000644,

a 2 = 0,019018 e a 3 = −0,01258

logo, P ( 22) = 484.0,0000644 + 22.0,019018 − 0,01258 ∴ P ( 22) = 0,3746464 Outro sistema pode ser obtido, assumindo a parábola como:

P ( x ) = a1 + a 2 ( x − 20) + a 3 ( x − 20)( x − 25) Substituindo os valores tabulados, temos:

0,34202 = a1 0,42262 = a1 + 5a 2 0,50000 = a1 + 10a2 + (10)(5)a 3 Resolvendo por substituição sucessiva vem:

a1 = 0,34202 a 2 = 0,01612 a 3 = −0,0000644 P ( 22) = 0,34202 + 0,1612( 22) + ( −0,0000644( 2)( −3)) = 0,37429864 Prof. José Augusto Lucas Matos

Página 11 de 31

Material sujeito a correções

2 Integração Numérica 2.1 Métodos de Newton-Cotes Surge da idéia de calcular f ( x ) em certo número finito de pontos, fazendo uma interpolação polinomial e integrando o polinômio ao invés de f ( x ) . Isto é: Se são conhecidos n + 1 pontos de f ( x )

( x0 , f ( x0 ), ( x1 , f ( x1 ), ( x2 , f ( x2 ), ....... , ( x n , f ( xn + 1 ) b

e como o objetivo é calcular a integral de f ( x ) , τ = ∫ f ( x )dx e só dispomos de pontos da a

função, obtemos a expressão do polinômio P ( x ) de grau máximo n que passe por estes pontos e aproximamos o valor da integral para:

b

τ = ∫a P ( x )dx Equação 3

A obtenção de P ( x ) pode ser feita pelos métodos de Newton ou pelo método de Lagrange. Quando os pontos constituem uma tabela de pontos equivalentes, este processo é conhecido como método de Newton-Cotes, no caso ainda de mais particular em que x0 − a e x n − b , obtemos as fórmulas fechadas de Newton-Cotes. No caso geral, o problema pode ser resolvido com a Equação 4, onde P ( x ) é obtido através de uma das fórmulas citadas no parágrafo anterior.

Prof. José Augusto Lucas Matos

Página 12 de 31

Material sujeito a correções

Nos métodos fechados de Newton-Cotes, temos:

x0 = a ,

xn = b

x i + 1 − x i = h com h > 0 e i = 0,1,2,3,......, n e podemos dizer que aos pontos conhecidos provêm da partição do intervalo de integração em n partes iguais, ou seja:

xx − x n

h= A integral aproximada seria x

τ ' = ∫x n Pn ( x )dx 0

Equação 5

onde Pn ( x ) é um polinômio de interpolação obtido por diferenças finitas. Finalmente observamos que, como se trata de uma tabela de pontos eqüidistantes, podemos fazer uma mudança de variáveis de x para u , como método de Newton por diferenças finitas. Teremos então; n

τ ' = ∫0 Pn ( u)du Equação 6

sendo:

u=

x − x0 h

Esta última fórmula Equação 4 resolve o problema dos métodos “fechados” de Newton-Cotes. Estes métodos usuais, se entanto, se baseiam na partição do intervalo de integração em n m subintervalos que contém n + 1 pontos da tabela, sendo dois nas extremidades. Aplicando-se um método “fechado” de Newton-Cotes a cada subintervalo, o valor de τ ' será a soma de todas as integrais ao do intervalo [a , b] , Obtêm-se assim as fórmulas compostas cuja precisão se revela plenamente satisfatória em casos ordinários.

Prof. José Augusto Lucas Matos

Página 13 de 31

Material sujeito a correções

Por exemplo:

Onde (P2= P2i ) e:

P21 é o polinômioi nterpolante de grau 2 do 1º intervalo, .. .. P2i é o polinômioi nterpolante de grau 2 do i - ésimo intervalo. .. Assim:

τ ' = τ 1' + τ 2' + ......... + τ n' 2 Donde: x

x

x

τ ' = ∫x 2 P21 ( x )dx + ∫ x 4 P21 ( x )dx + ...... + ∫ x n P21 ( x )dx 0

Prof. José Augusto Lucas Matos

2

n− 2

Página 14 de 31

Material sujeito a correções

2.1.1 Método Dos Trapézios No caso de m = 1 , teremos o método dos trapézios, onde cada integral parcial é aproximada pela área de um trapézio.

Calculando as expressões τ i' , (1,2,............., n) pela Equação 4, vem:

 u P1i ( u) = f i − 1 +   ∆f i − 1 ou P1i ( u) = f i − 1 + u∆f i − 1 1 Por outro lado, para cada valor τ i' temos: u =

x − xi −1 e usando a Equação 4 vem: h

1

τ i' = h∫0 ( f i −1 + u∆f i −1 )du = 1

u2  = h( f i −1u + ∆f i −1  = 2 0

h = hfi −1 + ∆f i −1 2 h 2

' Como: ∆f i − 1 = f i − f i − 1 , vem: τ i = ( f i + f i −1 )

O valor de τ ' que chamaremos de τ t , é:

h 2

τ t = ( f 0 + 2 f1 + 2 f 2 + ...... + 2 f n −1 + f n ) Equação 7

Que é a fórmula dos trapézios para a aproximação de integrais numericamente. Prof. José Augusto Lucas Matos

Página 15 de 31

Material sujeito a correções

Quando n − 1 temos a fórmula fechada de Newton-Cotes para um subintervalo.

τ '=

h [ f 0 + f1 ] 2

onde podemos notar que quanto maior o valor de n , mais exata será a aproximação, considerando-se apenas o erro intrínseco ao processo.

Prof. José Augusto Lucas Matos

Página 16 de 31

Material sujeito a correções

2.1.2 Método de Simpson (1ª regra) Uma aproximação, bem melhor, da integral pode ser obtida quando usamos parábolas do 2º grau interpoladas em cada subintervalo.

Neste caso m = 2 e n deve ser par. Seja a Equação 5 τ ' = ∫

x2i x2i−2

P2' ( x )dx onde temos u =

x − x2i − 2 , então: h

 u  u P2i ( u) = f 2 i − 2 +   ∆f 2 i − 2 +   ∆f 2 i − 2  1  2 Integrando, teremos: 2

τ i' = h∫0 ( f 2i − 2 + u∆f2i − 2 +

u(u − 1) 2 )∆ f 2i − 2 )du = 2

   u3 u2  u2 = h f 2i − 2u + ∆f 2i − 2 +  − ∆2 f 2i − 2  2  6 4  

2

= 0

1   = h f 2i − 2 + 2∆f 2i − 2 + ∆2 f2i − 2  3   como:

∆f 2 i − 2 = f 2 i − 1 − f 2 i − 2 e ∆2 f 2 i − 2 = f 2 i − 2 f 2 i − 1 + f 2 i − 2

Prof. José Augusto Lucas Matos

Página 17 de 31

Material sujeito a correções

vem:

τ i' = h( 2 f 2 i − 2 + 2 f 2 i − 1 − 2 f 2 i − 2 +

f 2i 2 f − f 2i −1 + 2i − 2 ) 3 3 3

donde:

τ i' =

h [ f 2i + 4 f 2i −1 + f 2i − 2 ] 3

O valor de τ i' que chamaremos de τ s , será: n

2

h [ f 2 i + 4 f 2i −1 + f 2 i − 2 ] 3 i =1

τs = ∑ isto é:

τ i' =

h [ f 0 + 4 f1 + 2 f 2 + 4 f 3 + ......... + 2 f n − 2 + 4 f n −1 + f n ] 3

Que é conhecida como fórmula da 1ª regra de Simpson e só pode ser aplicada para para um número de pontos ímpar.

2.1.2.1 Outra dedução da 1ª regra de Simpson Consideremos um polinômio P2 ( x ) que passe pelos pontos x0 = 0, x1 = 1 e x1 = 2 .

que é um polinômio de segundo Grau.

Prof. José Augusto Lucas Matos

Página 18 de 31

Material sujeito a correções

O polinômio tem a seguinte expressão:

P ( x ) = L0 ( x ) y0 + L1 ( x ) y1 + L2 ( x ) y2 onde:

L0 ( x ) =

( x − 1)( x − 2) 1 2 = ( x − 3 x + 2) (0 − 1)(0 − 2) 2

( x − 0)( x − 2) = x(2 − x ) = − x 2 + x (1 − 0)(1 − 2) ( x − 0)( x − 1) 1 2 L2 ( x ) = = ( x − x) ( 2 − 0)( 2 − 1) 2 L1 ( x ) =

daí: 2

∫0

2

P ( x )dx = ∫ [L0 ( x ) y0 + L1 ( x ) y1 + L2 ( x ) y2 ]dx 0

2

2

2

2

∫0 P ( x )dx = ∫0 L0 ( x ) y0dx + ∫0 L1 ( x ) y1dx + ∫0 L2 ( x ) y2dx 2

2

2

2

∫0 P ( x )dx = y0 ∫0 L0 ( x )dx + y1 ∫0 L1 ( x )dx + y2 ∫0 L2 ( x )dx 2

∫0

P ( x )dx = y0 ∫

21

0

2

2

21

0

0

( x 2 + 3 x + 2)dx + y1 ∫ ( − t 2 + 2t )dx + y2 ∫ 2

2

2

( t 2 − t )dx 2

 1  t 3 3t 2  1  t 3 t 2   − t 3 2t 2      −  P ( x ) dx y + 2 t + y + y = − +  0  1 2  ∫0   2 3 2 3 2 2   3 2   0   0  0   2

2

1  8

∫0 P ( x )dx = y0  2  3 − 2

∫0 P ( x )dx =

1  8 12    8  + 4   + y1  − + 4 + y2   − 2   2   3   2  3

y0 4 y1 y2 1 + + = [ y0 + 4 y1 + y2 ] 3 3 3 3

No caso geral em que x0 ≠ 0, h ≠ 1, x1 = x0 + h e x 2 = x1 + h . Por meio de uma mudança de variável recaímos no caso anterior ( x0 = 0 e h = 1) . 2

1

∫0 P ( x )dx = 3 [ y0 + 4 y1 + y2 ] Prof. José Augusto Lucas Matos

Página 19 de 31

Material sujeito a correções

Exemplo: Seja calcular a integral numérica da curva dada pela tabela abaixo ente os pontos x = 2 e x = 10 pelo método se Simpson 1ª regra: x 2 4 6 8 10

y 28 54 38 27 12

Os intervalos entre os pontos x são constantes. Logo h = 2

τ i' =

2 [28 + 4(54) + 2(38) + 4(27) + 12] = 440 3

2.1.3 Segunda Regra de Simpson De forma análoga a dedução da 1ª regra a expressão da segunda regra de Simpson é obtida aproximando-se a função f ( x ) por um polinômio interpolador de Gregory-Newton do 3º grau

P3 ( x ) . b

b

Ι = ∫ f ( x )dx = ∫ P3 ( x )dx a a b u( u − 1) 2 u( u − 1)( u − 2) 3  Ι = ∫  y0 + uy0 + ∆ y0 + ∆ y0 dx a 2! 3!   3

  u3 u2  2  u4 u3 u2  3  u2    ∆ y0  Ι = h uy0 + ∆ y0 +  − −  ∆ y0 +  +  2 6 4 24 6 4   0     Que resulta em:

Ι=

3h [ y0 + 3 y1 + 3 y2 + y3 ] 8

Que também é conhecida como a regra dos

3 . 8

A mesma só pode ser aplicada para os seguintes números de pontos: 4, 7, 10, 13, 16, 19,...... .

Prof. José Augusto Lucas Matos

Página 20 de 31

Material sujeito a correções

3 Diferenciação Numérica Equações diferenciais ordinárias ocorrem com muita freqüência na descrição de fenômenos da Natureza. Ex.: Podemos supor que sob determinadas condições ambientais, a taxa de crescimento de uma colônia de bactérias seja proporcional ao número de indivíduos num dado tempo; se y(t ) for o número de indivíduos num tempo t , temos a equação y( t ) = k ( y( t )) . Nem sempre é possível obter uma solução analítica para uma Equação Diferencial Ordinária (EDO) e nestes casos os métodos numéricos são a saída para encontrar uma solução aproximada. Ex.: As equações y' = x 2 + y 2 e y' ' = 6 y 2 + x não podem ser resolvidas em termos de funções elementares. Consideremos a equação y = x 2 , que é a equação de uma parábola. Consideremos, também,

dy = 2x . dx Esta última equação, à sua maneira descreve uma parábola. Ou seja: ao invés de termos, explicitamente o valor de y para cada x , ela descreve a inclinação (ou direção) da parábola em cada ponto da mesma.

dy = 2x dx é chamada de uma EDO. Diferencial, devido a relação entre derivadas e ordinária pois a derivada é comum (não uma derivada parcial

∂y ). ∂x

Além disso, como y = x 2 satisfaz a exigência de que

y = x 2 é uma solução da EDO dy

Prof. José Augusto Lucas Matos

dy

= 2 x para cada x, dizemos que dx

= 2x . dx

Página 21 de 31

Material sujeito a correções

Do cálculo conhecemos a forma como se apresenta uma equação diferencial de ordem n.

(

y n = f x , y , y (1) , y ( 2 ) ,..........., y ( n )

)

onde:

dl y y = l dx l

l = 1,2,3,...., n

 x ∈ [a, b ]   y[a , b] → R

1. y ( 2 ) = 3 y (1) − 2 y é uma EDO de segunda ordem com f ( x , y , y (1) ) = 3 y (1) − 2 y . Associadas a esta equação podem existir condições cujo número coincide com a ordem da EDO. Se tais condições se referem a um único valor de x , temos um problema de valor inicial (PVI), caso contrário temos um problema de valores de contorno. 2. O PVI seguinte é de 2ª ordem:

 y ( 2 ) = 3 y (1 ) − 2 y   y ( 0) = −1  (1 )  y ( 0) = 0

Prof. José Augusto Lucas Matos

Equação 8

Página 22 de 31

Material sujeito a correções

Trataremos de métodos numéricos para conseguir os valores de y ( x ) , em pontos distintos daqueles das condições iniciais associadas aos PVI’s. O tipo de PVI será o mais simples: (1ª ordem).

 y (1 ) = f ( x , y )   y ( x0 ) = y0 = η

sendo η um número dado.

PVI’s de ordem superior podem ser reduzidos a PVI’s de 1ª ordem através variáveis auxiliares.

 y ( 2 ) = 3 y (1 ) − 2 y  EX.:  y(0) = −1 , para reduzir à 1ª ordem, basta fazer:  (1 )  y ( 0) = 0

 y' = z e teremos :  ( 2) (1 ) , logo: y = z e  (0) z = 0

 y (1 ) = z  ( 0)  y = −1  (1 )  z = 3z − 2 y  (0) z = 0 Este problema tem solução única se a função f ( x , y ) satisfaz: 1. É definida e contínua na faixa: a ≤ x ≤ b, com ( −∞ < y < ∞ ) , onde a e b são finitos. 2. Existe uma constante L tal que para todo x ∈ [a , b] e todo par de números y e y * temos f ( x , y ) − f ( x , y * ≤ L y − y * . Então existe uma função y ( x ) , satisfazendo: •

y ( x ) é contínua e diferenciável para todo x ∈ [a , b] ;

• •

y (1) ( x ) = f ( x , y( x )), para x ∈ [a , b] e y(a ) = η , onde η é um número dado.

Observemos que a condição 2 acima será satisfeita se f ( x , y ) tem derivada contínua em relação a y e limitada na faixa em questão, pois então, do teorema do valor médio: _ ∂f f ( x, y) − f ( x, y ) = ( x , y )( y − y * ) ∂y *

onde y é um valor entre y e y * . Prof. José Augusto Lucas Matos

Página 23 de 31

Material sujeito a correções

A existência de

∂f não é necessária para que a equação do item 2 seja satisfeita. É bom ∂y

notar que L na é em geral, possível de ser computado.

3.1 Solução numérica de um PVI de 1ª ordem  y ( 2 ) = 3 y (1 ) − 2 y  Suponhamos que o PVI  y(0) = −1 satisfaça as condições de existência e unicidade.  (1 )  y ( 0) = 0 Para obtermos sua solução numérica, tomemos “m” subintervalos de [a , b], com ( m ≥ 1) , e (b − a ) façamos x j = x0 + jh onde h = para j = 0,1,2,............, m e x j ∈ [a , b] . m

A solução numérica y m ( x ) é a função linear por partes, cujo gráfico é uma poligonal com vértice nos pontos ( x j , y j ) onde y j foi calculado usando algum método numérico

b−a para n = 0,1,2,........ , teremos uma seqüência de 2n funções poligonais { yn ( x )} que convergem uniformemente para a solução y ( x ) do PVI. Se por exemplo m = 2 n , então hn =

Convém observar que os métodos numéricos vistos em seguida têm por objetivo calcular os vértices { yo , y1 , y2 ,........., yn ,} .

Prof. José Augusto Lucas Matos

Página 24 de 31

Material sujeito a correções

Convencionou-se usar a notação y( x j ), j = 0,1,2,...., m para indicar a solução exata o PVI nos pontos x j ∈ Ιh . A notação y( x j ) ≅ y j significa que y j é aproximação para y ( x j ) , com x j ∈ Ιh .

3.2 Método de Euler Seja o PVI:

 y (1 ) = f ( x , y )   y ( x0 ) = y0 = η Desejamos

as

aproximações

y1 , y2 , y3 ,.........., yn ,

para

as

soluções

exatas

y( x1 ), y( x2 ), ,.........., yn . Observemos a figura abaixo:

Prof. José Augusto Lucas Matos

Página 25 de 31

Material sujeito a correções

Como desconhecemos o valor de y( x1 ) , tomemos y1 como aproximação de y( x1 ) . Para isto tracemos a tangente T á curva y( x ) no ponto ( x0 , y( x0 )) , cuja equação é:

y ( x ) − y ( x 0 ) = ( x − x 0 ) y' ( x 0 ) Fazemos x = x1 e lembrando que y' ( x0 ) = f ( x0 , y( x0 )) e y1 ≅ y ( x1 ) , temos:

y1 = y0 + hf ( x0 , y( x0 ) ) O erro cometido na aproximação de y( x1 ) por y1 é e1 = y1 − y( x1 ) (diferença entre a solução numérica e a solução exata. Generalizando, teremos:

(

y j + 1 = y j + hf x j − y j

)

com j = 0,1,2,......., m − 1

Equação 9

Cujo erro é:

e j + 1 = y j + 1 − y( x j + 1 ) O método de EULER consiste em calcular recursivamente as fórmulas:

( A )  y0 = y ( a ) = η  (B )  y j + 1 = y j + hf ( x j , y j )

j = 0,1,2,...........m − 1

Estas fórmulas admitem várias interpretações analíticas:

(

1. Aproximando-se a derivada que aparece no PVI no ponto x j , y j

) por uma diferença

dividida vem:

y j +1 − y j h

= f (xj, yj )

E resolvendo-se para y j + 1 , obtemos a Equação 2. 2. Integrando-se y' ( t ) = f (t , y( t )) entre x e x + k obtemos:

y ( x + k ) − y( x ) = ∫ Prof. José Augusto Lucas Matos

x+r x

f ( t , ( y( t ))dt Página 26 de 31

Material sujeito a correções

Fazendo-se x = x j e k = h ,

y( x j + 1 ) − y( x j ) = ∫

x j +1 xj

f ( t , ( y( t ))dt

Aproximando-se a integral de forma bem grosseira: tamanho do intervalo vezes o valor do integrando à esquerda e identificando-se y ( x j ) com y j obtém-se a Equação 2. 3. Supondo uma expansão da solução y ( x ) pela série de Taylor em torno de x j ,

1 y( x j + h ) = y( x j ) + hf x j , y( x j ) + h2 y ( 2 ) ( x j ) + .... 2

(

)

Truncando-se a série após o termo em h e identificando-se y = x j com x j , teremos Equação 2. Ex.: Achar as aproximações para a solução do PVI:

 y' = x − y + 2 de [0;1] com h = 0,1   y ( 0) = 2 onde: x0 = 0 , y0 = 2 , a = 0 , b = 1 , m =

1−0 ⇒ m = 10 . 0,1

(

)

Usando-se y( x j + h ) = y( x j ) + hf x j , y( x j ) ,

j = 0,1,2,...,9

teremos: j = 0,

y1 = y0 + hf ( x0 , y0 )) = y0 + h( x0 - y 0 + 2) = 2 + 0,1 f (0;2) = 2 + 0,1(0 − 2 + 2) = 2 x1 = x0 + h ⇒ x1 = 0 + 0,1 ⇒ x1 = 0,1 j = 1, y2 = y1 + hf ( x1 , y1 )) = y1 + h( x1 - y 1 + 2) = 2 + 0,1 f (0;2) = 2 + 0,1(0,1 − 2 + 2) = 2,01 x 2 = x0 + 2h ⇒ x 2 = 0 + 2(0,1) ⇒ x2 = 0,2 j = 2, y3 = y2 + hf ( x 2 , y2 )) = y2 + h( x 2 - y 2 + 2) = 2 + 0,2 f (0,2;2) = 2 + 0,1(0,2 − 2 + 2) = 2,029 Até o valor j = 9 , teremos a tabela a seguir: Prof. José Augusto Lucas Matos

Página 27 de 31

Material sujeito a correções

j 0 1 2 3 4 5 6 7 8 9 10

xj 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

yj 2 2 2,01 2,029 2,0561 2,09049 2,131441 2,1782969 2,2304672 2,2874205 2,348684

y( x j )

y j − y( x j )

2 2,004837 2,018731 2,040818 2,070320 2,106531 2,148812 2,196585 2,249329 2,306570 2,367879

-0,004837 -0,0083731 -0,011818 -0,014220 -0,016041 -0,017371 -0,018288 -0,118862 -0,019149 -0,0192201

Olhando-se a tabela notamos que o erro cresce em valor absoluto devido à propagação. Para determinarmos uma expressão matemática para o erro usamos a fórmula de Taylor e desenvolvendo y ( x ) em torno de x0 , vem:

y ( x ) = y ( x0 ) +

x − x 0 (1 ) y ( x0 ) 1!

lembrando que: x1 − x0 = h , y (1) ( x0 ) = f ( x0 , y( x0 )) e y ( x1 ) = y1 , podemos escrever:

y1 = y0 + hf ( x0 , y0 ) Equação 10

O erro é então obtido a partir do resto da fórmula de Taylor:

( x − x0 ) 2 ( 2 ) y (ε ); x0 < ε < x1 2! como x1 − x0 = h e usando a notação do erro como ε 1 , vem:

h2 ( 2 ) ε1 = y (ε ) 2! Numa etapa j de cálculo teremos:

h2 ( 2 ) εj = y (ε ); x j − 1 < ε < x j 2! Que é chamado Erro local de truncamento. Prof. José Augusto Lucas Matos

Página 28 de 31

Material sujeito a correções

3.3 Método de Runge-Kutta Como vimos anteriormente, aumentando-se o número de pontos da malha, aumentamos o erro acumulado. Todos os métodos de passo simples são escritos na forma:

y j + 1 = y j + hφ ( x j ; y j ; h),

j = 0,1,2,..........m − 1

onde: φ é a função incremento e h p comprimento do passo. O método de EULER y j + 1 = y j + hf ( x j , y j ) é um método de passo simples (ordem um). Supondo que y ( x ) possua derivadas sucessivas suficientes, pode-se desenvolver y( x + h) em torno de x , através de Taylor:

y ( x + h) = y ( x ) +

h h y' ( x ) + y" (ξ ), x < ξ < ( x + h) 1! 2!

como:

φ ( x j ; y j ; h) = f ( x , y ) vem:

y( x + h) − y( x ) − hφ ( x , y( x ); h) = y( x ) +

h h y' ( x ) + y" (ξ ) − y( x ) − hf ( x , y ) = 1! 2!

h h2 = [ y( x ) + hy' ( x )] − [ y( x ) − hf ( x , y )] + y" (ξ ) = y" (ξ ) 2! 2! Logo Taylor fornece tantos métodos quantos se queira. Assim:

y j + 1 = y j + hy' ( x j ) +

h2 y" ( x j ) com j = 0,1,2,......m − 1 2!

Ex.: Achar as aproximações para a solução do PVI abaixo.

 y' = x − y + 2 de [0;1] com h = 0,1   y ( 0) = 2 onde: x0 = 0 , y0 = 2 , a = 0 , b = 1 , m =

Prof. José Augusto Lucas Matos

1−0 ⇒ m = 10 . 0,1 Página 29 de 31

Material sujeito a correções

Para j = 0

h2 y1 = y0 + h( y( x0 )) + y" ( x0 ) 2! y1 = y0 + h( x0 − y0 + 2) +

h2 ( y0 − x0 − 1) 2!

(0,1) 2 y1 = 2 + 0,1(0 − 2 + 2) + ( 2 − 0 − 1) 2! y1 = 2,005 x1 = x0 + h ⇒ x1 = 0 + 0,1 ⇒ x1 = 0,1 Para j = 1

h2 y2 = y1 + h( y' ( x1 )) + y" ( x1 ) 2! y2 = y1 + h( x1 − y1 + 2) +

h2 ( y1 − x1 − 1) 2!

(0,1) 2 y1 = 2,005 + 0,1(0,1 − 2,005 + 2) + ( 2,005 − 0,1 − 1) 2! y1 = 2,019025 x 2 = x0 + 2h ⇒ x 2 = 0 + 2(0,1) ⇒ x2 = 0,2 O que faremos até j = 9 .

j 0 1 2 3 4 5 6 7 8 9 10

xj 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

yj 2 2,005 2,019025 2,0412176 2,0708020 2,1070758 2,1494036 2,1972102 2,2499753 2,3072276 2,3685410

y( x j ) 2 2,004837 2,018731 2,040818 2,070320 2,106531 2,148812 2,196585 2,249329 2,306570 2,367879

y j − y( x j ) 0,000163 0,000264 0,000399 0,000482 0,000544 0,000591 0,000625 0,000646 0,000657 0,000662

Se compararmos com a tabela do método anterior vemos que o erro e bem menor.

Prof. José Augusto Lucas Matos

Página 30 de 31

Material sujeito a correções

∂f ∂f ( x0 , y0 ) + f ( x0 , y0 ). ( x0 , y0 ) ∂x ∂y y" ( x0 ) = 1 + ( x0 − y0 + 2).( −1) = 1 − x0 + y0 − 2 = ( y0 − x0 − 1) y" ( x0 ) =

∂f ( x 0 , y0 ) + ∂x + f ( x0 , y0 ) + ∂f + ( x0 , y0 ) ∂y

y" ( x0 ) =

Prof. José Augusto Lucas Matos

(y' em relação a x) (y' ) (y' em relação a y)

Página 31 de 31