Equações Diferenciais Ordinárias - Cálculo Numérico - Instituto Tecnológico de Aeronáutica

Prévia do material em texto

CCI-22
á lMatemática Computacional
Carlos Alberto Alonso Sanches
Juliana de Melo Bezerra
CCI-22
8) E õ Dif i i 8) Equações Diferenciais 
OrdináriasOrdinárias
Método de Euler Séries de Taylor Runge-KuttaMétodo de Euler, Séries de Taylor, Runge Kutta,
Adams-Bashforth, Adams-Moulton, Diferenças Finitas
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
DefiniçõesDefiniçõesf çf ç
ƒ Grande parte dos fenômenos físicos é modelada com 
õ dif i i i t é l f ã equações diferenciais, isto é, envolve uma função 
desconhecida e algumas de suas derivadas
ƒ Forma geral de uma equação diferencial com derivadas ƒ Forma geral de uma equação diferencial com derivadas 
até a ordem n:
y(n)(x) = f(x y(x) y’(x) y(n-1)(x)) onde a ≤ x ≤ by(n)(x) = f(x, y(x), y (x), ..., y(n )(x)), onde a ≤ x ≤ b
ƒ A solução desta equação diferencial é qualquer função 
y(x) que a satisfaça definida em [a b] e com n y(x) que a satisfaça, definida em [a,b] e com n 
derivadas nesse intervalo
ƒ Quando y é função de uma única variável x, a equação Quando y é função de uma única variável x, a equação 
acima é chamada de Equação Diferencial Ordinária ; 
caso contrário, seria uma Equação Diferencial Parcial
ƒ Neste curso, abordaremos métodos numéricos de 
resolução de Equações Diferenciais Ordinárias (EDOs)
Condições iniciaisCondições iniciaisçç
ƒ A resolução de uma equação diferencial geralmente tem como 
t f íli d resposta uma família de curvas
ƒ Exemplo:
ƒ y’ = 2x + 3y = 2x + 3
ƒ ∫y’dx = ∫(2x + 3)dx ⇒ y = x2 + 3x + c
ƒ Para especificar uma dessas curvas, é preciso impor condições 
iniciais à função y:
ƒ y(t1) = k1
ƒ y’(t ) = kƒ y (t2) = k2
ƒ ... 
ƒ y(n-1)(tn-1) = kn-1
ƒ Exemplo:
ƒ y’’ = -(1 - y2)y’ - y
(0) 1ƒ y(0) = 1
ƒ y’(0) = 2 
PVI e PVCPVI e PVC
ƒ A ordem de uma equação diferencial é a mais alta ordem de 
derivação que aparece neladerivação que aparece nela
ƒ De modo geral, para individualizar a solução de uma equação 
diferencial de ordem m, são necessárias m condições adicionais
ƒ Dada uma Equação Diferencial Ordinária de ordem m ≥ 1, se a 
função e suas derivadas até a ordem m-1 são especificadas em um 
mesmo ponto então temos um Problema de Valor Inicial (PVI)mesmo ponto, então temos um Problema de Valor Inicial (PVI)
ƒ Exemplo onde m = 3:
ƒ y’’’ + (x + 1)y’’ + cos(xy’) – (x2 - 1)y = x2 + y2sen(x + y)
(0) 1 1 ’(0) 2 2 ’’(0) 3 3ƒ y(0)=1,1; y’(0)=2,2; y’’(0)=3,3
ƒ Se as m condições adicionais não são dadas em um mesmo ponto, 
então temos um Problema de Valor de Contorno (PVC)
ƒ Exemplo (barra de comprimento L sujeita a uma carga uniforme q):
ƒ y(4)(x) + ky(x) = q
ƒ y(0) = y’(0) = 0; y(L) = y’’(L) = 0
k é uma constante que 
depende do material da barray(0) y (0) 0; y(L) y (L) 0
ƒ Ao contrário de um PVI, é comum que um PVC não tenha unicidade 
de solução
p m
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Problemas de Valor InicialProblemas de Valor Inicialmm
ƒ Embora haja garantia teórica da resolução analítica de 
 PVI ss s l ã st s d difí il bt ã : um PVI, essa solução costuma ser de difícil obtenção: 
por isso, utilizam-se métodos numéricos 
ƒ Dado o PVI y’ = f(x y) onde y(x ) = y construímos ƒ Dado o PVI y = f(x,y), onde y(x0) = y0, construímos 
x1, x2, ..., xn igualmente espaçados (embora não seja uma 
condição necessária), e calculamos as aproximações ç ), p ç
yi ≈ y(xi) nesses pontos
ƒ Se no cálculo de yi+1 usarmos apenas yi, teremos então 
um método de passo simples (ou passo um); se usarmos 
outros valores yj, j ≤ i, teremos um método de passo 
múltiplomúltiplo
ƒ Características dos métodos de passo simples:
ƒ Geralmente é preciso calcular f(x y) e suas derivadas em vários ƒ Geralmente, é preciso calcular f(x,y) e suas derivadas em vários 
pontos
ƒ Temos dificuldades em estimar o erro do resultado
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Método de EulerMétodo de Euler
ƒ Vamos resolver a equação diferencial ordinária de 
 d ’ f( ) à d l primeira ordem y’ = f(x,y), sujeita à condição inicial 
y(x0) = y0:
)(fdy )(0 fyy −
y(x)
y 
ŷ
),(
),(
00
yx
yxf
dx
y
00
= ),( 00
0
0 yxf
xx
yy =−⇒
ƒ Sejam y = y(x ) e h = x - x
y1
y0
ƒ Sejam y1 = y(x1) e h = x1 - x0
ƒ Equação da reta tangente a (x0,y0):
ƒ ŷ = y0 + h.f(x0,y0)
x1 xx0 
ƒ Quando h tende a zero, ŷ tende a y1:
ƒ y1 ≈ y0 + h.f(x0,y0) 
ƒ Importante: h deve ser pequeno
ƒ Generalizando temos a expressão do Método de Euler:
1 0 ƒ Importante: h deve ser pequeno...
h
ƒ Generalizando, temos a expressão do Método de Euler:
yi+1 ≈ yi + h.f(xi, yi)
Método de EulerMétodo de Euler
ƒ A expressão do Método de Euler pode ser deduzida de 
outro modo
ƒ Sabemos que y’(x) ≈ [y(x+h) – y(x)]/h, onde h é algum 
l ã fivalor pequeno, mas não fixo
ƒ Dividamos [a,b], onde a=x0 e b=xn, em n subintervalos 
d t h hde tamanho h:
ƒ xi = x0 + h.i, com 0 ≤ i ≤ n
S j 0 i i ã ( ) d ( ) é ƒ Seja yi, 0≤i≤n, uma aproximação para y(xi), onde y(x) é 
uma solução de y’(x) = f(x,y)
P t tƒ Portanto:
ƒ y’(xi) ≈ (yi+1 – yi)/h
 h ’( )ƒ yi+1 ≈ yi + h.y’(xi)
ƒ yi+1 ≈ yi + h.f(xi,yi)
ExemploExemplompmp
ƒ Considerando y como função de x, vamos resolver y’ = 2x + 3 no 
intervalo 1 ≤ x ≤ 1 5 quando y(1) 1intervalo 1 ≤ x ≤ 1,5, quando y(1) = 1
ƒ Pelo Método de Euler, temos:
 h f( ) ƒ yi+1 ≈ yi + h.f(xi, yi) 
ƒ yi+1 ≈ yi + h.(2xi + 3) 
C id d h 0 1
Repare que somente yi é usado no 
cálculo de yi+1 (passo simples)
ƒ Considerando h = 0,1:
x0 = 1,0 x1 = 1,1 x2 = 1,2 x3 = 1,3 x4 = 1,4 x5 = 1,5
y0 = 1 0 y1 = 1 5 y2 = 2 02 y3 = 2 56 y4 = 3 12 y5 = 3 70y0 = 1,0 y1 = 1,5 y2 = 2,02 y3 = 2,56 y4 = 3,12 y5 = 3,70
ƒ Considerando h = 0,01:
x0 = 1,00 x10 = 1,10 x20 = 1,20 x30 = 1,30 x40 = 1,40 x50 = 1,50
y0 = 1,0 y10 = 1,509 y20 = 2,038 y30 = 2,587 y40 = 3,156 y50 = 3,747
ƒ Os resultados não se alteraram muito. Mais adiante, veremos uma 
estimativa para os erros cometidosCCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Métodos de série de TaylorMétodos de série de Tayloryy
ƒ Suponhamos que, de alguma maneira, estejam disponíveis as 
i õ ( ) aproximações y1, y2, ..., yn para y(x) em x1, x2, ..., xn, 
respectivamente
ƒ A série de Taylor de k-ésima ordem de y(x) em torno de x = xi éA série de Taylor de k ésima ordem de y(x) em torno de x = xi é
T
k
i
i
)k(
2
i
i
iii E)xx(!k
)x(y)xx(
!2
)x(''y)xx)(x('y)x(y)x(y +−++−+−+= L
onde 1ki
)1k(
T )xx()!1k(
)ξ(yE +
+
−+= , ξ entre xi e x
ƒ Considerando xi+1 = xi + h, podemos obter a seguinte aproximação 
para yi+1 = y(xi+1):
hh k2
!k
hy
2
h''yh'yyy
k
)k(
i
2
iii1i ++++≈+ L
É fá il ifi é i d T l d 1ª d é i l ƒ É fácil verificar que a série de Taylor de 1ª ordem é equivalente 
ao Método de Euler: yi+1 ≈ yi + y’ih
Métodos de série de TaylorMétodos de série de Tayloryy
ƒ Para se encontrar as séries de Taylor de ordens mais altas, será 
i l l l d ’’( ) ’’’( ) (k)( )preciso calcular os valores de y’’(x), y’’’(x), ..., y(k)(x)
ƒ Considerando y’(x) = f(x,y(x)), vamos calcular y’’(x):
ƒ y’’(x) = f’(x y(x))ƒ y (x) = f (x,y(x))
ƒ y’’ = f’ = fx + fy.f, onde fx e fy são, respectivamente, as derivadas 
parciais de f em relação a x e y
ƒ Desse modo, a série de Taylor de 2ª ordem é 
yi+1 ≈ yi + h.f(xi,yi) + h2[fx(xi,yi) + fy(xi,yi).f(xi,yi)]/2
V m s l l ’’’(x):ƒ Vamos calcular agora y (x):
ƒ y’’’ = fxx + fxy.f + (fyx + fyy.f).f + f’.fy
ƒ y’’’ = f + f f + f f + f f2 + (f + f f) fy = fxx + fxy.f + fyx.f + fyy.f + (fx + fy.f).fy
ƒ y’’’ = fxx + 2fxy.f + fyy.f2 + fx.fy + fy2.f
ƒ É possível perceber como se torna difícil o cálculo de derivadas p p
mais altas. Isso somente será viável quando y’(x) tiver uma 
expressão simples...
ExemploExemplompmp
ƒ Usando a série de Taylor de 2ª ordem, vamos calcular 
( ) d ’ ( ) y(2,1) onde xy’ = x – y e y(2) = 2
ƒ xy’ = x – y ⇔ y’ = (x - y)/x ⇔ y’ = 1 – y/x
ƒ y’’ = -y’/x + y/x2
ƒ Série de Taylor de 2ª ordem em torno de x = 2:ƒ Série de Taylor de 2 ordem em torno de x = 2:
ƒ y(x) ≈ y(2) + (x - 2)y’(2) + (x - 2)2y’’(2)/2
’(2) 1 2/2 0
Expressão simples para 
y’: é viável utilizar ƒ y (2) = 1 – 2/2 = 0
ƒ y’’(2) = 0/2 + 2/22 = 1/2
( ) 2 ( 2)2/4
y : é viável utilizar 
método de série de 
Taylor de 2ª ordem
ƒ y(x) ≈ 2 + (x - 2)2/4
ƒ y(2,1) ≈ 2 + (0,1)2/4 
ƒ y(2,1) ≈ 2,0025
ExemploExemplompmp
ƒ Dado que y’ = x – y e y(0) = 2, determinar y(0,2) e y(0,4) 
l d é d l d ª dutilizando série de Taylor de 4ª ordem
ƒ Sabemos que x0 = 0 e y0 = 2, e consideraremos h = 0,2
ƒ Cálculos das derivadas (viáveis porque são simples):ƒ Cálculos das derivadas (viáveis porque são simples):
ƒ y’ = x – y ⇒ y0’ = y’(0) = 0 – 2 = -2
ƒ y’’ = 1 - y’⇒ y0’’ = y’’(0) = 1 – (-2) = 3
’’’ ’’ ’’’ ’’’( ) ƒ y’’’ = -y’’ ⇒ y0’’’ = y’’’(0) = -3
ƒ y(4) = -y’’’⇒ y0(4) = y(4)(0) = 3
ƒ Série de Taylor de 4ª ordem (em torno de x0 = 0 e x1 = 0,2):y 0 1
ƒ y1 = y(0,2) ≈ y(0) + h.y’(0) + h2y’’(0)/2 + h3y’’’(0)/6 + h4y(4)(0)/24 
ƒ y1 ≈ y0 + h.y0’ + h2y0’’/2 + h3y0’’’/6 + h4y0(4)/24 ≈ 1,6552
ƒ y2 = y(0,4) ≈ y1 + h.y1’ + h2y1’’/2 + h3y1’’’/6 + h4y1(4)/24 y2 y( , ) y1 y1 y1 y1 y1
• y1’ = x1 – y1 ≈ 0,2 – 1,6552 = -1,4562
• y1’’ = 1 - y1’ ≈ 1 – (-1,4562) ≈ 2,4562
• y1’’’ = -y1’’ ≈ -2,4562 Solução exata: y(x) = x +3.e-x -1
(0 2) 1 6562• y1(4) = -y1’’’ ≈ 2,4562 
ƒ Portanto, y2 ≈ 1,40995
y(0,2) ≈ 1,6562
y(0,4) ≈ 1,4110
Método de Euler AperfeiçoadoMétodo de Euler Aperfeiçoadop f çp f ç
ƒ Vejamos agora o Método de Euler Aperfeiçoado (também 
h d d Mét d d H )
L2
chamado de Método de Heun):
ƒ A reta L1, com coeficiente angular 
y’i = f(xi,yi), une os pontos Q = (xi,yi) 
L0
ŷi+1 = yi + hy’i
y L1
Py y p y
e P = (xi+1, ŷi+1):
ƒ L1: y = yi + (x-xi).f(xi,yi)
ƒ Por P traça-se a reta L2 com 
L
yi+1 y (x)
yi
ŷi+1 yi y i
QPor P, traça se a reta L2 com 
coeficiente angular f(xi+1, ŷi+1):
ƒ L2: y = ŷi+1 + (x-xi+1).f(xi+1,ŷi+1)
P r P tr ç s r t L c m inclin çã médi xxi 
yi
xi 1 
Q
ƒ Por P, traça-se a reta L0 com inclinação média 
entre L1 e L2: [f(xi,yi) + f(xi+1,ŷi+1)]/2
ƒ Por Q, traça-se a reta L paralela a L0:
É l ó 
xxi xi+1 
h
ƒ L: y = yi + (x-xi).[f(xi,yi) + f(xi+1,ŷi+1)]/2
ƒ A partir de L e de xi+1, obtém-se o valor de yi+1:
ƒ yi 1 = yi + (xi 1-xi) [f(xi yi) + f(xi 1 ŷi 1)]/2
É passo simples: só usa yi
Só calcula f(x,y)
yi+1 yi + (xi+1 xi).[f(xi,yi) + f(xi+1,ŷi+1)]/2
ƒ yi+1 = yi + h[f(xi,yi) + f(xi+1,yi + h.y’i))]/2
ƒ yi+1 = yi + h[f(xi,yi) + f(xi + h,yi + h.f(xi,yi))]/2
Coincide com um 
Método de Runge-
Kutta de 2ª ordem
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Métodos de RungeMétodos de Runge--KuttaKuttagg
ƒ A ideia básica destes métodos é aproveitar as qualidades 
d é d d é i d T l dos métodos de série de Taylor e, ao mesmo tempo, 
eliminar sua maior dificuldade de implementação: o cálculo 
das derivadas de f(x y)das derivadas de f(x,y)
ƒ Características dos Métodos de Runge-Kutta de ordem n :
1) São métodos de passo simples1) São métodos de passo simples
2) Não exigem o cálculo de qualquer derivada de f(x,y); por esse 
motivo, calculam f(x,y) em vários pontos
3) Após expandir f(x,y) por Taylor para função de duas variáveis em 
torno de (xi,yi) e agrupar os termos semelhantes, sua expressão 
coincide com a do método de série de Taylor de ordem nm m y m
ƒ O Método de Euler (equivalente ao método de série de 
Taylor de 1ª ordem) é um Método de Runge-Kutta de 1ª y g
ordem, e o Método de Euler Aperfeiçoado é um Método de 
Runge-Kutta de 2ª ordem
RungeRunge--Kutta de ordem nKutta de ordem ngg mm
ƒ Fórmula geral dos Métodos de Runge-Kutta:
 Φ( h)hƒ yi+1 = yi + Φ(xi, yi, h)h
ƒ Φ(xi, yi, h) é chamada função incremento, e pode ser 
interpretada como a inclinação no intervalo consideradop ç
ƒ Fórmula geral da função incremento de ordem n : 
ƒ Φ(xi, yi, h) = a1k1 + a2k2 + ... + ankn
k f( )ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + p1h, yi + q11k1h)
ƒ k3 = f(xi + p2h, yi + q21k1h + q22k2h)3 ( i p2 , yi q21 1 q22 2 )
ƒ ...
ƒ kn = f(xi + pn-1h, yi + q(n-1)1k1h + ... + q(n-1)(n-1)kn-1h)
ƒ a p e q : constantes obtidas igualando se a fórmula geral de ƒ ai, pij e qij: constantes obtidas igualando-se a fórmula geral de 
Runge-Kutta com os termos da expansão em série de Taylor
ƒ ki: relações de recorrência (cálculo computacional eficiente)ki relações de recorrênc a (cálculo computac onal ef c ente)
ƒ Os termos desprezados são de ordem O(hn+1), o que acarreta 
um erro global de ordem O(hn), pois h<1 Veremos isso depois...
RungeRunge--Kutta de 2ª ordemKutta de 2ª ordemgg mm
ƒ Série de Taylor de 2ª ordem em torno de xi:
ƒ y = y + f(x y )h + f’(x y )h2/2! + O(h3)ƒ yi+1 = yi + f(xi,yi)h + f (xi,yi)h2/2! + O(h3)
ƒ yi+1 = yi + f(xi,yi)h + [fx(xi,yi) + fy(xi,yi)f(xi,yi)]h2/2 + O(h3)
ƒ yi+1 = yi + hf(xi,yi) + h2fx(xi,yi)/2 + h2fy(xi,yi)f(xi,yi)/2+ O(h3)
é ªƒ Consideremos a definição do Método de Runge-Kutta de 2ª ordem:
yi+1 = yi + (a1k1 + a2k2)h, onde k1 = f(xi,yi) e k2 = f(xi+p1h, yi+q11k1h)
ƒ Expandindo k2 por Taylor (função de 2 variáveis) em torno de (xi yi):Expandindo k2 por Taylor (função de 2 variáveis) em torno de (xi,yi):
ƒ f(xi+p1h, yi+q11k1h) = f(xi,yi) + p1hfx(xi,yi) + q11k1hfy(xi,yi) + O(h2)
ƒ Substituindo na fórmula de Runge-Kutta:
2ƒ yi+1 = yi + a1k1h + a2[f(xi,yi) + p1hfx(xi,yi) + q11k1hfy(xi,yi) + O(h2)]h
ƒ yi+1 = yi + a1hf(xi,yi) + a2hf(xi,yi) + a2p1h2fx(xi,yi) + a2q11h2fy(xi,yi)f(xi,yi) + O(h3)
ƒ Desprezando os termos de O(h3), para que ambas expressões de yi+1p ( ), p q p yi+1
sejam iguais, é preciso que:
ƒ a1 + a2 = 1
ƒ a2p1 = ½a2p1 = ½
ƒ a2q11 = ½
ƒ 3 equações e 4 incógnitas: há infinitas soluções
RungeRunge--Kutta de 2ª ordemKutta de 2ª ordemgg mm
ƒ Há três versões mais utilizadas: a2 = ½, a2 = 1 ou a2 = 2/3
½ ½ƒ Método de Heun (a2 = ½, a1 = ½, p1 = q11 = 1):
ƒ yi+1 = yi + (½k1 + ½k2)h
ƒ k = f(x y )ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + h, yi + k1h)
ƒ Método do Ponto Médio (a2 = 1, a1 = 0, p1 = q11 = ½):Método do Ponto Médio (a2 1, a1 0, p1 q11 ½):
ƒ yi+1 = yi + k2h
ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + ½h, yi + ½k1h)
ƒ Método de Ralston (a2 = 2/3, a1 = 1/3, p1 = q11 = 3/4):
 (k /3 2k /3)hƒ yi+1 = yi + (k1/3 + 2k2/3)h
ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + 3h/4, yi + 3k1h/4)k2 f(xi 3h/4, yi 3k1h/4)
ƒ Este método fornece um limitante mínimo para o erro de 
aproximação nos Métodos de Runge-Kutta de 2ª ordem
RungeRunge--Kutta de 3ª e 4ª ordensKutta de 3ª e 4ª ordensgg
ƒ De modo semelhante, podem ser deduzidas as fórmulas 
d R K tt d d ide Runge-Kutta de ordens superiores
ƒ Em cada ordem, também haverá infinitas versões
Mé d d R K i h idƒ Métodos de Runge-Kutta mais conhecidos:
ƒ 3ª ordem:
ƒ y = y + (k + 4k + k )h/6ƒ yi+1 = yi + (k1 + 4k2 + k3)h/6
ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + ½h, yi + ½k1h)2 i yi 1
ƒ k3 = f(xi + h, yi - k1h + 2k2h)
ƒ 4ª ordem:
 (k k k k )h/ƒ yi+1 = yi + (k1 + 2k2 + 2k3 + k4)h/6
ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + ½h yi + ½k1h)k2 f(xi ½h, yi ½k1h)
ƒ k3 = f(xi + ½h, yi + ½k2h)
ƒ k4 = f(xi + h, yi + k3h)
ExemploExemplompmp
ƒ Usando o Método de Runge-Kutta de 2ª ordem (Método de Heun), 
resolva y’ x y tal que y(0) 2resolva y = x – y, tal que y(0) = 2
ƒ Consideraremos h = 0,2
ƒ f(x,y) = x - y
ƒ x0 = 0, xi = x0 + 0,2i
ƒ y0 = 2
ƒ k1 = f(x y )ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + h, yi + k1h) 
ƒ yi+1 = yi + (½k1 + ½k2)h
i xi yi k1 k2
0 0,0 2 -2 -1,4
1 0 2 1 66 1 46 0 9681 0,2 1,66 -1,46 -0,968
2 0,4 1,4172 -1,0172 -0,61376
3 0,6 1,254104 -0,654104 -0,323283
4 0,8 1,1563652 -0,356369 -0,0850914
5 1,0 1,1122192
Erros e ordens superioresErros e ordens superiorespp
ƒ Há um conhecido Método de Runge-Kutta de 5ª ordem, chamado 
Método de Butcher:Método de Butcher:
ƒ yi+1 = yi + (7k1 + 32k3 + 12k4 + 32k5 + 7k6)h/90
ƒ k1 = f(xi, yi)
ƒ k2 = f(xi + h/4, yi + k1h/4)
ƒ k3 = f(xi + h/4, yi + k1h/8 + k2h/8)
ƒ k4 = f(xi + h/2, yi – k2h/2 + k3h)4 f( i , yi 2 3 )
ƒ k5 = f(xi + 3h/4, yi + 3k1h/16 + 9k4h/16)
ƒ k6 = f(xi + h, yi - 3k1h/7 + 2k2h/7 + 12k3h/7 - 12k4h/7 + 8k5h/7)
E id t t é ssí l bt fó l s d R K tt d ƒ Evidentemente, é possível obter fórmulas de Runge-Kutta de 
ordens superiores, mas, de modo geral, o ganho em precisão não 
compensa o esforço computacional exigido nesses cálculos
ƒ Além disso, vimos que o Método de Runge-Kutta de ordem n se 
baseia na série de Taylor de mesma ordem: despreza termos de 
O(hn+1) mas com multiplicadores de valor desconhecido (y(n+1)(ξ))O(h ), mas com multiplicadores de valor desconhecido (y (ξ))...
ƒ Por esses motivos, costuma-se utilizar passo adaptativo, de modo 
análogo ao cálculo da quadratura (vide Capítulo 7)
Resultados experimentaisResultados experimentaisp mp m
ƒ Dado um PVI com solução analítica conhecida, podemos 
ê é ª ªresolvê-lo com métodos de Runge-Kutta de 1ª a 5ª 
ordens, com diversos tamanhos do passo h
S s s s lt d s btid s s l 㠃 Se compararmos os resultados obtidos com a solução 
exata, teremos um gráfico semelhante ao abaixo:
é ú d h d d f ã 
100
Erro 
relativo (%)
ƒ nf é o número de chamadas da função 
f(x,y) em cada iteração do método
ƒ (b-a)/h: número total de iterações
10-2
1
Euler
ƒ nf(b-a)/h é uma estimativa do tempo total 
gasto na execução do método
ƒ Conclusões:
10-6
10-4
10 2 Heun
RK de 3ª ordem
RK de 4ª ordem
ƒ Métodos de ordem superior alcançam uma 
precisão maior com o mesmo esforço 
computacional
Total de chamadas = nf(b-a)/h
Butcher ƒ Depois de um certo passo h, sua 
diminuição passa a proporcionar um ganho 
muito pequeno na precisão 
MatLabMatLab
ƒ ode23(fun, tspan, y0)
ƒ Através do método de Runge-Kutta adaptativo no tamanho do passo, com fórmulas 
de 2ª e 3ª ordens, retorna uma matriz de 2 colunas com pares (x,y)
ƒ fun é a função f(x,y)y
ƒ tspan é um vetor que descreve o intervalo desejado
ƒ y0 é valor de y no início desse intervalo
E lƒ Exemplo:
ƒ Dada a equação y’ = x - y, calcular os pares (x,y) no intervalo [0,1], onde y(0) = 2
ƒ [x,y] = ode23(inline('x - y'), [0 1], [2])
ƒ O passo h é encontrado automaticamente para que cada yi tenha erro relativo 10-3
ƒ ode45(fun, tspan, y0)
ƒ Idem utilizando fórmulas de Runge Kutta de 4ª e 5ª ordensƒ Idem, utilizando fórmulas de Runge-Kutta de 4 e 5 ordens
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Métodos de passo múltiploMétodos de passo múltiplop m pp m p
ƒ Vimos que, para encontrar uma aproximação de y(xi+1), 
os métodos de passo simples precisam apenas de y(xi), 
além de cálculos de y’ = f(x,y) em outros pontos
ƒ Por outro lado, suponhamos que, além de y(xi), também 
sejam conhecidas aproximações y(xi-1), ..., y(xi-m) em m 
pontos equidistantes isto é x x h i≤j≤i m+1pontos equidistantes, isto é, xj – xj-1 = h, i≤j≤i-m+1
ƒ Os métodos que utilizam o valor de y em mais de um 
t ã h d ét d d últi lponto são chamados métodos de passo múltiplo
ƒ Esses métodos baseiam-se na percepção de que, uma 
 ál l t h d i f ã li vez que o cálculo tenha começado, informação valiosa 
já está à disposição: a curvatura formada pelos valores 
anteriores permite uma melhor aproximação da soluçãoanteriores permite uma melhor aproximação da solução
Métodos de AdamsMétodos de Adamsmm
ƒ Entre os métodos de passo múltiplo, há uma classe p p ,
conhecida como Métodos de Adams, que se baseiam na 
integração numérica de y’ = f(x,y) de xi até xi+1:
∫ ∫∫ + ++ +=⇔= +1i 1i1i
x x
i1i
x
dx))x(y,x(f)x(y)x(ydx))x(y,x(fdx)x('y ∫ ∫∫
i ii x xx
ƒ Por sua vez, isso pode ser feito através de dois modos:
ƒ Adams–Bashforth (métodos explícitos ou fórmulas abertas) : 
sem usar o ponto xi+1
ƒ Adams-Moulton (métodos implícitos ou fórmulas fechadas) : 
usando o ponto xusando o ponto xi+1
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Métodos explícitosMétodos explícitosppƒ Na aproximação dessa integral, os Métodos de Adams-
Bashfort utilizam m+1 pontos: xi, xi-1, ..., xi-m
ƒ Por isso, são chamados métodos de ordem m+1
 é f é d ã d l ô ƒ Isso é feito através da integração do polinômio 
interpolador pm(x):
1ix +
dx)x(p)x(y)x(y
1i
ix
mi1i ∫++≈+
ƒ Seja y’k = f(xk,yk) = fk, k ≥ 0
ƒ A função y’ = f(x,y) é aproximada pelo polinômio pm(x), 
que interpola os pontos (xi, fi), (xi-1, fi-1), ..., (xi-m, fi-m). 
Basta escolher o valor de m
P d ( ) é d f d Lƒ Podemos expressar pm(x) através da forma de Lagrange:
ƒ pm(x) = L-m(x)fi-m + ... + L-1(x)fi-1 + L0(x)fi
Ordem 4: caso com pOrdem 4: caso com p33(x)(x)m m pm m p33( )( )
ƒ Pontos de interpolação: (xi, fi), (xi-1, fi-1), (xi-2, fi-2), (xi-3, fi-3) 
ƒ f(x y(x)) = y’(x) ≈ p (x) = L (x)f + L (x)f + L (x)f + L (x)fƒ f(x,y(x)) = y (x) ≈ p3(x) = L-3(x)fi-3 + L-2(x)fi-2 + L-1(x)fi-1 + L0(x)fi
ƒ L-3(x) = [(x - xi-2)(x - xi-1)(x - xi)]/(-h)(-2h)(-3h)
ƒ L-2(x) = [(x - xi-3)(x - xi-1)(x - xi)]/(h)(-h)(-2h)
/ƒ L-1(x) = [(x - xi-3)(x - xi-2)(x - xi)]/(2h)(h)(-2h)
ƒ L0(x) = [(x - xi-3)(x - xi-2)(x - xi-1)]/(3h)(2h)(h)
ƒ Sejam s = (x - xi)/h, dx = h.ds e x = hs + xi. Então:j i i
ƒ L-3(s) = -(s + 2)(s + 1)s/6 = -(s3 + 3s2 + 2s)/6
ƒ L-2(s) = (s + 3)(s + 1)s/2 = (s3 + 4s2 + 3s)/2
ƒ L-1(s) = -(s + 3)(s + 2)s/2 = -(s3 + 5s2 + 6s)/2-1( ) ( )( ) ( )
ƒ L0(s) = (s + 3)(s + 2)(s + 1)/6 = (s3 + 6s2 + 11s + 6)/6
ƒ Substituindo na integral:
++ 1111xx hhhh1i1i ∫∫∫∫∫∫ +−+−=≈ −−−−−−++ 1
0
0i
1
0
11i
1
0
22i
1
0
33i
x
3
x
ds)s(Lf
6
hds)s(Lf
2
hds)s(Lf
2
hds)s(Lf
6
hdx)x(pdx))x(y,x(f
1i
i
1i
i
i1i2i3i
x
3 f24
h55f
24
h59f
24
h37f
24
h9dx)x(p
1i
+−+−= −−−∫+ i1i2i3i
x
3 24242424
)(p
i
∫
]f9f37f59f55[
24
hyy 3i2i1iii1i −−−+ −+−+=
Ordem 4: estimativa de erroOrdem 4: estimativa de errom mm m
ƒ Pontos de interpolação: (xi, fi), (xi-1, fi-1), (xi-2, fi-2), (xi-3, fi-3) 
ƒ Vimos anteriormente que o erro na interpolação com p (x) é ƒ Vimos anteriormente que o erro na interpolação com p3(x) é 
E3(x) = (x – xi-3)(x – xi-2)(x – xi-1)(x – xi)f(4)(ξ)/4!, onde ξ ∈ (xi,xi-3)
ƒ Portanto, o erro cometido é:
dx))ξ(y,ξ(f)xx)(xx)(xx()xx(
!4
1)x(e )4(i1i2i
x
x
3i1i
1i
i
−−−−= −−−+ ∫+
ƒ Com s = (x - xi)/h, dx = h.ds e x = hs + xi:
ds))ξ(yξ(sf)1s)(2s()3s(h)x(e )4(
15
1i +++= ∫ ds))ξ(y,ξ(sf)1s)(2s()3s(!4)x(e 01i +++∫+
ƒ Como g(s) = s(s+1)(s+2)(s+3) não muda de sinal em [0;1], o Teorema 
d V l Médi i t i t i t (0 1) t l do Valor Médio para integrais garante que existe η ∈ (0;1) tal que:
∫∫ ==+++ 1
0
)4(
5
)4(
5
)4(
1
0
5
30
251))η(y,η(f
24
hds)s(g))η(y,η(f
!4
hds))ξ(y,ξ(sf)1s)(2s()3s(
!4
h
00
ƒ Portanto: 720
251)η(yh
720
251))η(y,η(fh)x(e )5(5)4(51i ==+
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Métodos implícitosMétodos implícitosmpmp
ƒ Na aproximação da integral, os Métodos de Adams-p g
Moulton utilizam m+2 pontos: xi+1, xi, ..., xi-m
ƒ Neste caso, o método tem ordem m+2, e a integração é 
f i é d ( )feita através de pm+1(x):
dx)x(p)x(y)x(y
1ix∫++ dx)x(p)x(y)x(y
ix
1mi1i ∫ ++ +≈
ƒ O polinômio p (x) interpola os pontos (x f ) ƒ O polinômio pm+1(x) interpola os pontos (xi+1, fi+1), 
(xi, fi), ..., (xi-m, fi-m)
ƒ De modo análogo aos métodos explícitos basta escolher ƒ De modo análogo aos métodos explícitos, basta escolher 
o valor de m e calcular a integração da forma de 
Lagrange:g g
ƒ pm+1(x) = L-m(x)fi-m + ... + L-1(x)fi-1 + L0(x)fi + L1(x)fi+1
Ordem 4: caso com pOrdem 4: caso com p33(x)(x)m m pm m p33( )( )
ƒ Pontos de interpolação: (xi+1, fi+1), (xi, fi), (xi-1, fi-1), (xi-2, fi-2) 
ƒ f(x y(x)) = y’(x) ≈ p (x) = L (x)f + L (x)f + L (x)f + L (x)fƒ f(x,y(x)) = y (x) ≈ p3(x) = L-2(x)fi-2 + L-1(x)fi-1 + L0(x)fi + L1(x)fi+1
ƒ L-2(x) = [(x - xi-1)(x - xi)(x - xi+1)]/(-3h)(-2h)(-h)
ƒ L-1(x) = [(x - xi-2)(x - xi)(x - xi+1)]/(h)(-h)(-2h)
/ƒ L0(x) = [(x - xi-2)(x - xi-1)(x - xi+1)]/(2h)(h)(-h)
ƒ L1(x) = [(x - xi-2)(x - xi-1)(x - xi)]/(3h)(2h)(h)
ƒ Sejam s = (x - xi)/h, dx = h.ds e x = hs + xi. Então:j i i
ƒ L-2(s) = -(s + 1)s(s - 1)/6 = -(s3 - s)/6
ƒ L-1(s) = (s + 2)s(s - 1)/2 = (s3 + s2 - 2s)/2
ƒ L0(s) = -(s + 2)(s + 1)(s - 1)/2 = -(s3 + 2s2 – s - 2)/20( ) ( )( )( ) ( )
ƒ L1(s) = (s + 2)(s + 1)s/6 = (s3 + 3s2 + 2s)/6
ƒ Substituindo na integral:
1111xx 1i1i ∫∫∫∫∫∫ +−−−− +−+−=≈ ++ 1
0
11i
1
0
0i
1
0
11i
1
0
22i
x
x
3
x
x
ds)s(Lf
6
hds)s(Lf
2
hds)s(Lf
2
hds)s(Lf
6
hdx)x(pdx))x(y,x(f
1i
i
1i
i
]ff5f19f9[
24
hyy 2i1ii1ii1i −−++ +−++= 24
yi+1 está presente em fi+1 = f(xi+1,yi+1): formulação implícita
Ordem 4: estimativa de erroOrdem 4: estimativa de errom mm m
ƒ Pontos de interpolação: (xi+1, fi+1), (xi, fi), (xi-1, fi-1), 
(xi-2, fi-2) 
ƒ De forma análoga, com s = (x - xi)/h, dx = h.ds e g ( i)
x = hs + xi:
d))ξ(ξ(f)1()1()2(h)( )4(
15 ∫ ds))ξ(y,ξ(f)1s(s)1s()2s(!4h)x(e )4(01i −++= ∫+
ƒ Como g(s) = (s+2)(s+1)s(s-1) é sempre menor ou igual a 
zero em [0;1], então existe η ∈ (0;1) tal que:
720
19)(yh)x(e )5(51i η−=+ 720
ExemploExemplompmp
ƒ Seja o PVI y’ = 0,04y, onde y(0) = 1000 
ƒ Através do Método de Adams-Bashforth de ordem 4, vamos 
calcular uma aproximação de y(1) usando h = 0,2
ƒ x0 = 0 e y0 = 1000x0 0 e y0 000
ƒ Apenas para testar o método, vamos utilizar y1, y2 e y3 obtidos através 
da solução exata do PVI: y(x) = 1000.e0,04x
ƒ Em seguida utilizamos a sua fórmula de cálculo:ƒ Em seguida, utilizamos a sua fórmula de cálculo:
ƒ yi+1 = yi + h(55fi - 59fi-1 + 37fi-2 – 9fi-3)/24
i f f( ) l ã ti xi yi fi = f(xi,yi) solução exata
0 0,0 1000 40 1000
1 0,2 1008,0321 40,321284 1008,0321, , , ,
2 0,4 1016,1287 40,645148 1016,1287
3 0,6 1024,2903 40,971612 1024,2903
4 0 8 1032 517487 41 30069948 1032 51754 0,8 1032,517487 41,30069948 1032,5175
5 1,0 1040,810756 1040,810774
Alguns casosAlguns casosgg
ƒ Métodos explícitos (Adams-Bashforth):p ( )
Ordem Fórmula Erro
2 yi+1 = yi + h(3fi – fi-1)/2 5h3f’’(ξ)/12
3 yi+1 = yi + h(23fi – 16fi-1 + 5fi-2)/12 9h4f(3)(ξ)/24
4 yi+1 = yi + h(55fi – 59fi-1 + 37fi-2 – 9fi-3)/24 251h5f(4)(ξ)/720
 h( f f f f f )/ h6f(5)(ξ)/5 yi+1 = yi + h(1901fi – 2774fi-1 + 2616fi-2 – 1274fi-3 + 251fi-4)/720 475h6f(5)(ξ)/1440
ƒ Métodos implícitos (Adams-Moulton):
Ordem Fórmula Erro
2 yi+1 = yi + h(fi+1 + fi)/2 -h3f’’(ξ)/12
3 yi+1 = yi + h(5fi+1 + 8fi - fi-1)/12 -h4f(3)(ξ)/24
4 yi+1 = yi + h(9fi+1 + 19fi - 5fi 1 + fi 2)/24 -19h5f(4)(ξ)/7204 yi+1 yi h(9fi+1 19fi 5fi-1 fi-2)/24 19h f (ξ)/720
5 yi+1 = yi + h(251fi+1 + 646fi - 264fi-1 + 106fi-2 - 19fi-3)/720 -27h6f(5)(ξ)/1440
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Métodos de previsãoMétodos de previsão--correçãocorreçãopp çç
ƒ Uma das principais desvantagens dos métodos de passo 
últi l é ã t i i i i d d d múltiplo é que não se auto-iniciam: precisam de dadosque 
geralmente são fornecidos por algum método de passo 
simples (Runge-Kutta ou série de Taylor por exemplo)simples (Runge Kutta ou série de Taylor, por exemplo)
ƒ Por outro lado, parece difícil utilizar métodos implícitos, 
pois na expressão de yi+1 aparece fi+1...p p yi+1 p fi+1
ƒ Na verdade, eles são usados em pares previsor-corretor :
1) Através de um método explícito (chamado previsor), encontra-
se a primeira aproximação y0i+1 para yi+1
2) Calcula-se então fi+1 = f(xi+1, y0i+1)
3) Com um método implícito (chamado corretor) utiliza-se o valor 3) Com um método implícito (chamado corretor), utiliza-se o valor 
acima para calcular uma nova aproximação y1i+1 para yi+1
4) Volta-se ao passo 2, e o processo continua até que um 
d i d l i d j l ddeterminado erro relativo de yi+1 seja alcançado
5) Caso se deseje calcular yi+2, calcula-se fi+1 e volta-se ao passo 1
ExemploExemplompmp
ƒ Seja o PVI y’ = -y2, onde y(1) = 1. Deseja-se obter valores de y com 
 l ti 10 4erros relativos menores que 10-4
ƒ Consideremos, por exemplo, h = 0,1
ƒ Novamente, vamos utilizar a solução exata y(x) = 1/x para calcular y1, y p y1
y2 e y3, pois usaremos métodos de ordem 4:
x0 = 1 y0 = 1/x0 = 1 f0 = -y02 = -1
x 1 1 y 1/x 0 9090909 f y 2 0 8264462x1 = 1,1 y1 = 1/x1 = 0,9090909 f1 = -y12 = -0,8264462
x2 = 1,2 y2 = 1/x2 = 0,8333333 f2 = -y22 = -0,6944443
x3 = 1,3 y3 = 1/x3 = 0,7692307 f3 = -y32 = -0,5917158
ƒ Previsor: y04 = y3 + h(55f3 – 59f2 + 37f1 – 9f0)/24 = 0,7144362
ƒ f04 = f(x4,y04) = -(y04)2 = -0,510419
ƒ Corretor: y1 = y + h(9f0 + 19f 5f + f )/24 = 0 7142698ƒ Corretor: y14 = y3 + h(9f04 + 19f3 - 5f2 + f1)/24 = 0,7142698
ƒ f14 = f(x4,y14) = -(y14)2 = -0,5101814
ƒ Corretor: y24 = y3 + h(9f14 + 19f3 - 5f2 + f1)/24 = 0,7142787
ƒ |y24 – y14|/|y24| = 1,2591374.10-5 < 10-4
ƒ Calcular f24, usar o previsor no cálculo de y05, e continuar o processo...
Convergência Convergência gg
ƒ Questões sobre os métodos de previsão-correção:Q p ç
ƒ Em que condições há garantia de convergência para yi+1?
ƒ Quantas iterações do corretor são necessárias para se atingir Q ç p g
essa convergência na precisão desejada?
ƒ Teorema: Se f(x,y) e ∂f/∂y são contínuas em x e y em 
todo o intervalo [a,b], as iterações do corretor vão 
convergir desde que h.|∂f/∂y| < 2, para x = xi e todo y 
tal que |y y | ≤ |y0 y |tal que |y – yi+1| ≤ |y0i+1 – yi+1|
ƒ Na prática, basta escolher h suficientemente pequeno...
ƒ Além disso, resultados experimentais mostram que, se o 
par previsor-corretor for da mesma ordem e h 
satisfizer as condições do teorema bastam apenas uma satisfizer as condições do teorema, bastam apenas uma 
ou duas iterações do corretor
Voltando ao exemplo anteriorVoltando ao exemplo anteriormpmp
ƒ Seja o PVI y’ = -y2, onde y(1) = 1
ƒ ∂f/∂y = -2y
P t d ê i j ti f it ƒ Para que o teorema da convergência seja satisfeito, 
h.|2y| < 2, ou seja, h < 1/|y| garante a convergência
ƒ No exemplo anterior, todos os valores obtidos para y 
são menores que 1, ou seja, 1/|y| > 1
ƒ Portanto, o espaçamento h = 0,1 satisfaz a condição 
exigida para a convergência
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Sistema de equações de 1ª ordemSistema de equações de 1ª ordemm q ç mm q ç m
ƒ Consideremos inicialmente um sistema de n EDOs de 
i i dprimeira ordem:
ƒ Este sistema pode ser escrito de forma vetorial:
y’ = f(x, y), onde y’, f e y têm componentes y’i, fi e yi, 
respectivamente, com 1 ≤ i ≤ n
ƒ Como no caso das equações diferenciais de primeira 
ordem, para que este sistema possua solução única, 
ã á i di õ i i i i d serão necessárias n condições iniciais, que podem ser 
expressas por yi(x0) = yi0
Exemplo: 2 equações de 1ª ordemExemplo: 2 equações de 1ª ordemmp q ç mmp q ç m
ƒ Resolva o sistema abaixo, onde n = 2:
’ ƒ y’ = z 
ƒ z’ = y + ex
ƒ y(0) = 1 z(0) = 0 x ∈ [0; 0 2] h = 0 1y(0) = 1, z(0) = 0, x ∈ [0; 0,2], h = 0,1
ƒ Consideraremos:
ƒ y’ = f(x, y, z) = z y f(x, y, z) z 
ƒ z’ = g(x, y, z) = y + ex
ƒ x0 = 0, y0 = 1, z0 = 0
ƒ Podemos utilizar, por exemplo, o Método de Heun:
⎥⎤⎢⎡
+⎤⎡⎤⎡ + 2f1fi1i kkhyy R t ⎥⎦⎢⎣ +
+⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+
+
2g1g
2f1f
i
i
1i
1i
kk2
h
z
y
z
y
⎤⎡⎤⎡ )zyx(fk ⎤⎡⎤⎡ )hkhkhx(fk
Repete-se nas 
demais variáveis
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
)z,y,x(g
)z,y,x(f
k
k
iii
iii
1g
1f ⎥⎦
⎤⎢⎣
⎡
+++
+++=⎥⎦
⎤⎢⎣
⎡
)hkz,hky,hx(g
)hkz,hky,hx(f
k
k
1gi1fii
1gi1fii
2g
2f
Continuação do exemploContinuação do exemploç mpç mp
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
2
0
)010(g
)0,1,0(f
)zyx(g
)z,y,x(f
k
k 000
1
1f
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+++
+++=⎥⎦
⎤⎢⎣
⎡
1052,2
2,0
)2,0;1;1,0(g
)2,0;1;1,0(f
)hkz,hky,hx(g
)hkz,hky,hx(f
k
k
1g01f00
1g01f00
2g
2f
⎦⎣⎦⎣⎦⎣⎦⎣ 2)0,1,0(g)z,y,x(gk 0001g
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+
++⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+
++⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
2053,0
01,1
1052,22
2,00
2
1,0
0
1
kk
kk
2
h
z
y
z
y
2g1g
2f1f
0
0
1
1
⎦⎣⎦⎣
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
1152,2
2053,0
)2053,0;01,1;1,0(g
)2053,0;01,1;1,0(f
)z,y,x(g
)z,y,x(f
k
k
111
111
1g
1f
⎤⎡⎤⎡ +⎤⎡⎤⎡ +⎤⎡⎤⎡ 04111416802053010011kkhyy
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+++
+++=⎥⎦
⎤⎢⎣
⎡
2519,2
4168,0
)4168,0;0305,1;2,0(g
)4168,0;0305,1;2,0(f
)hkz,hky,hx(g
)hkz,hky,hx(f
k
k
1g11f11
1g11f11
2g
2f
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+
++⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+
++⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
4237,0
0411,1
2519,21152,2
4168,02053,0
2
1,0
2053,0
01,1
kk
kk
2
h
z
y
z
y
2g1g
2f1f
1
1
2
2
xi yi zi
0 1 0
0,1 1,01 0,2053
0,2 1,0411 0,4237
Resultados:
Equações de ordem superiorEquações de ordem superiorq ç m pq ç m p
ƒ Uma EDO de ordem m > 1, y(m) = f(x, y, y’, y’’, ..., y(m-1)), 
pode ser facilmente transformada em um sistema com 
m EDOs de primeira ordem. Considerando y = z1, temos 
as seguintes EDOs:as seguintes EDOs:
ƒ y’ = z1’ = z2
ƒ y’’ = z2’ = z3
ƒ ...
ƒ y(m-1) = zm-1’ = zm
ƒ y(m) = zm’ = f(x, y, y’, y’’, ..., y(m-1)) = f(x, z1, z2, z3, ..., zm)
ƒ Este sistema pode ser resolvido de modo análogo ao p g
caso anterior
Exemplo: EDO 2ª ordem, passo simplesExemplo: EDO 2ª ordem, passo simplesmp m, p mpmp m, p mp
ƒ Dado o PVI y’’ = ex - 2y2, com y(0) = 0 e y’(0) = 0, calcule y(0,5) e 
y’’(0 5) usando o Método de Runge-Kutta de 4ª ordem com h = 0 5y (0,5) usando o Método de Runge Kutta de 4 ordem com h = 0,5
ƒ Sejam: y’ = f(x, y, z) = z, y’’ = z’ = g(x, y, z) = ex – 2y2
ƒ Dados iniciais: x0 = 0, y0 = 0, z0 = 0
ê d ál lƒ Sequência de cálculos:
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
1
0
)0,0,0(g
)0,0,0(f
)z,y,x(g
)z,y,x(f
k
k
000
000
1g
1f
⎦⎣⎦⎣⎦⎣⎦⎣ gyg 000g
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+++
+++=⎥⎦
⎤⎢⎣
⎡
2840,1
25,0
)25,0;0;25,0(g
)25,0;0;25,0(f
)2/hkz,2/hky,2/hx(g
)2/hkz,2/hky,2/hx(f
k
k
1g01f00
1g01f00
2g
2f
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+++
+++=⎥⎦
⎤⎢⎣
⎡
2762,1
321,0
)321,0;0625,0;25,0(g
)321,0;0625,0;25,0(f
)2/hkz,2/hky,2/hx(g
)2/hkz,2/hky,2/hx(f
k
k
2g02f00
2g02f00
3g
3f
⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
+++
+++=⎥⎦
⎤⎢⎣
⎡
5972,1
6381,0)6381,0;1605,0;5,0(g
)6381,0;1605,0;5,0(f
)hkz,hky,hx(g
)hkz,hky,hx(f
k
k
3g03f00
3g03f00
4g
4f
⎤⎡⎤⎡ +++⎤⎡⎤⎡ 1483,0kk2k2khyy 4f3f2f1f01 y(0,5) ≈ 0,1483⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
++++⎥⎦
⎤⎢⎣
⎡=⎥⎦
⎤⎢⎣
⎡
6431,0
483,0
kk2k2k6
h
z
y
z
y
4g3g2g1g
4f3f2f1f
0
0
1
1
y’’(0,5) = e0,5 - 2y(0,5)2 ≈ 1,6047
y’(0,5) ≈ 0,6431
y( , ) ,⇒
Exemplo: EDO 2ª ordem, passo múltiploExemplo: EDO 2ª ordem, passo múltiplomp m, p m pmp m, p m p
ƒ Agora, já seria possível continuar o exemplo anterior 
com o par previsor corretor de 2ª ordemcom o par previsor-corretor de 2 ordem
ƒ Lembrando:
ƒ h = 0,5, x0 = 0, y0 = 0, z0 = 0, x1 = 0,5, y1 = 0,1483, z1 = 0,6431
ƒ y’ = f(x, y, z) = z
ƒ y’’ = z’ = g(x, y, z) = ex – 2y2
ƒ Considerando fi = f(xi, yi, zi) = zi:n ran fi f( i, yi, zi) zi
ƒ Previsor de y: yi+1 = yi + h(3fi – fi-1)/2
ƒ Corretor de y: yi+1 = yi + h(fi+1 + fi)/2
ƒ Considerando g = g(x y z ) = exi – 2y 2:ƒ Considerando gi = g(xi, yi, zi) = exi – 2yi2:
ƒ Previsor de z: zi+1 = zi + h(3gi – gi-1)/2
ƒ Corretor de z: zi+1 = zi + h(gi+1 + gi)/2 
S ê i d ál l d b d idƒ Sequência de cálculos que deve ser obedecida:
ƒ Previsor de y: y2 = y1 + h(3f1 – f0)/2 = 0,6306
ƒ Previsor de z: z2 = z1 + h(3g1 – g0)/2 = 1,5966
ƒ Corretor de y: y2 = y1 + h(f2 + f1)/2 = 0,7082 (pode ser iterado)
ƒ Corretor de z: z2 = z1 + h(g2 + g1)/2 = 1,5250 (pode ser iterado)
ƒ E assim por diante... 
Exemplo: EDO 3ª ordemExemplo: EDO 3ª ordemmp mmp m
ƒ Dado o PVI y’’’ = 1 – y’, onde y(0) = 0, y’(0) = 2 e y’’(0) =0, y y y y y
calcule y(1), y’(1) e y’’(1) com o Método de Heun
ƒ Sejam:
’ f( ) ƒ y’ = f(x, y, z, w) = z
ƒ y’’ = z’ = g(x, y, z, w) = w
ƒ y’’’ = w’ = t(x, y, z, w) = 1 – z
ƒ Dados iniciais: x0 = 0, y0 = 0, z0 = 2, w0 = 0
ƒ Fórmulas do Método de Heun (Runge-Kutta de 2ª ordem):
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
=
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
)w,z,y,x(t
)w,z,y,x(g
)w,z,y,x(f
k
k
k
iiii
iiii
iiii
1t
1g
1f
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
++++
++++
++++
=
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
)hkw,hkz,hky,hx(t
)hkw,hkz,hky,hx(g
)hkw,hkz,hky,hx(f
k
k
k
1ti1gi1fii
1ti1gi1fii
1ti1gi1fii
2t
2g
2f
⎥⎦⎢⎣⎥⎦⎢⎣ ),,y,( iiii1t ⎦⎣⎦⎣ ),,y,( 1ti1gi1fii2t
⎥⎤⎢⎡ +
+
+⎥
⎤⎢⎡⎥⎤⎢⎡ + 2f1fi1i kk
kk
hz
y
z
y
⎥⎥
⎥
⎦⎢
⎢⎢
⎣ +
++⎥⎥⎦⎢
⎢
⎣
=⎥⎥⎦⎢
⎢
⎣ +
+
2t1t
2g1g
i
i
1i
1i
kk
kk
2w
z
w
z
Continuação do exemploContinuação do exemploç mpç mp
ƒ Utilizando h = 0,5:
i xi yi zi wi kf1 kg1 kt1 kf2 kg2 kt2
0 0 0 2 0 2 0 -1 2 -0 5 -10 0 0 2 0 2 0 1 2 0,5 1
1 0,5 1 1,875 -0,5 1,875 -0,5 -0,875 1,625 -0,938 -0,625
2 1 1,875 1,5156 -0,875
S l ã t ( ) ƒ Solução exata: y(x) = x + sen x
ƒ y(1) ≈ 1,8415
’(1) (1) 1 5403ƒ y’(1) = z(1) ≈ 1,5403
ƒ y’’(1) = w(1) ≈ -0,8415
E d é (h2) 0 25 d d é ƒ Estimativa de erro é O(h2) ≈ 0,25: na verdade, é menor...
CCICCI--2222
ƒ Definições
ƒ Problemas de Valor Inicial (PVI)
ƒ Métodos de passo simplesMétodos de passo simples
ƒ Método de Euler
ƒ Métodos de série de TaylorMétodos de série de Taylor
ƒ Métodos de Runge-Kutta
ƒ Métodos de passo múltiploMétodos de passo múltiplo
ƒ Métodos explícitos (Adams-Bashforth)
ƒ Métodos implícitos (Adams-Moulton)M t mp c t ( am M u t n)
ƒ Métodos de previsão-correção
ƒ Equações de ordem superiorEquações de ordem super or
ƒ Problemas de Valor de Contorno (PVC)
Problemas de Valor de ContornoProblemas de Valor de Contornomm
ƒ Como vimos anteriormente, dada uma Equação Diferencial Ordinária 
de ordem m > 1 se as condições iniciais envolvendo a função e suas de ordem m > 1, se as condições iniciais, envolvendo a função e suas 
derivadas até a ordem m-1, não são especificadas em um mesmo 
ponto, então temos um Problema de Valor de Contorno (PVC)
ª éƒ Concretamente, a forma geral dos PVC de 2ª ordem é:
ƒ y’’ = f(x, y, y’)
ƒ a1y(w) + b1y’(w) = c1 a1, a2, b1, b2, c1 e c2: constantes reais conhecidas b ã d l i lt ta1y(w) b1y (w) c1ƒ a2y(z) + b2y’(z) = c2
ƒ Se f(x,y,y’) = 0 e c1 = c2 = 0, o PVC é homogêneo: tem solução y(x) = 0
V l ã d PVC d 2ª d é d Mé d d 
ai e bi não podem ser nulos simultaneamente
ƒ Veremos a resolução de um PVC de 2ª ordem através do Método das 
Diferenças Finitas :
ƒ As derivadas são aproximadas por diferenças finitasp p ç
ƒ A equação diferencial transforma-se em um sistema de equações 
algébricas
ƒ Se esse sistema for linear pode ser resolvido com os métodos Se esse sistema for linear, pode ser resolvido com os métodos 
estudados no Capítulo 3; caso contrário, utilizam-se os métodos do 
Capítulo 4
Aproximações das derivadasAproximações das derivadasp m çp m ç
ƒ Considere o intervalo [a,b] dividido em n partes iguais 
d h h d b ã ê de tamanho h, onde x0 = a e xn = b. São três as 
aproximações mais usadas para y’(xi):
Diferença avançada
y’(xi) ≈ (yi+1 - yi)/hyi+1
Diferença centrada
y’(xi) ≈ (yi+1 – yi-1)/2h
yi-1
yi
Diferença atrasada
y’(xi) ≈ (yi – yi-1)/h
xxi-1 xi+1xi
hh
ƒ Nessas aproximações, podemos estimar os erros 
cometidos através da fórmula de Taylor de ordem k de 
( ) t d d ξ tá t y(x) em torno de xi, onde ξ está entre x e xi:
ƒ y(x) = y(xi) + y’(xi).(x-xi) + ... + y(k)(xi).(x-xi)k/k! + y(k+1)(ξ).(x-xi)k+1/(k+1)! 
Estimativa do erroEstimativa do erromm
ƒ O erro cometido no cálculo de y’(xi) através da 
diferença avançada pode ser estimado com a fórmula 
de Taylor de y(x) em torno de xi, considerando k = 1:
( ) ( ) ’( ) ( ) ’’(ξ) ( )2/2ƒ y(x) = y(xi) + y’(xi).(x - xi) + y’’(ξ).(x - xi)2/2
ƒ No ponto xi+1 = xi + h:
( ) ( ) ’( ) ( ) ’’(ξ ) ( )2/ƒ y(xi+1) = y(xi) + y’(xi).(xi+1 - xi) + y’’(ξi+1).(xi+1 - xi)2/2
ƒ y(xi+1) = y(xi) + y’(xi).h + y’’(ξi+1).h2/2
’( ) [ ( ) ( )]/h ’’(ξ ) h/2ƒ y’(xi) = [y(xi+1) - y(xi)]/h - y’’(ξi+1).h/2
ƒ Se y’’(x) for limitada em [a,b], então:
’( ) ( )/h O(h)ƒ y’(xi) = (yi+1 - yi)/h + O(h)
ƒ Um resultado análogo pode ser obtido em relação à 
diferença atrasada :diferença atrasada :
ƒ y’(xi) = (yi – yi-1)/h + O(h)
Estimativa do erroEstimativa do erromm
ƒ O erro cometido no cálculo de y’(xi) através da 
dif d d i d fó l diferença centrada pode ser estimado com a fórmula 
de Taylor de y(x) em torno de xi, considerando k = 2:
( ) ( ) ’( ) ( ) ’’( ) ( )2/2 ’’’(ξ) ( )3/6ƒ y(x) = y(xi) + y’(xi).(x - xi) + y’’(xi).(x - xi)2/2 + y’’’(ξ).(x - xi)3/6
ƒ Nos pontos xi+1 e xi-1:
( ) ( ) ’( ) h ’’( ) h2/2 ’’’(ξ ) h3/6ƒ y(xi+1) = y(xi) + y’(xi).h + y’’(xi).h2/2 + y’’’(ξi+1).h3/6
ƒ y(xi-1) = y(xi) - y’(xi).h + y’’(xi).h2/2 - y’’’(ξi-1).h3/6
S bt i d s õ s:ƒ Subtraindo as equações:
ƒ y(xi+1) - y(xi-1) = 2y’(xi).h + [y’’’(ξi+1) + y’’’(ξi-1)].h3/6
’(x ) [ (x ) (x )]/2h [ ’’’(ξ ) ’’’(ξ )] h2/12ƒ y (xi) = [y(xi+1) - y(xi-1)]/2h - [y (ξi+1) + y (ξi-1)].h2/12
ƒ Se y’’’(x) for limitada em [a,b], então:
’( ) ( )/2h O(h2)ƒ y’(xi) = (yi+1 – yi-1)/2h + O(h2)
ƒ Como geralmente h<1, esta fórmula é mais precisa
Aproximação da derivada segundaAproximação da derivada segundap m ç gp m ç g
ƒ Com a fórmula de Taylor de y(x) em torno de xi, agora 
 k 3 é í l i id ál l com k = 3, é possível estimar o erro cometido no cálculo 
de y’’(xi)
ƒ Nos pontos xi+1 e xi-1:
ƒ y(xi+1) = y(xi) + y’(xi).h + y’’(xi).h2/2! + y’’’(xi).h3/3! + y(4)(ξi+1).h4/4!y( i+1) y( i) y ( i) y ( i) y ( i) y (ξi+1)
ƒ y(xi-1) = y(xi) - y’(xi).h + y’’(xi).h2/2! - y’’’(xi).h3/3! + y(4)(ξi-1).h4/4!
S d s õ sƒ Somando as equações:
ƒ y(xi+1) + y(xi-1) = 2y(xi) + y’’(xi).h2 + [y(4)(ξi+1) + y(4)(ξi-1)].h4/24
ƒ y’’(xi) = [y(xi+1) – 2y(xi) + y(xi-1)]/h2 - [y(4)(ξi+1) + y(4)(ξi-1)].h2/24
ƒ Se y(4)(x) for limitada em [a b] então:Se y (x) for limitada em [a,b], então:
ƒ y’’(xi) = (yi+1 – 2yi + yi-1)/h2 + O(h2)
Exemplo: PVC linear de 2ª ordemExemplo:PVC linear de 2ª ordemmp mmp m
ƒ y’’(x) + 2y’(x) + y(x) = x, onde y(0) = 0 e y(1) = -1
Us m s s xim õ s m O(h2):ƒ Usaremos as aproximações com erro O(h2):
ƒ y’(xi) ≈ (yi+1 – yi-1)/2h
ƒ y’’(xi) ≈ (yi+1 – 2yi + yi-1)/h2
ƒ Substituindo-as na equação, e considerando x = xi:
ƒ (yi+1 – 2yi + yi-1)/h2 + 2(yi+1 – yi-1)/2h + yi = xi
ƒ yi+1 – 2yi + yi-1 + hyi+1 – hyi-1 + h2yi = h2xi
ƒ (1 - h)yi-1 + (h2 – 2)yi + (1 + h)yi+1 = ih3, pois xi = ih, 0<i<n
ƒ Como x0 = 0, y0 = 0, xn = 1 e yn = -1, chegamos ao sistema linear abaixo:
y(x) = 2.e-x(1-x) + x - 2
ƒ Soluções com h = 0,1
(a tabela ao lado não 
xi yi solução exata erro
0,1 -0,2720 -0,2713 0,0007(
está completa): 0,2 -0,4911 -0,4900 0,0011
0,3 -0,6641 -0,6629 0,0013
0,4 -0,7969 -0,7956 0,0013
Exemplo: PVC não linear de 2ª ordemExemplo: PVC não linear de 2ª ordemmp mmp m
ƒ y’’ = y.sen y + x.y, onde y(0) = 1 e y(1) = 5
U i ã ’’( ) ( 2 )/h2 O(h2)ƒ Usaremos a aproximação y’’(xi) ≈ (yi+1 – 2yi + yi-1)/h2 com erro O(h2):
ƒ Substituindo-a na equação e considerando x = xi:
ƒ (yi+1 – 2yi + yi-1)/h2 = yi.sen yi + xi.yiyi 1 yi yi 1 yi yi i yi
ƒ yi-1 – yi.[2 + h2(sen yi + ih)] + yi+1 = 0, pois xi = ih
ƒ Como x0 = 0, y0 = 1, xn = 1 e yn = 5, chegamos ao sistema não linear abaixo:
ƒ 1 y [2 + h2(sen y + h)] + y = 0 quando i=1ƒ 1 – y1.[2 + h2(sen y1 + h)] + y2 = 0, quando i=1
ƒ yi-1 - yi.[2 + h2(sen yi + ih)] + yi+1 = 0, para 1<i<n-1
ƒ yn-2 – yn-1.[2 + h2(sen yn-1 + (n-1)h)] + 5 = 0, quando i=n-1
ƒ Soluções com h = 0,1:
yi Resultado
y1 1 3186
yi Resultado
y 3 2091y1 1,3186
y2 1,6513
y3 2,0037
y6 3,2091
y7 3,6525
y8 4,1035
y4 2,3803
y5 2,7829
y9 4,5538