EDO PVI

Prévia do material em texto

Capítulo 1
Solução numérica de problemas de valor inicial
1.1 Introdução
Começamos pela seguinte pergunta: o que é uma equação diferencial ordinária (EDO)? Resposta: É uma equação que envolve
uma ou mais derivadas de uma função incógnita que depende de uma única variável independente. Vejamos alguns exemplos.
Exemplo 1. A seguinte equação:
du
dt
(t) = 2,
é uma equação diferencial ordinária.
• Neste caso, a EDO acima contém somente a derivada primeira da função u;
• Procura-se saber a expressão u = u(t) que satisfaz a EDO acima.
Exemplo 2. Determinar a função y = y(x), tal que:
d2y(x)
dx2
= y(x), ∀ x ∈ R (1.1.1)
A eq. (1.1.1) é uma equação diferencial ordinária, onde aparecem a função y = y(x) e derivada segunda desta função.
Observação 3. A ordem de uma equação diferencial ordinária é a ordem da maior derivada que aparece na EDO. Logo, a ordem
da EDO (1.1.1) é 2.
Seja u = u(x) uma função de uma variável, as seguintes notações são usadas na literatura:
u(0)(x) := u(x)
u(1)(x) :=
du(x)
dx
...
u(n)(x) :=
dnu(x)
dxn
Usando-se a notação acima, a seguinte EDO:
x2
d2u(x)
dx2
+ sen(x)
du(x)
dx
− exln(x) = 0
pode ser expressa da forma a seguir:
x2u(2)(x) + sen(x)u(1)(x)− exln(x) = 0
Observação 4. Também, são usadas as seguintes notações para as derivadas da função u = u(x):
u
′
(x) := u(x)
u
′′
(x) :=
d2u(x)
dx2
u
′′′
(x) :=
d3u(x)
dx3
uIV (x) :=
d4u(x)
dx4
1
Por que resolver EDOs?
• As equações diferenciais ordinárias são instrumentos essenciais para a modelação de fenômenos de várias áreas como a
física, química, engenharia, economia, ecologia etc.
• A solução de uma EDO permite conhecer o comportamento quantitativo e qualitativo da função procurada.
• O conhecimento do comportamento da função procurada permite compreender melhor o fenômeno estudado e assim propor
medidas, caso sejam necessárias, para atingir objetivos da pesquisa em andamento.
A seguir apresentamos alguns problemas físicos modelados por EDOs.
Exemplo 5. (Temperatura de um corpo imerso em um meio de temperatura constante Te) Consideremos um corpo com
temperatura T imerso em um ambiente com temperatura Te. Suponhamos que se considere que massa do corpo concentra-se
em um ponto do corpo. Então, a taxa de transferência de calor entre o corpo e o meio, v(t), pode ser descrita pela lei de
Stefan-Boltzmann, isto é:
v(t) = �γS(T 4(t)− T 4e ).
Onde:
• t é a variável tempo;
• � é constante de Boltzmann, isto é, � = 5.6x10−8 J/m2K4s (as letras: J, k, m e s, representam unidades de Joule, kelvin,
metro e segundo, respectivamente);
• γ é a constante de emissibilidade do corpo;
• S é a área superficial do corpo.
E a taxa de variação da temperatura do corpo pode ser modelada pela seguinte EDO:
dT
dt
(t) = −�γS(T
4(t)− T 4e )
MC
(1.1.2)
T (0) = T0
onde M , T0 e C são a massa do corpo, a temperatura inicial do corpo e o calor específico do material que consitui o corpo,
respectivamente.
Neste caso, procura-se a função T=T(t), solução da EDO (1.1.2)
Exemplo 6. (Dinâmica de população) Considere uma população de bactérias confinadas em um ambiente que pode ter, no
máximo, B elementos como população. Assuma que no tempo inicial a população de bactérias é y0 � B (o símbolo � significa
muito menor que) e que a taxa de crescimento das bactérias seja constante C > 0. Nesse caso, a taxa de mudança de população
é proporcional ao número de bactérias existentes.
A EDO que modela a taxa de crescimento da população de bactérias, y(t), que considera no máximo uma população com B
elementos, é dada por:
dy
dt
= Cy
(
1− y
B
)
Considerando-se duas populações de bactérias y1 e y2 em competição de sobrevivência. Uma proposta para modelagem da
dinâmica de crescimento de ambas as populações é dada pelas seguintes EDOs:
dy1
dt
=C1y1 (1− b1y1 − d2y2) ,
dy2
dt
=C2y2 (1− b2y2 − d1y1) .
Onde: Ci e bi são a taxa de crescimento e a quantidade de nutrientes associadas à população yi. E as constantes d1e d2 indicam
o tipo de interação entre as duas populações.
Neste caso, procuram-se as funções y1 = y1(t) e y2 = y2(t).
Exemplo 7. ( Trajetória de uma bola de beisebol) Deseja-se modelar a trajetória de uma bola lançada pelo Pitcher para o
Catcher (Veja Figura abaixo)
2
As equações que descrevem o movimento da bola são:
dx
dt
=v (1.1.3)
dv
dt
=F (1.1.4)
onde: x = (x(t), y(t), z(t))t designa a posição da bola no tempo t, v = (vx, vy, vz)t é a velocidade da bola e F um vetor do tipo
força, causador da trajetória descrita pela bola. As componentes da força F são:
Fx =− f(v)vvx +Bw(vxsen(φ)− vycos(φ)).
Fy =− f(v)vvy +Bwvxcos(φ).
Fz =− g − f(v)vvx −Bwvxsin(φ).
onde as constantes φ e w são o ângulo e o módulo da velocidade angular impressa na bola pelo Pitcher no momento do chute.
A constante B = 4.1 10−4 e f(v) é o coeficiente de frição definido por:
f(v) = 0.0039 +
0.0058
1 + e
v−35
5
e v = ||v|| (módulo do vetor velocidade).
Resolvendo-se o sistema (1.1.3) e (1.1.4), obtemos a posição e a velocidade da bola para todo instante de tempo t.
Exemplo 8. Considere o circuito elétrico linear mostrado na figura abaixo:
A fonte de entrada para o circuito é vs(t) e, considerando-se, a saída vC(t) como a diferença de voltagem através do capacitor
C e, analisando-se as malhas do circuito pode-se determinar uma equação diferencial ordinária que modela a relação entre vs(t)
e vC(t).
A equação diferencial ordinária que relaciona a voltagem de entrada vs(t) e a voltagem de saída vC(t) é dada por:
L1C
d³vC(t)
dt3
+R1C
d2vC(t)
dt2
+
(
L1
L2
+ 1
)
dvC(t)
dt
+
R1
L2
vC(t) =
dvs(t)
dt
(1.1.5)
• Os dados para a EDO acima são: L1, L2, R1, R2, C e a voltagem vs(t).
• Procura-se a voltagem vC(t).
• A ordem da EDO (1.1.5) é 3.
3
1.1.1 Solução analítica de uma equação diferencial ordinária
Toda equação diferencial ordinária de ordem n pode ser expressa como:
F (x, u(x), u(1)(x), u(2)(x), . . . , u(n)(x)) = 0,
onde F é uma função de (n+ 2) variáveis independentes com valor em R. Por exemplo:
u
′′
= u′ − x2u3 → u′′ − u′ + x2u3 = 0→ F (x, y, z, w) = w − z + x2y3
Onde: y = u, z = u′, w = u
′′
Definição 9. Uma função u = u(x), com x ∈ I (intervalo), é solução da equação diferencial ordinária de ordem n:
F (x, u(x), u(1)(x), u(2)(x), . . . , u(n)(x)) = 0, (1.1.6)
se, e somente se, a função u satisfaz a equação (1.1.6).
Isto implica que a candidata a solução, u = u(x). deve possuir todas as derivadas que aparecem na Eq. (1.1.6), isto é, devem
existir u(i)(x), i = 1, 2, 3, ...n.
Observação 10. Em alguns casos usaremos u = u(x) ou u = u(t), sendo u a variável depedente e x ou t a variável independente,
respectivamente.
Exemplo 11. Verifique que a função u(t) = k sen(t) , com k constante, é solução da EDO:
u
′′
(t) = −u(t), ∀t ∈ R (1.1.7)
Solução. A função candidata a solução é:
u(t) = k sen(t).
Logo,
u
′
(t) = k cos(t) e u
′′
(t) = −k sen(t) = −u(t).
Então,
u
′′
(t) = −u(t), ∀t ∈ R.
Portanto, u(t) = k sen(t) é solução da EDO (1.1.7).
Verifique que a função y(t) = t2cos(t) é solução da seguinte EDO:
tu
′
(t)− 2u(t) + t3sen(t) = 0 (1.1.8)
Solução. Derivando-se a função y(t) com relação a t, tem-se:
y
′
(t) = 2tcos(t)− t2sen(t) (1.1.9)
Multiplicando a Eq. (1.1.9) por t, obtém-se:
ty
′
(t) = 2t2cos(t)− t3sen(t). (1.1.10)
Mas t2cos(t) = y(t), então a relação (1.1.10) pode ser expressa por:
ty
′
(t) = 2y(t)− t3sen(t).
Logo,
ty
′
(t)− 2y(t) + t3sen(t) = 0
Portanto, a função y(t) = t2cos(t) satisfaz a EDO (1.1.8).
1.1.2 Como determinar soluções analíticas de EDOs?
Cada EDO precisa de um método específico para encontrar sua solução analítica, caso seja possível determinar a solução. Vejamos
um exemplo.
Exemplo 12. Determinar a solução da seguinte EDO:
4d2u
dx2
= 0 (1.1.11)
Solução. Integrando-se a equação (1.1.11) em relação a x, isto é:
ˆ
d
dx
(
du
dx
)
dx =
ˆ
0dx
Tem-se:
du
dx
+ C = 0 ou
du
dx
= C1 (C1 = −C) (1.1.12)
Integrando-se a Eq. (1.1.12) em relação a x, tem-se:
u = C1x+ C2 (1.1.13)
Observação 13. Do visto acima, temos:
• A solução para a EDO (1.1.11) dada por (1.1.13) define uma família de soluções.
• Para determinar uma solução específica da EDO, precisamos conhecer os valores das constantes C1e C2. Para isto,
precisamos de duas informações que permitem determinar o valor dessas constantes.
Figura 1.1.1: Família de soluções da EDO (1.1.11)
Exemplo 14. Determinar a solução da seguinte EDO:
du
dx
= −λu, λ > 0 (Constante dada), (1.1.14)
e esboçar o gráfico da solução.
Solução. A Eq. (1.1.14) pode ser expressa por:
du
u
= −λdx (1.1.15)
Integrando-se a Eq. (1.1.15), isto é: ˆ
du
u
=
ˆ
(−λ)dx,
obtém-se:
ln(|u(x)|) = −λx+ C
Então,
|u(x)| = e(−λx+C)
ou
|u(x)| = Ae−λx, A = eC > 0 (1.1.16)
Caso tenhamos a informação que a função u = u(x) possui valores positivos, então da Eq. (1.1.16) obtém-se a seguinte
expressão para u:
u(x) = Ae−λx (A > 0). (1.1.17)
No caso da função u = u(x) seja negativa para todo x, então a expressão para u é:
u(x) = −Ae−λx (A > 0) (1.1.18)
A expressão (1.1.17) ou (1.1.18) define uma família de soluções para a EDO. A escolha de um integrante dessa família de
soluções depende da informação adicional fornecida. Supondo que a informação adicional seja do tipo:
u(x0) = u0,
onde x0 e u0 são dados, então a solução analítica é:
u(x) = Ae−λx, se u0 > 0.
u(x) = −Ae−λx, se u0 < 0.
Onde o valor da constante A deve ser determinando usando-se a informação adicional.
Por exemplo, para u(0) = 1, a solução analítica é: u(x) = Ae−λx. Substituindo a informação adicional, temos que: 1 =
Ae−λ 0 = A. Logo, a solução procurada que passa pelo ponto (0,1) é u(x) = e−λx.
5
Gráfico das soluções dadas pelas Eqs. (1.1.17) ou (1.1.18), dependendo da informação adicional.
u0 = −6,−4,−2, 0, 2, 4, 6. x ∈ [−1, 1].
1.2 Problema de valor inicial
É uma EDO da forma:
{
y
′
= f(t, y), t ∈ (a, b) ⊂ R
y(a) = y0
(1.2.1)
Os dados do problema de valor inicial (PVI) são a função f(t, y) e o ponto (a, y0) ∈ R2
• Dom y = [a, b] (Intervalo) e Im y ∈ R;
• Domf ⊂ R2
Observe que neste caso a informação adicional é dada no extremo esquerdo do intervalo, daí o nome da informação adicional de
condição inicial.
Exemplo 15. Considere a seguinte EDO: {
y
′
= −2y x ∈ [a, b].
y(a) = 1
(1.2.2)
Ela possui a forma dada em (1.2.1) logo é um problema de valor inicial com: f(x, y) = −2y e y0 = 1.
Exemplo 16. Considere o seguinte problema de valor inicial:{
dy
dt =
2−y
t , t ∈ (1, 2, 5)
y(1) = y0.
(1.2.3)
A solução analítica dessa equação é:
y(t) =
(y0 − 2)
t
+ 2. (1.2.4)
Algumas soluções do PVI anterior são:
y0 = 2 y(t) = 2
y0 = 3 y(t) =
1
t
+ 2
6
y0 = 6 y(t) =
4
t
+ 2
y0 = −3 y(t) = −5
t
+ 2
y0 = −6 y(t) = −4
t
+ 2.
Graficamente
Todo PVI possui uma única solução? Resposta. Não.
Exercício 17. Verificar que a expressão (1.2.4) é a solução analítica do PVI (1.2.3).
Exemplo 18. Calcular a solução do seguinte PVI:
{
x′ = x2/3 t > 0
x(0) = 0
(1.2.5)
Solução:
Pode-se verificar que as seguintes funções são solução do PVI (1.2.5):
x = x(t) = 0 ∀ t ≥ 0.
x = x(t) =
1
27
t3 ∀ t ≥ 0.
xa(t) =
{
0 0 ≤ t ≤ a
1
27 (t− a)3 a ≤ t
Então, não temos uma única solução do PVI (1.2.5). Existe algum teorema que nos indique sob que condições o PVI tem
uma única solução? A resposta é sim!!!
Teorema 19. (Existência) Seja f uma função contínua em um retângulo B centrado em (x0, y0), isto é:
B = {(x, y) / |x− x0| ≤ α, |y − y0| ≤ β} , (1.2.6)
então o PVI (1.2.1) tem uma solução y = y(x) no intervalo x0 − r ≤ x ≤ x0 + r, onde r = min {α, β/M} , sendo M =
max(x,y)∈B |f(x, y)|.
Teorema 20. (Unicidade) Se f e ∂f∂x são continuas no retângulo B definido em (1.2.6), então o PVI (1.2.1) tem uma única
solução em
{
x /|x− x0| ≤ min{α, β/M}, M = max(x,y)∈B |f(x, y)|
}
.
Observação 21. Temos:
•
{
x /|x− x0| ≤ min{α, β/M}, M = max(x,y)∈B |f(x, y)|
}
= [x0 − r, x0 + r], r = min{α, β/M}, M = max(x,y)∈B |f(x, y)|
• o intervalo [x0−r, x0+r] dado pelo teorema de existência ou pelo teorema de existência e unicidade pode ser muito menor
que o intervalo para t que define o retângulo B;
• ambos os teoremas garantem a existência de solução localmente.
7
1.3 Métodos numéricos para PVI
Vejamos um exemplo.
Exemplo 22. (Solução analítica e solução numérica) Seja o problema de valor inicial:
{
y
′
= ex
y(0) = 1
(1.3.1)
Neste caso:
• a função f(x, y) = ex;
• a função f e a ∂f∂x (x, y) = e
x são funções contínuas. Logo, temos existência e unicidade de solução em torno do ponto
(x0, y0) = (0, 2);
• a solução analítica do PVI acima é y(x) = ex. (Verificar)
A Figura abaixo mostra o esboço das funções y(x) = ex e y = ex + 2, sendo que a primeira função é a solução do PVI acima.
Figura 1.3.1:
Para obter uma solução numérica em algum intervalo [a, b] deve-se:
• escolher abscisas do intervalo [a, b] onde deseja-se conhecer valores aproximados para a solução do PVI;
• usar alguma expresão discreta que seja consistente com o PVI e com os pontos onde deseja-se determinar a solução.
Da figura abaixo, podemos observar que yk é uma aproximação para y(xk).
x x
y(x )
y
y
k−1
k
k
k−1 k
y(x )
k−1
8
O valor desta aproximação, yk, pode ser calculado usando-se o conceito de derivada da função y no ponto (xk−1, y(xk−1)) e
do triângulo determinado pelos pontos (xk−1, yk−1), (xk,yk−1) e (xk, yk) (Veja Figura acima). Isto é:
yk − yk−1
xk − xk−1 = y
′
(xk−1)
Como y é solução de um PVI, então y
′
(xk−1) = f(xk−1, yk−1), logo:
yk − yk−1
xk − xk−1 = f(xk−1.yk−1) (1.3.2)
ou
yk = yk−1 + (xk − xk−1)f(xk−1.yk−1). (1.3.3)
Onde:
• as notações yk−1 e yk representam o valor aproximado para y(xk−1) e y(xk), respectivamente;
• a fórmula (1.3.3) indica como calcular uma aproximação para y em xk, yk, em função das informações conhecidas em xk−1.
Neste caso, em função de yk−1 e f(xk−1.yk−1);
• o lado esquerdo da relação (1.3.2) é um aproximação para a primeira derivada de y em xk.
Para a função f do PVI (1.3.1), a expressão (1.3.3) fica:
yk = yk−1 + (xk − xk−1)exk−1 (1.3.4)
Suponhamos que desejamos obter aproximações para y nos pontos x0 = 0, x1 = 0.25, x2 = 0.5, x3 = 0.75 e x4 = 1.0,
Lembremos que em x0 = 0 conhecemos o valor exato da solução, isto é, conhecemos y(x0) = y0 = 1, então vamos aplicar a regra
discreta (1.3.3) para k = 1, 2, e 3. Então:
• k=1:
y1=y0 + (x1 − x0)ex0 = 1 + (0.25− 0.0)e0.0 = 1.25
• k=2:
y2=y1 + (x2 − x1)ex1 = 1.25 + (0.5− 0.25)e0.25 = 1.5710064
• k=3:
y3=y2 + (x3 − x2)ex2 = 1.5710064 + (0.75− 0.5)e0.5 = 1.9831867
• k=4:
y4=y3 + (x3 − x2)ex3 = 1.9831867 + (1− 0.75)e0.75 = 2.5124367
Suponhamos que dividimos em 8 partes iguais o intervalo e deseja-se calcular aproximações para y nos pontos:
x1 = 0.125, x2 = 0.250, x3 = 0.375, x4 = 0, 5, x5 = 0.625, x6 = 0.75, x7 = 0.875, x9 = 1.0
Vamos calcular aproximações para y(xi), i = 1, 2, 3, 4, 5, 6, 7, 8 usando a Eq. (1.3.4) do mesmo PVI acima.
k xk yk ek = y(xk)− yk
0 0 1 0.
1 0.125 1.125 0.0081485
2 0.250 1.2666436 0.0173819
3 0.375 1.4271467 0.0278447
4 0.5 1.6090207 0.0397006
5 0.625 1.8151108 0.0531351
6 0.75 2.0486416 0.0683585
7 0.875 2.3132666 0.0856087
8 1.0 2.613126 0.1051559
9
Figura 1.3.2: Solução exata (curva vermelha) e solução numérica (pontos discretos).
Figura 1.3.3: Erro: e(xi) = y(xi)− yi, onde y(xi) e yi,são o valor exato e aproximado no ponto xi, respectivamente.
10
Figura 1.3.4: Solução exata (curva) e solução aproximada (pontos) para n=8.
Figura 1.3.5:Erro: e = y(xi)− yi.
11
Então, o que são métodos numéricos para PVIs de primeira ordem? São fórmulas discretas e recursivas associadas
ao PVI que permitem obter soluções aproximadas para a solução exata do PVI .
O método apresentado acima é conhecido como o método de Euler e é reescrito abaixo:
yk = yk−1 + hf(xk−1, yk−1),
onde h = xk − xk−1, denominado tamanho de passo.
Como foram obtidas algumas das fórmulas discretas (métodos numéricos) para resolver o PVI? Foram obtidos
usando-se:
1. Fórmula de Taylor;
2. Integração numérica. Aqui veremos o método de obtenção das fórmulas pela fórmula de Taylor.
1.3.1 Método explícito e implícito
Um método numérico para o PVI é explícito quando o cálculo da aproximação yk é direto, isto é, não é necessário resolver
nenhuma equação para o cálculo desta. Caso contrário, o método numérico é implícito.
Exemplo:
• Método explícito: yk = yk−1 + hf(xk−1, yk−1)
• Método implícito: yk = yk−1 + hf(xk, yk) —> Para calcular yk precisa-se resolver uma equação não linear.
1.4 Métodos de passo simples
Se para calcular yk usa-se apenas a informação em yk−1 tem-se um método de passo simples ou passo um. Se além das informações
em yk−1, o método precisa de outros valores anteriores, isto é, yk−2, yk−3 etc, então o método é denominado de passo múltiplo.
Exemplos
• Método de Euler: yk = yk−1 + hf(xk−1, yk−1) é um método de um passo;
• Método de dois passos: yk+2 = yk+1 + h3 (3f(xk+1, yk+1)− 2f(xk, yk)).
1.4.1 Métodos numéricos baseados na fórmula de Taylor
Teorema 23. Seja g ∈ Cn+1[a, b] (isto é, g tem derivadas contínuas de ordem 1, 2, . . . ,n+ 1). Suponha que x0 ∈ [a, b] é uma
valor fixo, então para todo x ∈ [a, b] existe um número c = c(x) entre x0 e x, tal que:
g(x) = Pn(x) +Rn(x),
sendo:
Pn(x) =
n∑
k=0
1
k!
g(k)(x0)(x− x0)k, (1.4.1)
e
Rn(x) =
1
(n+ 1)!
g(n+1)(c)(x− x0)n+1. (1.4.2)
• O polinômio de grau n, Pn(x) =
∑n
k=0
1
k!g
(k)(x0)(x− x0)k é denominado o polinômio de Taylor de ordem n em torno do
ponto x0 da função g;
• O polinômio de Taylor de ordem n da função g não é um polinômio interpolador associado a g;
• O polinômio de Taylor de g pode ser usado para aproximar a função g, isto é:
g(x) ≈ Pn(x),
o erro dessa aproximação, En(x) = Rn(x) = 1(n+1)!g
(n+1)(c)(x− x0)n+1, onde c é um valor entre x0 e x.
Exemplo 24. Usando-se o polinômio de Taylor para y = y(x) em torno do ponto x0, determinar aproximações para y(x0 + h).
12
Solução: Sabe-se que o polinômio de Taylor de ordem n em torno de x0 para a função y(x) é dado por:
Pn(x) =
n∑
k=0
1
k!
y(k)(x0)(x− x0)k
ou
Pn(x) = y(x0) + y
(1)(x0)(x− x0) + 1
2
y(2)(x0)(x− x0)2 + · · ·+ 1
n!
y(n)(x0)(x− x0)n
Também sabe-se que a função y = y(x) pode ser aproximada por seu polinômio de Taylor , isto é:
y(x) ≈ Pn(x)
ou seja:
y(x) ≈ y(x0) + y(1)(x0)(x− x0) + 1
2
y(2)(x0)(x− x0)2 + · · ·+ 1
n!
y(n)(x0)(x− x0)n (1.4.3)
Aplicando-se a relação (1.4.4) para x = x0 + h, tem-se:
y(x0 + h) ≈ y(x0) + y(1)(x0)(h) + 1
2
y(2)(x0)(h)
2 + · · ·+ 1
n!
y(n)(x0)(h)
n (1.4.4)
A fórmula (1.4.4) fornece um caminho para determinar um valor aproximado para y(x0+h). Por exemplo, as seguintes expressões
são aproximações para o valor de y(x0 + h):
Aproximação de Taylor para y(x0 + h) por um:
y(x0 + h) ≈ y(x0) + y′(x0)h polinômio de grau 1
y(x0 + h) ≈ y(x0) + y′(x0)h+ 12y
′′
(x0)h
2 polinômio de grau 2
y(x0 + h) ≈ y(x0) + y′(x0)h+ 12y
′′
(x0)h
2 + 16y
′′′
(x0)h
3 polinômio de grau 3
...
...
y(x0 + h) ≈ y(x0) + y(1)(x0)(h) + 12y(2)(x0)(h)2 + · · ·+ 1n!y(n)(x0)(h)n polinômio de ordem n
Observa-se que para determinar um valor aproximado de y(x0 + h) usando-se o polinômio de Taylor de ordem n, precisa-se
conhecer o valor da função y e das derivadas desta função no ponto x0.
O que tem a ver o polinômio interpolador com o seguinte Problema de Valor Inicial ?
(∗)
{
y
′
(x) = f(x, y(x)),
y(x0) = y0
Resposta: Num PVI é conhecida a função f(x, y), isso implica que em um PVI conhece-se a derivada da função incógnita,
isto é:
y
′
(x) = f(x, y(x)) (1.4.5)
Logo, derivando sucessivamente a informação (1.4.5), podemos determinar as derivadas da função y(x). Com essas derivadas e
usando o polinômio de Taylor de ordem n em torno de x0 pode-se determinar uma aproximações para y(x0 + h).
x0 x1=x0+h x2=x1+h x3=x2+h
y0
y(x1+h)=?
y(x0+h)=?
y(x2+h)=?
x4
y(x4)
Figura 1.4.1:
• Seguindo o processo realizado acima para o cálculo aproximado de y(x0 + h), podemos obter um valor aproximado para
y(x0 + 2h). Neste caso, usam-se o polinômio de Taylor de ordem n em torno do ponto (x0 + h) e os valores das derivadas
da função y em tal ponto.
• Portanto, uma vez escolhido a ordem do polinomio de Taylor e o cálculo das expressões para as derivadas necessárias,
podem-se determinar aproximações para y nos pontos x = x0 + h, x0 + 2h, x0 + 3h, x0 + 4h, etc.
13
• O método de obtenção de aproximações para a função procurada usando o polinômio de Taylor de ordem n, chama-se o
método de Taylor de ordem n para resolver um PVI.
• Observe que o método de Taylor de ordem n é um método explícito, cujo cálculo do valor aproximado de y em x0+ h está
em função das informações de y e das derivadas de y em x0.
• A fórmula do método de Taylor de ordem n para resolver o PVI é dada por:
yk = yk−1 + hy
(1)
k−1 +
1
2!
y
(2)
k−1h
2 + . . .
1
n!
y
(n)
k−1h
n (1.4.6)
onde:
y
(j)
k−1 =
djy
dxj
(xk−1) =
d
dx
(
dj−1y
dxj−1
)
=
dj−1
dxj−1
(
dy
dx
)
dj−1
dxj−1
f(x, y(x))
∣∣∣∣
x=xk−1
j = 2, 3, 4, . . .
Como calcular as derivadas de ordem maior igual a 2 da função y de um PVI? Resposta: Usando a regra da cadeia.
Sabe-se que num PVI temos:
y
′
(x) = f(x, y(x)),
então:
y
′′
(x) =
d
dx
f(x, y(x)) = fx + fy
dy
dx
= fx + fyf,
onde: fx = ∂f∂x (x, y(x)) e fy =
∂f
∂y (x, y(x)).
E
y
′′′
=
d
dx
y
′′
=
d
dx
(fx + fyf) = fxx + 2fxyf + fyyf
2 + fxfy + (fy)
2f
Exemplo 25. Use o método de Taylor de segunda ordem para encontrar uma aproximação do seguinte PVI:
(∗)
{
y
′
= xe3x − 2y x ∈ [0, 1]
y(0) = −0.08
com tamanho de passo h= 0.25.
Solução:
Para h = 0.25, a Figura abaixo mostra os pontos onde a solução será aproximada.
0 0,25 0.5 0.75 1.0
y1
y2
y3
 
x1 x2 x3
|| || || || ||
y0
?
?
?
y4 ?
x0 x4
Da fórmula (1.4.6), o método de Taylor de segunda ordem é:
yk = yk−1 + hy
′
k−1 +
h2
2
y
′′
k−1 (1.4.7)
E do PVI dado acima, (*), tem-se f(x, y(x)) = xe3x − 2y(x). Então:
y
′
= xe3x − 2y
y
′′
=
d
dx
(
y
′)
=
d
dx
(
xe3x − 2y) = e3x + 3xe3x − 2dy
dx
14
y
′′
= e3x + 3xe3x − 2 (xe3x − 2y) = e3x + xe3x + 4y
Então,
y
′
k−1 = xk−1e
3xk−1 − 2yk−1 (1.4.8)
e
y
′′
k−1 = e
3xk−1 + xk−1e3xk−1 + 4yk−1 (1.4.9)
Substituindo (1.4.8) e (1.4.9) em (1.4.7), tem-se:
yk = yk−1 + h
(
xk−1e3xk−1 − 2yk−1
)
+
h2
2
(
e3xk−1 + xk−1e3xk−1 + 4yk−1
)
Os dados iniciais para determinar as aproximaçõesy1, y2, y3 e y4 para y(x1), y(x2), y(x3) e y(x4), respectivamente, são:
x0 = 0 e y0 = −0.08
Para k=1:
y1 = y0 + h(x0e
3x0 − 2y0) + h
2
2
(e3x0 + x0e
3x0 + 4y0)
y1 = −0.08 + 0.25(−2(−0.08)) + (0.25)
2
2
(1 + 4(−0.08))
y1 = −0.0187 ≈ y(x1) = y(0.25)
Para k=2:
y2 = y1 + h(x1e
3x1 − 2y1) + h
2
2
(e3x1 + x1e
3x1 + 4y1)
Substituindo os valores de x1 = 0.25 e y1 = −0.0187, tem-se:
y2 = 0.2035 ≈ y(x2) = y(0.5)
Para k=3:
y3 = y2 + h(x2e
3x2 − 2y2) + h
2
2
(e3x2 + x2e
3x2 + 4y2)
Substituindo os valores de x2 = 0.5 e y2 = 0.2035, tem-se:
y3 = 0.8974 ≈ y(x3) = y(0.75)
Para k=4:
y4 = y3 + h(x3e
3x3 − 2y3) + h
2
2
(e3x3 + x3e
3x3 + 4y3)
Substituindo os valores de x23 = 0.75 e y3 = 0.8974, tem-se:
y4 = 2, 8597 ≈ y(x4)= y(1.0)
A solução analítica do PVI acima:
y(x) =
1
5
xe3x − 1
25
e3x − 1
25
e−2x,
então:
k xk yk y(xk) e(xk) = y(xk)− yk
0 0 - 0.08 - 0.08 0.
1 0.25 -0.01875 - 0.0030912 0.0156588
2 0.5 0.2035 0.2541862 0.0508971
3 0.75 0.8974 1.0347257 0.1373798
4 1.0 2,8597 3.2082725 0.3496202
B.
O esboço do gráfico da solução exata e da solução aproximada são mostradas na figura abaixo:
15
Figura 1.4.2: Solução exata (curva sólida) e solução aproximada (pontos discretos).
1.4.2 Métodos de Runge-Kutta explícito
São métodos que se caracterizam pelas seguintes propriedades:
1. são de um passo um;
2. não exigem o cálculo de qualquer derivada de f(x, y); pagando, por isso, o preço de termos de cálcular f(x, y) em vários
pontos;
3. após expandir f(x, y) para funções de duas variáveis em torno de (xn, yn) usando a série de Taylor e agrupar os termos
semelhantes, suas expressões coincidem com a do método de Taylor de mesma ordem.
1.4.2.1 Método de Runge-Kutta de segunda ordem
Os métodos de Runge-Kutta procuram fórmulas de alta ordem que não precisem calcular derivadas, como é o caso das fórmulas
de Taylor.
Por exemplo, o método de Taylor de segunda ordem é:
yk = yk−1 + hy
′
k−1 +
h2
2
y
(2)
k−1 = yk−1 + hf(x, y(x))|xk−1 +
h2
2
df(x, y(x))
dx
|xk−1 ,
que precisa do cálculo da primeira e segunda derivada da função f com relação a x. Já o método explícito de Runge-Kutta de
segunda ordem é uma expressão da forma:
yk = yk−1 + h(ak1 + bk2),
onde a e b são constantes e as funções k1 e k2 são valores da função f(x, y) calculadas em pontos convenientes. A ordem desse
método deve ser igual à ordem do método de Taylor de segunda ordem, e para isso pede-se que:
h(ak1 + bk2) = hf(x, y(x))|xk−1 +
h2
2
df(x, y(x))
dx
|xk−1 + . . . . (1.4.10)
As expressões consideradas para k1 e k2, são, respectivamente:
k1 = f(xk−1, yk−1) (1.4.11)
k2 = f(xk−1 + ph, yk−1 + q h f(xk−1, yk−1)). (1.4.12)
Logo, para determinar a fórmula do método explícito de Runge-Kutta de segunda ordem precisamos conhecer os valores de a, b,
p e q.
Aplicando expansão em série de Taylor para funções de duas variáveis na função k2, temos:
k2 = f(xk−1 + ph, yk−1 + q h f(xk−1, yk−1)).
= f(xk−1, yk−1) + phfx(xk−1, yk−1) + qhfy(xk−1, yk−1) + . . . (1.4.13)
16
Logo, substituindo (1.4.13) no lado esquerdo da relação (1.4.10), temos:
h(ak1 + bk2) = haf + bh [f + phfx + qhfy] + . . .
= h(a+ b)f + h2b (pfx + qfy) + . . . (1.4.14)
Usando o fato que f
′
= fx + fyf , o lado direito da relação (1.4.10) fica:
hf +
h2
2
df
dx
+ · · · = hf + h
2
2
[fx + fyf ] + . . . . (1.4.15)
Igualando (1.4.14) e (1.4.15), temos:
h(a+ b)f + h2b (pfx + qfy) + · · · = hf + h
2
2
[fx + fyf ] + . . . (1.4.16)
Então:
a+ b = 1
bp =
1
2
bq =
1
2
Temos três equações e quatro incógnitas, portanto temos infinitas soluções. Para a = 0, temos que b = 1 e p = q = 1/2.
Já para a = 1/2, temos que b = 1/2 e p = q = 1. Escolhendo outros valores para a, tal que b 6= 0, temos outros métodos de
Runge-Kutta de segunda ordem.
Dois dos métodos de Runge-Kutta bastante usados, são:
1. Método de explícito de Runge-Kutta de segunda ordem (Método de Euler melhorado ou conhecido como o método de
Heun)
yk = yk−1 +
h
2
(k1 + k2) , k = 1, 2, 3, . . . (1.4.17)
k1 = f(xk−1, yk−1)
k2 = f(xk−1 + h, yk−1 + hk1)
2. Método explícito de Runge-Kutta de segunda ordem (Método de Euler modificado):
yk = yk−1 + hk2 (1.4.18)
k1 = k1(xk−1, yk−1)
k2 = f(xk−1 +
h
2
, yk−1 +
h
2
k1)
1.4.2.2 Método de Runge-Kutta de terceira ordem
São métodos do tipo:
yk = yk−1 + h(ak1 + bk2 + ck3), (1.4.19)
onde k1 e k2 são expressões da forma (1.4.11) e (1.4.13), respectivamente. E k3 possui a seguinte expressão:
k3 = f(xk−1 + ph, yk−1 + hqk1 + hrk2). (1.4.20)
Das expressões acima, para determinar o método de Runge-Kutta de terceira ordem precisamos conhecer o valor das seguintes
seis incógnitas: a. b, c, p, q e r. Lembrando que desejamos que o método resultante seja um método de terceira ordem.
Alguns dos métodos de Runge-Kutta de terceira ordem mais usados são:
1. Método de Heun de terceira ordem:
yk = yk−1 +
h
4
(k1 + 3k3)
k1 = f(xk−1, yk−1)
k2 = f((xk−1 +
1
3
h, yk−1 +
1
3
hk1)
k3 = f((xk−1 +
2
3
h, yk−1 +
2
3
hk2)
17
2. Método de Nystrom:
yk = yk−1 +
h
4
(k1 +
3
2
[k2 + k3])
k1 = f(xk−1, yk−1)
k2 = f((xk−1 +
2
3
h, yk−1 +
2
3
hk1)
k3 = f((xk−1 +
2
3
h, yk−1 +
2
3
hk2)
1.4.2.3 Método de Runge-Kutta de quarta ordem
yk = yk−1 +
h
6
(k1 + 2(k2 + k3) + k4) ,
k1 = f(xk−1, yk−1)
k2 = f(xk−1 +
h
2
, yk−1 +
h
2
k1)
k3 = f(xk−1 +
h
2
, yk−1 +
h
2
k2)
k4 = f(xk−1 + h, yk−1 + hk3)
1.5 Método de passo múltiplo
Podem ser obtidos por várias técnicas. Algumas dessas técnicas são baseadas na série de Taylor ou nas fórmulas de integração
numérica.
O PVI a ser resolvido é:
y
′
(x) = f(x, y(x)), x ∈ [a, b] (1.5.1)
y(a) = y0 (1.5.2)
1.5.1 Métodos numéricos para PVIs obtidos por Série de Taylor
Por série de Taylor, temos que:
y(xk) = y(xk−1) + hy
′
(xk−1) +
h
2!
y(2)(xk−1) + . . . ,
considerando os dois primeiros termos da série de Taylor e sabendo que y
′
(xk−1) = f(xk−1, y(xk−1)), temos que:
y(xk) ≈ y(xk−1) + hf(xk−1, y(xk−1)).
Daí, temos a seguinte fórmula discreta:
yk = yk−1 + hf(xk−1, yk−1).
Também, sabemos que:
y(xk+1) = y(xk) + hy
′
(xk) +
h
2!
y(2)(xk) +
h3
3!
y(3)(xk) + . . . (1.5.3)
y(xk−1) = y(xk)− hy′(xk) + h
2!
y(2)(xk) +
h3
3!
y(3)(xk) + . . . (1.5.4)
Fazendo (1.5.3)-(1.5.4), temos:
y(xk+1)− y(xk−1) = 2hy′(xk) + hy(2)(xk) + h
3
3
y
′′′
(xk) + . . . (1.5.5)
Considerando unicamento o primeiro termo do lado direito da expressão acima, temos:
y(xk+1)− y(xk−1) ≈ 2hy′(xk) (1.5.6)
ou:
y(xk+1) ≈ y(xk−1) + 2hy′(xk) (1.5.7)
18
Daí, temos a seguinte fórmula discreta:
yk+1 = yk−1 + 2hfk. (1.5.8)
Tal fórmula é explícita e de dois passos. Para calcular com essa fórmula y2 precisamos y0 e y1. Mas, não temos y1 e nem podemos
usar essa fórmula para determinar ele, logo precisamos de uma fórmula de um passo para poder calcular ele.
1.5.2 Métodos numéricos para PVIs obtidos por Integração Numérica
1.5.2.1 Integrando a EDO (1.5.1) em [xk−1, xk]
Sejam xk e xk−1 dois pontos onde a solução y será aproximada (Veja Figura abaixo).
x x
y(x )k−1
k
k−1 k
y(x )
Integrando a Eq. (1.5.1) no intervalo [xk−1, xk]., tem-se:
ˆ xk
xk−1
y
′
(x)dx =
ˆ xk
xk−1
f(x, y(x)) dx (1.5.9)
y(xk)− y(xk−1) =
ˆ xk
xk−1
f(x, y(x)) dx (1.5.10)
y(xk) = y(xk−1) +
ˆ xk
xk−1
f(x, y(x)) dx (1.5.11)
Podemos observar que o integrando da Eq. (1.5.9) é uma função de x, isto é, f(x, y(x)) = F (x), e a integral definida da
função F (x) em [xk−1, xk] pode ser aproximada de diferentes formas.
Por exemplo:
1. F (x) ≈ p0(x) = F (xk−1), então:
´ xk
xk−1
f(x, y(x)) dx = f(xk−1, y(xk−1))(xk−xk−1). Daí, temos: yk = yk−1+h f(xk−1, yk−1).
x xk−1 k
F(x )
F(x )
k−1
k
19
Tal fórmula é a fórmula do método explícito de Euler.
2. F (x) ≈ p0(x) = F (xk), então
´ xk
xk−1
f(x, y(x)) dx = f(xk, y(xk))(xk − xk−1). Logo, yk = yk−1 + h f(xk, yk).
x x
k
k
k−1F(x )
F(x )
k−1
Tal fórmula é a fórmula do método implícito de Euler.
3. F (x) ≈ p1(x) = a + bx, onde p1 é o polinômio interpolador de grau 1 de F que passa pelos pontos (xk−1, Fk−1) e (xk, Fk),
logo
´ xk
xk−1
f(x, y(x)) dx = f(xk,yk)+f(xk−1,yk−1)2 (xk − xk−1). Então, yk = yk−1 + h
(
f(xk−1,yk−1)+f(xk,yk)
2
)
.
x x
k
F(x )
F(x )
k−1
k
k−1
Tal fórmula é conhecida como o método do trapézio para PVI. Essa fórmula é implícita.
4. F (x) ≈ p1(x), onde p1(x) é o polinômio interpolador que passa pelospontos (xk−2, Fk−2) e (xk−1, Fk−1) e integrar esse
polinômio no intervalo [xk−1, xk], Fazendo isso, temos:
yk = yk−1 +
h
3
[3fk − 2fk−1] . (1.5.12)
Que é uma fórmula explícita.
1.5.2.2 Integrando a EDO (1.5.1) em [xk−2, xk]
O intervalo [xk−2, xk] contém dois sub-intervalos, cada um de comprimento h (para pontos igualmente espaçados). Integrando
a EDO, temos:
ˆ xk
xk−2
y
′
(x)dx =
ˆ xk
xk−2
f(x, y(x)) dx (1.5.13)
y(xk)− y(xk−2) =
ˆ xk
xk−2
f(x, y(x)) dx
ou:
y(xk) = y(xk−2) +
ˆ xk
xk−2
f(x, y(x)) dx (1.5.14)
Aplicando o método de integração 1/3 de Simpson na integral acima, temos:
y(xk) ≈ y(xk−2) + h
3
[fk−2 + 4fk−1 + fk] . (1.5.15)
Daí, temos a seguinte fórmula discreta:
yk = yk−2 +
h
3
[fk−2 + 4fk−1 + fk] . (1.5.16)
Esse método é um método implícito de dois passos, e portanto precisamos além da condição inicial y0 a informação y1. Essa
última informação deve ser obtida usando-se algum método de um passo.
20
1.5.2.3 Integrando a EDO (1.5.1) em [xk−3, xk]
Neste caso, o intervalo [xk−3, xk] contém três sub-intervalos, cada um de comprimento h. Integrando a EDO, temos:
ˆ xk
xk−3
y
′
(x)dx =
ˆ xk
xk−3
f(x, y(x)) dx (1.5.17)
y(xk)− y(xk−3) =
ˆ xk
xk−3
f(x, y(x)) dx
ou:
y(xk) = y(xk−3) +
ˆ xk
xk−3
f(x, y(x)) dx (1.5.18)
Aplicando o método de integração 3/8 de Simpson na integral acima, temos:
y(xk) ≈ y(xk−3) + 3h
8
[fk−3 + 3 (fk−2 + fk−1) + fk] (1.5.19)
Logo, temos a seguinte versão discreta:
yk = yk−3 +
3h
8
[fk−3 + 3 (fk−2 + fk−1) + fk] . (1.5.20)
Esse método é um método implícito de três passos. Para o cálculo de y3, precisamos y0, y1 e y2, logo devemos usar algum método
que permita calcular os valores y1 e y2.
21