EDO

Prévia do material em texto

16/12/2020
1
Métodos Numéricos em Engª Química
Soluções Numéricas de Equações Diferenciais
Profs. Emanuel Negão Macêdo
Josiel Lobato Ferreira
 
 
Faculdade de Engª Química – FEQ
1
2
16/12/2020
2
Problemas de Valor Inicial (PVI) e Problemas de Valor de Contorno (PVC)
3
4
16/12/2020
3
Problemas de Valor Inicial (PVI) e Problemas de Valor de Contorno (PVC)
6
Método de Euler ou Euler-Cauchy
Método de Euler melhorado (2a ordem)
Métodos da Série de taylor
Métodos de Runge-Kutta
Métodos para eq. dif. de segunda ordem
Métodos de Multiplos passos:
Adams Basforth e Adams Mouton
Problemas de Valor Inicial (PVI) 
Método de Shooting
Método de Diferenças Finitas
Problemas de Valor de Contorno (PVC) 
5
6
16/12/2020
4
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 ,
e em geral podemos escrevê-la como:
),(
0),,(
yxfy
yyxF


Problema do valor inicial
- uma equação diferencial
- uma condição que deve 
ser satisfeita pela solução
00 )(
),(
yxy
yxfy


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
)(x- 22
e),( yxyxf 
O valor da função em cada ponto da rede é
calculado a partir de expansões em série de
Taylor.
x0 x1= x0+h x2= x1+h.......
h é o passo .
Método de Euler ou Euler-Cauchy
2        y( x h) y( x) hy ( x) O(h ) y( x ) hy ( x) y( x) hf ( x ,y )
O valor de y para um passo h é dado pela expansão:
Como em geral h é pequeno, suprimimos os termos de ordem O(h2): h2, h3, 
.....
Resultando na aproximação


)(
!2
)()()(
)(
0
2
000
00
xy
h
xyhxyhxy
yxy
O que resulta no processo iterativo
),(
),(
),(
1
1112
0001
nnnn yxhfyy
yxhfyy
yxhfyy





A omissão dos termos de ordem
superior a 2 causa erros de
truncamento (que podem
ocorrer junto a erros de
arredondamento).
00 )(
),(
yxy
yxfy


7
8
16/12/2020
5
n xn yn 0,2(xn+yn) Exato erro
0 0,000 0,000 0,000 0,000 0,000
1 0,200 0,000 0,040 0,021 0,021
2 0,400 0,040 0,088 0,092 0,052
3 0,600 0,128 0,146 0,222 0,094
4 0,800 0,274 0,215 0,426 0,152
5 1,000 0,488 0,718 0,230
Exemplo:
passo h=0,2
O erro não é (em geral) conhecido. Podemos estimá-lo
utilizando um passo h´=2h
0)0( 

y
yxy
)(2,01 nnnn yxyy 
n xn yn 0,4(xn+yn) y n para h=0,2erro = yn(h=0,4) - y n(h=0,2) 
0 0,000 0,000 0,000 0,000 0,000
1 0,400 0,000 0,160 0,040 0,040
2 0,800 0,160 0,384 0,274 0,114
3 1,200 0,544
9
10
16/12/2020
6
Método de Euler melhorado (2a ordem)
Método chamado de preditor-corretor.
)],(),([
2
1
),(
*
111
*
1




nnnnnn
nnnn
yxfyxfhyy
yxhfyy
)(
2
1
),(
),(
21!
112
1
1
kkyy
kyxhfk
yxhfk
hxx
nn
nn
nn
nn







Exemplo:
o mesmo visto anteriormente
 
  
 
1
2
1 1 2
0 2
0 2
0 2 0 2 0 2
1
0 22 0 02
2
n n
n n n n
n n n n n
y x y
h ,
k , x y
k , x , y , x y
y y ( k k ) y , x y ,
  

 
    
      
n xn yn 0,22(x n+yn) + 0,02 Exato erro
0 0,000 0,0000 0,0200 0,0000 0,0000
1 0,200 0,0200 0,0684 0,0214 0,0014
2 0,400 0,0884 0,1274 0,0918 0,0034
3 0,600 0,2158 0,1995 0,2221 0,0063
4 0,800 0,4153 0,2874 0,4255 0,0102
5 1,000 0,7027 0,3946 0,7183 0,0156
11
12
16/12/2020
7
00 )(
),(
yxy
yxfy


13
14
16/12/2020
8
00 )(
),(
yxy
yxfy


15
16
16/12/2020
9
17
18
16/12/2020
10
Para RK4 - Se f(x,y) não depender de y, o método reduz-se à regra de integração de
Simpson
n xn yn 0,2214(xn+yn) + 0,0214 Exato erro
0 0,000 0,000000 0,021400 0,000000 0,000000
1 0,200 0,021400 0,070418 0,021403 0,000003
2 0,400 0,091818 0,130288 0,091825 0,000007
3 0,600 0,222106 0,203414 0,222119 0,000012
4 0,800 0,425521 0,292730 0,425541 0,000020
5 1,000 0,718251 0,401821 0,718282 0,000031
Comparação entre os métodos
xn y = ex - x - 1 Erro
Euler Euler melhorado Runge-Kutta
0,000 0,000000 0,000 0,0000 0,000000
0,200 0,021403 0,021 0,0014 0,000003
0,400 0,091825 0,052 0,0034 0,000007
0,600 0,222119 0,094 0,0063 0,000012
0,800 0,425541 0,152 0,0102 0,000020
1,000 0,718282 0,229 0,0156 0,000031
0)0( 

y
yxy
19
20
16/12/2020
11
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
12
232
05,001,0
kk
kk
K
y
f
hK





 h  h/2 se K  0,05
h  2h se 0,01  K
h não muda se
0,05  K  0,01
Estimativa de erro:  
hy
hy
yy
2passooparaobtidoé
~~onde
passooparaobtidoé~onde
~~~
15
1

21
22
16/12/2020
12
Métodos para eq. dif. de segunda ordem
P.V.I.
00
00
)(
)(
),,(
yxy
yxy
yyxfy


 Novamente o problema é obter os valores de yn e 
yn´ para a seqüência x1 = x0 + h; x2 = x0 + 2h; ...




)(
!3
)(
!2
)()()(
)(
!3
)(
!2
)()()(
32
32
xy
h
xy
h
xyhxyhxy
xy
h
xy
h
xyhxyhxy
Começamos mais uma vez pelas expansões em
série de Taylor da função e de sua derivada:
O método mais simples consiste em desprezar
os termos em derivadas de ordem y´´´ ou
superiores
)()()(
)(
2
)()()(
2
xyhxyhxy
xy
h
xyhxyhxy


001
0
2
001
0000
2
),,(
yhyy
y
h
yhyy
yyxfy



1o passo:
2o passo:
112
1
2
112
1111
2
),,(
yhyy
y
h
yhyy
yyxfy



Runge-Kutta-Nyström
)22(
3
1
))(
3
1
(
)2,,(
2
1
)(
),,
2
1
(
2
1
),,
2
1
(
2
1
)
2
1
(
2
1
),,(
2
1
43211
3211
1
34
3
23
12
1
1
kkkkyy
kkkyhyy
hxx
kyLyhxhfk
kyhL
kyKyhxhfk
kyKyhxhfk
kyhK
yyxhfk
nn
nnn
nn
nnn
n
nnn
nnn
n
nnn












Valores iniciais: x0, y0, y0´
passo h
Saída
00
00
)(
)(
),,(
yxy
yxy
yyxfy



23
24
16/12/2020
13
25
26
16/12/2020
14
Problema de Valor de Contorno
27
28
16/12/2020
15
Problema de Valor de Contorno
Iremos aproximar a solução do PVC usando o Método Shooting (MS) e o
Método de Diferenças Finitas (MDF)
PVC – Método de Shooting
29
30
16/12/2020
16
31
32
16/12/2020
17
PVC – Método de Diferenças Finitas
33
34
16/12/2020
18
Primeiramente, dividimos o intervalo [a; b] em N subintervalos de
tamanho h. Onde o passo da malha é dado por h = (b-a)/N e os pontos
da malha são dados por xk = a + kh (k = 0; 1; : : : ; N).
PVC – Discretização
x=a x=b
x
x0 x1 x2 x3 x4.................xk-1 xk xk+1 ....................xN-1 xN
Para k = 1; 2; ...;N-1
35
36
16/12/2020
19
v0 = 1 vN = 2
Para k = 1; 2; ...;N-1
37
38
16/12/2020
20
Exemplo
39
40
16/12/2020
21
Diferenças Finitas para EDOs
41
42
16/12/2020
22
De forma análoga ao que foi feito
no exemplo anterior obteve-se as
seguintes equações:
43
44
16/12/2020
23
45
46
16/12/2020
24
Equações diferenciais parciais
Método das Linhas
0:parabólicaCalor)(
0:ahiperbólicOnda)(
0:elípticaLaplace)(0
222
2
2
2
2
2
2
22










bacTc
t
T
bac
x
c
t
bac
Uma equação é dita quasilinear se for
linear nas derivadas mais altas:
;;;onde
),,,,(2
2
2
2
x
u
u
yx
u
u
x
u
u
uuuyxFcubuau
xxyxx
yxyyxyxx










Método das Linhas
MOL
O método das linhas é um método semi-discreto para a resolução de
EDP's que consiste em discretizar as variáveis espaciais e manter uma
das varáveis contínua (usualmente o tempo), de modo a transformar
a EDP em um sistema de EDO's que pode então ser resolvido através
dos métodos para a resolução de PVI's (como os métodos de
RungeKutta).
A abordagem utilizada para a discretização das variáveis espaciais
usualmente é o método de diferenças finitas, por isso o método das
linhas é muitas vezes chamado de método de diferenças finitas semi-
discreto.
Este método é aplicado principalmente para equações parabólicas,
pois sua aplicação em equaçõeselípticas origina um conjunto de
PVC's, o que por sua vez também precisam ser resolvido por métodos
de discretização. Quando aplicado em equações parabólicas, a
variável onde a informação se propaga somente em uma das direções
(como por exemplo o tempo) é mantida contínua enquanto as demais
são discretizadas.
47
48
16/12/2020
25
Para ilustrar a aplicação do método das linhas, considere que é necessário obter a
variação de temperatura ao longo de uma barra metálica com uma extremidade
isolada e outra mantida a uma temperatura Text e que perde calor para o meio
externo por convecção. Considere que se deseja obter como a temperatura varia ao
longo do tempo a partir de um estado inicial T(x; 0) = Tini.
Para resolver esta equação com o método das linhas, é preciso discretizar a equação
na direção espacial e manter a função contínua no tempo. Assim, deve-se definir um
domínio discreto na direção x. Neste caso, será considerado um domínio somente
com N = 6 pontos (ou Ni = 5 intervalos) para ilustração :
Neste caso, a discretização da EDP na direção x irá
resultar um conjunto de valores (T1; T2; T3; T4; T5;
T6) que representam a temperatura nos pontos
respectivos (x1; x2; x3; x4; x5; x6).
49
50
16/12/2020
26
A EDP avaliada neste exemplo envolve a derivada
segunda em relação a direção x, então deve-se
discretizar esta derivada. Considerando um esquema
central, a derivada segunda nos pontos xi são dadas por:
Para fechar o sistema de equações, é preciso ainda denir equações para T0 e T6,
que são obtidas através da aplicação das condições de contorno. A condição T(0; t)
= Text resulta diretamente em: T1 = Text
Assim, obtém-se uma equação para a temperatura na extremidade x = 0. Na
extremidade x = L, a condição é de derivada nula. Neste caso pode-se aplicar um
esquema de discretização para trás, de onde se obtém que T6 = T5.
Considere, por exemplo, um caso onde L = 1 cm (Δx = 0.2 cm), Text = 100 C,
Ta=25 C e h′=0.1 cm^-2. Além disso, considere que
alfa = 0.01 cm^2/s e Tini = 25 C. Neste caso, a barra metálica está inicialmente na
mesma temperatura que o ambiente externo. Em um dado instante, a temperatura na
extremidade x = 0 é aumentada para 100 C.
51
52
16/12/2020
27
MOL - Exemplo
53
54
16/12/2020
28
55
56
16/12/2020
29
57
58
16/12/2020
30
Método de Diferenças Finitas
1 – Método Simples Explícito
2 – Método Simples Implícito
3 – Método de Crank - Nicolson
59
60
16/12/2020
31
2
1
0 Estável
2
Logo se definirmos a malha x(Nint) o limite para
x
t , para um valor fixo de .
2
   


  
2
f
1
int
int
3 f
step
2
int
int
5 f
step
Exemplo : 0.25; t 1; t x
1
a)N 10 x 10
N
t
t 2.5 10 N 400
t
1
b)N 100 x 10
N
t
t 2.5 10 N 40000
t




     
   
     

    
     

61
62
16/12/2020
32
63
64
16/12/2020
33
65
66
16/12/2020
34
67
68
16/12/2020
35
69
70
16/12/2020
36
f
71
72
16/12/2020
37
Lista de Exercícios
Resolva as questões usando a metodologia adequada (no Fortran) e compare com a solução do
Mathematica. Discuta os resuldados.
Parte 1 – Problemas de Valor Inicial (PVI)
1)
2) 3)
Parte 1 – Problemas de Valor Inicial (PVI)
4) 5)
6)
73
74
16/12/2020
38
7)
8)
RK-4
RK-4 – Runge Kutta de 4ª ordem
9)
10)
75
76
16/12/2020
39
11)
12) Resolver o seguinte PVI usando a rotina IVPAG do IMSL – Fortran e comparar com a solução do
software de manipulação simbólica Mathematica.
13)
Parte 2 – Problemas de Valor de Contorno (PVC)
1)
2)
77
78
16/12/2020
40
Parte 2 – Problemas de Valor de Contorno (PVC)
3)
4)
5 = a); 6 = b); 7 = c); 8 = d)
With a) h = /4 and b) h = /4
Parte 2 – Problemas de Valor de Contorno (PVC)
9)
10)
79
80
16/12/2020
41
Parte 3 – EDPs
Resolva os problemas abaixo usando o Mathematica e por esquema numérico no Fortran.
1) Resolva o problema abaixo para as malhas com 10, 25 e 50 intervalos. Use: a) Método das Linhas; b)
Método simples explícito, c) Método simples implícito e c) Método de Crank-Nicolson. Compare com a
solução do Matemática.
1) Use K = 1, W=0.1 2) Use K = 1, Co = 0.1; Bim = 1 3) Use go = 10, qo = 0.1; Bi = 10 e
e Csat = 1 e  = 0.1
4) Use A0 = 1; A1 = 0.1 e  = 1 e 10 5) Use U0 = 1;  = 1 e 10
2
2
C C
K C
t X
t 0 C 0
C
X 0 W
X
X 1 C 0
 
    
  

    
 
   
2
2
im im sat
C C
K C
t X
t 0 C Co
C
X 0 B C B C
X
C
X 1 0
X
     
  
 
       
    
2
b(X 1)
02
0
g e
t X
t 0 1
X 0 q
X
X 1 Bi Bi
X
 

  
    
  
 
    
       
0 1
2
2
A A A Sin( t)
U U
A
t Y
t 0 U 0
U
Y 0 0
Y
Y 1 U 0
   
 
   
  

   
 
   
2
2
0
U U
t Y
t 0 U 0
Y 0 U U Sin( t)
Y 1 U 0
 
     
     

  
PPEQ - ITEC / UFPA 82
Equações diferenciais parciais
0:parabólicaCalor)(
0:ahiperbólicOnda)(
0:elípticaLaplace)(0
222
2
2
2
2
2
2
22










bacTc
t
T
bac
x
c
t
bac
Uma equação é dita quasilinear se for
linear nas derivadas mais altas:
;;;onde
),,,,(2
2
2
2
x
u
u
yx
u
u
x
u
u
uuuyxFcubuau
xxyxx
yxyyxyxx










81
82
16/12/2020
42
PPEQ - ITEC / UFPA 83
Equações de diferenças para Eq. de 
Laplace e Poisson
),(
0
2
2
2
2
2
22
yxfu
y
u
x
u
uuuu yyxx






Laplace
Poisson
Vamos ver o caso mais simples em duas dimensões (x e y):
(x-h,y) (x,y) (x+h,y)
h h
k
k
(x,y-k)
(x,y+k)
PPEQ - ITEC / UFPA 84
 
 ),(),(
2
1
),(
),(),(
2
1
),(
)(emtermosodesprezand,EE
),(
6
1
),(
2
1
),(),(),(:E
),(
6
1
),(
2
1
),(),(),(:E
3
21
2
2
2
1
kyxukyxu
k
yxu
yhxuyhxu
h
yxu
hO
yxuyxuhyxhuyxuyhxu
yxuyxuhyxhuyxuyhxu
y
x
xxxxxx
xxxxxx







(x-h,y) (x,y) (x+h,y)
h h
k
k
(x,y-k)
(x,y+k)
83
84
16/12/2020
43
PPEQ - ITEC / UFPA 85
 
 
)],(),(
),(),([
4
1
),(
),(),(2),(
1
),(
),(),(2),(
1
),(
),(),(2),(),(
2
2
2
kyhxukyhxu
kyhxukyhxu
hk
yxu
kyxuyxukyxu
k
yxu
yhxuyxuyhxu
h
yxu
yxuhyxuyhxuyhxu
xy
yy
xx
xx











Para as derivadas segundas, desprezando os termos O(h4),
temos
Juntando as aproximações das derivadas primeiras e segundas,
fazendo h=k, obtemos a equação de diferenças correspondente à
equação de Poisson:
),(),(4),(),(),(),( 2 yxfhyxuhyxuyhxuhyxuyhxu 
PPEQ - ITEC / UFPA 86
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:
•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)
Fronteira C
85
86
16/12/2020
44
PPEQ - ITEC / UFPA 87
Exemplo 1
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.
12
x
y
12
u=0
u=100
u=100
u=100 R
u=0
u=100
u=100
P02
P10 P20
P01 P11 P21
P12
PPEQ - ITEC / UFPA 88
A equação de transferência de calor é
ut = c2(uxx+uyy)
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+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
ui+1,jui-1,j
ui,j+1
ui,j-1
ui,j
87
88
16/12/2020
45
PPEQ - ITEC / UFPA 89
- 4u11 + u21 + u12 = -200
u11 - 4u21 + u22 = -200
u11 - 4u12 + u22 = -100
u21 +u12 - 4u22 = -100
Dando como resultados
u11 = u21 = 87,5 (88,1)
u12 = u22 = 62,5 (61,9)
PPEQ - ITEC / UFPA 90
Exemplo 2
Resolver a equação de Poisson apresentadaabaixo com
Δx=Δy=0,25cm e usando a eliminação gaussiana. O problema
apresenta uma chapa de w=1cm de largura por h=1,5 cm de
comprimento, k = 0,4 J/cm.s.C e Q = 400 J/cm3.s e com
temperatura igual a 0°C em todas as bordas.
2 2
2 2
T T Q
x y k
 
  
 

89
90
16/12/2020
46
PPEQ - ITEC / UFPA 91
Esquema do problema proposto:
PPEQ - ITEC / UFPA 92
Aplicando diferenças finitas na equação acima, considerando “i” e
“j” os índices nas direções x e y, respectivamente, temos:
 2 2 2 2
1 1 1 1
2 1 0i , j
i , j i , j i , j i , j i , j
Q
T T T T T x
k
   
   
 
       
 

Onde:
x
y



Como Δx=Δy tem-se que β=1 e substituindo os valores de Q e 
k resulta:
1 1 1 1
4 62 5
i , j i , j i , j i , j i , j
T T T T T ,
   
      
Aplicando-se a equação acima em cada ponto interno da malha
montou-se o seguinte sistema linear:
91
92
16/12/2020
47
PPEQ - ITEC / UFPA 93
4 1 0 0 0 1 0 0 0 0 0 0 0 0 0
1 4 1 0 0 0 1 0 0 0 0 0 0 0 0
0 1 4 1 0 0 0 1 0 0 0 0 0 0 0
0 0 1 4 1 0 0 0 1 0 0 0 0 0 0
0 0 0 1 4 0 0 0 0 1 0 0 0 0 0
1 0 0 0 0 4 1 0 0 0 1 0 0 0 0
0 1 0 0 0 1 4 1 0 0 0 1 0 0 0
0 0 1 0 0 0 1 4 1 0 0 0 1 0 0
0 0 0 1 0 0 0 1 4 1 0 0 0 1 0
0 0 0 0 1 0 0 0 1 4 0 0 0 0 1
0 0 0 0 0 1 0 0 0 0 4 1 0 0 0
0 0 0 0 0 0 1 0 0 0 1 4 1 0 0
0 0 0 0 0 0 0 1












1 1
2 1
3 1
4 1
5 1
3 2
2 2
3 2
4 2
5 2
1 3
2 3
3 3
4 3
5 3
62 5
0 0 0 1 4 1 0
0 0 0 0 0 0 0 0 1 0 0 0 1 4 1
0 0 0 0 0 0 0 0 0 1 0 0 0 1 4
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
T ,
T
T
T
T
T
T
T
T
T
T
T
T
T
T
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
     
  
  
     
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
62 5
,
,
,
,
,
,
,
,
,
,
,
,
,
,
 
  
 
  
 
 
 
 
 
 
 
 
 
  
 
  
 
 
 
PPEQ - ITEC / UFPA 94
0°C
0°C
0°C
0°C
48,46 62,41 48,46
68,94 90,21 68,94
74,60 90,03 74,60
68,94 90,21 68,94
48,46 62,41 48,46
y
x
w
h
Resolvendo o sistema acima por eliminação gaussiana temos as
seguintes temperaturas:
93
94
16/12/2020
48
PPEQ - ITEC / UFPA 95
A figura abaixo representa o perfil de temperatura na placa na qual
a temperatura é igual a zero nas bordas e cresce a medida que se
caminha para o interior da placa, devido à fonte interna de calor.
PPEQ - ITEC / UFPA 96
Abaixo tem-se as isolinhas de temperatura e nota-se o perfil de
temperatura não é suave devido malha ser bem grosseira o que
também pode observado na figura anterior.
95
96