Upload
vuonganh
View
217
Download
1
Embed Size (px)
Citation preview
FORTRAN E QTIPLOT:
SOLUÇÃO NUMÉRICA DE EQUAÇÕES DE MOVIMENTO
Edio Cunha da Costa
Conteúdo Programático
1. Equações diferenciais ordinárias
(a) Equações diferenciais ordinárias na Mecânica Clássica: 2a
Lei de Newton
i. MRU
ii. MRUV
iii. Movimento plano
iv. Força resistiva proporcional a velocidade
v. MHS
vi. MHS amortecido
vii. MHS forçado
viii. MHS acoplado
ix. Gravitação
x. Partículas carregadas em campos eletromagnéticos var-
iáveis
(b) Soluções numéricas
i. Método de Euler
ii. Método de Runge Kutta (4a ordem)
2. Introdução ao FORTRAN
(a) Histórico
(b) FORTRAN 77 e FORTRAN 90/95
(c) Programação em FORTRAN
i. Declaração de variáveis (integer, real, double precision)
ii. Operações matemáticas (+ - * / **)
iii. Funções (sqrt, log, exp, sin, cos, tan)
iv. Entrada e saída de dados (read, write, open, close, for-
mat)
v. Comparações (maior que, menor que, ...)
vi. Loop (do, do while, enddo)
vii. Condições (if-then-else)
3. Algoritmo e programas para resolver os problemas de mecânica
clássica usando os métodos de Euler e Runge Kutta
4. Grá�cos e análise dos resultados: QtiPlot
(a) Barra de menu/ícones
(b) Importar dados numéricos
(c) Opções da tabela
(d) Opções do grá�co
(e) Análise numérica
1. Equações Diferenciais Ordinárias
• Equações que envolvem uma função e suas derivadas or-dinárias (função de uma variável)
• Podem ser linear, não-linear, com coe�cientes constantes,coe�cientes variáveis, homogêneas, etc. A solução é umafunção que satisfaz a equação diferencial
dx
dt= Ax2 x2dy
dx= Ax+By
d2u
dt2=du
dt+ tu
(dv
dt
)2= v
A
x
d2y
dx2=dy
dxAd2x
dt2+B
dx
dt+ Cx = 0
• As constantes de integração estão dadas pelas condições ini-ciais ou pelas condições de contorno
• Um equação diferencial ordinária de ordem n pode ser escrita
como um sistema de n equações diferenciais de primeira or-
dem. Por exemplo seja a equação linear de segunda ordem
Ad2x
dt2+B
dx
dt+ Cx = 0 (1)
Se de�nirmos a variável u tal que u =dx
dt, a equação (1)
poderá ser escrita como o sistema de equações de primeira
ordem
Adu
dt+Bu+ Cx = 0
dx
dt= u
(a) Equações diferenciais ordinárias na Mecânica Clássica:
2a Lei de Newton
• O movimento de uma partícula é determinado pelas forçasque atuam na mesma:
d~p
dt=∑
~F (2)
em que m é a massa da partícula, ~p = m~v é o momento linear e∑ ~F é a soma (vetorial) de todas as forças que atuam na referidapartícula. Da equação (2) obtemos três equações diferenciais,em geral, acopladas
dpx
dt=∑Fx m
dvx
dt=∑Fx
dpy
dt=∑Fy ou m
dvy
dt=∑Fy
dpz
dt=∑Fz m
dvz
dt=∑Fz
e nem sempre de fácil solução.
• Para um sistema físico constituído por várias partículas, um
sistema de equações diferenciais deste tipo deve ser escrito
para cada uma delas.
• As equações de movimento obtidas a partir da 2a Lei de
Newton são equações diferenciais de segunda ordem para a
posição da partícula. Assim, esta equação pode também
ser escrita como um sistema de duas equações de primeira
ordem, uma para a velocidade e uma para a posição:
d~v
dt=
∑ ~F
me
d~r
dt= ~v
i. MRU
Se∑Fx = 0, então
mdvx
dt= 0 e
dx
dt= vx (3)
ii. MRUV
Se∑Fx = cte então
mdvx
dt= cte e
dx
dt= vx (4)
iii. Movimento no plano
Seja ~F a força (ou resultante das forças) que atua(m) na partícula.
Então∑Fx(t) = f(t),
∑Fy(t) = g(t) e
dvx
dt=
f(t)
me
dx
dt= vx (5)
dvy
dt=
g(t)
me
dy
dt= vy (6)
iv. Força resistiva proporcional a velocidade
Seja uma partícula em queda, sujeita a ação das forças peso e
resistência do ar, dadas respectivamente por P = mg e Far = αvy.
Com um referencial orientado positivamente para cima, a 2a Lei
de Newton �ca
dvy
dt= −g +
α
mvy e
dy
dt= vy (7)
v. MHS
Seja uma partícula de massa m sujeita a ação de uma força
elástica, dada pela Lei de Hooke. A 2a Lei de Newton �ca
dvx
dt= −
k
mx e
dx
dt= vx (8)
vi. MHS amortecido
Seja o MHS anterior, agora, sujeito a uma força de atrito, pro-
porcional a velocidade, dada por αv. A equação (8) �ca escrita
comodvx
dt= −
k
mx−
α
mvx e
dx
dt= vx (9)
vii. MHS forçado
O mesmo MHS anterior pode estar sujeito a uma força externa
(Fext), em geral variável. Temos
dvx
dt= −
k
mx−
α
mvx +
1
mFext(t) (10)
viii. MHS acoplado
Sejam dois MHS's acoplados por uma mola. Neste caso teremos
um sistema de equações acopladas.
mAdvAdt
= k2xB − (k1 + k2)xA (11)
mBdvBdt
= k2xA − (k3 + k2)xB (12)
dxAdt
= vA (13)
dxBdt
= vB (14)
ix. Gravitação
Consideremos duas partículas interagindo gravitacionalmente. Asequações de movimento para o sistema são
mAd~vAdt
=GmAmB
r3~r e mB
d~vBdt
= −GmAmB
r3~r (15)
comd ~rAdt
= ~vA ed~rBdt
= ~vB (16)
Observe que o sistema de equações (15) nos dá um sistema comseis equações escalares, a saber
mAdvAxdt
=GmAmB
r3 (xB − xA)
mAdvAy
dt=
GmAmB
r3 (yB − yA)
mAdvAzdt
=GmAmB
r3 (zB − zA)
mBdvBxdt
= −GmAmB
r3 (xB − xA)
mBdvBy
dt= −
GmAmB
r3 (yB − yA)
mBdvBzdt
= −GmAmB
r3 (zB − zA)
Nestas equações
r =√
(xB − xA)2 + (yB − yA)2 + (zB − zA)2
é a distância entre as partículas.
x. Partículas carregadas em campos eletromagnéticos variáveis
Seja uma partícula de massa m e carga elétrica q movendo-se
num campo eletromagnético, dado pelos campos ~E e ~B. Então
md~v
dt= q~v × ~B + q ~E (17)
ou seja
mdvx
dt= q (vyBz − vzBy) + qEx (18)
mdvy
dt= q (vzBx − vxBz) + qEy (19)
mdvz
dt= q (vxBy − vyBx) + qEz (20)
comdx
dt= vx,
dy
dt= vy e
dz
dt= vz. O sistema acima é um problema
mais complexo porque os campos ~E e ~B, em geral, dependem
da posição e do tempo.
(b) Soluções numéricas
i. Método de Euler
Consideremos y = y(x),dy
dx= f(x, y) com y(xo) = yo . A
de�nição de derivada é
dy
dx= lim
∆x→0
y(x+ ∆x)− y(x)
∆x
Assim, podemos aproximar f(x, y) ≈y(x+ ∆x)− y(x)
∆x, de modo
que, considerando valores pequenos de ∆x, escrevemos
y(x+ ∆x) = y(x) + f(x, y)∆x (21)
Assim, tendo o valor inicial y(xo) para ∆x = 0, obtemos os
valores de y(x) para cada valor subsequente de x, incrementado
por ∆x, ou seja
x1 = xo + ∆x y(x1) = y(xo) + f(xo, yo)∆xx2 = x1 + ∆x y(x2) = y(x1) + f(x1, y1)∆xx3 = x2 + ∆x y(x3) = y(x2) + f(x2, y2)∆x
· · · · · ·xn+1 = xn + ∆x y(xn+1) = y(xn) + f(xn, yn)∆x
ii. Método de Runge Kutta (4a ordem)
O método de Runge Kutta é um método mais preciso e e�-ciente em relação ao método de Euler. Consideremos a equaçãodiferencial
dy
dx= f(x, y), y(xo) = yo
A solução desta equação é aproximada por
y(xn+1) = y(xn) +1
6(k1 + 2k2 + 2k3 + k4) (22)
sendo os coe�cientes dados por
k1 = ∆x f(xn, yn)
k2 = ∆x f(xn +∆x
2, yn +
k1
2)
k3 = ∆x f(xn +∆x
2, yn +
k2
2)
k4 = ∆x f(xn + ∆x, yn + k3)
e
xn+1 = xn + ∆x
2. Introdução ao FORTRAN
(a) Histórico
- O FORTRAN foi a primeira linguagem de programação de alto
nível, criada em 1956.
- Já sofreu várias modi�cações para permitir a incorporação de
recursos de programação cada vez mais modernos.
- Em 1966 o American National Standards Institute (ANSI)
padronizou o FORTRAN IV em duas versões (básico e avançado)
com recursos de operações lógicas, operações básicas sobre números
complexos, comandos de entrada e saída, entre outros.
- Em 1977, o ANSI padronizou o FORTRAN 77, que inclui recur-
sos de programação estruturada, de operações sobre caracteres,
manipulação e uso de arquivos, entre outros.
- De lá pra cá foram criadas outras versões, dentre elas o FOR-
TRAN 90 e o FORTRAN 95.
(b) FORTRAN 77 e FORTRAN 90/95
Existem várias diferenças entre estas versões. Por exemplo no
FORTRAN 77 as linhas do programa fonte devem começar na
coluna 7 e terminar na coluna 72 enquanto o FORTRAN 95 não
precisa; o FORTRAN 95 possui recursos para otimizar operações
com matrizes.
(c) Programação em FORTRAN
Para programar em FORTRAN:
• escrevemos o �código fonte� em qualquer editor de texto,
salvando o arquivo como .for ou .f90 ou .f95, por exemplo.
• numa janela de terminal compilamos o código fonte para
criar o arquivo executável. Isto é feito através do comando
gfortran arquivo.f95
que criará, no diretório em que foi salvo o código fonte, um
arquivo a.out executável. Podemos alterar o nome do arquivo
executável com o comando
gfortran arquivo.f95 -o executavel
e podemos ainda usar a opção -W (warning, cuidado) para que o
compilador informe algum possível erro ou inconsistência. Neste
caso, usamos
gfortran -W arquivo.f95 -o executavel
• para executar o programa basta digitar na linha de comando
o nome do arquvi executável, precedido de ./
./a.out ./executavel
• se após o nome do arquivo executável digitarmos & o pro-
grama será executado em �segundo plano�
i. Declaração de variáveis:
• inteiras (integer a,b(10),c1,d(3,5)) (integer :: a,b(10),c1,d(3,5))
• reais (real a, b(10), c1, d(3,5),var) (real :: a, b(10), c1,
d(3,5),var)
• dupla precisão (double precision)
• complexas (complex)
• lógicas (logical)
• caractere (character)
Vetores e matrizes podem também ser declarados pelo comando
dimension
• real, dimension (3,2) :: a
• real, dimension (1:8) :: b
ii. Operações matemáticas
• adição +
• subtração -
• produto *
• divisão /
• potência **
iii. Funções
• sqrt(x) (raiz quadrada)
• cos(x) (cosseno-em radiano)
• cosh(x) (cosseno hiperbólico)
• sin(x) (seno-em radiano)
• sinh(x) (seno hiperbólico)
• tan(x) (tangente)
• tanh(x) (tangente hiperbólica)
• exp(x) (função exponencial)
• log(x) (logarítmo natural)
• log10(x) (logarítmo na base 10)
iv. Entrada e saída de dados (read, write, open, close, format)
read(stream, label [, end=end][, err=err]) list
write(stream, label) list
• stream é um número que referencia um arquivo, ou uma
variável caracter ou * para entrada via teclado
• label é o número de um formato de leitura dos dados ou *
pra livre de formato
• list lista dos itens a serem lidos, separados por vírgula, pos-
sivelmente textos colocados entre aspas
• [, end=end][, err=err] é opcional
label format (format descriptors)
• label é o inteiro referenciado pelo read
• format descriptions é a lista de itens separados por vírgula
que descrevem como os dados serão lidos/gravados. Os for-
matos são:
� nIw para saída de inteiros
� nFw.d saída de real ou complexo na forma de ponto �xo
� nEw.d saída de real ou complexo na forma de ponto �u-
tuante
� n contador opcional (quantos itens de entrada e saída)
� w número de caracter por número, incluindo sinais e es-paço
� d número de digitos decimais de w
open([unit=]stream, err=escape, action=action,
�le=name)
• stream identi�cador vinculado ao read ou write
• action umas das opções read, write ou readwrite
• �le nome do arquivo, entre aspas, no qual os dados serãolidos e/ou gravados
close(stream)
para fechar o arquivo indicado em read, write, open
v. Comparações (maior que, menor que, ...)
FORTRAN 77: .gt. .ge. .lt. .le. .eq. .ne.
FORTRAN 95: > >= < <= == /=
vi. Loop (do, do while, enddo)
Podemos repetir uma sequencia de comandos e instruções us-
ando loop que tem a estrutura
do var=inicio, �nal, passo
xxx
end do
var é uma variável inteira
passo é opcional; indica o incremento da variável var
[nome:] do var = inicio, �nal, passo
xxx
end do [nome]
ou
[nome:] do while (expressão lógica)
xxx
end do [nome]
vi. Condições (if-then-else)
if (expressão lógica) then
...
else if (expressão lógica) then
...
else if (expressão lógica) then
...
else
...
end if
3. Algoritmo e programas para resolver os problemas de
mecânica clássica usando os métodos de Euler e Runge
Kutta
O algoritmo geral para escrevermos programas para a solução
numérica de equações de movimento pode ser assim escrito:
• início do programa
• declarar e fornecer os valores de m, xo, vo, to, ∆t e parâmetros
(g, µ, α, . . . )
• escrever as forças F1, F2, . . .
• para t = to, tmax, ∆t faça
� calcule os valores de t, v, x, F1, F2, . . .
� escreva os valores num arquivo de dados
• �m para
• �m do programa
Podemos inserir comentários no arquivo fonte (programa) us-
ando o símbolo ! na primeira coluna ou ao �nal de uma linha de
comando
A linguagem apresenta muitos outros recursos que exigiriam
curso mais extenso. No entanto, existem vários livros, man-
uais e refências disponíveis em livrarias e na rede mundial de
computadores, que podem ser consultados para implementar o
código fonte.
Por exemplo, ao invéz da estrutura if-then-else podemos utilizar
a seguinte sequencia
select case (valor da expressão)
case (valor 1)
...
case (valor 2, valor 3)
...
case (valor 4)
...
case default
...
end select
Podem ser escritas subrotinas
programas no editor kate
• O editor kate é um editor de texto que sinaliza diferentes
formas de código fonte como FORTRAN, C, C++, entre
outros.
• Se o terminal konsole estiver instalado no sistema opera-
cional, uma janela do konsole poderá ser aberta interior-
mente ao kate, na qual poderemos compilar o programa fonte
e executar o arquivo executável.
4. Grá�cos e análise dos resultados: QtiPlot
O QtiPlot é um software que, através da leitura dos dados obti-
dos pela integração numérica, permite a construção de grá�-
cos, ajuste de curvas, cálculos sobre os dados como derivadas,
integrais, médias, construção de novas tabelas, exportação degrá�cos em formatos espcí�cos (eps, pdf)etc..., etc..., etc...
Ao longo do minicurso, veremos algumas destas ferramentas dis-postas na
1. Barra de menu/ícones
2. Importação e exportação de dados numéricos
3. Opções da tabela
4. Opções do grá�co
5. Análise numérica no grá�co e na tabela
Referências
1. Linguagem de Programação Estruturada FORTRAN 77.
MAXIMILIAN EMIL HEHL. McGraw-Hill.
2. Self-study guide 2 Programming in Fortran 95. RACHAEL
PADMAN MICHAELMAS. University of Cambridge - Depart-
ment of Physics. Computational Physics.