NotasAulaENG1714 3 MarcioCarvalho

Prévia do material em texto

EQUAÇÃO DIFERENCIAL ORDINÁRIA
PROBLEMA DE VALOR INICIAL
cayyxf
dx
dy  )( :inicial cond.,),(
 Geralmente a variável x representa o tempo e a equação diferencialrepresenta a lei da natureza que descreve a taxa de variação de uma grandeza
Exemplo: A taxa de crescimento de uma população é proporcional à população:
habitantes de numero :, NKN
dt
dN 
 É comum uma situação ser descrito por um sistema de equações diferenciais
siyyyxf
dx
dy
si
i ,,2,1,),,,,( 21  
Em notação vetorial: cyyfy  )(,),( ax
dx
d
 Equações diferenciais de ordem superior podem ser escritas comoum sistema de equações diferenciais de primeira ordem 















333213
2232
1121
2
2
321
3212
2
3
3
)0(,),,,(
)0(,
)0(,
,,
)0(,)0(,)0( :iniciais cond.,,,,





xg
dx
d
dx
d
dx
d
dx
yd
dx
dyy
yyy
dx
yd
dx
dyyxg
dx
yd
MÉTODO DE EULER
Método de Passo a Passo Explícito: O valor da função no instante k+1é calculado somente em função do valor da função no instante k
 Também conhecido como Método da Tangente
 Solução numérica: A função y(x) não será obtida para todos os valores de x
 Definir pontos xi onde a função y(x) será calculada: Malha
x
y
0x 1x 2x 3x 4x
1y
0y
2y
Série de Taylor de y(x) em xk:
),(
)(2)()()(
1
2
1
kkkk
kkkk
yxfhyy
xyhxyhxyxy



 
Exemplo:
1)0(; :Resolver  yy
dx
dy Solução exata: xexy )(
Método de Euler: kkkkk yhyyyxfhyy   11 ),(
xk y(xk) yk Erro0.0 1.000 1.000 0.0000.2 1.221 1.200 -0.0210.4 1.492 1.440 -0.0520.6 1.822 1.728 -0.094
h = 0.2
xk y(xk) yk Erro0.0 1.000 1.000 0.0000.1 1.105 1.100 -0.0050.2 1.221 1.210 -0.0110.3 1.350 1.331 -0.0190.4 1.492 1.464 -0.0280.5 1.649 1.610 -0.0390.6 1.822 1.771 -0.051
h = 0.1
A solução apresenta erros.
Aparentemente o erro é proporcional ao passo h
Dois tipos de Erros: Erro de Truncamento (Discretização)Erro de Arredondamento
Erro de Truncamento (Discretização) )(),( kkk xyyxhE 
Erro que surge devido à aproximação da curva y(x)por uma reta no intervalo h
Propagação do erro: Solução numérica tende a aproximar a solução exata que passa por (xk ,yk) 
Erro pode ser decrescido diminuindo-se o passo h x
y
Família de curvas
Aumento do número de intervalos leva a umaumento do número de cálculos necessários e consequentementeum aumento do tempo computacional e dos erros de arredondamento
hn
1
Erro
Truncamento
Total
Arredondamento
Existe um tamanho de passo ótimo para produzir o erro mínimo
Análise do Erro de Truncamento para o Método de Euler
)(),( kkk xyyxhE ERRO GLOBAL
ERRO LOCAL: Diferença entre o valor yk e o valor da solução exataque passa pelo ponto (xk-1, yk-1)
x
y
Erro local
),( 11  kk yx
Erro global
Deseja-se obter uma estimativa do erro e saber se o erro cresce ou não
x
y
),( 00 yx
))(,( 11 xyx
),( 11 yx
))(,( 22 xyx
0 Se yf
 ),(
1 ),(),(
kk yxy
kkkk yxfhyyyxfdx
dy

 
FAMÍLIA DE CURVAS DIGERVE
))(,(),( 1111 xyxyyxy 
ERRO AUMENTA EXPONENCIALMENTE
x
y
),( 00 yx
))(,( 11 xyx
),( 11 yx
))(,( 22 xyx
0 Se yf
FAMÍLIA DE CURVAS DIGERVE
))(,(),( 1111 xyxyyxy 
ERRO LOCAL COMPENSA 
ERRO DA DISCRETIZAÇÃO
SE O ERRO LOCAL EM TODOS OS PONTOS FOR MENOR DO QUE 
O ERRO GLOBAL SATISFAZ A SEGUINTE CONDIÇÃO:
Lh
exyy
Lh
y
fL
e
exyy
Lnh
nn
xxxLh
Lnh
nn
n
1)(
:pequenofor Se
max onde ;11)( 0







Erro pode crescer exponencialmente se L > 0
Erro pode ficar em uma faixa aceitável se h for muito pequeno
Análise do Estabilidade
ESTÁVEL: Produz solução limitada 
INSTÁVEL: Produz solução que tente ao infinito 
)( 
O Método de Euler é condicionalmente estável
y
dx
dy  :Exemplo
Solução exata:
Solução pelo método de Euler:
xeyxy 0)( 
n
nnnn hyyyhyy )1(011   
Observe que, por serie de Taylor:  2)(1
2hhe h 
A solução pelo método de Euler produz os dois primeiros termos da serie de Taylor 
0 Para  Solução exata decai com x
Solução aproximada pode crescer com x se 1h1 
:)1(2 Para 0  y
Para <0, o Método de Euler só é estável para 2h
Metodo de Euler:
 -2.0000 h = h = h =0.2000 0.4000 1.0000Xk Exata Yk Yk Yk0.00 1.0000 1.0000 1.0000 1.00000.20 0.6703 0.60000.40 0.4493 0.3600 0.20000.60 0.3012 0.21600.80 0.2019 0.1296 0.04001.00 0.1353 0.0778 -1.00001.20 0.0907 0.0467 0.00801.40 0.0608 0.02801.60 0.0408 0.0168 0.00161.80 0.0273 0.01012.00 0.0183 0.0060 0.0003 1.00002.20 0.0123 0.00362.40 0.0082 0.0022 0.00012.60 0.0055 0.00132.80 0.0037 0.0008 0.00003.00 0.0025 0.0005 -1.00003.20 0.0017 0.0003 0.00003.40 0.0011 0.00023.60 0.0007 0.0001 0.00003.80 0.0005 0.00014.00 0.0003 0.0000 0.0000 1.00004.20 0.0002 0.00004.40 0.0002 0.0000 0.00004.60 0.0001 0.00004.80 0.0001 0.0000 0.00005.00 0.0000 0.0000 -1.0000
-1.5000
-1.0000
-0.5000
0.0000
0.5000
1.0000
1.5000
0.00 1.00 2.00 3.00 4.00 5.00
Exata
h=0.2
h=0.4
h=1.0
Microsoft Excel 
Worksheet
1)0(, problema oResolver  yyy
Comparar a solução exata com a obtida pelo método de Euler com
h = 0.1, h = 0.5, h = 1.5 e h = 2. Metodo de Euler:
 -1.0000 h=0.1 h=0.5 h=1.50.00 1.0000 1.0000 1.0000 1.00000.10 0.9048 0.90000.20 0.8187 0.81000.30 0.7408 0.72900.40 0.6703 0.65610.50 0.6065 0.5905 0.50000.60 0.5488 0.53140.70 0.4966 0.47830.80 0.4493 0.43050.90 0.4066 0.38741.00 0.3679 0.3487 0.25001.10 0.3329 0.31381.20 0.3012 0.28241.30 0.2725 0.25421.40 0.2466 0.22881.50 0.2231 0.2059 0.1250 -0.50001.60 0.2019 0.18531.70 0.1827 0.16681.80 0.1653 0.15011.90 0.1496 0.13512.00 0.1353 0.1216 0.06252.10 0.1225 0.10942.20 0.1108 0.09852.30 0.1003 0.08862.40 0.0907 0.07982.50 0.0821 0.0718 0.03132.60 0.0743 0.06462.70 0.0672 0.05812.80 0.0608 0.05232.90 0.0550 0.04713.00 0.0498 0.0424 0.0156 0.25003.10 0.0450 0.03823.20 0.0408 0.03433.30 0.0369 0.03093.40 0.0334 0.02783.50 0.0302 0.0250 0.00783.60 0.0273 0.02253.70 0.0247 0.02033.80 0.0224 0.01823.90 0.0202 0.01644.00 0.0183 0.0148 0.0039
-0.6000
-0.4000
-0.2000
0.0000
0.2000
0.4000
0.6000
0.8000
1.0000
1.2000
0.00 1.00 2.00 3.00 4.00 5.00
exata
h=0.1
h=0.5
h=1.5
Microsoft Excel 
Worksheet
MÉTODO DE EULER DE ORDEM ELEVADA












),(2),()()(
),()( Onde
)(2)()()(
2
1
2
1
kkkkkk
k
kkkk
yxf
y
f
x
fhyxfhxyxy
dx
dy
y
f
x
fyxf
dx
dy
dx
dxy
xyhxyhxyxy 
1)( :Exata Sol. 2)0(, :Resolver  xexyyxy
dx
dy x
Usar método de Euler de segunda ordem com h = 0.1
MÉTODO DE RUNGE-KUTTA
Método de Passo a Passo Explícito: O valor da função no instante k+1é calculado somente em função do valor da função no instante k
 O coeficiente que multiplica o passo h no cálculo da função no instante k+1é calculado de forma que a expansão coincida com o desenvolvimentoem série de Taylor até os termos de ordem N (ordem do método)
 Derivação para Runge-Kutta de segunda ordem:
),(Aprox.
2)()(
:ordem) (2Taylor de Serie
),( c.i.,),(
21211
2
1
a
00
kkkkkk
kk
kkk
FhyhxFhFhyy
dx
dy
y
F
x
FhFhxyxy
yxyxF
dx
dy
 








  ()
As constantes são determinadas de tal forma que asexpressões anteriores coincidam2121 ,,, 
Por série de Taylor:
  










 y
FF
x
FhFhyy
y
FFh
x
FhyxFFhyhxF
k
k
k
kkk
k
k
k
kkkkk
2122211
2121 ),(),(


A solução aproximada fica:
()Comparando com
1 e 2121;21;1 2121221221  





  ),( ),(;22121211 hfyhxFf yxFfffhyy kk kkkk
RUNGE KUTTA DE QUARTA ORDEM
 34
23
12
1
43211
, 2
,2
2,2
),( 6336
hfyhxFf
fhyhxFf
fhyhxFf
yxFf
ffffhyy
kk
kk
kk
kk
kk



 


 



 
Método de Runge-Kutta mais utilizado
Boa combinação entre precisão e simplicidade de programação
)5.0 ordem, quarta Kutta-(Runge 1)0(, problema oResolver  hyyy
 
 
 

304
203
102
01
432101
5.02
5.02
5.0
2265.0
),(
fyf
fyf
fyf
yf
ffffyy
yyxF



 


 



Metodo de Runge-Kutta (4 ordem):
h = 0.5000 f1 f2 f3 f4 y - approx0.00 1.0000 1.00000.50 0.6065 -1.0000 -0.7500 -0.8125 -0.5938 0.60681.00 0.3679 -0.6068 -0.4551 -0.4930 -0.3603 0.36821.50 0.2231 -0.3682 -0.2761 -0.2991 -0.2186 0.22342.00 0.1353 -0.2234 -0.1675 -0.1815 -0.1326 0.13552.50 0.0821 -0.1355 -0.1017 -0.1101 -0.0805 0.08223.00 0.0498 -0.0822 -0.0617 -0.0668 -0.0488 0.04993.50 0.0302 -0.0499 -0.0374 -0.0405 -0.0296 0.03034.00 0.0183 -0.0303 -0.0227 -0.0246 -0.0180 0.01844.50 0.0111 -0.0184 -0.0138 -0.0149 -0.0109 0.01115.00 0.0067 -0.0111 -0.0084 -0.0091 -0.0066 0.00685.50 0.0041 -0.0068 -0.0051 -0.0055 -0.0040 0.00416.00 0.0025 -0.0041 -0.0031 -0.0033 -0.0024 0.00256.50 0.0015 -0.0025 -0.0019 -0.0020 -0.0015 0.00157.00 0.0009 -0.0015 -0.0011 -0.0012 -0.0009 0.00097.50 0.0006 -0.0009 -0.0007 -0.0007 -0.0005 0.00068.00 0.0003 -0.0006 -0.0004 -0.0005 -0.0003 0.00038.50 0.0002 -0.0003 -0.0003 -0.0003 -0.0002 0.00029.00 0.0001 -0.0002 -0.0002 -0.0002 -0.0001 0.00019.50 0.0001 -0.0001 -0.0001 -0.0001 -0.0001 0.000110.00 0.0000 -0.0001 -0.0001 -0.0001 0.0000 0.0000
5.1 ordem, quarta Kutta-Runge h
Metodo de Runge-Kutta (4 ordem):
h = 1.5000 f1 f2 f3 f4 RK EULER0.00 1.0000 1.0000 1.00001.50 0.2231 -1.0000 -0.2500 -0.8125 0.2188 0.2734 -0.50003.00 0.0498 -0.2734 -0.0684 -0.2222 0.0598 0.0748 0.25004.50 0.0111 -0.0748 -0.0187 -0.0607 0.0164 0.0204 -0.12506.00 0.0025 -0.0204 -0.0051 -0.0166 0.0045 0.0056 0.06257.50 0.0006 -0.0056 -0.0014 -0.0045 0.0012 0.0015 -0.03139.00 0.0001 -0.0015 -0.0004 -0.0012 0.0003 0.0004 0.015610.50 0.0000 -0.0004 -0.0001 -0.0003 0.0001 0.0001 -0.0078
-0.6000
-0.4000
-0.2000
0.0000
0.2000
0.4000
0.6000
0.8000
1.0000
1.2000
0.00 2.00 4.00 6.00 8.00 10.00 12.00
Exata
RK
Euler
Exercício
Escreva uma rotina MatLab para solução de um problema de valor inicial usando os métodos de Euler explícito de primeira ordem e Runge-Kutta de quarta ordem. Utilize a rotina desenvolvida para resolver o problema:
Determine o valor de para diferentes passos de tempo. 
RUNGE KUTTA PARA SISTEMA DE EDOS
 
 
 334
223
112
1
43211
0
0
,, 2
,2,2
2,2,2
,,
226
)(;),,(
)(;),,(
hgzhfyhxFf
ghzfhyhxFf
ghzfhyhxFf
zyxFf
ffffhyy
BxzzyxG
dx
dz
AxyzyxF
dx
dy
kkk
kkk
kkk
kkk
kk



 


 




  
 
 334
223
112
1
43211
,, 2
,2,2
2,2,2
,,
226
hgzhfyhxGg
ghzfhyhxGg
ghzfhyhxGg
zyxGg
gggghzz
kkk
kkk
kkk
kkk
kk



 


 


Exercício
As equações de Lotka-Volterra descrevem a evolução da população de um sistema predador-presa, onde x e y são os números de presas e predadores, a é a taxa de crescimento da presa, c é a taxa de morte do predador e b e d são taxas caracterizando o efeito da interação predador-presa na morte da presa e no crescimento do predador, respectivamente.
a) Determine a evolução da população no intervalo 0 < t < 30 utilizando o método de RK-4ª ordem. Compare os resultados em termos do valor de passo de tempo necessário para obter a solução do problema. Utilize os seguintes dados:
dxycy
dt
dy
bxyax
dt
dx


1)0( 2)0(
4,08,0
6,0 2,1






ty
tx
d
c
b
a b) Construa o gráfico da população da presa em função da população do predador. Comente o resultado obtido.
c) O que acontece quando ?2)0( 2)0(  ty tx
MÉTODOS IMPLÍCITOS
 Nos métodos vistos até agora, o valor da função no instante k+1é calculado somente em função de informações no instante k.
 Como visto anteriormente, os métodos explícitos somente são estáveiscomo o tamanho do passo muito pequeno.
x
y
0x 1x 2x
1y
2y
),( ),(
)(2)()()(
111
111
1
2
111







kkkk
kkkk
kkk
x
k
yxfhyy
yxfhyy
xyhxyhxyhxy
hx
k

Série de Taylor para trás ao redor de xk+1:
Lado direito da equação não é conhecido !!
MÉTODO DE EULER IMPLÍCITO
1)0(, problema oResolver  yyy
)1(
)(),(),(
1
11111
h
yy
yhyyyxfhyy
yyxF
k
k
kkkkkkk





Metodo de Euler Implicito
 -1.0000 h=1.5 (EXP) h=0.5 (IMP)0.00 1.0000 1.0000 1.00000.50 0.60651.00 0.36791.50 0.2231 -0.5000 0.40002.00 0.13532.50 0.08213.00 0.0498 0.2500 0.16003.50 0.03024.00 0.01834.50 0.0111 -0.1250 0.0640 -0.6000
-0.4000
-0.2000
0.0000
0.2000
0.4000
0.6000
0.8000
1.0000
1.2000
0.00 1.00 2.00 3.00 4.00 5.00
Exata
Euler Explicito
Euler Implicito
Em casos simples, a equação pode ser facilmente rearrumada
para determinar a função no passo k+1
Se F(x,y) for uma função não linear, a função no passo k+1
será calculada através da solução de uma equação 
não linear (pelo Método de Newton, por exemplo).
O método de Euler Implícito é incondicionalmente estável
1)0(, problema oResolver  yyy x
k
x
kk
x
kkkkkkk
x
yyhy
yhyyyxfhyy
yyxF
k
k







1
1
11
11111 ),(
),(
O valor de yk+1 não pode ser explicitada em função de yk
Para cada passo, deve-se determinar a raiz da equação não linear 
conhecidos ,,;0 cbacayy b 
Cálculo pelo Método de Newton, por exemplo
   
)(1
)()()()1(
)(
)0(
ofim_equant1
repetir ,)( Enquanto0
0)(
j
k
jjjj
j
k
b
yy
jj
yfyfyy
yf
j
yy
cayyyf










MÉTODO IMPLÍCITO DE SEGUNDA ORDEM
(MÉTODO DO TRAPÉZIO)
x
y
0x 1x 2x
1y
2yTeorema do Valor Médio:
 
   ),(),(21)()(21)(Tomar 
)(
)()()( que tal,
111
1
1
11









kkkkkk
kk
xk
kk
kk
yxfyxfxyxyy
yhyy
xx
xyxyyxx



    111 ,,2   kkkkkk yxfyxfhyy
Método de segunda ordem (erro decai com h2)
Método Estável
Mesmas dificuldades do Método de Primeira Ordem (Euler)
Exercício
Escreva uma rotina SciLab para solução de um problema de valor inicial usando o método de Euler implícito de primeira ordem. Utilize a rotina desenvolvida para resolver o problema:
Determine o valor de para diferentes passos de tempo. 
MÉTODOS PREDITOR-CORRETOR Na resolução de uma equação diferencial deve-se decidir entre o uso deum método explícito, mais fácil e não estável, e o uso de ummétodo implícito, mais difícil, porém estáveis.
Quando a equação diferencial é não linear, deve-se usar técnicas iterativaspara resolver a equação não linear resultante em cado passo.
Para o Método de Newton, o chute inicial deve ser bom para o processoconvergir em poucas iterações.
Uma opção é usar um método explícito para obter o chute inicial doprocesso iterativo (Preditor) e um método implícito para obtera solução (Corretor).
Método de Heun: Preditor: Euler e Corretor: Trapézio
    111* 1
* 1
,,2 :Inicial Chute
),(




kkkkkkk
kkkk
yxfyxfhyyy
yxfhyy
MÉTODOS DE MÚLTIPLOS PASSOS 
 Nos métodos vistos até agora, o valor yk+1 depende somente de informações no instante anterior k no instante presente k+1.
Métodos de Múltiplos Passos:Usar informações em vários pontos anteriorespara obter maior precisão.
Passo Único: ),()( 11 kkkkkk yxfh yyhOh yydxdy  
Dois Passos:
),(2
),(2)(2
11
11311
kkkk
kk
kkkk
yxfhyy
yxf
h
yyhO
h
yy
dx
dy




)()(2)()(
)(6)(2)()()(
)(6)(2)()()(
311
32
1
32
1
hOxyhxyxy
xyhxyhxyhxyxy
xyhxyhxyhxyxy
kkk
kkkkk
kkkkk








MÉTODO LEAPFROG
 Idéia geral dos Métodos de Múltiplos Passos: Usar informações em vários pontos anteriores a xk para descrever como a função se comportaentre xk e xk+1.
 








1
11
)(
)(,)()()(
1
1
k
k
k
k
k
k
x
x
kk
x
x
x
x
kk
dxxpyy
dxxyxfdxxyxyxy
 p(x) é um polinômio interpolador de grau N que aproxima f(x,y) e passapelo conjundo de dados Nkkkifx ii  ,,1,,),( 
Para N = 0: método de passo único
),(
)(
1
1
kkkk
x
x
kkk
yxfhyy
hfdxffxp
k
k



Método de Euler
Para N = 1: método de Dois Passos
Método de Adams-Bashforth de 2a ordem
p(x) é um polinômio linear que interpola ),( e ),( 11 kkkk fxfx 
 
   22 )()(
)(
1
2
1
11
111
1
1
hffhf
h
xxffhfdxxp
ff
h
xxff
xx
xxf
xx
xxxp
kkk
x
x
k
kkk
x
x
kk
k
kk
kk
k
k
kk
k
k
k
k
k













   1111 322   kkkkkkkkk ffhyyffhfhyy
Para N = 2: método de Três Passos — Método de Adams-Bashforth 3a ordem
p(x) é um polinômio quadrático que interpola ),( e ),(,),( 1122 kkkkkk fxfxfx 
 211 5162324   kkkkk fffhyy
Métodos de Adams-Bashforth são Métodos Explícitos
 Pode-se formar os polinômios interpoladores utilizando-se pontos para frente.Métodos Implícitos
 Situação mais comum é formar um polinômio de ordem N com os pontos Nkkkk xxxx  ,,, 11
Método de Adams-Moulton de ordem N
Para N = 0: método de passo único
Regra do Trapézio    111 ,,2   kkkkkk yxfyxfhyy
Para N = 3: Método de Adams-Moulton de 4a ordem
        2211111 ,,5,19,924   kkkkkkkkkk yxfyxfyxfyxfhyy
SOLUÇÃO DE EQUAÇÕES RÍGIDAS (STIFF)
EQUAÇÕES INSTÁVEIS





1)0( 1)0(c.c.
01110
y
y
yyy Solução exata: )1()( xexy 





1)0( 1)0(c.c.
01110
y
y
yyy

Solução exata: )2(1212111)( 11xx eexy    (2)
(1)
0
2
4
6
8
10
0 0.2 0.4 0.6 0.8 1 1.2
Microsoft Excel 
Worksheet
 O problema (1) é instável.
 Pequenas alterações na condição inicial podem produzir grandesalterações na solução do problema para x grande.
 Esses problemas são chamados de mal condicionados.
A solução numérica é extremamente difícil, pois erros de arredondamentoe da discretização podem causar o mesmo efeito que a pequenamudança na condição inicial do problema e a solução tenderá adivergir para o infinito.
 Estes problemas requerem método numéricos estáveis com passos bemmenores do que o usual.
As instabilidades são mais pronunciadas em problemas não-lineares.
2y(0)c.c.
)2( :Exemplo

 yxyy Solução exata: )1(2)( xy
Este problema é instável. A solução para condição inicial é:0)0( yy 
)2()2(
2)( 200
0
xeyy
yxy 
1
1.5
2
2.5
3
0 0.2 0.4 0.6 0.8 1 1.2
Microsoft Excel 
Worksheet
MÉTODOS ESTÁVEIS E INSTÁVEIS



 1)0( 12 :Problema y yy Solução exata: )1(2
121)( 2   xexy
Este problema é estável, pois a solução não muda muito alterando-se a c.c.
21 21)(1)0( Se 2    xexyy 
Aplicando-se o método de Leapfrog (segunda ordem, dois passos, explícito):
  hyhyyhyy kkkkk 24122 111  

 
ky
eyy
k
h
 quando 2
121;1 210
-10
-5
0
5
10
0 2 4 6 8
Método implícitos são estáveis e portanto devem ser usados para problemas stiff
Resolver o problema:   0)0(;)sin(100  yyxy
Solução exata: 0001.1 01.0)cos(01.0)sin()(
100xexxxy

 Solução por Runge-Kutta:
Metodo de Runge-Kutta de quarta ordem:
h = 0.0300 y(0) = 0.0000
X f1 f2 f3 f4 Yk Yexato
0.00 0.0000 00.03 0.0000 1.4999 -0.7500 5.2495 0.033747047 0.0204960.06 -0.3752 1.6865 -1.4060 6.8397 0.068874772 0.0500020.09 -0.8911 1.9421 -2.3077 9.0234 0.105880707 0.0799120.12 -1.6002 2.2930 -3.5468 12.0236 0.14545912 0.1097730.15 -2.5747 2.7752 -5.2496 16.1467 0.188574804 0.1395360.18 -3.9137 3.4383 -7.5896 21.8144 0.236564525 0.1691740.21 -5.7535 4.3504 -10.8055 29.6059 0.291276497 0.198660.24 -8.2817 5.6055 -15.2252 40.3183 0.355262162 0.2279660.27 -11.7560 7.3323 -21.3001 55.0471 0.43203987 0.2570680.30 -16.5308 9.7080 -29.6503 75.2989 0.526457436 0.2859380.33 -23.0937 12.9765 -41.1288 103.1450 0.645190639 0.3145510.36 -32.1148 17.4727 -56.9085 141.4339 0.797428656 0.3428810.39 -44.5154 23.6576 -78.6019 194.0818 0.995816526 0.3709020.42 -61.5628 32.1644 -108.4264 266.4737 1.257751018 0.398590.45 -84.9991 43.8645 -149.4308 366.0140 1.607162455 0.425918
M eto do de R unge-Kutta de quarta o
h = 0.03 y(0) = 0
X f1 f2 f3 f4 Yk Yexato
0 =$E$3 =(SIN(A7)-0.01*COS(A7)+0.01*
=A7+$B$3 =100*(SIN(A7)-F7) =100*(SIN(A7+$B$3/2)-(F7+$B$3/2*B8)) =100*(SIN(A7+$B$3/2)-(F7+$B$3/2*C8)) =100*(SIN(A7+$B$3)-(F7+$B$3*D8)) =F7+$B$3/6*(B8+2*C8+2*D8+E8) =(SIN(A8)-0.01*COS(A8)+0.01*
=A8+$B$3 =100*(SIN(A8)-F8) =100*(SIN(A8+$B$3/2)-(F8+$B$3/2*B9)) =100*(SIN(A8+$B$3/2)-(F8+$B$3/2*C9)) =100*(SIN(A8+$B$3)-(F8+$B$3*D9)) =F8+$B$3/6*(B9+2*C9+2*D9+E9) =(SIN(A9)-0.01*COS(A9)+0.01*
=A9+$B$3 =100*(SIN(A9)-F9) =100*(SIN(A9+$B$3/2)-(F9+$B$3/2*B10)) =100*(SIN(A9+$B$3/2)-(F9+$B$3/2*C10)) =100*(SIN(A9+$B$3)-(F9+$B$3*D10)) =F9+$B$3/6*(B10+2*C10+2*D10+E10) =(SIN(A10)-0.01*COS(A10)+0.0
=A10+$B$3 =100*(SIN(A10)-F10) =100*(SIN(A10+$B$3/2)-(F10+$B$3/2*B11=100*(SIN(A10+$B$3/2)-(F10+$B$3/2*C11=100*(SIN(A10+$B$3)-(F10+$B$3*D11)) =F10+$B$3/6*(B11+2*C11+2*D11+E11) =(SIN(A11)-0.01*COS(A11)+0.0
=A11+$B$3 =100*(SIN(A11)-F11) =100*(SIN(A11+$B$3/2)-(F11+$B$3/2*B12=100*(SIN(A11+$B$3/2)-(F11+$B$3/2*C12=100*(SIN(A11+$B$3)-(F11+$B$3*D12)) =F11+$B$3/6*(B12+2*C12+2*D12+E12) =(SIN(A12)-0.01*COS(A12)+0.0
=A12+$B$3 =100*(SIN(A12)-F12) =100*(SIN(A12+$B$3/2)-(F12+$B$3/2*B13=100*(SIN(A12+$B$3/2)-(F12+$B$3/2*C13=100*(SIN(A12+$B$3)-(F12+$B$3*D13)) =F12+$B$3/6*(B13+2*C13+2*D13+E13) =(SIN(A13)-0.01*COS(A13)+0.0
=A13+$B$3 =100*(SIN(A13)-F13) =100*(SIN(A13+$B$3/2)-(F13+$B$3/2*B14=100*(SIN(A13+$B$3/2)-(F13+$B$3/2*C14=100*(SIN(A13+$B$3)-(F13+$B$3*D14)) =F13+$B$3/6*(B14+2*C14+2*D14+E14) =(SIN(A14)-0.01*COS(A14)+0.0
=A14+$B$3 =100*(SIN(A14)-F14) =100*(SIN(A14+$B$3/2)-(F14+$B$3/2*B15=100*(SIN(A14+$B$3/2)-(F14+$B$3/2*C15=100*(SIN(A14+$B$3)-(F14+$B$3*D15)) =F14+$B$3/6*(B15+2*C15+2*D15+E15) =(SIN(A15)-0.01*COS(A15)+0.0
=A15+$B$3 =100*(SIN(A15)-F15) =100*(SIN(A15+$B$3/2)-(F15+$B$3/2*B16=100*(SIN(A15+$B$3/2)-(F15+$B$3/2*C16=100*(SIN(A15+$B$3)-(F15+$B$3*D16)) =F15+$B$3/6*(B16+2*C16+2*D16+E16) =(SIN(A16)-0.01*COS(A16)+0.0
=A16+$B$3 =100*(SIN(A16)-F16) =100*(SIN(A16+$B$3/2)-(F16+$B$3/2*B17=100*(SIN(A16+$B$3/2)-(F16+$B$3/2*C17=100*(SIN(A16+$B$3)-(F16+$B$3*D17)) =F16+$B$3/6*(B17+2*C17+2*D17+E17) =(SIN(A17)-0.01*COS(A17)+0.0
Microsoft Excel 
Worksheet
O método só é estável para passos muito pequenos h < 0.030
 Solução por Trapézio:
    
    
 
10021
)sin()sin(100221001
)sin(100)sin(1002
,,2
1
1
111
111
h
xxhhy
y
yxyxhyy
yxfyxfhyy
kkk
k
kkkkkk
kkkkkk



 







Metodo do Trapezio
h = 0.5000 y(0) = 0.0000
X Yk Yexato
0.00 0.0000 00.50 0.460986095 0.4706026531.00 0.844567185 0.8359843631.50 0.988636033 0.9966879462.00 0.920867137 0.9133675582.50 0.59974723 0.6064229383.00 0.157533472 0.1510048333.50 -0.347014762 -0.3413845224.00 -0.744664953 -0.750191044.50 -0.980244479 -0.9753246275.00 -0.95713432-0.96166473
M eto do do T rapezio
h = 0.5
X Yk
0 =$E$3
=A7+$B$3 =(B7*(1-100*$B$3/2)+$B$3/2*100*(SIN(A7)+SIN(A8)))/(1+$B$3/2*100)
=A8+$B$3 =(B8*(1-100*$B$3/2)+$B$3/2*100*(SIN(A8)+SIN(A9)))/(1+$B$3/2*100)
=A9+$B$3 =(B9*(1-100*$B$3/2)+$B$3/2*100*(SIN(A9)+SIN(A10)))/(1+$B$3/2*100)
=A10+$B$3 =(B10*(1-100*$B$3/2)+$B$3/2*100*(SIN(A10)+SIN(A11)))/(1+$B$3/2*100)
=A11+$B$3 =(B11*(1-100*$B$3/2)+$B$3/2*100*(SIN(A11)+SIN(A12)))/(1+$B$3/2*100)
=A12+$B$3 =(B12*(1-100*$B$3/2)+$B$3/2*100*(SIN(A12)+SIN(A13)))/(1+$B$3/2*100)
=A13+$B$3 =(B13*(1-100*$B$3/2)+$B$3/2*100*(SIN(A13)+SIN(A14)))/(1+$B$3/2*100)
=A14+$B$3 =(B14*(1-100*$B$3/2)+$B$3/2*100*(SIN(A14)+SIN(A15)))/(1+$B$3/2*100)
=A15+$B$3 =(B15*(1-100*$B$3/2)+$B$3/2*100*(SIN(A15)+SIN(A16)))/(1+$B$3/2*100)
=A16+$B$3 =(B16*(1-100*$B$3/2)+$B$3/2*100*(SIN(A16)+SIN(A17)))/(1+$B$3/2*100)
Microsoft Excel 
Worksheet
PROBLEMA DE VALOR DE CONTORNO
ByAyyyxf
dx
yd  )1(;)0( :cont. cond.,),,(22
MÉTODO DE TIRO (SHOOTING METHOD) 
 Baseado nos Problemas de Valor Inicial.
 Estima-se a derivada em x = 0:
 Utilizando-se os métodos de valor inicial, determina-se a solução em x = 1.
 Verifica-se se o valor obtido em x = 1 é igual a B.
 Repetir até convergir.
)0(y
A
B
Problema Linear
 Usar o Princípio da Superposição.
)( e )( 21 xyxy : Funções que satisfazem a eq. Diferencial e a cond. de contorno em x = 0. 
)0()0( )0()0( 21 21 yy Ayy  
)()()( 2211 xycxycxy  satisfaz a eq. Diferencial 
Aycc
ycycy

 )0(1 Se )0()0()0( 21 2211
)(
)1()1()1( 2211
xy
Bycycy 
é solução do problema se:
1221
212211
21 1)1()1( )1()1()1(
1
cce
yy
yBc
Bycyc
cc 





Problema Não-Linear
 O Princípio da Superposição não pode ser usado.
 Cada valor de corresponde um valor de .

 Utilizar métodos iterativos para determinar o valor de que fornece a condição de contorno correta em x = 1, isto édeterminar a raiz da equação:
 )0(y )1(y
)()1( gy 
Bgf  )()( 
 Por exemplo: Método da Secante:
 

)()()(
)()( Chute
)()( Chute
01
01112
111
000




gg
Bg
Bgf
Bgf




 O Problema de Valor Inicial pode ser mal-condicionado mesmo queo Problema de Valor de Contorno seja bem-condicionado.
MÉTODO DE DIFERENÇAS FINITAS
A equação diferencial é transformada em um conjunto de equações algébricas
 A função desconhecida é calculada apenas em N pontos – NÓS
 Derivadas são aproximadas por diferença dos valores nodais
 A equação aproximada é escrita em cada ponto nodal, gerando N equações
Exemplo:
LYLy
Yy
Lx
dx
dyyxf
dx
yd





)( )0(
0;,,
0
2
2
Obter y(x).
i=1 i=2 i=Ni i+1i-1
1x 1ix ix 1ix Nxhx 
Aproximação das derivadas por uso de série de Taylor truncada:
)(2)(2
2
1
2
1 ByhyhyyAyhyhyy iiiiiiii   
Diferentes aproximações para primeira derivada:
 DIFERENÇA PARA FRENTE
 DIFERENÇA PARA TRÁS
 DIFERENÇA CENTRAL
h
yyyA iii
 1)(
h
yyyB iii 1)( 
h
yyyBA iii 2)()( 11  
Problema Linear
Aproximação para segunda derivada (Deve-se eliminar o termo de 1a ordem):
 DIFERENÇA CENTRAL
2 11
2)()(
h
yyyyBA iiii 

i i+1i-1
1ix ix 1ix
hx 
e d
Interpretação geométrica 
2 112
2
11
2
2
2
h
yyy
dx
yd
h
h
yy
h
yy
h
dx
dy
dx
dy
dx
dy
dx
d
dx
yd
iii
iiii
ed









EXEMPLO: Condução de calor em uma barra com convecção natural
AT BT
Thc ,
Lx
Equação diferencial que descrevea variação da temperatura na barra.
 
B
A
c
TLxT
TxT
TT
kA
Ph
dx
Td


 
)( )0(
022
i=1 i=2 i=Ni i+1i-1
1x 1ix ix 1ix Nx
Incógnitas do problema: Ni TTTT ,,,,, 21 
No ponto i:  
 








 




T
kA
PhT
h
T
kA
Ph
h
T
h
TT
kA
Ph
h
TTT
TT
kA
Ph
dx
Td
c
ii
c
i
i
ciii
i
c
i
12212
2 11
2
2
121
02
0
Para se obter as N equações necessárias para determinar as N incógnitas,a aproximação deve ser satisfeita em todos os nós:
(c.c.) 
1213
1212
(c.c.) 1
423222
322212
1
BN
cc
cc
A
TTNi
T
kA
PhT
h
T
kA
Ph
h
T
h
i
T
kA
PhT
h
T
kA
Ph
h
T
h
i
TTi





 






 


















































 


 



B
c
c
c
A
N
c
c
T
TkA
Ph
TkA
Ph
TkA
Ph
T
T
T
T
T
T
hkA
Ph
hh
hkA
Ph
hh





4
3
2
1
222
222
10000
0121
00121 00001
O sistema pode ser escrito em forma matricial como:
Exemplo: Desenvolver o problema se os nós não são uniformemente distribuidos.
i=1 i=2 i=Ni i+1i-1
1x 1ix ix 1ix Nx
1ix ix
    


TLxT
K
hLx
dx
dT
TxT
dx
Td
e
A)0(
022
EXEMPLO: Condição de contorno na derivada da função
AT Thc ,
L
x
No ponto i:
0121
02
0
12212
2 11
2
2











iii
iii
i
T
h
T
h
T
h
h
TTT
dx
Td
Para se obter as N equações necessárias para determinar as N incógnitas,a aproximação deve ser satisfeita em todos os nós:
   

 

















T
K
h
K
h
h
T
h
TT
K
h
h
TTNi
T
h
T
h
T
h
i
T
h
T
h
T
h
i
TTi
ee
NN
eNN
A
11
01213
01212
(c.c.) 1
11
423222
322212
1

   









































TK
h
T
T
T
T
T
T
K
h
hh
hhh
hhh
e
A
N
e






00
0
110000
00121
000121 000001
4
3
2
1
222
222
O sistema pode ser escrito em forma matricial como:
Exercício
Escreva uma rotina SciLab para solução de um problema de valor de contorno pelo método de diferenças finitas (diferença central). Utilize a rotina desenvolvida para resolver o problema:
L = 5 m
Problema Não-Linear
EXEMPLO: Problema de Convecção e Difusão 
AV BV
Lx
Equação diferencial que descrevea velocidade em cada ponto.
B
A
VLxu
Vxu
dx
udK
dx
duu



)( )0(
022
i=1 i=2 i=Ni i+1i-1
1x 1ix ix 1ix Nx
Incógnitas do problema: Ni uuuu ,,,,, 21 
Aproximação por diferença central:
2 112
2
11
2
2
h
uuu
dx
ud
h
uu
dx
du
iii
i
ii
i




No ponto i:
022
022
0
1211212
2 1111
2
2




 

 




 

 





i
Ai
i
Bi
ii
i
Ai
iiiii
i
ii
i
u
h
Ku
h
uu
h
Ku
h
K
h
uuuK
h
uuu
dx
udK
dx
duu
  
Para se obter as N equações necessárias para determinar as N incógnitas,a aproximação deve ser satisfeita em todos os nós:
(c.c.) 
0223
0222
(c.c.) 1
42324222
32213212
1
BN
A
VuNi
u
h
Ku
h
uu
h
Ku
h
Ki
u
h
Ku
h
uu
h
Ku
h
Ki
Vui





 

 





 

 







































B
A
N T
V
u
u
u
u
u
ABA
ABA





00
0
10000
000
00001
4
3
2
1
333
222
Função de 
u1 e u3
Os coeficientes da matriz dependem da solução do problema.
Sistema de equações Não-Linear
Solução de Sistema de Equações Não-Linear
Método de Picard:
1. Chute inicial;
2. Calcular coeficientes da matriz usando o valor atual das incógnitas;
3. Resolver o sistema de equações e determinar o novo valor dasincógnitas;
• Comparar solução atual com anterior;
• Se não covergiu, voltar para 2.
 )0()0(3)0(2)0(1)0( ,,,, Nuuuuc 
 )(kcAA 
   fcAc kk 1)()1(  
Convergência Ruim
Método de Newton:
Generalização do Método de Newton para 1 equação não-linear
       
     
 
 xf
xfx
xfxxfxxf
xfxxfxxfxxf





 0
2 
PROCEDIMENTO ITERATIVO
)1(
)()1(
)(
)(
)(
)0(
 :Raiz 1
)( )(
do ,)( While0
 :inicial Chute

 




i
ii
i
i
i
x
ii
xxx
xf
xf
x
xf
i
x












0),,,,(
0),,,,( 0),,,,(
0),,,,(
321
3213
3212
3211
NN
N
N
N
xxxxf
xxxxf
xxxxf
xxxxf




Sistema a ser resolvido:
Expansão por série de Taylor até termos de primeira ordem de cada equação:
N
N
NNN
NN
NNN
N
N
N
NN
N
N
N
NN
x
x
fx
x
fx
x
fxxxf
xxxxxxf
x
x
fx
x
fx
x
fxxxf
xxxxxxf
x
x
fx
x
fx
x
fxxxf
xxxxxxf




























221121
2211
222
211
2212
22112
122
111
1211
22111
),,,(
0),,,(
),,,(
0),,,(
),,,(
0),,,(
N
N
NNN
NN
N
N
N
N
N
N
x
x
fx
x
fx
x
fxxxf
x
x
fx
x
fx
x
fxxxf
x
x
fx
x
fx
x
fxxxf



















221121
222
211
2212
122
111
1211
),,,(
),,,(
),,,(

fJx
f
f
f
x
x
x
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
f
N
x
N
J
N
NNN
N
N
1
2
1
2
1
21
2
2
2
1
2
1
2
1
1
1
























































  



Sistema em
Forma matricial
Matrix Jacobiana
j
i
ij x
fJ 

PROCEDIMENTO ITERATIVO
 
)1(
)()1(
1
)(
)0(
 :Raiz
1
do , While
0 :inicial Chute








i
ii
i
c
ii
xcc
fJx
cf
i
c



Solução de um sistema linear
Convergência Quadrática
Voltando ao Problema Não-Linear:
B
A
VLxu
Vxu
dx
udK
dx
duu



)( )0(
022
 
 
  0,,,,,, ,,3,2
022,,,,,,
0,,,,,,
111
2 1111111
11111




 

 





BNNiiiN
iiiii
iNiiii
ANiii
Vuuuuuuf
Ni
h
uuuK
h
uuuuuuuuf
Vuuuuuuf




Sistema de equações algébricas não-linear resultante:
Solução pelo Método de Newton…
Cálculo da matriz Jacobiana:
1;0;;0
;0;2;22;2;0
0;;0;2;22;2
0;;0;1
11
5
3324
32423
3322
3
1
3
2
4
2223
21322
2221
2
1
2
1
1
1




















 N
N
N
NN
N
N
u
f
u
f
u
f
u
f
h
u
h
K
u
f
h
uu
h
K
u
f
h
u
h
K
u
f
u
f
u
f
u
f
h
u
h
K
u
f
h
uu
h
K
u
f
h
u
h
K
u
f
u
f
u
f
u
f





Exercício