Métodos Computacionais - Nilo Maurício

Prévia do material em texto

N
.M
.S
ot
om
ay
orME´TODOS COMPUTACIONAIS DA FI´SICA
Apostila
Araguaina 2014
N. M. Sotomayor
N
.M
.S
ot
om
ay
or
Suma´rio
Resumo 1
1 Integrac¸a˜o Nume´rica 2
1.1 Interpolac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Diferenc¸as Finitas . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Polinoˆmio Interpolador de Gregory-Newton . . . . . . 3
1.2 Integrac¸a˜o Nume´rica . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Exerc´ıcio . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Implementac¸a˜o computacional . . . . . . . . . . . . . . . . . . 8
1.3.1 Exemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.2 Exemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Soluc¸a˜o Nume´rica de Equac¸o˜es Diferenciais Ordina´rias EDOs 15
2.1 Equac¸o˜es diferenciais ordina´rias nume´ricas . . . . . . . . . . . 15
2.1.1 Equac¸a˜oes diferenciais ordina´rias . . . . . . . . . . . . 15
2.1.2 Problema de valor inicial . . . . . . . . . . . . . . . . . 16
2.2 Algoritmos nume´ricos para equac¸o˜es diferenciais ordina´rias . . 16
2.3 Me´todo de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Me´todo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4.1 Soluc¸a˜o Exa´ta do PVI do exemplo . . . . . . . . . . . 21
2.4.2 Soluc¸a˜o do PVI atrave´s do Me´todo de Euler . . . . . . 22
3 Os Me´todos de Runge-Kutta de primeira e segunda ordem 27
3.1 Definic¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.1 Me´todo de Runge-Kutta de primeira ordem s = 1 . . . 28
3.1.2 Me´todos de Runge-Kutta de Segunda Ordem . . . . . . 29
3.1.3 Me´todo de Heun (Karl Heun) . . . . . . . . . . . . . . 31
3.1.4 Me´todo de Runge-Kutta de segunda ordem (Ponto
Me´dio) . . . . . . . . . . . . . . . . . . . . . . . . . . 33
i
N
.M
.S
ot
om
ay
or
4 O Me´todo de Runge-Kutta de Terceira ordem 42
4.1 Resumo de teoria ba´sica . . . . . . . . . . . . . . . . . . . . . 42
5 O Me´todo de Runge-Kutta de Quarta Ordem 52
5.1 Teoria ba´sica . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.1.1 Interpretac¸a˜o geome´trica . . . . . . . . . . . . . . . . . 53
6 Integrac¸a˜o nume´rica das equac¸a˜oes de Newton 57
6.1 Forc¸a sobre um objeto em queda livre . . . . . . . . . . . . . . 57
6.1.1 Implementac¸a˜o computacional . . . . . . . . . . . . . . 58
6.2 Trajeto´rias em duas dimenso´es . . . . . . . . . . . . . . . . . . 61
7 Integrac¸a˜o nume´rica das equac¸a˜oes diferenciais de segunda
ordem 68
7.1 Me´todo de Numerov . . . . . . . . . . . . . . . . . . . . . . . 68
7.2 Derivac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
8 Integrac¸a˜o Nume´rica de Equac¸o˜es Diferenciais Parciais 72
8.1 Introduc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
8.2 Equac¸o˜es diferenciais parciais . . . . . . . . . . . . . . . . . . 72
8.3 Me´todo das diferenc¸as finitas . . . . . . . . . . . . . . . . . . 72
8.4 Soluc¸a˜o nume´rica da equac¸a˜o de Schro¨dinger dependente do
tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
8.5 Fundamentos . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
ii
N
.M
.S
ot
om
ay
or
Lista de Figuras
1.1 Considere o gra´fico de um conjunto de valores da grandeza y
em func¸a˜o da grandeza x. Se na˜o conhecemos o verdadeiro
valor funcional y = f(x) podemos aproximar a dependeˆncia
de y em func¸a˜o de x por um polinoˆmio interpolador Pn(x). . . 5
1.2 A curva em vermelho representa a forma funcional exata de
f(x) no intervalo [0.4 : 2.0], e´ necessa´rio calcular o valor da
integral nesse intervalo. Os pontos em preto representam nove
valores funcionais yi para correspondentes nove valores de xi,
no intervalo escolhido, esses valores foram calculados com o
co´digo computacional mostrado na Figura 1.4. . . . . . . . . . 10
1.3 Fluxograma empregado para descrever o algoritmo do ca´lculo
dos nove valores funcionais (xi, yi) da func¸a˜o f(x). . . . . . . 11
1.4 Co´digo computacional ou programa fonte em linguagem for-
tran 77 empregado para o o ca´lculo dos nove valores funcionais
(xi, yi) da func¸a˜o f(x). . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Tabela que mostra os nove valores funcionais (xi, yi) da func¸a˜o
f(x) obtidos com o co´digo mostrado na Figura 1.4. . . . . . . 12
1.6 Algoritmo em fluxograma para o ca´lculo nume´rico da a´rea do
exerc´ıcio proposto. . . . . . . . . . . . . . . . . . . . . . . . . 13
1.7 Ca´lculo da integral. . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1 Interpretac¸a˜o geome´trica do Me´todo de Euler. O esquema
mostra uma func¸a˜o y = f(x) com valores conhecidos: (x0, y0)
e dy/dx avaliada em x = x0. Considerando-se que a forma
funcional exa´ta e´ desconhecida, deseja-se encontrar uma apro-
ximac¸a˜o para o valor y(x1) com x1 muito pro´ximo de x0.
Empregando-se a definic¸a˜o da derivada de uma func¸a˜o pode se
obter uma aproximac¸a˜o para y1 sendo y1 = y0+htanθ, esta ex-
pressa˜o e´ a mesma dada pela fo´rmula de Euler y1 = y0+hdy/dx. 20
iii
N
.M
.S
ot
om
ay
or
2.2 Representac¸a˜o esquema´tica da aplicac¸a˜o recurssiva do Me´todo
de Euler para a soluc¸a˜o nume´rica de uma EDO. Dado o valor
da func¸a˜o e o da sua derivada no ponto (x0, y0) e´ poss´ıvel
fazer a previssa˜o do valor da func¸a˜o no ponto x1. Apo´s a
determinac¸a˜o de (x1, y1 e dy/dx em x1 podem ser empregados
estes valores para a determinac¸a˜o de y2 e assim sucessivamente. 21
2.3 Fluxograma do Algoritmo para a soluc¸a˜o exa´ta do PIV 2.37 . . 22
2.4 Programa fonte em FORTRAN 77 para o ca´lculo exa´to da
func¸a˜o 2.37. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Fluxograma do Algoritmo para o Me´todo de Euler . . . . . . . 24
2.6 Programa fonte em FORTRAN 77 que implementa o algoritmo
para o Me´todo de Euler. . . . . . . . . . . . . . . . . . . . . . 25
2.7 Gra´fico da soluc¸a˜o do PVI do exemplo, sao˜ comparadas a
soluc¸a˜o exata e as soluc¸o˜es aproximadas com diferentes valores
do passo h = (xf = x0)/n. A curva em vermelho representa
a soluc¸a˜o exa´ta com n = 20, a curva em verde representa a
soluc¸a˜o aproximada obtida pela aplicac¸a˜o do Me´todo de Euler
com n = 20, a curva em azul representa a soluc¸a˜o aproximada
obtida pela aplicac¸a˜o do Me´todo de Euler com n = 40 e a
curva em rouxo representa a soluc¸a˜o aproximada obtida pela
aplicac¸a˜o do Me´todo de Euler com n = 80. Observa-se que
quanto menor o passo a soluc¸a˜o nume´rica aproxima-se mais
da soluc¸a˜o exa´ta. . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1 Esquema que mostra a interpretac¸a˜o geome´trica do Me´todo de
Heun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2 Fluxograma que explica a implementac¸a˜o computacional do
Me´todo de Heun . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3 Programa fonte em Fortran 77, que implementa o Me´todo de
integrac¸a˜o nume´rica de Heun. . . . . . . . . . . . . . . . . . . 36
3.4 Comparac¸a˜o dos ca´lculos da soluc¸a˜o exa´ta do PIV e dos Me´todos
de Euler e Heun (RK-2). . . . . . . . . . . . . . . . . . . . . . 37
3.5 Comparac¸a˜o dos valores nume´ricos da soluc¸a˜o exa´ta do PIV
e dos Me´todos de Euler e Heun (RK-2). . . . . . . . . . . . . 38
3.6 Programa fonte em Fortran 77, que implementa o Me´todo de
integrac¸a˜o nume´rica RK-2 Ponto Me´dio. . . . . . . . . . . . . 39
3.7 Interpretac¸a˜o geome´trica do Me´todo de integrac¸a˜o nume´rica
RK-2 Ponto Me´dio. . . . . . . . . . . . . . . . . . . . . . . . . 40
3.8 Comparac¸a˜o dos valores nume´ricos da soluc¸a˜o exa´ta do PIV
e dos Me´todos de Euler, Heun e Midpoint (RK-2). . . . . . . . 41
iv
N
.M
.S
ot
om
ay
or
4.1 Programa fonte em Fortran 77, que implementa o Me´todo de
integrac¸a˜o nume´rica de Runge-Kutta de ordem 3 RK-3. . .. . 49
4.2 Comparac¸a˜o de trajeto´rias da soluc¸a˜o exa´ta e do me´todo de
Runge-Kutta de ordem 3 para o PIV em estudo. . . . . . . . . 50
4.3 Comparac¸a˜o dos valores nume´ricos das trajeto´rias da soluc¸a˜o
exa´ta e do me´todo de Runge-Kutta de ordem 3 para o PIV em
estudo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.1 Co´digo fortran 77 mostrando a implementac¸a˜o do Me´todo de
integrac¸a˜o nume´rica de Runge-Kutta quarta Ordem . . . . . . 54
5.2 k1 e´ a inclinac¸a˜o em (x0, y0). Avanc¸a-se meio passo e calcula-
se yM1 com k1 atrave´s do me´todo de Euler. . . . . . . . . . . 55
5.3 k2 e´ a inclinac¸a˜o em (xM1, yM1). Volta-se para x0 e calcula-se
yM2, com k2. Usando-se (xM , yM2) calcula-se k3. . . . . . . . 55
5.4 Usando-se k3 e o me´todo de Euler encontra-se y1T . Em se-
guida, empregando-se (x1, y1T ) encontras-se k4. . . . . . . . . 56
5.5 Finalmente, volta-se para x0 e calcula-se a aproximac¸a˜o para
y1 empregando-se a inclina˜c¸a˜o me´dia, usando o me´todo de Eu-
ler. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.1 Co´digo fortran 77 mostrando a implementac¸a˜o do Me´todo de
Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.2 Queda de um corpo de massa m sem resisteˆncia do a´r. . . . . 60
6.3 Co´digo fortran 77 mostrando a implementac¸a˜o do Me´todo de
Euler-Richardson . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.4 Queda de um corpo de massa m com resisteˆncia do a´r propor-
cional a v2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.5 Movimento de um proje´til em 2D. . . . . . . . . . . . . . . . . 63
6.6 Programa fonte para o ca´lculo das trajeto´rias de um proje´til
em 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
6.7 Gra´fico das trajeto´rias nume´ricas de um proje´til em 2D lanc¸ado
do n´ıvel do solo, comaˆngulo de 45o e para diferentes valores do
paraˆmetro k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7.1 Implementac¸a˜o computacional do Me´todo de Numerov. . . . . 69
8.1 Programa fonte em fortran 77 para o ca´lculo nume´rico da evoluc¸a˜o
temporal de um pacote de onda gaussiano atrave´s da equac¸a˜o de
Schro¨dinger empregando-se diferenc¸as finitas e o me´todo de passo
me´dio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.2 A subrotina INITIALDATA e´ invocada pelo programa principal
para gerar os dados inicias que sera˜o empregados em todo o processo. 76
v
N
.M
.S
ot
om
ay
or
8.3 A subrotina CLEANMATRIX e´ invocada pelo programa principal
para zerar inicialmente os vetores a serem empregados no programa. 76
8.4 A subrotina WAVEPACKET e´ invocada pelo programa prin-
cipal para construir o pacote de onda gaussiano inicial (em
t = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.5 A subrotina EVOLVE e´ invocada pelo programa principal para
integrar nume´ricamente as equac¸o˜es diferenciais para as partes
real e imagina´ria do pacote de onda. . . . . . . . . . . . . . . 77
8.6 A func¸ao˜ V e´ empregada para definir o potencial. . . . . . . . 78
8.7 Gra´fico das partes real, imagina´ria e densidade de probabi-
lidade do pacote de onda gaussiano evoluindo em func¸a˜o do
tempo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
vi
N
.M
.S
ot
om
ay
or
Resumo
Sa˜o apresentados os me´todos ba´sicos para integrac¸a˜o nume´rica de func¸o˜es,
equac¸o˜es diferenciais ordina´rias (problemas de valor inicial), equac¸o˜es dife-
renciais de segunda ordem e, equac¸o˜es diferenciais parciais nume´ricas.
1
1Palavras chave: Me´todos Nume´ricos, Equac¸o˜es Diferenciais Ordina´rias, Linguagens e
Aplicativos Computacionais.
1
N
.M
.S
ot
om
ay
or
Cap´ıtulo 1
Integrac¸a˜o Nume´rica
1.1 Interpolac¸a˜o
Denomina-se interpolac¸a˜o o me´todo que permite construir um novo con-
junto de dados a partir de um conjunto discreto de dados pontuais previa-
mente conhecidos. Atrave´s da interpolac¸a˜o, pode-se construir uma func¸a˜o
que se ajuste aos dados pontuais existentes permitindo assim fazer previso˜es
em todo o intervalo do domı´nio. Dados discretos podem ser interpolados por
polinoˆmios um exemplos deste tipo de func¸o˜es sa˜o o polinoˆmio de Newton e
o polinoˆmio de Gregory-Newton.
1.1.1 Diferenc¸as Finitas
Uma diferenc¸a finita e´ uma expresa˜o matema´tica da forma:
f(x+ b)− f(x+ a). (1.1)
Uma diferenc¸a finita progressiva e´ definida por:
∆f(x) = f(x+ h)− f(x), (1.2)
Segundo a ordem as principais diferenc¸as ordina´rias ou finitas sa˜o:
∆0f(x) = f(x) (1.3)
∆1f(x) = f(x+ h)− f(x)
∆2f(x) = ∆1f(x+ h)−∆1f(x)
∆3f(x) = ∆2f(x+ h)−∆2f(x)
−−−−−−
2
N
.M
.S
ot
om
ay
or
∆0f(x) = f(x) = y0 (1.4)
∆1f(x) = f(x+ h)− f(x)
= y1 − y0
∆2f(x) = ∆1f(x+ h)−∆1f(x)
= [f(x+ h+ h)− f(x+ h)]− [f(x+ h)− f(x)]
= f(x+ 2h)− 2f(x+ h) + f(x)
= y2 − 2y1 + y0
∆3f(x) = ∆2f(x+ h)−∆2f(x)
= [∆1f(x+ 2h)−∆1f(x+ h)]− [∆1f(x+ h)−∆1f(x)]
= ∆1f(x+ 2h)− 2∆1f(x+ h) + ∆1f(x)
= [f(x+ 3h)− f(x+ 2h)]− 2[f(x+ 2h)− f(x+ h)]
+ [f(x+ h)− f(x)]
= f(x+ 3h)− 3f(x+ 2h) + 3f(x+ h)− f(x)
= y3 − 3y2 + 3y1 − y0
1.1.2 Polinoˆmio Interpolador de Gregory-Newton
Seja o polinoˆmio Pn(x),
Pn(x) = y0+N1y0(x−x0)+N2y0(x−x0)(x−x1)+...+Nny0(x−x0)(x−x1)...(x−xn−1)
(1.5)
onde,
• N1y0: Operador diferenc¸a dividida de ordem n;
• Nny0 = ∆
ny0
n!hn
;
• ∆ny0: Operador diferenc¸a finita.
Pn(x) = y0+
∆1y0
1!h1
(x−x0)+∆
2y0
2!h2
(x−x0)(x−x1)+...+∆
ny0
n!hn
(x−x0)(x−x1)...(x−xn−1)
(1.6)
Fazemos as seguintes aproximac¸o˜es:
P1(x) = y0 +
∆1y0
h
(x− x0) (1.7)
3
N
.M
.S
ot
om
ay
or
P2(x) = y0 +
∆1y0
h
(x− x0) + ∆
2y0
2h2
(x− x0)(x− x1) (1.8)
P3(x) = y0+
∆1y0
h
(x−x0)+∆
2y0
2h2
(x−x0)(x−x1)+∆
3y0
6h3
(x−x0)(x−x1)(x−x2).
(1.9)
Fazemos a mudanc¸a de varia´vel,
x− x0 = hu (1.10)
x− x1 = h(u− 1)
x− x2 = h(u− 2)
x− x3 = h(u− 3)
Com as novas varia´veis as aproximac¸o˜es ficam na forma,
P1(x) = y0 +∆
1y0u (1.11)
P2(x) = y0 +∆
1y0u+
∆2y0
2
u(u− 1) (1.12)
P3(x) = y0 +∆
1y0u+
∆2y0
2
u(u− 1) + ∆
3y0
6
u(u− 1)(u− 2). (1.13)
1.2 Integrac¸a˜o Nume´rica
Seja uma func¸a˜o f(x) integra´vel no intervalo [a, b], enta˜o∫ b
a
f(x) dx = F (b)− F (a) (1.14)
onde,
F
′
(x) = f(x). (1.15)
Quando a forma anal´ıtica de F (x) for de dif´ıcil obtenc¸a˜o ou se forem
conhecidos somente valores discretos de f(x), se faz necessa´rio o emprego de
me´todos nume´ricos para avaliar a integral de f(x). Esses me´todos consistem
em aproximar a func¸a˜o f(x) por um polinoˆmio interpolador e determinar
4
N
.M
.S
ot
om
ay
or
x a
0
=
y
0
x
1
y
1
x
2
y
2
x
3
y
3
x
4
y
4
x
5
x
y
5
x
6
y
6
x b
n
=
y
n
y
p(x)nf(x)
Figura 1.1: Considere o gra´fico de um conjunto de valores da grandeza y
em func¸a˜o da grandeza x. Se na˜o conhecemos o verdadeiro valor funcional
y = f(x) podemos aproximar a dependeˆncia de y em func¸a˜o de x por um
polinoˆmio interpolador Pn(x).
anal´ıticamente a integral desse polinoˆmio no intervalo [a, b]. Empregaremos
para o propo´sito os polinoˆmios de Gregory-Newton.
Inicialmente, vamos aproximar a func¸a˜o inco´gnita por um polinoˆmio inter-
polador de grau n = 1 de forma que,∫ b
a
f(x) dx ≈
∫ b=x1
a=x0
P1(x)dx (1.16)
∫ b=x1
a=x0
P1(x)dx =
∫ b=x1
a=x0
(
y0 +∆
1y0u
)
dx (1.17)
lembrando que,
x− x0 = hu, (1.18)
enta˜o,
dx = hdu (1.19)
fazemos a mudanc¸a de varia´vel nos limites de integrac¸a˜o. Assim si,
5
N
.M
.S
ot
om
ay
or
x −→ x0 ⇒ u→ 0, (1.20)
x −→ x1 ⇒ u→ 1 (1.21)
de forma que,∫ b=x1
a=x0
(
y0 +∆
1y0u
)
dx =
∫ 1
0
(
y0 +∆
1y0u
)
hdu (1.22)
=
∫ 1
0
y0hdu+
∫ 1
0
∆1y0hudu (1.23)= hy0u
∣∣∣∣1
0
+ h∆1y0
u2
2
∣∣∣∣1
0
(1.24)
= hy0 + h∆
1y0
1
2
= hy0 +
h
2
(y1 − y0) (1.25)
finalmente,
∫ b=x1
a=x0
P1(x)dx =
h
2
(
y0 + y1
)
(1.26)
conhecida como a Regra do Trape´zio.
A seguir vamos aproximar a func¸a˜o inco´gnita por um polinoˆmio interpolador
de grau n = 2, seguindo o procedimento anterior:∫ b
a
f(x) dx ≈
∫ b=x2
a=x0
P1(x)dx (1.27)
∫ b=x2
a=x0
P2(x)dx =
∫ b=x2
a=x0
[
y0 +∆
1y0u+
∆2y0
2
u(u− 1)
]
dx (1.28)
lembrando que,
x− x2 = h(u− 2), (1.29)
enta˜o,
dx = hdu (1.30)
assim,
6
N
.M
.S
ot
om
ay
or
∫ b=x2
a=x0
P2(x)dx =
∫ b=x2
a=x0
[
y0 +∆
1y0u+
∆2y0
2
u(u− 1)
]
hdu (1.31)
fazemos a mudanc¸a de varia´vel nos limites de integrac¸a˜o,
x− x0 = hu (1.32)
x− x2 = h(u− 2) (1.33)
Assim si,
x −→ x0 ⇒ u→ 0, (1.34)
x −→ x2 ⇒ u→ 2 (1.35)
de forma que,
∫ b=x2
a=x0
P2(x)dx =
∫ 2
0
[
y0 +∆
1y0u+
∆2y0
2
u(u− 1)
]
hdu (1.36)
integrando,
∫ 2
0
[
y0+∆
1y0u+
∆2y0
2
u(u−1)
]
hdu = h
∫ 2
0
[
y0+∆
1y0u+
1
2
∆2y0u
2−1
2
∆2y0u
]
du
(1.37)
h
∫ 2
0
[
y0+∆
1y0u+
1
2
∆2y0u
2−1
2
∆2y0u
]
du = h
[
y0u+∆
1y0
u2
2
+∆2y0
u3
6
−∆2y0u
2
4
]∣∣∣∣∣
2
0
(1.38)
fazendo a avaliac¸a˜o obtemos,
∫ b=x2
a=x0
P2(x)dx =
h
3
[
y0 + 4y1 + y2
]
(1.39)
Expressa˜o conhecida como a Regra 1/3 de Simpson.
7
N
.M
.S
ot
om
ay
or
1.2.1 Exerc´ıcio
Aproximar a func¸a˜o inco´gnita f(x) por um polinoˆmio interpolador de grau
n = 3 e encontrar a expressa˜o para
∫ b
a
f(x)dx.
Resposta
∫ b=x3
a=x0
P3(x)dx =
3
8
h
[
y0 + 3y1 + 3y2 + y3
]
(1.40)
Expressa˜o conhecida como a Regra 3/8 de Simpson.
1.3 Implementac¸a˜o computacional
1.3.1 Exemplo 1
Seja a func¸a˜o f(x) = 1/x utilize o me´todo do trape´zio para aproximar:
I =
∫ 6
1
1
x
dx. (1.41)
Segundo o me´todo do trape´zio a a´rea sob a curva no intervalo de x = 1 a
x = 6 pode ser aproximada pela expressa˜o:
I =
h
2
(y0 + y1) (1.42)
Empregamos h = 5 (distaˆncia entre ponto final e ponto inicial h = 6 − 1),
como passo de integrac¸a˜o, consideramos tambe´m y0 = 1/x0 = 1/1 e y1 =
1/x1 = 1/6, de forma que,
I =
5
2
(
1
1
+
1
6
)
=
35
12
= 2, 916 (1.43)
Devido ao erro de se aproximar uma func¸a˜o por um polinoˆmio de grau 1 e´
que se faz necessa´rio dividir o intervalo de integrac¸a˜o em va´rios subintervalos
menores para aplicar repetitivamente o me´todo do trape´zio a cada um destes
subintervalos aumentando considera´velmente a precisa˜o do resultado.
Consideremos,
I =
∫ b=xm
a=x0
f(x)dx (1.44)
8
N
.M
.S
ot
om
ay
or
Divida-se o intervalo [a, b] em m subintervalos iguais, em seguida aplica-se o
me´todo do trape´zio a cada par de pontos consecutivos da abscissa, de forma
que,
I =
∫ b=xm
a=x0
f(x)dx ≈ h
2
(y0 + y1) +
h
2
(y1 + y2) + ...+
h
2
(ym−1 + ym) (1.45)
I =
∫ b=xm
a=x0
f(x)dx ≈ h
2
(
y0 + 2y1 + 2y2 + ...+ 2ym−1 + ym
)
(1.46)
I =
∫ b=xm
a=x0
f(x)dx ≈ h
2
m∑
i=0
ciyi (1.47)
com c0 = cm = 1, e i = 1, 2, ...,m− 1.
1.3.2 Exemplo 2
Seja a func¸a˜o y = exp(x) sin(10x)+ 8, encontre o valor da integral da func¸a˜o
no intervalo de x = 0.4 a x = 2.0 pelo me´todo exa´to (por integrac¸a˜o direta),
e empregando a aproximac¸a˜o do me´todo do trape´zio sucessivas vezes.
Ca´lculo exa´to
I =
∫ 2.0
0.4
[exp(x) sin(10x) + 8]dx (1.48)
O intervalo sera´ dividido em m = 08 partic¸o˜es iguais, para calcular a integral
de f(x) no intervalo desejado sera´ aplicado o me´todo do trape´zio a cada
partic¸a˜o, assim a integral desejada sera´ a soma das a´reas dos oito trape´zios.
A figura 1.5 mostra os nove valores de x e y para a func¸a˜o f(x) obtidos pelo
co´digo computacional em fortran 77.
A seguir e´ apresentado um fluxograma que ilustra o algoritmo para o ca´lculo
da a´rea do exemplo proposto.
A seguir e´ apresentado na figura 1.7 o programa fonte redigido em fortran
77 para o ca´lculo aproximado da integral pedida. O programa usa como
dados de entrada os nove valores funcionais mostrados na tabela da Figura
1.5 esses dados esta˜o em um arquivo de texto denominado dados.dat. Apo´s
leitura dos dados nume´ricos o programa calcula a a´rea dos oito trape´zios e
as soma sequencialmente armazenando esses valores em uma varia´vel do tipo
acumulador (chamada integral no programa).
9
N
.M
.S
ot
om
ay
or
2
4
6
8
10
12
14
16
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
y x x=exp( )*sin(10* )+8
y0
y1
y2
y3
y4
y5
y6
y7
y8
x
y
Figura 1.2: A curva em vermelho representa a forma funcional exata de
f(x) no intervalo [0.4 : 2.0], e´ necessa´rio calcular o valor da integral nesse
intervalo. Os pontos em preto representam nove valores funcionais yi para
correspondentes nove valores de xi, no intervalo escolhido, esses valores fo-
ram calculados com o co´digo computacional mostrado na Figura 1.4.
10
N
.M
.S
ot
om
ay
or
x : xF
<
>=
Início
Fim
Ler
x , y , x , m0 0 f
h=(x -x )/mf 0
x=x , y=y0 0
i=1
y=8+exp(x)sin(10x)
x = x + i h0
i=i+1
print
x, y
Figura 1.3: Fluxograma empregado para descrever o algoritmo do ca´lculo dos
nove valores funcionais (xi, yi) da func¸a˜o f(x).
11
N
.M
.S
ot
om
ay
or
Figura 1.4: Co´digo computacional ou programa fonte em linguagem fortran
77 empregado para o o ca´lculo dos nove valores funcionais (xi, yi) da func¸a˜o
f(x).
m x y
0 0.40000000000000002 6.8709833462630909
1 0.60000000000000009 7.4908717676112584
2 0.80000000000000004 10.201857270801881
3 1.0000000000000000 6.5211972999713232
4 1.2000000000000002 6.2185151746646303
5 1.4000000000000001 12.017110915969930
6 1.6000000000000001 6.5740055374670368
7 1.8000000000000000 3.4567919067612700
8 2.0000000000000000 14.745803672878749
Figura 1.5: Tabela que mostra os nove valores funcionais (xi, yi) da func¸a˜o
f(x) obtidos com o co´digo mostrado na Figura 1.4.
12
N
.M
.S
ot
om
ay
or
i : m-1
<
>=
Fim
i=1
Início
Ler
«dados.dat»
m
x = x(1)0
x = x(9)f
h = (x - x )/(m-1)f 0
integral=0.0
integral=integral+Ai
i=i+1
Figura 1.6: Algoritmo em fluxograma para o ca´lculo nume´rico da a´rea do
exerc´ıcio proposto.
13
N
.M
.S
ot
om
ay
or
Figura 1.7: Ca´lculo da integral.
14
N
.M
.S
ot
om
ay
or
Cap´ıtulo 2
Soluc¸a˜o Nume´rica de Equac¸o˜es
Diferenciais Ordina´rias EDOs
2.1 Equac¸o˜es diferenciais ordina´rias nume´ricas
Equac¸o˜es diferenciais ordina´rias nume´ricas corresponde a` parte da ana´lise
nume´rica que estuda a soluc¸a˜o nume´rica de equac¸o˜es diferenciais ordina´rias
(ODEs). E´ conhecido que muitas equac¸o˜es diferenciais na˜o podem ser resol-
vidas de forma anal´ıtica, neste caso e´ necessa´rio a procura por uma apro-
ximac¸a˜o da soluc¸a˜o exata. A seguir sera˜o apresentados alguns algoritmos
empregados para obter essas soluc¸o˜es aproximadas. Sa˜o exemplos deste al-
goritmos os me´todos de Taylor, Euler, Runge-Kutta, etc.
2.1.1 Equac¸a˜oes diferenciais ordina´rias
Uma equac¸a˜o diferencial ordina´ria (frequentemente chamada de ODE ou
EDO) e´ uma relac¸a˜o que envolve uma varia´vel independente, uma func¸a˜o
inco´gnita e as derivadas desta func¸a˜o desconhecida. Uma equac¸a˜o diferencial
ordina´ria de ordem n e´ uma relac¸a˜o do tipo,
F (x; y, y
′
, ..., yn) = 0 (2.1)
onde,
• y e´ uma func¸a˜o de x;
• y′ = dy
dx
e´ a primeira derivada em relac¸a˜o a x;
• yn = d
ny
dxn
e´ a n−e´sima derivada em relac¸a˜o a x [1].
15
N
.M
.S
ot
om
ay
or
Em geral, uma ODE de n−e´sima ordempossui n soluc¸o˜es linearmente
independentes. Mais ainda, qualquer combinac¸a˜o linear de soluc¸o˜es linear-
mente independentes e´ tambe´m uma soluc¸a˜o.
2.1.2 Problema de valor inicial
Denomina-se problema de valor inicial (PVI) a uma equac¸a˜o diferencial
ordina´ria associada a um valor espec´ıfico da func¸a˜o inco´gnita em um dado
ponto do domı´nio da soluc¸a˜o (condic¸o˜es iniciais).
2.2 Algoritmos nume´ricos para equac¸o˜es di-
ferenciais ordina´rias
2.3 Me´todo de Taylor
E´ poss´ıvel de aproximar o valor de uma func¸a˜o f(x), entorno de um certo
ponto particular x0, se e´ conhecido, tanto o valor da func¸a˜o f(x0), quanto os
valores das suas derivadas, nesse ponto. A aproximac¸a˜o e´ conhecida como
expansa˜o de Taylor da func¸a˜o entorno do ponto x0.
f(x) = f(x0)+
(x− x0)
1!
df
dx
∣∣∣∣
x=x0
+
(x− x0)2
2!
d2f
dx2
∣∣∣∣
x=x0
+...+
(x− x0)n
n!
dnf
dxn
∣∣∣∣
x=x0
(2.2)
considerando x− x0 = h a equac¸a˜o anterior pode se escrever:
y = f(x0 + h) = f(x0) +
h
1!
df
dx
∣∣∣∣
x=x0
+
h2
2!
d2f
dx2
∣∣∣∣
x=x0
+ ...+
hn
n!
dnf
dxn
∣∣∣∣
x=x0
(2.3)
ou ainda,
y = y0 +
h
1!
df
dx
∣∣∣∣
x=x0
+
h2
2!
d2f
dx2
∣∣∣∣
x=x0
+ ...+
hn
n!
dnf
dxn
∣∣∣∣
x=x0
(2.4)
onde h e´ o valor do domı´nio para o qual sera´ encontrada a aproximac¸a˜o da
func¸a˜o.
Para ilustrar a obtenc¸a˜o da aproximac¸a˜o atrave´s da expansa˜o de Taylor con-
sideremos apenas os treˆs primeiros termos da se´rie:
y = f(x0 + h) = f(x0) +
h
1!
df
dx
∣∣∣∣
x=x0
+
h2
2!
d2f
dx2
∣∣∣∣
x=x0
(2.5)
16
N
.M
.S
ot
om
ay
or
a derivada
d2f
dx2
pode-se obter sem dificuldade derivando
df
dx
mais uma vez.
Exemplo
Consideramos o seguinte problema de valor inicial:
dy
dx
= 2xy E. D.
y(x = 1) = 1 C. I.
}
x0 = 1, y0 = 1 (2.6)
Aplicamos o me´todo de Taylor para obter uma aproximac¸a˜o de y para x =
1, 1.
Para isto, tomamos h = 0, 1.
dy
dx
= 2xy, (2.7)
d2y
dx2
= 2x
dy
dx
+ 2y, (2.8)
quando h = 0, 1
dy
dx
∣∣∣∣
x=x0
= 2x0y0 = 2(1)(1) = 2, (2.9)
como
d2y
dx2
∣∣∣∣
x=x0
= 2x0
dy
dx
∣∣∣∣
x=x0
+ 2y0 (2.10)
d2y
dx2
∣∣∣∣
x=x0
= 2(1)(2) + 2(1) = 6 (2.11)
Assim, a equac¸a˜o 2.14 podemos escrever na forma
y1 = y0 + h
dy
dx
∣∣∣∣
x=x0
+
h2
2
d2y
dx2
∣∣∣∣
x=x0
(2.12)
y1 = (1) + (0, 1)(2) +
(0, 1)2
2
(6) (2.13)
de onde y1 = 1, 23.
17
N
.M
.S
ot
om
ay
or
2.4 Me´todo de Euler
O me´todo de Euler e´ o me´todo nume´rico mais elementar de resoluc¸a˜o de
equac¸o˜es diferenciais ordina´rias. O me´todo de Euler tambe´m e´ conhecido
como me´todo do Trape´zio ou me´todo de Heun.
Seja uma expansa˜o da soluc¸a˜o exata y(x), em se´rie de Taylor, em torno do
valor inicial x0
y = f(x0 + h) = f(x0) +
h
1!
df
dx
∣∣∣∣
x=x0
+
h2
2!
d2f
dx2
∣∣∣∣
x=x0
+
h3
3!
d3f
dx3
∣∣∣∣
x=x0
+ ... (2.14)
Truncando a se´rie apo´s o termo de derivada primeira, sendo x1 = x0+h e y1
uma aproximac¸a˜o de f(x1), tem-se:
f(x0 + h) = f(x0) + h
dy
dx
∣∣∣∣
x=x0
(2.15)
y1 = y0 + hf(x0, y0) (2.16)
onde,
f(x0, y0) =
dy
dx
∣∣∣∣
x=x0
(2.17)
Para fazer a previsa˜o de y2 = f(x2), com x2 = x1 + h,
f(x1 + h) = f(x1) + h
dy
dx
∣∣∣∣
x=x1
(2.18)
y2 = y1 + hf(x1, y1) (2.19)
onde,
f(x1, y1) =
dy
dx
∣∣∣∣
x=x1
(2.20)
Da mesma forma,
y3 = y2 + hf(x2, y2) (2.21)
onde,
f(x2, y2) =
dy
dx
∣∣∣∣
x=x2
(2.22)
18
N
.M
.S
ot
om
ay
or
As sucessivas aproximac¸o˜es para y podem, enta˜o, ser obtidas pela fo´rmula
de recorreˆncia.
yi+1 = yi + hf(xi, yi) (2.23)
com,
f(xi, yi) =
dy
dx
∣∣∣∣
x=xi
(2.24)
que e´ conhecido como Me´todo de Euler cuja interpretac¸a˜o geome´trica e´ mos-
trada na figura 2.1. O esquema mostra uma func¸a˜o y = f(x) com valores
conhecidos: (x0, y0) e dy/dx avaliada em x = x0. Considerando-se que a
forma funcional exa´ta e´ desconhecida, deseja-se encontrar uma aproximac¸a˜o
para o valor y(x1) com x1 muito pro´ximo de x0. Empregando-se a definic¸a˜o
da derivada de uma func¸a˜o pode se obter uma aproximac¸a˜o para y1 sendo
y1 = y0 + htanθ, esta expressa˜o e´ a mesma dada pela fo´rmula de Euler
y1 = y0+ hdy/dx. A equac¸a˜o 2.15 e´ usada para encontrar uma aproximac¸a˜o
de y para x = x0+ h, para encontrar aproximac¸o˜es consecutivas e´ necessa´rio
repetir a aplicac¸a˜o da fo´rmula de Euler como mostrado na figura 2.2.
Exemplo
A seguir vamos analisar a soluc¸a˜o exa´ta e a soluc¸a˜o nume´rica do problema
de valor inicial (PVI) seguinte:
dy
dx
= x− 2y + 1 (2.25)
Primeiramente, mostramos a soluc¸a˜o exata da equac¸a˜o diferencial,
dy
dx
+ 2y = x+ 1 (2.26)
e2x
[dy
dx
+ 2y] = e2x(x+ 1) (2.27)
d
dx
[
e2xy] = e2x(x+ 1) (2.28)
d
[
e2xy] = e2x(x+ 1)dx (2.29)
e2xy =
∫
(x+ 1)e2xdx (2.30)
19
N
.M
.S
ot
om
ay
or
Figura 2.1: Interpretac¸a˜o geome´trica do Me´todo de Euler. O esquema mostra
uma func¸a˜o y = f(x) com valores conhecidos: (x0, y0) e dy/dx avaliada em
x = x0. Considerando-se que a forma funcional exa´ta e´ desconhecida, deseja-
se encontrar uma aproximac¸a˜o para o valor y(x1) com x1 muito pro´ximo de
x0. Empregando-se a definic¸a˜o da derivada de uma func¸a˜o pode se obter uma
aproximac¸a˜o para y1 sendo y1 = y0 + htanθ, esta expressa˜o e´ a mesma dada
pela fo´rmula de Euler y1 = y0 + hdy/dx.
Integrando por partes o lado direito,
e2xy =
1
2
(x+ 1)e2x − 1
2
∫
e2xdx (2.31)
e2xy =
1
2
(x+ 1)e2x − 1
4
e2x + C (2.32)
y =
1
2
(x+ 1)− 1
4
+ C e−2x (2.33)
20
N
.M
.S
ot
om
ay
or
Figura 2.2: Representac¸a˜o esquema´tica da aplicac¸a˜o recurssiva do Me´todo
de Euler para a soluc¸a˜o nume´rica de uma EDO. Dado o valor da func¸a˜o e
o da sua derivada no ponto (x0, y0) e´ poss´ıvel fazer a previssa˜o do valor da
func¸a˜o no ponto x1. Apo´s a determinac¸a˜o de (x1, y1 e dy/dx em x1 podem ser
empregados estes valores para a determinac¸a˜o de y2 e assim sucessivamente.
y =
1
2
x+
1
4
+ C e−2x (2.34)
As condic¸o˜es iniciais po PVI sa˜o:
(x0, y0) = (0, 1) (2.35)
Substitu´ındo esses valores na EDO 2.34 obtemos,
C =
3
4
(2.36)
de forma que,
y(x) =
1
2
x+
1
4
+
3
4
e−2x (2.37)
e´ a func¸a˜o ingo´gnita que satisfaz o problema de valor inicial 2.25.
2.4.1 Soluc¸a˜o Exa´ta do PVI do exemplo
A seguir apresentamos o algoritmo em fluxograma e o co´digo computacional
em fortran 77 para a avaliac¸a˜o da soluc¸a˜o exata do problema de valor incial
21
N
.M
.S
ot
om
ay
or
proposto. A figura 2.3 mostra um fluxograma para a representac¸a˜o da soluc¸a˜o
exa´ta do problema 2.37.
Figura 2.3: Fluxograma do Algoritmo para a soluc¸a˜o exa´ta do PIV 2.37
A implementac¸a˜o computacional do algoritmo 2.3 e´ mostrado na figura . O
co´digo fonte foi redigido em linguagem fortran 77.
2.4.2 Soluc¸a˜o do PVI atrave´s do Me´todo de Euler
O algoritmo para a resoluc¸a˜o do problema proposto no item anterior e´ mos-
trado na forma de fluxograma na figura 2.5,
22
N
.M
.S
ot
om
ay
or
Figura 2.4: Programa fonte em FORTRAN 77 para o ca´lculo exa´to da func¸a˜o
2.37.
A implementac¸a˜o computacional do algoritmo mostrado na figura 2.5 foi
realizada em linguagem de programac¸a˜o Fortran 77. O co´digo do programa
fonte e´ mostrado na figura 2.6 a seguir:
O gra´fico mostrado na figura 2.7 representa a soluc¸a˜o do problema 2.37, ava-
liado entre os pontos x0 = 0 e xf = 1, o passo e´ h = 0, 01. A curva na cor
vermelha representa a soluc¸a˜o exata da equac¸a˜o e a curva verde represen-
tam a soluc¸a˜o nume´rica. Pode-se determinar que quanto maior for o passo,
maior sera´ oerro encontrado. A diferenc¸a entre as curvas vermelha e verde
representa o erro entre a soluc¸a˜o exata e a soluc¸a˜o nume´rica.
O me´todo de Euler e´ de passo simples por que a aproximac¸a˜o yi+1 e calculada
a partir somente do valor yi do passo anterior, sendo φ a func¸a˜o incremento.
um Me´todo de passo simples e´ definido na forma:
23
N
.M
.S
ot
om
ay
or
Figura 2.5: Fluxograma do Algoritmo para o Me´todo de Euler
yi+1 = yi + hφ(xi, yi;h) (2.38)
24
N
.M
.S
ot
om
ay
or
Figura 2.6: Programa fonte em FORTRAN 77 que implementa o algoritmo
para o Me´todo de Euler.
25
N
.M
.S
ot
om
ay
or
Figura 2.7: Gra´fico da soluc¸a˜o do PVI do exemplo, sao˜ comparadas a soluc¸a˜o
exata e as soluc¸o˜es aproximadas com diferentes valores do passo h = (xf =
x0)/n. A curva em vermelho representa a soluc¸a˜o exa´ta com n = 20, a curva
em verde representa a soluc¸a˜o aproximada obtida pela aplicac¸a˜o do Me´todo de
Euler com n = 20, a curva em azul representa a soluc¸a˜o aproximada obtida
pela aplicac¸a˜o do Me´todo de Euler com n = 40 e a curva em rouxo representa
a soluc¸a˜o aproximada obtida pela aplicac¸a˜o do Me´todo de Euler com n = 80.
Observa-se que quanto menor o passo a soluc¸a˜o nume´rica aproxima-se mais
da soluc¸a˜o exa´ta.
26
N
.M
.S
ot
om
ay
or
Cap´ıtulo 3
Os Me´todos de Runge-Kutta
de primeira e segunda ordem
3.1 Definic¸a˜o
Em ana´lise nume´rica os me´todos de Runge-Kutta sa˜o uma famı´lia im-
portante de me´todos iterativos impl´ıcitos e expl´ıcitos para a aproximac¸a˜o de
soluc¸o˜es de equac¸o˜es diferenciais ordina´rias. Estas te´cnicas foram desenvol-
vidas por volta de 1900 pelos matema´ticos alema˜es Carl David Tolme´ Runge
e Martin Wilhelm Kutta.
Consideremos o problema de valor inicial,
dy
dx
= f(x, y) (3.1)
Eles sa˜o me´todos de passo simples e apresentam a forma geral.
yi+1 = yi + hφ(xi , yi ; h) (3.2)
onde
φ(xi , yi ; h) = b1k1 + b2k2 + ...+ bsks (3.3)
ou,
yi+1 = yi + h(b1k1 + b2k2 + ...+ bsks) (3.4)
de forma compacta,
27
N
.M
.S
ot
om
ay
or
yi+1 = yi + h
s∑
i=1
biki (3.5)
com,
k1 = f(xi, yi), (3.6)
k2 = f(xi + c2h , yi + a21hk1), (3.7)
k3 = f(xi + c3h , yi + h(a31k1 + a32k2)), (3.8)
... (3.9)
ks = f(xi + csh , yi + h(as1k1 + as2k2 + ...+ as,s−1ks−1)). (3.10)
Sendo a e c constantes definidas para cada me´todo particular.
Para especificar um me´todo part´ıcular e´ necessa´rio escolher:
• o inteiro s (nu´mero de esta´gios, ordem do me´todo);
• os coeficientes aij com i > j, para i = 2, 3, ..., s e j = 1, 2, ..., s− 1;
• os coeficiente bi (para i = 1, 2, ..., s)
• os coeficientes ci (para i = 2, 3, ..., s).
3.1.1 Me´todo de Runge-Kutta de primeira ordem s = 1
Para este tipo particular o nu´mero de esta´gios s = 1, de forma que a equac¸a˜o
4.2 fornece,
yi+1 = yi + hb1k1 (3.11)
yi+1 = yi + hb1f(xi, yi). (3.12)
Por outro lado a expansa˜o em se´rie de Taylor de ordem 1 fornece,
yi+1 = yi + h
dy
dx
∣∣∣∣
x=xi
(3.13)
mas em acordo com a definic¸a˜o 3.1,
dy
dx
∣∣∣∣
x=xi
= f(xi, yi) = (3.14)
28
N
.M
.S
ot
om
ay
or
de forma que,
yi+1 = yi + hf(xi, yi) (3.15)
da comparac¸a˜o das duas equac¸o˜es 3.12 e 3.15, determina-se que,
b1 = 1 (3.16)
Este e´ o chamado me´todo de Runge-Kutta de primeira ordem RK-1 e e´
equivalente ao Me´todo de Euler.
3.1.2 Me´todos de Runge-Kutta de Segunda Ordem
Para este tipo particular o nu´mero de esta´gios s = 2, de forma que a equac¸a˜o
4.2 fornece,
yi+1 = yi + h[b1k1 + b2k2] (3.17)
yi+1 = yi + h
[
b1f(xi , yi) + b2f(xi + c2h , yi + a21hk1)
]
(3.18)
yi+1 = yi + h b1f(xi , yi) + h b2f
[
xi + c2h , yi + a21hf(xi, yi)
]
(3.19)
para,
f
[
xi + c2h , yi + a21hf(xi , yi)
]
(3.20)
comparando com,
f
(
xi +∆x , yi +∆y
)
(3.21)
fazendo expansa˜o em se´rie de Taylor de func¸a˜o de duas varia´veis obtemos,
f(x+∆x, y +∆y) = f(xi, yi) + c2h
∂f(xi, yi)
∂x
+ a21hf(xi, yi)
∂f(xi, yi)
∂y
+O2
(3.22)
substitu´ındo em 3.19,
29
N
.M
.S
ot
om
ay
or
yi+1 = yi+hb1f(xi, yi)+hb2
[
f(xi, yi)+c2 h
∂f(xi, yi)
∂x
+a21hf(xi, yi)
∂f(xi, yi)
∂y
]
(3.23)
yi+1 = yi+h
[
b1f(xi, yi)+b2f(xi, yi)
]
+b2c2h
2∂f(xi, yi)
∂x
+b2a21h
2f(xi, yi)
∂f(xi, yi)
∂y
(3.24)
yi+1 = yi + h
[
b1 + b2
]
f(xi, yi) + b2c2h
2∂f(xi, yi)
∂x
+ b2a21h
2f(xi, yi)
∂f(xi, yi)
∂y
(3.25)
Vamos empregar agora a expansa˜o em se´rie de Taylor de uma func¸a˜o de uma
varia´vel pro´ximo do ponto x0,
y1 = y0 + h
dy
dx
+
h2
2
d2y
dx2
+ ...+O(h3) (3.26)
considerando que para o nosso problema de valor inicial PVI definimos,
dy
dx
= f(x, y) (3.27)
e,
d2y
dx2
=
d
dx
f(x, y) (3.28)
podemos re-escrever 3.26, ate´ segunda ordem, da forma seguinte,
yi+1 = yi + hf(xi, yi) +
h2
2
d
dx
f(xi, yi) (3.29)
entretanto,
d
dx
f(xi, yi) =
[
∂f(xi, yi)
∂x
+
∂f(xi, yi)
∂y
dy
dx
]
(3.30)
usando 4.12 e substitu´ındo 3.30 em 3.29 obtemos,
yi+1 = yi + hf(xi, yi) +
h2
2
[
∂f(xi, yi)
∂x
+ f(xi, yi)
∂f(xi, yi)
∂y
]
(3.31)
ou de forma equivalente,
30
N
.M
.S
ot
om
ay
or
yi+1 = yi + hf(xi, yi) +
h2
2
[
∂f(xi, yi)
∂x
]
+
h2
2
[
f(xi, yi)
∂f(xi, yi)
∂y
]
(3.32)
comparando 3.32 com 3.25 indicada novamente a seguir,
yi+1 = yi + h
[
b1 + b2
]
f(xi, yi) + b2c2h
2∂f(xi, yi)
∂x
+ b2a21h
2f(xi, yi)
∂f(xi, yi)
∂y
concluimos:
b1 + b2 = 1 (3.33)
b2c2 =
1
2
(3.34)
b2a2,1 =
1
2
(3.35)
Estas constantes correspondem ao Me´todo de Runge-Kutta de segunda
ordem. Temos um conjunto de treˆs equac¸o˜es lineares alge´bricas com quatro
inco´gnitas, de acordo com isto ha´ inumeras possibilidades para a escolha
destas constantes. Existira´ um me´todo de Runge-Kutta de segunda ordem
para cada escolha do conjunto de constantes.
A seguir sa˜o apresentados alguns exemplos de casos part´ıculares dos me´todos
de Runge-Kutta de segunda ordem.
3.1.3 Me´todo de Heun (Karl Heun)
Chamado tambe´m Me´todo de Euler Melhorado ou Me´todo de Euler Modifi-
cado, e´ o exemplo expl´ıcito da aplicac¸a˜o da regra do trape´zio. Para o caso
part´ıcular da seguinte escolha das constantes: c2 = 1, b1 =
1
2
, b2 =
1
2
e
a2,1 = 1, substitu´ımos estes valores na equac¸a˜o 3.19 obtendo-se,
yi+1 = yi + b1hf(xi , yi) + b2hf
[
xi + c2h , yi + a21hf(xi, yi)
]
(3.36)
temos enta˜o,
yi+1 = yi +
h
2
f(xi , yi) +
h
2
f
[
xi + h , yi + hf(xi, yi)
]
(3.37)
31
N
.M
.S
ot
om
ay
or
yi+1 = yi +
h
2
[
f(xi , yi) + f
[
xi + h , yi + hf(xi, yi)
]]
(3.38)
Na expressa˜o 3.38 fazemos a identificac¸a˜o e interpretac¸a˜o dos seguintes ter-
mos,
f(xi, yi) =
dy
dx
∣∣∣∣
x=xi
(3.39)
e´ a inclinac¸a˜o no ponto (xi, yi)
f
[
xi + h, yi + hf(xi, yi)
]
=
dy
dx
∣∣∣∣
x=xi+h
(3.40)
e´ a inclinac¸a˜o no ponto
[
xi + h, yi + hf(xi, yi)
]
A equac¸a˜o 3.38 pode re-escrever da forma:
yi+1 = yi +
h
2
[
dy
dx
∣∣∣∣
x=xi
+
dy
dx
∣∣∣∣
x=xi+h
]
(3.41)
Pode-se interpretar a equac¸a˜o 3.41 como a escolha de uma inclinac¸a˜o me´dia
para a melhor previsa˜o do ponto yi+1.
yi+1 = yi + h
[ dy
dx
∣∣∣∣
x=xi
+
dy
dx
∣∣∣∣
x=xi+h
2
]
(3.42)
yi+1 = yi + h Φ (3.43)
O algoritmo deste me´todo particular indica que e´ necessa´rio conhecer as
derivadas da func¸a˜o nos pontos (xi, yi) e (xi+1, yi+1), em seguida e´ calculada
a me´dia aritme´tica destas inclinac¸o˜es. Volta-se para o ponto (xi, yi) e a partir
dele e´ calculada uma melhor aproximac¸a˜o de yi+1 empregando-sea fo´rmula
3.43.
32
N
.M
.S
ot
om
ay
orq
x0
x1
y0
y1
y1x h0+
h
y f x= ( )
dx0y0
q
x0
x1
y0
y1
y1
x h0+
h
y f x= ( )
dx1y1
x0
x1
y0
y1
y1
x h0+
h
y f x= ( )
(dx0y0+dx1y1)/2
(a) (b) (c)
Figura 3.1: Esquema que mostra a interpretac¸a˜o geome´trica do Me´todo de
Heun
3.1.4 Me´todo de Runge-Kutta de segunda ordem (Ponto
Me´dio)
Para o caso part´ıcular da seguinte escolha das constantes: c2 =
1
2
, b1 = 0
, b2 = 1 e a2,1 =
1
2
, substitu´ımos estes valores na equac¸a˜o 3.19 obtendo-se,
yi+1 = yi + b1hf(xi , yi) + b2hf
[
xi + c2h , yi + a21hf(xi, yi)
]
(3.44)
temos enta˜o,
yi+1 = yi + h f
[
xi +
h
2
, yi +
h
2
f(xi, yi)
]
(3.45)
Na expressa˜o 3.45 fazemos a identificac¸a˜o e interpretac¸a˜o dos seguintes ter-
mos,
f(xi, yi) =
dy
dx
∣∣∣∣
x=xi
(3.46)
e´ a inclinac¸a˜o no ponto (xi , yi)
f
[
xi +
h
2
, yi +
h
2
f(xi, yi)
]
=
dy
dx
∣∣∣∣∣
x=xi+
h
2
(3.47)
33
N
.M
.S
ot
om
ay
or
e´ a inclinac¸a˜o no ponto me´dio
[
xi +
h
2
, yi +
h
2
f(xi, yi)
]
A equac¸a˜o 3.45 pode re-escrever da forma:
yi+1 = yi + h
dy
dx
∣∣∣∣
x=xi+
h
2
(3.48)
ou,
yi+1 = yi + h Φpm (3.49)
O algoritmo deste me´todo particular indica que e´ necessa´rio calcular a
derivada da func¸a˜o no ponto xi+ h/2, em seguida e´ realizado o avanc¸o ate´ o
ponto xi + h e calcula-se a aproximac¸a˜o para y1 empregando-se para isto a
inclinac¸a˜o computada para o ponto xi + h/2.
34
N
.M
.S
ot
om
ay
or
i=i+1
Início
x:xF
<
>=
Fim
Ler
x,y,x ,n0 0 f
h=(x -x)/nf 0
x=x,y=y0 0
i=1
dx0y0=x-2.0*y+1.0
x=x +i*h0
ylast=y
y=y+h*dx0y0
dx1y1=x-2.0*y+1.0
y=ylast+h*(dx0y0+dx1y1)/2
Print
x,y
Figura 3.2: Fluxograma que explica a implementac¸a˜o computacional do
Me´todo de Heun
35
N
.M
.S
ot
om
ay
or
Figura 3.3: Programa fonte em Fortran 77, que implementa o Me´todo de
integrac¸a˜o nume´rica de Heun.
36
N
.M
.S
ot
om
ay
or
0.75
0.8
0.85
0.9
0.95
1
0 0.1 0.2 0.3 0.4 0.5 0.6
"trayexato.dat"
"trayeuler.dat"
"trayheun.dat"
Figura 3.4: Comparac¸a˜o dos ca´lculos da soluc¸a˜o exa´ta do PIV e dos Me´todos
de Euler e Heun (RK-2).
37
N
.M
.S
ot
om
ay
or
x y (Exato) y (Euler) y (Heun)
0.0000000000000000 1.0000000000000000 1.0000000000000000 1.0000000000000000
2.99999999999999989E-002 0.97132340018818653 0.96999999999999997 0.97135000000000005
5.99999999999999978E-002 0.94519032753786814 0.94269999999999998 0.94524043000000002
8.99999999999999967E-002 0.92145265855845404 0.91793800000000003 0.92152343697400008
0.12000000000000000 0.89997089579991507 0.89556172000000001 0.90005977294211326
0.14999999999999999 0.88061366551128839 0.87542801680000004 0.88071829415688230
0.17999999999999999 0.86325724455327324 0.85740233579200009 0.86337548943695175
0.20999999999999999 0.84778511486129260 0.84135819564448011 0.84791503595172113
0.23999999999999999 0.83408754385460560 0.82717670390581133 0.83422738085933101
0.27000000000000002 0.82206118928049221 0.81474610167146266 0.82220934729331796
0.29999999999999999 0.81160872707051990 0.80396133557117488 0.81176376328084687
0.32999999999999996 0.80263850086877442 0.79472365543690437 0.80279911225790157
0.35999999999999999 0.79506419196997880 0.78694023611069008 0.79522920392449170
0.39000000000000001 0.78880450847891770 0.78052382194404868 0.78897286425608626
0.41999999999999998 0.78378289257180978 0.77539239262740578 0.78395364355638208
0.44999999999999996 0.77992724480544939 0.77146884906976143 0.78009954150140060
0.47999999999999998 0.77716966448133407 0.76868071812557570 0.77734274818601912
0.51000000000000001 0.77544620512980877 0.76695987503804119 0.77561940024159282
0.54000000000000004 0.77469664423370432 0.76624228253575877 0.77486935114753208
0.56999999999999995 0.77486426636222783 0.76646774558361319 0.77503595491074573
0.59999999999999998 0.77589565893415158 0.76757968084859640 0.77606586233494035
Figura 3.5: Comparac¸a˜o dos valores nume´ricos da soluc¸a˜o exa´ta do PIV e
dos Me´todos de Euler e Heun (RK-2).
38
N
.M
.S
ot
om
ay
or
Figura 3.6: Programa fonte em Fortran 77, que implementa o Me´todo de
integrac¸a˜o nume´rica RK-2 Ponto Me´dio.
39
N
.M
.S
ot
om
ay
or
q
x0
x1
y0
y1
y1
x h0+
h
y f x= ( )
dx dy0 0
x0
x1xM
y0
yM
y1
y1
x h0+x h/20+
y f x= ( )
dx dyM M
dx dyM M
x0
x1
y0
y1
y1
y1
x h0+
y f x= ( )
(a) (b) (c)
Figura 3.7: Interpretac¸a˜o geome´trica do Me´todo de integrac¸a˜o nume´rica RK-
2 Ponto Me´dio.
40
N
.M
.S
ot
om
ay
or
0.75
0.8
0.85
0.9
0.95
1
0 0.1 0.2 0.3 0.4 0.5 0.6
"trayeuler06.dat"
"trayexato06.dat"
"trayheun06.dat"
"traymid06.dat"
Figura 3.8: Comparac¸a˜o dos valores nume´ricos da soluc¸a˜o exa´ta do PIV e
dos Me´todos de Euler, Heun e Midpoint (RK-2).
41
N
.M
.S
ot
om
ay
or
Cap´ıtulo 4
O Me´todo de Runge-Kutta de
Terceira ordem
4.1 Resumo de teoria ba´sica
yi+1 = yi + h(b1k1 + b2k2 + ...+ bsks) (4.1)
yi+1 = yi + h
(
b1k1 + b2k2 + b3k3
)
(4.2)
com os coefientes dados por:
k1 = f(xi, yi), (4.3)
k2 = f(xi + c2h , yi + a21hk1), (4.4)
k3 = f
[
xi + c3h , yi + h(a31k1 + a32k2)
]
. (4.5)
de forma que a equac¸a˜o 4.2 pode-se escrever como,
yi+1 = yi + hb1k1 + hb2f
[
xi + c2h , yi + a21hk1
]
+
hb3f
[
xi + c3h , yi + h(a31k1 + a32k2)
] (4.6)
(4.7)
Lembramos a expansa˜o em se´rie de Taylor de uma func¸a˜o de duas varia´veis
ate´ terceira ordem,
42
N
.M
.S
ot
om
ay
or
f(x+∆x, y +∆y) = f(xi, yi) + (∆x)
∂f(xi, yi)
∂x
+ (∆y)
∂f(xi, yi)
∂y
+
(∆x)2
2!
∂2f(xi, yi)
∂x2
+ (∆x)(∆y)
∂2f(xi, yi)
∂y∂x
+
(∆y)2
2!
∂2f(xi, yi)
∂y2
+
(∆x)3
3!
∂3f(xi, yi)
∂x3
+
(∆x)2(∆y)
2
∂3f(xi, yi)
∂x∂x∂y
+
(∆x)(∆y)2
2
∂3f(xi, yi)
∂x∂y∂y
+
(∆y)3
3!
∂3f(xi, yi)
∂y3
+O(∆4).
(4.8)
A seguir vamos expandir f [xi+ c2h , yi+a21hk1] e f [xi+ c3h , yi+h(a31k1+
a32k2)] ate´ segunda ordem usando a expressa˜o 4.8 para logo substituir esses
resultados na equac¸a˜o 4.6 (fica como exerc´ıcio para o aluno a verificac¸a˜o).
yi+1 = yi + hb1f + hb2f + h
2b2c2
df
dx
+ h2b2a21f
df
dy
+
h3b2c
2
2
2
d2f
dx2
+
h3b2c2a21f
∂2f
∂x∂y
+
h3b2a
2
21f
2
2
∂2f
∂y2
+ hb3f + h
2b3c3
df
dx
+ h2b3a31f
df
dy
+
h3b3a32c2
df
dx
df
dy
+ h3b3a32a21f
(
df
dy
)2
+
h3b3c
2
3
2
∂2f
∂x2
+ h3b3c3a31f
∂2f
∂x∂y
+
h3b3c3a32f
∂2f
∂x∂y
+
h3b3a
2
31
2
f 2
∂2f
∂y2
+ h3b3a31a32f
2∂
2f
∂y2
+
h3b3a
2
32
2
f 2
∂2f
∂y2
+ ...
(4.9)
Organizando termos obtemos,
43
N
.M
.S
ot
om
ay
or
yi+1 = yi + h
(
b1 + b2 + b3
)
f + h2
(
b2c2 + b3c3
)
df
dx
+
h2
(
b2a21 + b3a31 + b3a32
)
f
df
dy
+ h3
(
b2c2a21 + b3c3a31 + b3c3a32
)
f
∂2f
∂x∂y
+
h3
2
(
b2c
2
2 + b3c
2
3
)
∂2f
∂x2
+ h3b3a32c2
∂f
∂x
∂f
∂y
+ h3b3a32a21f
(
∂f
∂y
)2
+
h3
2
(
b2a
2
21 + b3a
2
31 + 2b3a31a32 + b3a
2
32
)
f 2
∂2f
∂y2
+ ...
(4.10)
Guardemos a equac¸a˜o 4.10 por um momento, faremos agora a expansa˜o de
uma func¸a˜o de uma varia´vel ate´ terceira ordem,
y1 = y0 + h
dy
dx
+
h2
2
d2y
dx2
+
h3
6
d3y
dx3
+ ...+O(h4) (4.11)
considerando que para o nosso problema de valor inicial PVI definimos,
dy
dx
= f(x, y) (4.12)
d2y
dx2
=
d
dx
f(x, y) =
[
∂f(x, y)
∂x
+
∂f(x, y)
∂y
dy
dx
]
(4.13)
d2y
dx2
=
[
∂f(x, y)
∂x+
∂f(x, y)
∂y
dy
dx
]
(4.14)
d3y
dx3
=
d
dx
(
d2y
dx2
)
(4.15)
d3y
dx3
=
∂2
∂x2
f(x, y) + 2f(x, y)
∂2
∂x ∂y
f(x, y) +
∂
∂x
f(x, y)
∂
∂y
f(x, y)+[
f(x, y)
]2
∂2
∂y2
f(x, y) + f(x, y)
[
∂
∂y
f(x, y)
]2
(4.16)
44
N
.M
.S
ot
om
ay
or
Prova
d3y
dx3
=
d
dx
(
d2y
dx2
)
=
d
dx
[
∂f(x, y)
∂x
+
∂f(x, y)
∂y
dy
dx
]
(4.17)
d3y
dx3
=
d
dx
[
∂f(x, y)
∂x
]
+
d
dx
[
∂f(x, y)
∂y
dy
dx
]
(4.18)
d3y
dx3
=
∂2f
∂y∂x
dy
dx
+
∂2f
∂x2
+
df
dx
[
df
dy
]
+ f
d
dx
[
df
dy
]
(4.19)
d3y
dx3
=
∂2f
∂y∂x
dy
dx
+
∂2f
∂x2
+
[
∂f
∂x
+ f
∂f
∂y
]
∂f
∂y
+ f
[
d
dx
[
∂f
∂y
]]
(4.20)
d3y
dx3
=
∂2f
∂y∂x
dy
dx
+
∂2f
∂x2
+
∂f
∂x
∂f
∂y
+ f
∂f
∂y
∂f
∂y
+ f
[[
∂
∂y
(
df
dy
)]
dy
dx
+
(
∂
∂x
df
dy
)]
(4.21)
d3y
dx3
=
∂2f
∂y∂x
dy
dx
+
∂2f
∂x2
+
∂f
∂x
∂f
∂y
+ f
(
∂f
∂y
)2
+ f 2
∂2f
∂y2
+ f
∂2f
∂x∂y
(4.22)
d3y
dx3
= f
∂2f
∂y∂x
+
∂2f
∂x2
+
∂f
∂x
∂f
∂y
+ f
(
∂f
∂y
)2
+ f 2
∂2f
∂y2
+ f
∂2f
∂x∂y
(4.23)
d3y
dx3
= 2f
∂2f
∂y∂x
+
∂2f
∂x2
+
∂f
∂x
∂f
∂y
+ f
(
∂f
∂y
)2
+ f 2
∂2f
∂y2
(4.24)
que e´ igual a:
d3y
dx3
=
∂2
∂x2
f(x, y) + 2f(x, y)
∂2
∂x ∂y
f(x, y) +
∂
∂x
f(x, y)
∂
∂y
f(x, y)+[
f(x, y)
]2
∂2
∂y2
f(x, y) + f(x, y)
[
∂
∂y
f(x, y)
]2
(4.25)
45
N
.M
.S
ot
om
ay
or
Voltando para,
y1 = y0 + h
dy
dx
+
h2
2
d2y
dx2
+
h3
6
d3y
dx3
(4.26)
substu´ındo os valores encontrados obtemos,
y1 = y0+hf+
h2
2
[
∂f
∂x
+f
∂f
∂y
]
+
h3
6
[
2f
∂2f
∂y∂x
+
∂2f
∂x2
+
∂f
∂x
∂f
∂y
+f
(
∂f
∂y
)2
+f 2
∂2f
∂y2
]
(4.27)
ou de forma compacta,
y1 = y0 + hf +
h2
2
∂f
∂x
+
h2
2
f
∂f
∂y
+
h3
3
f
∂2f
∂y∂x
+
h3
6
∂2f
∂x2
+
h3
6
∂f
∂x
∂f
∂y
+
h3
6
f
(
∂f
∂y
)2
+
h3
6
f 2
∂2f
∂y2
OK
(4.28)
Comparamos agora as equac¸o˜es 4.10 e 4.28 e igualamos os coeficientes das
duas se´ries obtendo-se o seguinte sistema de equac¸o˜es alge´bricas,
b1 + b2 + b3 = 1 (1)
b2c2 + b3c3 =
1
2
(2)
b2a21 + b3a31 + b3a32 =
1
2
(3)
b2c
2
2 + b3c
2
3 =
1
3
(4)
b2c2a21 + b3c3a31 + b3c3a32 =
1
3
(5)
b2a
2
21 + b3a
2
31 + 2b3a31a32 + b3a
2
32 =
1
3
(6)
b3a32c2 =
1
6
(7)
b3a32a21 =
1
6
(8)
(4.29)
46
N
.M
.S
ot
om
ay
or
dividindo (7) e (8) obtemos,
c2 = a21 (4.30)
subtraindo (6) de (5) obtemos,
c3 = a31 + a32 (4.31)
por ser um sistema indeterminado escolhemos,
b3 =
1
6
(4.32)
e,
a21 =
1
2
(4.33)
substitu´ındo em (8) obtem-se,
a32 = 2 (4.34)
e,
c2 =
1
2
(4.35)
usando (2) e (4) encontra-se,
b2 =
2
3
(4.36)
e,
c3 = 1 (4.37)
em (1),
b1 =
1
6
(4.38)
e finalmente,
a31 = −1. (4.39)
Com esses resultados a equac¸a˜o 4.1 pode-se escrever como,
yi+1 = yi + h
(
1
6
k1 +
2
3
k2 +
1
6
k3
)
(4.40)
47
N
.M
.S
ot
om
ay
or
ou ainda,
yi+1 = yi +
h
6
(
k1 + 4k2 + k3
)
(4.41)
com,
k1 = f(xi, yi), (4.42)
k2 = f
(
xi +
h
2
, yi +
h
2
k1
)
, (4.43)
k3 = f
(
xi + h , yi − hk1 + 2hk2
)
. (4.44)
48
N
.M
.S
ot
om
ay
or
Figura 4.1: Programa fonte em Fortran 77, que implementa o Me´todo de
integrac¸a˜o nume´rica de Runge-Kutta de ordem 3 RK-3.
49
N
.M
.S
ot
om
ay
or
0.75
0.8
0.85
0.9
0.95
1
0 0.1 0.2 0.3 0.4 0.5 0.6
"trayexato.dat"
"trayrk3.dat "
Figura 4.2: Comparac¸a˜o de trajeto´rias da soluc¸a˜o exa´ta e do me´todo de
Runge-Kutta de ordem 3 para o PIV em estudo.
50
N
.M
.S
ot
om
ay
orPIV SOLUÇÃO EXATAx y
0.0000000000000000 1.0000000000000000
2.99999999999999989E-002 0.97132340018818653
5.99999999999999978E-002 0.94519032753786814
8.99999999999999967E-002 0.92145265855845404
0.12000000000000000 0.89997089579991507
0.14999999999999999 0.88061366551128839
0.17999999999999999 0.86325724455327324
0.20999999999999999 0.84778511486129260
0.23999999999999999 0.83408754385460560
0.27000000000000002 0.82206118928049221
0.29999999999999999 0.81160872707051990
0.32999999999999996 0.80263850086877442
0.35999999999999999 0.79506419196997880
0.39000000000000001 0.78880450847891770
0.41999999999999998 0.78378289257180978
0.44999999999999996 0.77992724480544939
0.47999999999999998 0.77716966448133407
0.51000000000000001 0.77544620512980877
0.54000000000000004 0.77469664423370432
0.56999999999999995 0.77486426636222783
0.59999999999999998 0.77589565893415158
MÉTODO DE RUNGE-KUTTA ORDEM 3
y y
0.0000000000000000 1.0000000000000000
2.99999999999999989E-002 0.97132300000000005
5.99999999999999978E-002 0.94518957377200008
8.99999999999999967E-002 0.92145159375381391
0.12000000000000000 0.89996955873996676
0.14999999999999999 0.88061209151718600
0.17999999999999999 0.86325546575559120
0.20999999999999999 0.84778316045184854
0.23999999999999999 0.83408544031977472
0.27000000000000002 0.82205896061731232
0.29999999999999999 0.81160639498680254
0.32999999999999996 0.80263608496835115
0.35999999999999999 0.79506170992413427
0.39000000000000001 0.78880197618499237
0.41999999999999998 0.78378032429988320
0.44999999999999996 0.77992465333395522
0.47999999999999998 0.77716706122239898
0.51000000000000001 0.77544360024505132
0.54000000000000004 0.77469404674118048
0.56999999999999995 0.77486168423516111
0.59999999999999998 0.77589309919204230
Figura 4.3: Comparac¸a˜o dos valores nume´ricos das trajeto´rias da soluc¸a˜o
exa´ta e do me´todo de Runge-Kutta de ordem 3 para o PIV em estudo.
51
N
.M
.S
ot
om
ay
or
Cap´ıtulo 5
O Me´todo de Runge-Kutta de
Quarta Ordem
5.1 Teoria ba´sica
yi+1 = yi + hφ(xi, yi, h) (5.1)
yi+1 = yi + h[a1k1 + a2k2 + a3k3 + a4k4] (5.2)
onde,
k1 = f(xi, yi) (5.3)
k2 = f
(
xi +
h
2
, yi +
h
2
k1
)
(5.4)
k3 = f
(
xi +
h
2
, yi +
h
2
k2
)
(5.5)
k4 = f(xi + h, yi + hk3) (5.6)
A expressa˜o para o valor aproximado da func¸a˜o e´,
yi+1 = yi +
h
6
[
k1 + 2k2 + 2k3 + k4
]
(5.7)
Se definimos a inclinac¸a˜o me´dia por,
φ =
k1 + 2k2 + 2k3 + k4
6
(5.8)
52
N
.M
.S
ot
om
ay
or
teremos,
yi+1 = yi + hφ (5.9)
5.1.1 Interpretac¸a˜o geome´trica
• k1 e´ a inclinac¸a˜o no in´ıcio do intervalo xi, yi;
• k2 e´ a inclinac¸a˜o na metade do intervalo, usando o valor da inclinac¸a˜o
k1 para determinar o valor de yM no ponto xM = xi +
h
2
(usando o
me´todo de Euler);
• k3 e´ a inclinac¸a˜o no ponto me´dio, mas agora usando a inclinac¸a˜o k2
para determinar o valor de y;
• k4 e´ a inclinac¸a˜o no final do intervalo, com o valor y determinado usando
k3.
O me´todo de 4a ordem e´ um dois mais utilizados. Pode ser implementado
utilizando o seguinte algoritmo:
53
N
.M
.S
ot
om
ay
or
Figura 5.1: Co´digo fortran 77 mostrando a implementac¸a˜o do Me´todo de
integrac¸a˜o nume´rica de Runge-Kutta quarta Ordem
54
N
.M
.S
ot
om
ay
or
x0
y0
x1
k1
y1
xM
yM1
yM
Figura 5.2: k1 e´ a inclinac¸a˜o em (x0, y0). Avanc¸a-se meio passo e calcula-se
yM1 com k1 atrave´s do me´todo de Euler.
x0
y0
x1
k1
k2
y1
xM
yM1
yM2
yM
Figura 5.3: k2 e´ a inclinac¸a˜o em (xM1, yM1). Volta-se para x0 e calcula-se
yM2,com k2. Usando-se (xM , yM2) calcula-se k3.
55
N
.M
.S
ot
om
ay
or
x0
y0
x1
k1
k2
k3
k4
y1
xM
yM1
y 1T
yM2
yM
Figura 5.4: Usando-se k3 e o me´todo de Euler encontra-se y1T . Em seguida,
empregando-se (x1, y1T ) encontras-se k4.
x0
y0
x1
k1
6
+ 2k2+ 2k3+ k4
y1
xM
yM
k
k =
e
Figura 5.5: Finalmente, volta-se para x0 e calcula-se a aproximac¸a˜o para y1
empregando-se a inclina˜c¸a˜o me´dia, usando o me´todo de Euler.
56
N
.M
.S
ot
om
ay
or
Cap´ıtulo 6
Integrac¸a˜o nume´rica das
equac¸a˜oes de Newton
6.1 Forc¸a sobre um objeto em queda livre
Consideremos a segunda lei de Newton para o movimento de queda livre de
uma part´ıcula de massa m sobre a superf´ıcie da terra,
F (y, v, t) = m
d2y
dt2
(6.1)
Podemos extrair duas equac¸o˜es,
a(t) =
F (y, v, t)
m
(6.2)
e,
v(t) =
dy(t)
dt
(6.3)
Primeiramente aplicamos o me´todo de Euler a`s equac¸o˜es de movimento de
Newton,
vn+1 = vn + an∆t (6.4)
yn+1 = yn + vn∆t (6.5)
onde ∆t e´ o passo de tempo an = an(yn, vn, tn). Temos enta˜o um sistema de
duas equac¸o˜es diferenciais acopladas de primeira ordem.
Uma variac¸a˜o do algoritmo de resoluc¸a˜o de Euler e´ o Me´todo de Euler-
Cromer,
57
N
.M
.S
ot
om
ay
or
vn+1 = vn + an∆t (6.6)
yn+1 = yn + vn+1∆t (6.7)
Calcula a posic¸a˜o avaliando v ao final do intervalo ao inve´s do in´ıcio, pore´m,
uma escolha poss´ıvel que poderia fornecer melhores resultados seria calcular
a velocidade na metade do intervalo. Este me´todo e´ chamado algoritmo de
Euler-Richardson,
vn+1 = vn + amid∆t (6.8)
yn+1 = yn + vmid∆t (6.9)
6.1.1 Implementac¸a˜o computacional
Primeiramente realizamos a simulac¸a˜o da queda livre de um objeto de massa
m usando o me´todo de Euler. A figura 6.1 mostra o co´digo fonte em fortran
77 que realiza este processo.
A seguir, apresentamos a simulac¸a˜o da queda de um objeto de massa m
quando ale´m da forc¸a da gravidade temos resisteˆncia do a´r.
F = −mg + Fa (6.10)
A forma de Fa pode ser obtida atrave´s de experimentos, ha´ duas apro-
ximac¸a˜oes mais usadas Fa ∼ v, e Fa ∼ v2,
Fa = αv (6.11)
e,
Fa = βv
2 (6.12)
A velocidade terminal ocorre quando Fa = mg, para o primeiro caso vt =
mg
α
e para o segundo caso vt =
(
mg
β
) 1
2
. Da equac¸a˜o 6.10 obtemos,
F1 = −mg
(
1− αv
mg
)
= −mg
(
1− v
mg/α
)
(6.13)
F2 = −mg
(
1− βv
2
mg
)
= F2 = −mg
(
1− v
2
mg/β
)
(6.14)
de forma que podemos escrever
58
N
.M
.S
ot
om
ay
or
Figura 6.1: Co´digo fortran 77 mostrando a implementac¸a˜o do Me´todo de
Euler
59
N
.M
.S
ot
om
ay
or
0
2
4
6
8
10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
y
(m
)
t (s)
"fall.dat"
Figura 6.2: Queda de um corpo de massa m sem resisteˆncia do a´r.
60
N
.M
.S
ot
om
ay
or
F1 = −mg
(
1− v
vt
)
(6.15)
e,
F2 = −mg
(
1− v
2
v2t
)
(6.16)
Assim, a acelerac¸a˜o dos sistemas e´ dada por,
a = −g
(
1− v
vt
)
(6.17)
e,
a = −g
(
1− v
2
v2t
)
(6.18)
Usamos enta˜o qualquer uma das duas equac¸o˜es em conjunto com as duas
seguintes,
vn+1 = vn + amid∆t (6.19)
yn+1 = yn + vmid∆t (6.20)
que sa˜o as equac¸o˜es correspondentes ao me´todo de Euler-Richardson. Pre-
viamente a acelerac¸a˜o, velocidade e posic¸a˜o tem que ser calculadas no ponto
me´dio do intervalo correspondente a um passo de integrac¸a˜o.
6.2 Trajeto´rias em duas dimenso´es
Consideremos o movimento de um proje´til em duas dimenso˜es. O proje´til e´
lanc¸ado do cha˜o com velocidade inicial fazendo um aˆngulo θ com a horizontal
(eixo x), sera´ levado em conta a forc¸a de atrito com o ar proporcional ao
quadrado da velocidade.
A seguna lei de Newton descreve a dinaˆmica do sistema,
~F = m~a (6.21)
ou,
61
N
.M
.S
ot
om
ay
or
Figura 6.3: Co´digo fortran 77 mostrando a implementac¸a˜o do Me´todo de
Euler-Richardson
62
N
.M
.S
ot
om
ay
or
0
2
4
6
8
10
0 0.2 0.4 0.6 0.8 1 1.2 1.4
y
(m
)
t (s)
"falld.dat"
Figura 6.4: Queda de um corpo de massa m com resisteˆncia do a´r proporci-
onal a v2.
q
v0
|v0 |cosq
|v0 |senq
mg
Fa
x
y
Figura 6.5: Movimento de um proje´til em 2D.
63
N
.M
.S
ot
om
ay
or
Fx = m
dvx
dt
(6.22)
Fy = m
dvy
dt
(6.23)
Fz = m
dvz
dt
(6.24)
vamos considerar que na˜o ha´ movimento na direc¸a˜o x enta˜o,
m
dvx
dt
= −Fa cos θ (6.25)
m
dvy
dt
= −Fa sin θ −mg (6.26)
consideremos tambe´m atrito proporcional ao quadrado da velocidade,
dvx
dt
= −Fa
m
cos θ (6.27)
dvy
dt
= −Fa
m
sin θ − g (6.28)
Lembrando que,
Fa = βv
2 (6.29)
re-escrevemos as equac¸o˜es da forma,
dvx
dt
= −βv
2
m
cos θ (6.30)
dvy
dt
= −βv
2
m
sin θ − g (6.31)
ou ainda,
dvx
dt
= −βv
m
v cos θ (6.32)
dvy
dt
= −βv
m
v sin θ − g (6.33)
Finalmente,
64
N
.M
.S
ot
om
ay
or
dvx
dt
= −βv
m
vx (6.34)
dvy
dt
= −βv
m
vy − g. (6.35)
Vamos simplificas as duas equac¸a˜oes introduzindo a constante k = β/m,
dvx
dt
= −k v vx
dvy
dt
= −k v vy − g.
(6.36)
A equac¸a˜o de acoplamento das duas equac¸o˜es e´,
v2 = v2x + v
2
y (6.37)
A seguir e´ apresentado o programa fonte para o ca´lculo das trajeto´rias
A seguir sa˜o apresentados os resultados obtidos para as trajeto´rias do proje´til
lanc¸ado do cha˜o com aˆngulo de 45o e paraˆmetro k variando de 0, 0.01, 0.02, e 0.04.
65
N
.M
.S
ot
om
ay
or
Figura 6.6: Programa fonte para o ca´lculo das trajeto´rias de um proje´til em
2D.
66
N
.M
.S
ot
om
ay
or
0
1
2
3
4
5
6
0 5 10 15 20
y
(m
)
x (m)
"tray2dk004.dat"
"tray2dk0.dat"
"tray2dk001.dat"
"tray2dk002.dat"
q=45
o
Figura 6.7: Gra´fico das trajeto´rias nume´ricas de um proje´til em 2D lanc¸ado
do n´ıvel do solo, comaˆngulo de 45o e para diferentes valores do paraˆmetro k.
67
N
.M
.S
ot
om
ay
or
Cap´ıtulo 7
Integrac¸a˜o nume´rica das
equac¸a˜oes diferenciais de
segunda ordem
7.1 Me´todo de Numerov
O me´todo de Numerov e´ um me´todo nume´rico para a resoluc¸a˜o de equac¸o˜es
diferenciais ordina´rias de segunda ordem nas quais na˜o aparece o termo de
derivada primeira ordem. Este me´todo foi desenvolvido por Boris Vasil’evich
Numerov. E´ um me´todo multipasso linear de quarta ordem.
A fo´rmula de iterac¸a˜o e´, [
d2
dx2
+ f(x)
]
y(x) = 0 (7.1)
A func¸a˜o y(x) e´ avaliada no intervalo [a, b] o qual e´ particionado em intervalos
iguais xn. Iniciando com valores iniciais de y(x) em dois valores consequtivos
xn−1 e xn a aproximac¸a˜o da func¸a˜o pode ser calculada pela equac¸a˜o:
yn+1 =
[
2− 5h
2
6
f(xn)
]
y(xn)−
[
1 +
h2
12
f(xn−1)
]
y(xn−1)
1 +
h2
12
f(xn+1)
(7.2)
onde h = xn − xn−1.
68
N
.M
.S
ot
om
ay
or
Figura 7.1: Implementac¸a˜o computacional do Me´todo de Numerov.
69
N
.M
.S
ot
om
ay
or
7.2 Derivac¸a˜o
Realizamos a expansa˜o em se´rie de Taylor de y(x) em torno do ponto x0 ate´
quinta ordem,
y(x) = y(x0) + (x− x0)y′(x0) + (x− x0)
2
2!
y′′(x0) +
(x− x0)3
3!
y′′′(x0) (7.3)
+
(x− x0)4
4!
y′′′′(x0) +
(x− x0)5
5!
y′′′′′(x0) +O[(x− x0)6]
Denotamos a distaˆncia x− x0 por h e, fazendo x = x0 + h,
y(x0 + h) = y(x0) + hy
′(x0) +
h2
2!
y′′(x0) +
h3
3!
y′′′(x0) (7.4)
+
h4
4!
y′′′′(x0) +
h5
5!
y′′′′′(x0) +O(h6)
Substitu´ındo h por −h obtemos tambe´m,
y(x0 − h) = y(x0)− hy′(x0) + h
2
2!
y′′(x0)− h
3
3!
y′′′(x0) (7.5)
+
h4
4!
y′′′′(x0)− h
5
5!
y′′′′′(x0) +O(h6)
Apenas aspoteˆncias ı´mpares mundam de sinal, assim, temos os pontos de
amostragem: (xn−1, yn−1) e (xn+1, yn+1). Tomando as equac¸o˜es para yn−1 e
yn+1 observamos que,
yn+1 = y(xn + h) = y(xn) + hy
′(xn) +
h2
2!
y′′(xn) +
h3
3!
y′′′(xn) (7.6)
+
h4
4!
y′′′′(xn) +
h5
5!
y′′′′′(xn) +O(h6)
yn−1 = y(xn − h) = y(xn)− hy′(xn) + h
2
2!
y′′(xn)− h
3
3!
y′′′(xn) (7.7)
+
h4
4!
y′′′′(xn)− h
5
5!
y′′′′′(xn) +O(h6)
A soma destas duas equac¸o˜es fornece,
70
N
.M
.S
ot
om
ay
or
yn−1 + yn+1 = 2yn + h2y′′n +
h4
12
y′′′′n +O(h6) (7.8)
Resolvemos esta equac¸a˜o para y′′n e substituimos ela pela expressa˜o y
′′
n =
−anyn que obtemos da equac¸a˜o diferencial.
h2anyn = 2yn − yn−1 − yn+1 + h
4
12
y′′′′n +O(h6) (7.9)
We take the second derivative of our defining differential equation and get,
y′′′′(x) = − d
2
dx2
[a(x)y(x)] (7.10)
We replace the second derivative d
2
dx2
with the second order difference quotient
and insert this into our equation for anyn (note that we take the mixed
forward and backward finite difference, not the double forward difference or
the double backward difference),
h2anyn = 2yn − yn−1 − yn+1 − h
4
12
an−1yn−1 − 2anyn + an+1yn+1
h2
+O(h6)
(7.11)
We solve for yn+1 to get,
yn+1 =
(
2− 5h2
6
an
)
yn −
(
1 + h
2
12
an−1
)
yn−1
1 + h
2
12
an+1
+O(h6) (7.12)
This yields Numerov’s method if we ignore the term of order h6. It follows
that the order of convergence (assuming stability) is 4.
71
N
.M
.S
ot
om
ay
or
Cap´ıtulo 8
Integrac¸a˜o Nume´rica de
Equac¸o˜es Diferenciais Parciais
8.1 Introduc¸a˜o
Equac¸o˜es diferenciais parciais nume´ricas e´ parte da ana´lise nume´rica que
estuda as soluc¸o˜es nume´ricas das equac¸o˜es diferenciais parciais. Dentre as
diferentes te´cnicas de resoluc¸a˜o de equac¸o˜es diferenciais parciais nume´ricas
empregaremos o me´todo de diferenc¸as finitas.
8.2 Equac¸o˜es diferenciais parciais
Uma equac¸a˜o diferencial parcial EDP e´ uma equac¸a˜o composta de func¸o˜es
multivaria´veis desconhecidas e as suas derivadas parciais. Uma EDP para a
func¸a˜o u(x1, x2, ..., xn), e´ uma equac¸a˜o da forma,
F (x1, x2, ..., xn, u, ux1 , ux2 , ..., uxn , ux1x1 , ux1x2 , ..., ux1xn , ...) = 0 (8.1)
Se F e´ uma func¸a˜o linear de u e de suas derivadas, logo a EDP e´ denominada
linear.
8.3 Me´todo das diferenc¸as finitas
Neste me´todo as func¸o˜es sa˜o representadas pelo seus valores em certos ponto
de uma malha e suas derivadas sa˜o aproximadas atrave´s de diferenc¸as nestes
valores.
72
N
.M
.S
ot
om
ay
or
8.4 Soluc¸a˜o nume´rica da equac¸a˜o de Schro¨din-
ger dependente do tempo
O objetivo deste cap´ıtulo e´ mostrar a evoluc¸a˜o temporal de um pacote de
onda gaussiano atrave´s da soluc¸a˜o nume´rica da equac¸a˜o de Schro¨dinger de-
pendente do tempo. Para isso, e´ empregada a aproximac¸a˜o de diferenc¸as fini-
tas separandose as partes real e imagina´ria da func¸a˜o de onda e empregando-
se o me´todo de passo me´dio. Finalmente, apresenta-se o aplicativo Gnuplot
do Projeto GNU como uma importante ferramenta para a gerac¸a˜o de visua-
lizac¸o˜es animadas quando usada em combinac¸a˜o com o gfortran e o shell do
sistema Linux.
8.5 Fundamentos
A equac¸a˜o de Schro¨dinger dependente do tempo, na˜o relativista, para uma
part´ıcula de massa m se movendo em uma dimensa˜o na presenc¸a de um
potencial V (x, t) e´:
− ~
2
2m
∂2
∂x2
Ψ(x, t) + V (x, t)Ψ(x, t) = i~
∂
∂t
Ψ(x, t) (8.2)
m e´ a massa da part´ıcula, ~ e´ a constante racionalizada de Planck, e V (x, t)
e´ a func¸a˜o potencial, aqual e´ considerada independente do tempo. A soluc¸a˜o
da equac¸a˜o nem sempre e´ direta, entretanto, uma aproximac¸a˜o nume´rica da
soluc¸a˜o pode ser obtida introduzindo uma malha retangular com partic¸o˜es
para a coordenada espacial e para a coordenada temporal. Denotaremos as
partic¸o˜es por,
tn = t0 + n4t (8.3)
xv = x0 + v4x (8.4)
A func¸a˜o de onda num ponto particular do domı´nio e´ denotada por Ψ(xv, tn).
Uma primeira aproximac¸a˜o para a soluc¸a˜o nume´rica da equac¸a˜o 8.2 e´ ob-
tida quando empregadas diferenc¸as finitas para relacionar Ψ(xv, tn) com
Ψ(xv, tn+1) para cada valor de xv (esquema expl´ıcito) [2].
1
∆t
[
Ψ(xv, tn+1)−Ψ(xv, tn)
]
=
1
(∆x)2
[
Ψ(xv+1, tn)−2Ψ(xv, tn)+Ψ(xv−1, tn)
]
(8.5)
73
N
.M
.S
ot
om
ay
or
Embora simples esse esquema pode conduzir a soluc¸o˜es na˜o esta´veis, a seguir
empregamos um algoritmo que trata separadamente as partes real e imagia´ria
da func¸a˜o de onda.
Ψ(x, t) = R(x, t) + iI(x, t), (8.6)
Substitu´ındo 8.6 na equac¸a˜o 8.2 e considerando } = 1 obtemos,
Ĥ(R + iI) = i
d
dt
(R + iI) (8.7)
Igualando as partes real e imagina´ria obtem-se,
ĤR(x, t) = − d
dt
I(x, t) (8.8)
ĤI(x, t) =
d
dt
R(x, t) (8.9)
com
Ĥ ≡ − 1
2m
∂2
∂x2
+ V (x). (8.10)
As equac¸o˜es 8.8 podem ser nume´ricamente resolvidas atrave´s do emprego do
me´todo de integrac¸a˜o de passo me´dio, neste caso as equac¸o˜es sa˜o:
R(x, t+∆t)−R(x, t)
∆t
= ĤI(x, t+
1
2
∆t) (8.11)
I(x, t+ 3
2
∆t)− I(x, t+ 1
2
∆t)
∆t
= −ĤR(x, t) (8.12)
onde, R(x, 0) e I(x, 1
2
∆t) sa˜o os valores iniciais. Para a obtenc¸a˜o da evoluc¸a˜o
temporal e´ escolhido um pacote de onda inicial (t = 0) com perfil gaussiano,
Ψ(x, 0) =
(
1
2piω2
) 1
4
eik0(x−x0)e−(x−x0)
2/4ω2 (8.13)
A figura 8.1 mostra o programa fonte para o ca´lculo da evoluc¸a˜o temporal
do pacote de onda gaussiano da equac¸a˜o 8.13. Mostra-se apenas a parte
do programa principal, observa-se que sa˜o empregadas quatro subrotinas, a
primeira (INITIALDATA) para a introduc¸a˜o dos dados inicias, a segunda
(CLEANMATRIX) para a depurac¸a˜o das varia´veis vetoriais empregadas, a
terceira (WAVEPACKET) para a construc¸a˜o do pacote de onda inicial e a
quarta (EVOLVE) para a integrac¸a˜o nume´rica das duas equac¸o˜es diferenciais.
74
N
.M
.S
ot
om
ay
or
Figura 8.1: Programa fonte em fortran 77 para o ca´lculo nume´rico da evoluc¸a˜o
temporal de um pacote de onda gaussiano atrave´s da equac¸a˜o de Schro¨dinger
empregando-se diferenc¸as finitas e o me´todo de passo me´dio.
O co´digo declara treˆs varia´veis vetoriais para armazenar os resultados nume´ricos
das partes real e imagina´ria assim como o mo´dulo da func¸a˜o de onda. Treˆs
arquivos de sa´ıda sao˜ criados “real.dat”, “imag.dat”e “prob.dat”que arma-
zenam as partes real, imagina´ria e a densidade de probabilidade respectiva-
mente.
A subrotina INITIALDATA ingressa os valores iniciais do programa. Valo-
res como o potencial (V0), intervalo de integrac¸a˜o (XMIN,XMAX), passo
espacial DX, passo temporal DT , sa˜o definidos nesta parte.
A subrotina CLEANMATRIX e´ chamada para zerar inicialmente os vetores
que devera˜o armazenar os resultados nume´ricos das partes real, imagina´ria e
densidade de probabilidade.
A subrotina WAVEPACKET e´ chamada para construir o pacote de onda
75
N
.M
.S
ot
om
ay
or
Figura 8.2: A subrotina INITIALDATA e´ invocada pelo programa principal para
gerar os dados inicias que sera˜o empregados em todo o processo.
Figura 8.3: A subrotina CLEANMATRIX e´ invocada pelo programa principal para
zerar inicialmente os vetores a serem empregados no programa.
inicial Ψ(x, t = 0) em acordo com a equac¸a˜o 8.13.
A subrotina EVOLVE e´ invocada pelo programa principal para realizar a
integrac¸a˜o nume´rica das equac¸o˜es 8.11. As partes real e imagina´ria sa˜o in-
tegradas separadamente de forma simples atrave´s de um processo iterativo
usando-se o me´todo de passo me´dio.
Finalmente e´ empregada tambe´m uma func¸a˜o fortran para definir o poten-
cial. O potencial e´ zero de XMIN ate´ x = 0 a partir de x > 0 o potencial
toma o valor V = V0 de forma que ha´ umprocesso de espalhamento por
uma barreira de potencial no ponto x = 0. O valor do potencial pode ser
alterado. Quando o programa e´ executado gera um conjunto de valores para
a parte real, para a parte imagina´ria e para a densidade de probabilidade
para cada passo de tempo. Cada func¸a˜o compreende um conjunto de valores
76
N
.M
.S
ot
om
ay
or
Figura 8.4: A subrotina WAVEPACKET e´ invocada pelo programa principal
para construir o pacote de onda gaussiano inicial (em t = 0).
Figura 8.5: A subrotina EVOLVE e´ invocada pelo programa principal para
integrar nume´ricamente as equac¸o˜es diferenciais para as partes real e ima-
gina´ria do pacote de onda.
nume´ricos no intervalo de [XMIN,XMAX], no caso do exemplo sa˜o 200
pares de dados para cada func¸a˜o. Isto pode aumentar ou diminuir depen-
dendo da presic¸a˜o desejada. Se sa˜o 1000 intervalos de tempo DT sa˜o 600000
dados para as treˆs func¸o˜es. Certamente, se tentamos realizar o gra´fico desses
dados sera´ imposs´ıvel visualizar alguma evoluc¸a˜o teremos uma superposic¸a˜o
muito grande de dados que impossibilitam qualquer ana´lise. O Gnuplot e´
um aplicativo do Projeto GNU da Free Software Foundation que permite o
gra´fico sincronizado e ainda a criac¸a˜o de uma animac¸a˜o em formato GIF.
77
N
.M
.S
ot
om
ay
or
Figura 8.6: A func¸ao˜ V e´ empregada para definir o potencial.
Para isto e´ necessa´rio criar um arquivo de texto com o conteu´do mostrado a
seguir,
set term gif animate delay 0.1 size 1450, 780
set output "onda.gif"
do for [n=1:1001] {
plot [-20:0][-1:1] "prob.dat" u 1:2 every :n::n::n w l t
sprintf("n=%i", n), "real.dat" u 1:2 every :n::n::n w l
smooth csplines,"imag.dat" u 1:2 every :n::n::n w l smooth
csplines }
O arquivo tem que ser salvo com a extensa˜o .gp, exempo ca.gp. Para
executar o arquivo e´ necessa´rio digitar no prompt de comando do Linux
gnuplot ca.gp. O gnuplot vai ler os dados das func¸o˜es e vai gerar uma
a˜presentac¸a˜o sequeˆncial dos dados salvando eles em um arquivo de animac¸a˜o
de nome “onda.gif”. Parte da evoluc¸a˜o pode ser vista na figura 8.7.
78
N
.M
.S
ot
om
ay
or
-20 -10 0 10 20
-0,6
-0,4
-0,2
0,0
0,2
0,4
0,6
x
-20 -10 0 10 20
-0,6
-0,4
-0,2
0,0
0,2
0,4
0,6
x
-20 -10 0 10 20
-0,6
-0,4
-0,2
0,0
0,2
0,4
0,6
x
Figura 8.7: Gra´fico das partes real, imagina´ria e densidade de probabilidade
do pacote de onda gaussiano evoluindo em func¸a˜o do tempo.
79
N
.M
.S
ot
om
ay
or
Refereˆncias Bibliogra´ficas
[1] Wolfram mathworld, web mathematics resource, 24/10/2007.
http://mathworld.wolfram.com/OrdinaryDifferentialEquation.html
[2] Harvey Gould, Jan Tobochnik, An Introduction to Computer Simu-
lation Methods, Application to Physical Systems, Second Edition,
Addison-Wesley Publishing Co. USA 1996.
80