Métodos Numéricos e Estatísticos
Prof. Marcone Jamilson Freitas Souza
Aula 7: Métodos numéricos para equações diferenciais
• 1a ordem
• Passos múltiplos
• 2a ordem
REDEMAT - Métodos Numéricos e
Estatísticos
1
Equações diferenciais de 1a ordem
Métodos numéricos são usados quando
não é possível obter uma solução geral,
ou a forma dela é tão complicada que
seu uso não é prático.
Uma equação diferencial de 1a ordem
tem a forma F ( x, y, y )  0 ,
e em geral podemos escrevê-la como:
y   f ( x, y )
REDEMAT - Métodos Numéricos e
Estatísticos
Problema do valor inicial
- uma equação diferencial
- uma condição que deve
ser satisfeita pela solução
y  f ( x, y)
y( x0 )  y0
2
Os métodos que estudaremos
partem da idéia de que o espaço da
variável independente (x) pode ser
discretizado, formando uma rede
x0 x1= x0+h x2= x1+h.......
h é o passo .
f ( x, y )  x e
- (x 2  y 2 )
O valor da função em cada
ponto da rede é calculado a
partir de expansões em série
de Taylor.
REDEMAT - Métodos Numéricos e
Estatísticos
3
Método de Euler ou Euler-Cauchy
O valor de y para um passo h é dado
pela expansão:
y ( x0 )  y 0
h2
y ( x0  h)  y ( x0 )  hy ( x0 ) 
y ( x0 )  
2!
Como em geral h é pequeno, suprimimos os
termos de ordem O(h2): h2, h3, .....
Resultando na aproximação
y ( x  h)  y ( x)  hy ( x)  O(h 2 )  y ( x)  hy ( x) 
y ( x)  hf ( x, y )
REDEMAT - Métodos Numéricos e
Estatísticos
4
O que resulta no processo iterativo
y1  y0  hf ( x0 , y0 )
y 2  y1  hf ( x1 , y1 )

y n 1  y n  hf ( xn , y n )
A omissão dos termos de ordem
superior a 2 causa erros de truncagem
(que podem ocorrer junto a erros de
arredondamento).
REDEMAT - Métodos Numéricos e
Estatísticos
5
Exemplo:
y  x  y
n
0
1
2
3
4
5
y (0)  0
passo h=0,2
xn
yn
0,000
0,200
0,400
0,600
0,800
1,000
0,000
0,000
0,040
0,128
0,274
0,488
0,2(x n+y n)
0,000
0,040
0,088
0,146
0,215
Exato
0,000
0,021
0,092
0,222
0,426
0,718
erro
0,000
0,021
0,052
0,094
0,152
0,230
yn 1  yn  0,2 ( xn  yn )
O erro não é (em geral) conhecido. Podemos estimá-lo
utilizando um passo h´=2h
n
0
1
2
3
xn
0,000
0,400
0,800
1,200
yn
0,000
0,000
0,160
0,544
REDEMAT - Métodos Numéricos e
Estatísticos
0,4(x n+y n) y n para h =0,2
erro = y n(h =0,4) - y n(h =0,2)
0,000
0,160
0,384
0,000
0,040
0,274
0,000
0,040
0,114
6
Método de Euler melhorado (2a ordem)
Método chamado de preditor-corretor.
y n*1  y n  hf ( xn , y n )
y n 1
1
 y n  h[ f ( xn , y n )  f ( xn 1 , y n*1 )]
2
xn 1  xn  h
k1  hf ( xn , y n )
k 2  hf ( xn 1 , y n  k1 )
y n !
1
 y n  (k1  k 2 )
2
REDEMAT - Métodos Numéricos e
Estatísticos
7
Exemplo:
o mesmo visto anteriormente
y  x  y
h  0,2
k1  0,2 xn  y n 
k 2  0,2 xn  0,2  y n  0,2 xn  y n 
y n 1  y n 
n
0
1
2
3
4
5
0,2
(k1  k 2 )  y n  0,22 xn  y n   0,02
2
xn
0,000
0,200
0,400
0,600
0,800
1,000
yn
0,0000
0,0200
0,0884
0,2158
0,4153
0,7027
REDEMAT - Métodos Numéricos e
Estatísticos
0,22(x n+y n) + 0,02
0,0200
0,0684
0,1274
0,1995
0,2874
0,3946
Exato
0,0000
0,0214
0,0918
0,2221
0,4255
0,7183
erro
0,0000
0,0014
0,0034
0,0063
0,0102
0,0156
8
Método de Runge-Kutta (4a ordem)
xn 1  xn  h
1
y n 1  y n  (k1  2k 2  2k 3  k 4 )
6
k1  hf ( xn , y n )
1
1
k 2  hf ( xn  h, y n  k1 )
2
2
1
1
k 3  hf ( xn  h, y n  k 2 )
2
2
k 4  hf ( xn  h, y n  k 3 )
REDEMAT - Métodos Numéricos e
Estatísticos
Se f(x,y) não depender
de y, o método reduz-se
à regra de integração
de Simpson
9
n
0
1
2
3
4
5
xn
0,000
0,200
0,400
0,600
0,800
1,000
yn
0,2214(x n+y n) + 0,0214
0,000000
0,021400
0,091818
0,222106
0,425521
0,718251
0,021400
0,070418
0,130288
0,203414
0,292730
0,401821
Exato
0,000000
0,021403
0,091825
0,222119
0,425541
0,718282
erro
0,000000
0,000003
0,000007
0,000012
0,000020
0,000031
Comparação entre os métodos
xn
0,000
0,200
0,400
0,600
0,800
1,000
y = ex - x - 1
0,000000
0,021403
0,091825
0,222119
0,425541
0,718282
REDEMAT - Métodos Numéricos e
Estatísticos
Euler
0,000
0,021
0,052
0,094
0,152
0,229
Erro
Euler melhorado
0,0000
0,0014
0,0034
0,0063
0,0102
0,0156
Runge-Kutta
0,000000
0,000003
0,000007
0,000012
0,000020
0,000031
10
Qual o valor mais adequado para o passo h?
Se a função f varia muito com y, então h deve ser pequeno, para
evitar erros de truncagem. Em geral, adota-se a “proposta” de
que
f
0,01  K  h
 0,05
h  h/2 se K  0,05
y
h  2h se 0,01  K
k3  k 2
h não muda se
K 2
0,05  K  0,01
k 2  k1
1 ~ ~
Estimativa de erro:    y  ~
y
15
onde ~
y é obtido para o passo h
~
onde ~
y é obtido para o passo 2h
REDEMAT - Métodos Numéricos e
Estatísticos
11
Métodos para eq. dif. de segunda ordem
P.V.I.
y   f ( x, y, y )
y ( x0 )  y 0
y ( x0 )  y0
Novamente o problema é obter os
valores de yn e yn´ para a seqüência
x1 = x0 + h; x2 = x0 + 2h; ...
Começamos mais uma vez pelas expansões
em série de Taylor da função e de sua derivada:
h2
h3
y ( x  h)  y ( x)  hy ( x) 
y ( x)  y ( x)  
2!
3!
h2
h3
y ( x  h)  y ( x)  hy ( x) 
y ( x)  y ( x)  
2!
3!
REDEMAT - Métodos Numéricos e
Estatísticos
12
O método mais simples consiste em desprezar
os termos em derivadas de ordem y´´´ ou
superiores
h2
y ( x  h)  y ( x)  hy ( x)  y ( x)
2
y ( x  h)  y ( x)  hy ( x)
1o passo:
2o passo:
y0  f ( x0 , y0 , y0 )
y1  f ( x1 , y1 , y1 )
h2
y1  y0  hy0 
y0
2
y1  y0  hy0
h2
y 2  y1  hy1 
y1
2
y 2  y1  hy1
REDEMAT - Métodos Numéricos e
Estatísticos
13
Runge-Kutta-Nyström
Valores iniciais:
x0 , y0 , y0 ´
passo h
1
k1  hf ( xn , y n , y n )
2
1
1
K  h( y n  k1 )
2
2
1
1
k 2  hf ( xn  h, y n  K , y n  k1 )
2
2
1
1
k 3  hf ( xn  h, y n  K , y n  k 2 )
2
2
L  h( y n  k 3 )
1
k 4  hf ( xn  h, y n  L, y n  2k 3 )
2
xn 1  xn  h
Saída
REDEMAT - Métodos Numéricos e
Estatísticos
1
y n 1  y n  h( y n  ( k1  k 2  k 3 ))
3
1
y n 1  y n  ( k1  2k 2  2k 3  k 4 )
3
14
Equações diferenciais parciais
Uma equação é dita quasilinear se for
linear nas derivadas mais altas:
au xx  2bu xy  cu yy  F ( x, y, u, u x , u y )
 2u
 2u
u
onde u xx  2 ; u xy 
; ux  ; 
xy
x
x
 2  0
(Laplace)
elíptica :
ac  b 2  0
2
 2
2  
c
2
t
 x2
(Onda)
hiperbólic a :
ac  b 2  0
T
 c 2  2T
t
(Calor)
parabólica :
ac  b 2  0
REDEMAT - Métodos Numéricos e
Estatísticos
15
Equações de diferenças para Eq. de
Laplace e Poisson
Vamos ver o caso mais simples em duas dimensões (x e y):
2
2

u

u
2
2
Laplace  u  0
 u  u xx  u yy  2  2
x y
Poisson
 2 u  f ( x, y )
(x,y+k)
k
h
h
(x-h,y)
(x,y)
(x+h,y)
k
(x,y-k)
REDEMAT - Métodos Numéricos e
Estatísticos
16
1 2
1
E1 : u ( x  h, y )  u ( x, y )  hu x ( x, y )  h u xx ( x, y )  u xxx ( x, y ) 
2
6
1 2
1
E 2 : u ( x  h, y )  u ( x, y )  hu x ( x, y )  h u xx ( x, y )  u xxx ( x, y ) 
2
6
E1  E 2 , desprezando termos em O(h 3 )
1
u x ( x, y )  u ( x  h, y )  u ( x  h, y )
2h
1
u y ( x, y )  u ( x, y  k )  u ( x, y  k )
2k
(x,y+k)
k
h
h
(x-h,y)
(x,y)
(x+h,y)
k
(x,y-k)
REDEMAT - Métodos Numéricos e
Estatísticos
17
Para as derivadas segundas, desprezando os termos O(h4),
temos
u ( x  h, y )  u ( x  h, y )  2u ( x, y )  h 2 u xx ( x, y )
1

u ( x  h, y )  2u ( x, y )  u ( x  h, y )
u
(
x
,
y
)

2
 xx
h

u yy ( x, y )  1 u ( x, y  k )  2u ( x, y )  u ( x, y  k )

k2
1
u xy ( x, y ) 
[u ( x  h, y  k )  u ( x  h, y  k )
4hk
u ( x  h, y  k )  u ( x  h, y  k )]
Juntando as aproximações das derivadas primeiras e segundas,
fazendo h=k, obtemos a equação de diferenças correspondente à
equação de Poisson:
u( x  h, y)  u( x, y  h)  u( x  h, y)  u( x, y  h)  4u( x, y)  h 2 f ( x, y)
REDEMAT - Métodos Numéricos e
Estatísticos
18
Para f(x,y) = 0 temos a equação de Laplace.
h é chamado de o comprimento da malha (mesh size).
Equações elípticas - em geral - devem levar em
conta problemas de contorno (condições
previamente definidas numa dada fronteira espacial, por exemplo). Casos mais comuns:
Fronteira C
•Dirichlet: se u é definido na fronteira C
•Neumann: se un=u/n (derivada na
direção normal) é definida na fronteira.
Para resolver o problema, é necessário
criar uma malha.: nós da rede ou da malha (Pij)
REDEMAT - Métodos Numéricos e
Estatísticos
19
Exemplo
Uma placa de 12 cm de lado tem suas bordas mantidas às
temperaturas mostradas na figura. Quais os valores das
temperaturas no interior da placa? Será escolhido um
comprimento h = 4 cm.
u=0
u=0
y
12
u=100
P02
P12
P01
P11
P21
P10
P20
u=100
u=100
R
x
12
u=100
REDEMAT - Métodos Numéricos e
Estatísticos
u=100
20
A equação de transferência de calor é
ut = c2(uxx+uyy)
ui-1,j
Para o regime estacionário ut = 0, a
equação se reduz à de Laplace
uxx+uyy = 0
Para cada ponto da malha, temos a seguinte equação:
ui,j+1
ui,j
ui+1,j
ui,j-1
ui+1,j + ui-1,j + ui,j+1 + ui,j-1 -4 ui,j = 0
P11:
- 4u11 + u21 + u01 + u12 + u10 = 0
- 4u11 + u21 + 100 + u12 + 100 = 0
- 4u11 + u21 + u12 = - 200
REDEMAT - Métodos Numéricos e
Estatísticos
21
- 4u11 + u21
+ u12
u11 - 4u21
u11
u21
= -200
+ u22
= -200
- 4u12 + u22
= -100
+u12
- 4u22 = -100
Dando como resultados
u11 = u21 = 87,5 (88,1)
u12 = u22 = 62,5 (61,9)
REDEMAT - Métodos Numéricos e
Estatísticos
22