Baixe metodos numéricos para edo's e outras Notas de estudo em PDF para Engenharia Mecânica, somente na Docsity!
M´etodos Num´ericos para
EDO’S
9.1 Introdu¸c˜ao
O estudo das equa¸c˜oes diferenciais foi motivado inicialmente por proble- mas da f´ısica, ou seja problemas de mecˆanica, eletricidade termodinˆamica, magnetismo etc.
Atualmente muitas outras ´areas do conhecimento tem a formula¸c˜ao te´orica de seus problemas utilizando essas equa¸c˜oes. Entre outras podemos as ´areas de Qu´ımica, Ecologia, Biologia, Economia e Sociologia.
Alguns Exemplos de Equa¸c˜oes Diferenciais Ordin´arias
- Seja x(t) o espa¸co percorrido por um corpo em queda livre num ins- tante de tempo t
a) Se desprezarmos a resistˆencia do ar teremos a trajet´oria do corpo regida pela equa¸c˜ao
x′′(t) + g = 0 ; g = ´e a acelera¸c˜ao da gravidade
b) Se considerarmos a resistˆencia do ar, um corpo caindo de para- quedas por exemplo, teremos mx′′(t) + bx′(t) + mg = 0 ; b = cte < 0
- Seja M (t) a massa de um material radioativo no instante t. Sabe-se que a taxa de desintegra¸c˜ao ´e proporcional a quantidade de
massa existente. Assim temos
M ′(t) = kM (t) ; k = cte < 0
Esta mesma equa¸c˜ao descreve a reprodu¸c˜ao de bact´erias numa cultura s´o que neste caso k = cte > 0
9.2 Considera¸c˜oes Gerais sobre EDO’s
Suponha f : Rn^ → R e y : R → R, com derivadas at´e ordem n. Dizemos que uma equa¸c˜ao diferencial ordin´aria tem ordem n se ´e uma equa¸c˜ao que pode ser escrita como
yn^ = f (x, y, y′, · · · , yn−^1 ) ; yn^ = derivada de ordem n (9.1)
Uma fun¸c˜ao y = φ(x) ´e dita uma solu¸c˜ao da equa¸c˜ao (9.1) se:
a) y = φ(x) possue derivadas at´e ordem n;
b) yn(x) = f (x, φ(x), · · · , φn−^1 (x)).
Exemplo 9.2.
y′^ = x ⇒ y = x^2 2
- k k = cte fam´ılia de solu¸c˜oes
Exemplo 9.2.
y′^ = −y ⇒ y = βe−x^ β = cte fam´ılia de solu¸c˜oes
9.3 EDO’s de Primeira Ordem
y′(x) = f (x, y(x))
No que segue sempre consideraremos EDO’s de primeira ordem. As EDO’s de ordem superior ser˜ao reduzidas a um sistema de EDO’s de primeira or- dem.
A demonstra¸c˜ao desse teorema foge aos nossos objetivos.
Voltando ao exemplo 9.4.3 podemos verificar que
f (x, y) = 4x
y ... ∂f ∂y
2 x √ y
Como lim y→ 0
2 x √ y
= ∞ temos que ∂f ∂y
torna-se ilimitada quando y ≈ 0
9.5 M´etodo de Picard
Considere o PVI{ y′^ = f (x, y) y(x 0 ) = y 0
D = {(x, y) ∈ R^2 : |x − x 0 | < a e |x − x 0 | < b}; f e
∂f ∂y cont´ınuas em D
Usando o Teorema Fundamental do C´alculo Integral temos
y(x) − y(x 0 ) =
∫ (^) x
x 0
y′(t) dt =
∫ (^) x
x 0
f (t, y(t) dt ⇒
y(x) = y(x 0 ) +
∫ (^) x
x 0
f (t, y(t) dt (9.2)
De acˆordo com (9.2) definimos o algoritmo de Picard
y 0 (x) := y 0 yk(x) := y 0 +
∫ (^) x x 0 f^ (t, yk−^1 (t))^ dt^ k^ = 1,^2 ,^ · · ·
Pode-se provar que
|φ(x) − yk(k)| ≤ M N k−^1 k!
hk^ ∀x ∈ [x 0 − h, x 0 + h], onde:
φ(x)solu¸c˜ao do PVI h = min{a,
b M
M = max{|f (x, y)|, (x, y) ∈ D} N = max{|
∂f ∂y (x, y)|, (x, y) ∈ D}
Observe que
lim k→∞ yk(x) = φ(x) ∀x ∈ [x 0 − h, x 0 + h]
Exemplo 9.5.
y′^ = x^2 + y^2 y(0) = 0
f (x, y) = x^2 + y^2
y 0 (x) = 0
y 1 (x) =
∫ (^) x
0
f (t, y 0 (t)) dt =
∫ (^) x
0
t^2 dt = x^3 3
y 2 (x) =
∫ (^) x
0
f (t, y 1 (t)) dt =
∫ (^) x
0
t^2 + ( t^3 3
)^2 dt = x^3 3
x^7 63
y 3 (x) =
∫ (^) x
0
f (t, y 2 (t)) dt =
∫ (^) x
0
t^2 + (
t^3 3
t^7 63 )^2 dt =
x^3 3
x^7 63
2 x^11 2079
x^15 59535
... φ(x) ≈
x^3 3
x^7 21
2 x^11 231
x^15 315
9.6 Solu¸c˜ao por S´erie de Taylor
y′^ = f (x, y(x) y(x 0 ) = y 0
Usando o desenvolvimento em s´erie de Taylor temos
yn+1 = y(xn+1) = y(xn + h) = y(xn) + hy′(xn) +
h^2 2 y′′(xn) +
h^3 6 y′′′(xn)+
hk k!
yk(xn) + hk+ (k + 1)!
yk+1(μ) μ ∈ (xn, xn + h)
y′(x) = f (x, y(x))
y′′(x) = d dx
f (x, y(x)) = ∂f ∂x
dx + ∂f ∂y
dy = ∂f ∂x
Assim fazendo k = 2 obtemos a seguinte f´ormula
yn+1 = yn + hf + h^2 2
∂f ∂x
) + y′′′(μ) h^3 6
Podemos ent˜ao definir
yn+1 = yn + h[f + h 2
∂f ∂x
)]
O erro cometido ser´a y′′′(μ)h^3 /6 e assim o erro ´e de ordem e(h) = O(h^3 )
E claro que podemos obter f´´ ormulas com erros de ordem superior utilizando mais termos na s´erie de Taylor. No entanto isso nos obrigaria a calcular y′′′(x), yiv(x),etc. E claro que o c´´ alculo dessas derivadas de ordem superior apesar de poss´ıvel ´e bastante penoso. Vamos nos contentar com o caso em que k = 2.
Algoritmo de Taylor de Ordem 3
Dado o PVI { y′^ = f (x, y) y(x 0 ) = y 0
Dado h > 0 considere xi = x 0 + hi, i = 0, · · · , n.
Podemos ent˜ao construir uma tabela de yi aproxima¸c˜ao para y(xi) onde y(x)
´e a solu¸c˜ao do PVI usando o algoritmo de Taylor
i = 0, · · · , n
T (xi, yi) = f (xi, yi) + h 2
[
∂f ∂x
(xi, yi) + f (xi, yi) ∂f ∂y
(xi, yi)] yi+1 = yi + hT (xi, yi)
Exemplo 9.6. Dado o PVI { y′^ = y y(0) = 1 Fa¸ca uma tabela para a solu¸c˜ao aproximada do PVI usando h = 0.1 e n = 5.
Solu¸c˜ao
f (x, y) = y;
∂f ∂x
∂f ∂y
... T (x, y) = y +
h 2 (f (x, y).1) = y(1 +
h 2
yi+1 = yi + hT (xi, yi) = yi[1 + h(1 +
h 2
)]
Os valores yi calculados atrav´es da equa¸c˜ao acima est˜ao listados na tabela abaixo bem com os valores da solu¸c˜ao exata do PVI que nesse caso ´e conhe- cida e ´e dada por y(x) = ex. Tamb´em listamos a diferen¸ca entre a a solu¸c˜ao exata e a aproximada.
xi 0.00 0.10 0.20 0.30 0.40 0. yi 1.00000 1.10500 1.22102 1.34923 1.49090 1. exi^ 1.00000 1.10517 1.22140 1.34986 1.49182 1. exi^ − yi 0.00000 0.00017 0.00038 0.00063 0.00092 0.
AT EN C¸ AO˜! Se um algoritmo tem ei = O(hp) ou seja se o erro em cada etapa ´e proporcional `a hp^ ent˜ao E = O(hp−^1 ) pois temos n itera¸c˜oes e n = (a − x 0 )/h.
9.8 M´etodos de Passo-Simples
Defini¸c˜ao 9.8.1 Um algoritmo para resolver um PVI ´e dito ser de passo simples se a aproxi¸c˜ao yi+1 for calculada utilizando apenas o resultados do passo anterior.
Assim podemos definir os algoritmos de passo simples como sendo aqueles do tipo
yi+1 = yi + φ(xi, yi)
9.9 M´etodo de Euler
Seja
y′^ = f (x, y) y(x 0 ) = y 0
Observando que
y(x + h) = y(x) + hy′(x) + h^2 2 y′′(μ) μ ∈ (x, x + h)
temos que
y(xi + h) ≈ y(xi) + hy′(xi, yi) = y(xi) + hf (xi, yi)(fig 9.1)
Definimos ent˜ao o Algoritmo de Euler
Dados x 0 , y 0 , h. Geramos aproxima¸c˜oes yi para y(xi) atrav´es de
i = 0, 1 , 2 , 3 , · · · { yi+1 = yi + hf (xi, yi) xi+1 = xi + h
AT EN C¸ AO˜! Note que o algoritmo de Euler ´e equivalente `a solu¸c˜ao por s´erie de Taylor com k = 1. Assim temos
ei = O(h^2 ) e ent˜ao E = O(h).
9.10 Interpreta¸c˜ao Geom´etrica do M´etodo de Eu-
ler
x
y
y(x)
xi
yi
xi+
y(xi+1)
r yi+
Figura 9.1: M´etodo de Euler
Seja r a reta tangente a curva y = y(x) no ponto (xn, yn) e s a reta vertical passando por (xi+1, 0).
Assim temos { r : yi + y′(xi)(x − xi) s : x = xi+
Observe agora que yi+1 ´e a itersec¸c˜ao das retas r e s pois a intersec¸c˜ao ´e dada por
yi + y′(xi)(xi+1 − xi) = yi + y′(xi)(h) = yi + hf (xi, yi)
Exemplo 9.10.
Dado o PVI
y′^ = −y y(0) = 1
Fa¸ca uma tabela da solu¸c˜ao aproximada, usando o m´etodo de Euler, com
Exemplo 9.10.2 Considere o PVI
y′^ = y y(0) = 1
Usando M´etodo de Euler Determine y(a).
(Observe que y(x) = ex^ ´e a solu¸c˜ao exata do PVI)
Solu¸c˜ao
h = a/n ⇒ a = nh
yi+1 = yi = hf (xi, yi) = yi + hyi = yi(1 + h)
y 0 = 1
y 1 = 1(1 + h)
y 2 = (1 + h)(1 + h) = (1 + h)^2
y 3 = (1 + h)^2 (1 + h) = (1 + h)^3 .. .
yn = (1 + h)n^ = (1 + a/n)n
Como y(xn) = y(x 0 + nh) = y(nh) = y(a) temos que
y(a) = y(xn) ≈ yn = (1 + a/n)n
Observe que lim n→∞ yn = lim n→∞ (1 + a/n)n^ = ea
9.11 M´etodo de Heun
Analizando o m´etodo de Euler observamos que ele considera a inclina¸c˜ao y′(x) constante no intervalo [xi, xi + h] e toma essa constante como a incli- na¸c˜ao no extremo esquerdo do intervalo.(Veja figura 9.1)
O m´etodo de Heun ao inv´es disso utiliza o valor dessa constante como a m´edia aritm´etica das inclina¸c˜oes nos pontos extremos do intervalo [xi, xi +h].
Podemos ent˜ao considerar
yi+1 = yi + h 2
(y′(xi) + y′(xi+1))
Usando o PVI podemos escrever as equa¸c˜oes
y′(xi) = f (xi, yi)
y′(xi+1) = f (xi+1, yi+1)
Como na ´ultima equa¸c˜ao n˜ao conhecemos o valor de yi+1 vamos utilizar
aproxima¸c˜ao dada por yi+1 = yi + hf (xi, yi).
Assim teremos
yi+1 = yi +
h 2 (f (xi, yi) + f (xi+1, hf (xi, yi))
Podemos ent˜ao considerar o algoritmo de Heun. Dados x 0 , y 0 , h. Geramos aproxima¸c˜oes yi para y(xi) atrav´es de
i = 0, 1 , 2 , 3 , · · ·
k 1 = f (xi, yi)
k 2 = f (xi + h, yi + hk 1 )
yi+1 = yi + h 2
(k 1 + k 2 )
xi+1 = xi + h
AT EN C¸ AO˜! Pode ser demonstrado que o m´etodo de Heun tem erro local ei = O(h^3 ) e erro global E = O(h^2 )
Exemplo 9.11.1 Vamos resolver o mesmo problema do exemplo 9.10.1 uti- lizando o m´etodo de Heun no sentido de comparar os resultados.
Dado o PVI
y′^ = −y y(0) = 1
Fa¸ca uma tabela da solu¸c˜ao aproximada, usando o m´etodo de Heun, com h = 0.1 e n = 10.
Solu¸c˜ao
x 0 = 0; y 0 = 1; f (x, y) = −y
k1 = f (xi, yi) = −yi
k2 = f (xi + h, yi + hk 1 ) = −(yi + hk 1 )
yi+1 = yi + h 2
(k 1 + k 2 )
A tabela abaixo mostra os resultados das itera¸c˜oes bem como as compara¸c˜oes com a solu¸c˜ao exata y(x) = e−x
O algoritmo de Hunge-Kutta ser´a ent˜ao dado por
Dados x 0 , y 0 , h. Geramos aproxima¸c˜oes yi para y(xi) atrav´es de
i = 0, 1 , 2 , 3 , · · ·
k 1 = f (xi, yi)
k 2 = f (xi +
h 2 , yi +
h 2 k 1 )
k 3 = f (xi + h 2
, yi + h 2
k 2 )
k 4 = f (xi + h, yi + hk 3 )
yi+1 = yi +
h 6 (k 1 + 2k 2 + 2k 3 + k 4 )
xi+1 = xi + h
Exemplo 9.12.
Dado
y′^ = 1 − x + 4y y(0) = 1
Determinar uma aproxima¸c˜ao para y(0.2) usando h = 0.2.
Solu¸c˜ao
f (x, y) = 1 − x + 4y
k 1 = f (0, 1) = 5
k 2 = f (0. 1 , 1 .5) = 6. 9
k 3 = f (0. 1 , 1 .69) = 7. 66
k 4 = f (0. 2 , 2 .532) = 10. 928
y 1 = 1 + 0. 2 /6[5 + 2(6.9 + 7.66) + 10.928] = 2. 5016
Exemplo 9.12.
Dado
y′^ = x^2 + y^2 y(0) = 0
Determine uma aproxima¸c˜ao para y(xi) usando h = 0.05 e n = 6.
Solu¸c˜ao
x 0 = 0; y 0 = 0; f (x, y) = x^2 + y^2
k 1 = f (xi, yi) = x^2 i + y i^2
k 2 = f (xi + h/ 2 , yi + h/ 2 k 1 ) = −(yi + hk 1 )
yi+1 = yi + h 2
(k 1 + k 2 )
A tabela abaixo mostra os resultados das itera¸c˜oes.
yi 0.00 0.05 0.10 0.15 0.20 0.25 0. yi 0.00 0.000813 0.003000 0.007315 0.014513 0.025364 0.
9.13 M´etodos de Predi¸c˜ao-Corre¸c˜ao
Os m´etodos do tipo predi¸c˜ao-corre¸c˜ao s˜ao baseados na seguinte t´ecnica
Vamos considerar o PVI
y′^ = f (x, y) y(x 0 ) = y 0
Pelo teorema fundmental do c´alculo integral temos
y(xk+1) = y(xk) +
∫ (^) xk+
xk
y′(t) dt =
∫ (^) xk+
xk
f (t, y(t)) dt (9.3)
Observando que f (t, y(t)) ´e uma fun¸c˜ao unicamente da vari´avel t podemos considerar g(t) = f (t, y(t)) e utilizar um m´etodo de integra¸c˜ao num´erica para calcular
∫ (^) xk+ xk g(t)^ dt.
Usando, por exemplo, a regra do Trap´ezio teremos: ∫ (^) xk+
xk
g(t) dt =
h 2 (g(xk) + g(xk+1)) =
h 2 (f (xk, yk) + f (xk+1, yk+1))
Usando a equa¸c˜ao (9.3) teremos
yk+1 = yk + h 2
(f (xk, yk) + f (xk+1, yk+1)
Exemplo 9.13.
Dado o PVI
y′^ = y y(0) = 1
Fa¸ca uma tabela para y(x) de x = 1 at´e x = 0.6 com h = 0.1 e usando o corretor 3 vezes em cada etapa.
Solu¸c˜ao
y 10 = y 0 + hf (0, 1) = 1 + (0.1)1 = 1. 1
y 11 = y 0 + h/2(f (0, 1) + f (0. 1 , 1 .1)) = 1 + (0.1)/2(1 + 1.1) = 1. 105
y 12 = y 0 + h/2(f (0, 1) + f (0. 1 , 1 .105)) = 1 + (0.1)/2(1 + 1.105) = 1. 10525 y 13 = y 0 + h/2(f (0, 1) + f (0. 1 , 1 .10525)) = 1 + (0.1)/2(1 + 1.10525) = 1. 1052625
As demais itera¸c˜oes est˜ao na tabela abaixo bem como a solu¸c˜ao exata como ´ultima coluna.
xi y^0 i y i^1 y^2 i y^3 i exi 0.10 1.1000000 1.1050000 1.1052500 1.1052625 1. 0.20 1.2157887 1.2213151 1.2215914 1.2216052 1. 0.30 1.3437657 1.3498737 1.3501791 1.3501944 1. 0.40 1.4852139 1.4919648 1.4923024 1.4923192 1.
Uma quest˜ao interessante a ser discutida ´e sobre o valor de m, ou seja, sobre quantidade de vezes que devemos utilizar o corretor. Uma t´ecnica para se determinar uma estimativa para m ´e a seguinte:
Vamos considerar a t´ıtulo de praticidade o exemplo anterior ou seja { y′^ = y y(0) = 1
y′(x) = y(x) ⇒ y′′(x) = y′(x) ⇒ y′′(0) = y′(0) = y(0) = 1
y′′′(x) = y′′(x) ⇒ y′′′(0) = y′′(0) = 1 ⇒
Como o erro da f´ormula do corretor ´e dada por ETi ≈
h^3 12 y′′′(ξ) podemos
considerar para efeitos pr´aticos que y′′′(ξ) ≈ y′′′(0) = 1
... EiT ≈
(0.1)^3
× 1 = 0. 00033
Assim, para o exemplo anterior, n˜ao se deve esperar mais do que 3 casas decimais exatas, bastando usar o corretor duas vezes.
Exemplo 9.13.
Dado o PVI
y′^ = x −
y y(0) = 1
Fa¸ca uma tabela para y(x) de x = 0 at´e x = 0.4 com h = 0.1. Determine uma estimativa para m.
Solu¸c˜ao
y′(x) = x −
y(x) ⇒ y′(0) = 0 −
y(0)
y′′(x) = 1 +
y′(x) (y(x))^2 ⇒ y′′(0) = 1 +
y′(0) (y(0))^2
y′′′(x) = −
y′′(x)(y(x))^2 − 2(y′(x))^2 y(x) y(x)^4
y′′′(0) = − y′′(0)(y(0))^2 − 2(y′(0))^2 y(0) (y(0))^4
... EiT ≈ 2
(0.01)^3
Assim n˜ao devemos esperar mais do que 3 casas decimais exatas.
xi y^0 i y^1 i y i^2 y^3 i 0.10 0.9000000 0.8994444 0. 0.20 0.7982261 0.7961792 0. 0.30 0.6903929 0.6857831 0. 0.40 0.5693739 0.5595193 0.5579727 0.
Podemos ent˜ao considerar as aproxima¸c˜oes
y(0.1) ≈ 0. 899 y(0.2) ≈ 0. 796 y(0.3) ≈ 0. 685 y(0.4) ≈ 0. 557
