MRK4 - O canivete suíço

Secção dedicada à linguagem de programação favorita dos quarkianos: Python!

MRK4 - O canivete suíço

Mensagempor jap em Domingo Jan 10, 2010 5:22 pm

Venho aqui apresentar-vos aquele que é conhecido pelo "canivete suíço" dos métodos numéricos para resolver computacionalmente problemas de mecânica - trata-se de um integrador de equações diferenciais acopladas de 1ª ordem pelo método de Runge-Kutta. Não vou explicar aqui o algoritmo deste método (há muitas referências on-line, vejam Runge-Kutta methods), mas como aplicá-lo para resolver alguns problemas em concreto. :wink:

A ideia de base é escrever as equações que determinam a dinâmica do vosso problema como um sistema de equações diferenciais acopladas de 1ª ordem. De uma forma geral, as vossas equações deverão ter a forma:

y^\prime = f(t,y)


onde y representa um conjunto de variáveis dinâmicas do sistema (tipicamente coordenadas e velocidades), y^\prime representa as primeiras derivadas em ordem ao tempo destas coordenadas e t a variável tempo. f representa aqui a função das coordenadas e do tempo de cada uma das primeiras derivadas dessas coordenadas. Reparem que a equação acima é uma forma compacta de representar não uma equação mas sim um conjunto de equações. Por exemplo, se o nosso sistema for descrito por quatro variáveis dinâmicas (digamos, por exemplo, x, y, v_x e v_y), então existem 4 y distintos estas quatro variáveis), e quatro "f" distintos, um para cada coordenada! Pode parecer confuso, mas na realidade é muito simples. :lol: Em linguagem pitónica, pensem e y como uma lista de coordenadas e em f como uma lista de funções...

Exemplo. Suponhamos que pretendemos fazer uma simulação de um oscilador harmónico amortecido a uma dimensão.

A força que actua sobre a partícula tem a forma

F  = - k x - c v


onde k é a constante da "mola" e c o coeficiente da força dissipativa, proporcional à velocidade.

Em termos da aceleração da partícula vem, dividindo pela massa

a = \frac{d^2x}{dt^2} = - \frac{k}{m} x - \frac{c}{m} v


Ora para utilizarmos o "canivete suíço" temos de transformar esta equação diferencial de 2ª ordem num sistema de equações diferenciais de 1ª ordem nas variáveis dinâmicas do sistema. Como proceder? :roll:

Pois bem, basta reparar que

v = \frac{dx}{dt},

pelo que a equação de 2ª ordem acima é equivalente ao seguinte sistema de 2 equações de 1ª ordem nas variáveis x e v:

\frac{dx}{dt} = v

\frac{dv}{dt} = - \frac{k}{m} x - \frac{c}{m} v

Reparem que estas equações têm, ambas a forma

y^\prime = f(y,t)

embora, neste exemplo simples, a variável tempo não aparece explicitamente no segundo membro das equações.

Trascrevendo em linguagem pitónica estas equações, e se colocarmos as variáveis dinâmicas x e v numa lista de nome y, tal que y[0] é a posição e y[1] a velocidade, então podemos definir a lista de funções f da seguinte forma:

Código: Seleccionar Todos
     f[0] = lambda t,y:    y[1]     # dxdt
     f[1] = lambda t,y: -(k/m)*y[0]-(c/m)*y[1]     # dvdt


Para integrar as equações com o canivete suíço basta chamar a função mrk4 do módulo pitónico que eu escrevi e que aqui vos deixo para vocês "brincarem" com ele. A chamada da função faz-se da seguinte forma:

Código: Seleccionar Todos
table = multirk4(f,t0,y,dt,tmax)


onde table vai guardar os valores de y para os tempos t numa grelha que vai de t=0 até t= tmax em passos no tempo de dt.

O ficheiro mrk4.py que aqui vos deixo implementa o exemplo acima.

O resultado é este:

Imagem

Aqui fica o programa:



Será que alguém consegue agora utilizar este programa para simular o movimento da Terra (ou de qualquer outro planeta do sistema solar) à volta do Sol? :D
José António Paixão
Departamento de Física da FCTUC
Avatar do utilizador
jap
Site Admin
Site Admin
 
Mensagens: 6790
Registado: Quinta Nov 09, 2006 9:34 pm
Localização: Univ. de Coimbra

Re: MRK4 - O canivete suíço

Mensagempor jap em Domingo Jan 10, 2010 6:56 pm

Obviamente, para simular o movimento planetário, é preciso substituir esta parte do código

Código: Seleccionar Todos
     f[0] = lambda t,y:    y[1]     # dxdt
     f[1] = lambda t,y: -(k/m)*y[0]-(c/m)*y[1]     # dvdt


pelas equações correspondentes (as leis de Newton) para o movimento do planeta....

Como escrever essas equações? :P
José António Paixão
Departamento de Física da FCTUC
Avatar do utilizador
jap
Site Admin
Site Admin
 
Mensagens: 6790
Registado: Quinta Nov 09, 2006 9:34 pm
Localização: Univ. de Coimbra


Voltar para Pitónica

Quem está ligado

Utilizadores a navegar neste fórum: Nenhum utilizador registado e 1 visitante

cron