FORTRAN E QTIPLOT: SOLUÇÃO NUMÉRICA DE EQUAÇÕES DE ... · clássica usando os métodos de...

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.