O pêndulo gravítico

Fernando Barão, IST - Dep de Física
Agosto 2018

Índice

  1. Experiência
  2. Resolução numérica da equação do pêndulo gravítico

O aparecimento de microcomputadores de baixo custo como é o caso do Raspberry Pi (RPI) e da disponibilização da funcionalidade de leitura fácil de sinais através da sua interface GPIO, permitiu a multiplicação de pequenos projectos experimentais e de análise dados alargado a um público alvo já não tão especializado. Adicionalmente, a utilização da linguagem de programação python e das múltiplas bibliotecas associadas, simplificou também muito os aspectos ligados ao desenvolvimento de aplicações para a leitura e a análise de sinais. Neste workshop faremos uma introdução à montagem de um pêndulo simples, à construção de um detector de passagem (photogate) constituído pelo par LED infravermelho/Fotodíodo e à leitura dos sinais com auxílio da linguagem python gerados pelo fotodíodo e colectados na interface GPIO do RPI. Os sinais de tempo de passagem registados serão analisados e apresentados graficamente com auxílio de bibliotecas do python.

1. Experiência

Introdução

Ir para o Índice

O pêndulo é um sistema mecânico simples que nos permite o estudo da conservação da energia mecânica e fazer a medição do valor da aceleração gravítica local (g). A medição de g é importante para os estudos geológicos (ver gravimetria).
Historicamente, o pêndulo teve um papel muito importante no desenvolvimento científico e foi objecto de estudo por grandes nomes da ciência como Galileu (1564-1642) e Newton (1643-1727). As experiências do plano inclinado e do pêndulo foram cruciais para o desenvolvimento das leis da dinâmica de Newton no século XVII.

O pêndulo simples consiste num fio inextensível de massa desprezável e comprimento L e um corpo de massa m na extremidade. Consideramos em primeira aproximação que a fricção com o ar não não tem influência no movimento do pêndulo.

A equação do movimento do pêndulo deriva-se a partir da segunda lei de Newton na direcção tangencial. \begin{equation} \vec F_t = m \; \vec a_t \end{equation} Assim, considerando a massa em movimento e sendo $\theta$ o ângulo de abertura do pêndulo e ainda o versor $\vec e_{\theta}$ que aponta no sentido crescente do ângulo de abertura, tem-se: \begin{array}{lcl} \vec F_t & = & -m g \sin\theta \; \vec e_{\theta} \\ \vec a_t & = & L \frac{d^2\theta}{dt^2} \; \vec e_{\theta} \end{array} Obtém-se então a equação do movimento, \begin{equation} \frac{d^2\theta}{dt^2} \; + \; \frac{g}{L} \; \sin\theta = 0 \end{equation}

Movimento oscilatório harmónico simples

Uma vez que,
$\sin\theta \simeq \theta \; ( 1 - \frac{\theta^2}{6} + \cdots)$
podemos fazer a aproximação, $\sin\theta \simeq \theta$, desde que o termo $\frac{\theta^2}{6}$ seja desprezável. Isto acontece tipicamente para,
$\frac{\theta^2}{6} < 10^{-3} \quad \Rightarrow \quad \theta \sim 0.1 \; rad$.

Assim, o pêndulo simples colocado a oscilar com um ângulo $\theta$ pequeno ($\theta < 0.1 rad \sim 10^{\circ}$), implica uma força de restituição proporcional ao deslocamento ($m g \sin\theta \simeq m g \theta$), dando origem a um movimento oscilatório harmónico simples, \begin{equation} \frac{d^2\theta}{dt^2} \; + \; \frac{g}{L} \; \theta = 0 \end{equation}

Frequência de oscilação

A solução da equação do oscilador harmónico é dada por, \begin{equation} \theta = A \; \sin( \omega t + B) \end{equation} onde $A, B$ são constantes que se obtêm a partir das condições iniciais do problema, isto é, velocidade e posição, e $\omega$ a frequência angular. \begin{equation} \omega = \sqrt{\frac{g}{L}} \quad \Rightarrow \quad Período = {2\pi} \sqrt{\frac{L}{g}} \end{equation}

A frequẽncia de oscilação na aproximação dos ângulos pequenos, só depende portanto do comprimento do fio.

A energia mecânica

No caso ideal, desprezando-se a dissipação de energia por atritos (ar, contacto do fio) o pêndulo é um sistema conservativo em que a energia mecânica se mantém constante, \begin{equation} E = \frac{1}{2} \; m \; v^2 + m \; g \; z = E_0 \end{equation}

Considerando por exemplo a energia inicial em $t=0$, quando se liberta o pêndulo de uma altura $h$,
$E_0 = m \; g \; h$
obtém-se facilmente a expressão da velocidade linear:
$\frac{1}{2} \; m \; v^2 + m \; g \; z = m \; g \; h \quad \Rightarrow \quad v^2 = 2 \; g \; ( h - z)$

Em particular, o ponto de velociade máxima do pêndulo obtém-se em $z=0$, \begin{equation} v = \sqrt{2 \; g \; h} \end{equation}

Como medir o período e a velocidade do pêndulo?

Ir para o Índice

A medição do tempo de passagem da massa $m$ faz-se com o auxílio de um detector de passagem também conhecido como photogate. Este detector é constituído por uma estrutura em U de plástico impressa numa impressora 3D, que alberga no seu interior um LED emissor de infravermelhos de um lado e no outro lado oposto, um fotodetector. A forma de U do detector permite colocar o emissor de luz infravermelha e o detector em faces opostas. A passagem da massa do pêndulo no interior do U provocará uma interrupção do feixe de luz que está a ser permanentemente colectado pelo fotodetector.

Lista de Componentes para a realização do photogate:

O circuito do emissor de infravermelhos (LED) inclui uma resistência eléctrica em série para limitar a corrente elétrica que atravessa o LED ($I < 20 mA$). O fotodetector de infravermelhos é alimentado com +5V e tem uma saída no estado HIGH (3.3V) quando não existe luz IV a ser colectada.

O fotodetector possui um sinal de saída que podemos ler usando o microcomputador Raspberry Pi (RPI). Vamos usar para tal o conjunto de pins chamado GPIO (General Purpose Input/Output) que o RPI possui.

O conjunto de pins existentes no Raspberry Pi (40 pins) constituem uma interface utilizável entre o computador e o utilizador. Poderemos desta forma, usar o computador para controlar circuitos eléctricos, motores, sensores, etc. que queiramos ligar ao RPi. Há diferentes tipos de pins. Uns providenciam diferença de potencial de forma permanente (+5V, +3.3V, GND) e outros podem ser usados para leitura de sinais (+3.3V) ou saída de um sinal de tensão 3.3V (3V3). Os pins GPIO possuem um número que corresponderá à sua identificação:

O Raspberry Pi disponibiliza a tensão necessária para que o LED infravermelho acenda e para o correcto funcionamento do fotodetector. E além do mais, recebe o sinal de saída do fotodetector que pode estar no estado HIGH ou LOW, dependendo se está a receber luz ou não.

Usaremos a linguagem de programação python para a leitura do sinal de saída do fotodetector.

Coloquemos o photogate com o feixe de luz na altura $z=0$. Assim, quando a massa inicia o atravessamento do feixe de luz infravermelha, o sinal de saída passa a LOW (0) e mantém-se assim enquanto a massa impede a luz de ser detectada. Podemos assim medir o tempo total $\Delta t$ que demora a massa a atravessar o feixe. Conhecendo o diâmetro da massa, $D$, podemos determinar a velocidade do corpo oscilante como sendo, \begin{equation} v = \frac{D}{\Delta t} \end{equation}

O corpo de massa $m$ que usámos são esferas de pesca com um diâmetro de aproximadamente $D=25.9 \pm 0.1$ mm.

Em relação ao período de tempo, que corresponde ao tempo que demora o corpo a realizar um ciclo completo do movimento, podemos escolher por exemplo como referência a posição em que o corpo interrompe o feixe de luz infravermelha. Assistiremos assim às seguintes transições (ver mais à frente detalhe sobre o fotodetector):

Podemos deta forma calcular o tempo de cada período de oscilação.

A utilização do microcomputador Raspebrry Pi na experiência do pêndulo

Ir para o Índice

O Raspberry Pi é um computador com um sistema operativo linux que permite a sua programação usando linguagens de alto nível. Entre estas, o python destaca-se como uma linguagem interpretada e de grande popularidade. Existem um conjunto vasto de bibliotecaS disponível que permite o seu uso de forma eficaz.

No nosso caso vamos querer ler o sinal de saída do fotodetector que para tal ligaremos a um dos pins disponíveis GPIO. Vejamos de seguida o mapa dos pins:
Esquema de dos pins disponíveis no Raspberry Pi.
A leitura do sinal do fotodetector será feita no pin GPIO20.

Vamos utilizar a numeração dos pins associadas ao nome GPIOn. Para podermos utilizar os pins do GPIO, necessitamos de importar a biblioteca RPi.GPIO e definir o pin GPIO20 como entrada.

Lance um novo Terminal e escreva python o que lança o interpretador de python,

  python

Agora vamos importar a biblioteca RPi.GPIO, definir o esquema de numeração dos pins que usaremos e definir ainda GPIO20 como entrada,

  >>> import RPi.GPIO as GPIO
  >>> GPIO.setwarnings(False)
  >>> GPIO.cleanup()
  >>> GPIO.setmode(GPIO.BCM)
  >>> GPIO.setup(20, GPIO.IN, pull_up_down=GPIO.PUD_UP)
Uma pequena nota: o argumento da função setup(),

pull_up_down = GPIO.PUD_UP
é importante porque permite a definição do estado do pin de entrada 20. Diz-se então que este não se encontra num estado flutuante! No nosso caso, quando a luz está a ser detectada pelo fotodetector, o sinal deste é LOW. O material resistivo sensível à luz que compôe o foto detector, possui uma resistência eléctrica que diminui com o aumento da luminosidade detectada. Portanto, quanto mais luz incidente menor é a resistência e daí que a diferença de potencial seja da ordem de zero quando há luz a atravessar. Por outrolado, quando não há luz no fotodetector este está no estado HIGH (grande resistência eléctrica e portanto grande diferença de potencial).
Ao termos colocado o argumento acima na função, estamos a definir a leitura do pin 20 como sendo HIGH (1), quando o fotodetector não tem sinal. O que nos parece lógico porque sabemos que a luz está a chegar a ele!

De seguida podemos colocar o pêndulo a oscilar e lermos o sinal de entrada do pin 20, imprimindo o seu estado no ecran de forma contínua (loop infinito):

  >>> (...)
  >>> while True:
  ...    print(GPIO.input(20))

Por último podemos testar a função da bibliota RPi.GPIO que detecta transições de estados. O que vai acontecer é que o programa ficará parado à espera que a transição ocorra, e somente nesse instante prosseguirá! Há três estados de transição que podemos testar:
GPIO.FALLING ($1 \rightarrow 0$)
GPIO.RISING ($0 \rightarrow 1$)
GPIO.BOTH ($1 \rightarrow 0$ ou $0 \rightarrow 1$)

Faríamos assim para detectar uma transição RISING:

  >>> (...)
  >>> GPIO.wait_for_edge(20, GPIO.RISING)
  >>> (...)


Podemos agora escrever as linhas de código python num ficheiro, cujo nome seja pendulum.py, e aí será mais fácil correr o programa, fazendo no terminal:

  python pendulum.py

Um programa python mínimo para detectar as oscilações do pêndulo

Façamos um programa mínimo que seja capaz de medir o período do pêndulo. Recorde as transições esperadas e enunciadas acima no sinal de entrada, e associaremos uma data temporal a cada transição luz-sombra. O character # é usado em cada início de linha para definir um comentário, com o objectivo de tornar o código mais legível:

ficheiro: pendulum.py

# importar a biblioteca RPI
import RPi.GPIO as GPIO
# importar a biblioteca time que permite obter o tempo decorrido com precisão de microsegundos desde 1 Jan 1970
import time

# --- definir o pin de entrada

pin = 20
GPIO.setwarnings(False)
GPIO.cleanup()
GPIO.setmode(GPIO.BCM)
GPIO.setup(pin, GPIO.IN, pull_up_down=GPIO.PUD_UP)

# --- detectar transições de sinal 1->0
#     o programa pára até a transição acontecer e depois disso registamos o tempo em segundos (com precisão de microsegundos) que decorreu desde 1 jan 1970

P  = 0. # período de tempo
while True:
    GPIO.wait_for_edge(pin, GPIO.FALLING) #blocking
    t0 = time.time()
    GPIO.wait_for_edge(pin, GPIO.FALLING) #blocking
    t1 = time.time()
    GPIO.wait_for_edge(pin, GPIO.FALLING) #blocking
    t2 = time.time()
    #diferencas de tempo
    dt0 = t1-t0
    dt1 = t2-t1
    # periodo
    P = t2-t0
    # terceira transicao 1->0 detectada: periodo
    print("tempo (segs):  t0={0:17.6f} t1={1:17.6f} t2={2:17.6f}  "
          "| diferença de tempo entre transições 1->0 (segs): {3:9.6f} {4:9.6f} "
          "| Período (segs): {5:9.6f}".format(t0, t1, t2, dt0, dt1, P))

O ficheiro com o código python pode ser descarregado aqui.

Exercício: desenvolvimento de um programa python para medição da velocidade

Ir para o Índice

A determinação da velocidade implica a deteção de transições $ 1 \rightarrow 0 $ e $ 0 \rightarrow 1 $ e ainda a medida dos intervalos de tempo entre estas transições. Inspirando-se agora do código anterior, desenvolva um pequeno programa que permita a medida da velocidade da massa.

2. Resolução numérica da equação do pêndulo gravítico

Ir para o Índice

Introdução

Ir para o Índice

A equação a resolver do pêndulo, tal como na generalidade das situações das equações do movimento, é uma equação diferencial ordinária de segunda ordem, \begin{equation} \frac{d^2\theta}{dt^2} \; + \; \frac{g}{L} \; \sin\theta = 0 \end{equation} Esta equação pode-se transformar facilmente num sistema de equações diferenciais de primeira ordem, introduzindo uma nova variável que é a velocidade angular ($\omega=d\theta/dt$), \begin{equation} \begin{cases} \omega = & \frac{d\theta}{dt} \\ \frac{d\omega}{dt} = & \frac{g}{L} \; \sin\theta \end{cases} \end{equation} cujas condições iniciais (em $t=0$) são: \begin{equation} \begin{cases} \theta(0) = & \theta_0 \\ \omega(0) = & 0 \end{cases} \end{equation}

Solução numérica

Ir para o Índice

Consideremos uma equação diferencial ordinária de primeira ordem (ODE) genérica, função do tempo $t$, \begin{equation} \frac{dy}{dt} = f(t, y(t)) \end{equation} A procura da sua solução numérica implica a discretização do tempo em pequenos intervalos de largura $h$ e a determinação da função $y$ em cada ponto discreto do tempo. Tem-se então entre dois instantes temporais $t_n$ e $t_{n+1}$ separados de $h$, \begin{equation} \int_{t_n}^{t_{n+1}} \frac{dy}{dt} \; dt = \int_{t_n}^{t_{n+1}} f(t, y(t)) \; dt \quad \Rightarrow \quad y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} f(t, y(t)) \; dt \end{equation} A forma mais simples e grosseira de calcular o integral à direita que corresponde à área subentendida pela função $f$ (ver na figura, rectângulo a vermelho) é calcular a área do rectângulo associado de lados $h$ e $f(t_n)$. Isto corresponde a fazer uma aproximação de que a função $f$ é constante no pequeno intervalo de integração. Vê-se facilmente da figura o erro que estamos a cometer no cálculo da área, tomando esta aproximação.

Vem então: \begin{equation} y_{n+1} = y_n + h \cdot f(t_n, y(t_n)) \end{equation} Este é o chamado método de Euler, que é muito simples de aplicar mas de precisão limitada.

Um método mais preciso resulta da aproximação da área a calcular como sendo a do trapézio definido na figura a azul. A área resulta assim da soma área do rectângulo já antes utilizada e do triâgulo acima deste, \begin{equation} \int_{t_n}^{t_{n+1}} f(t, y(t)) \; dt \simeq h \cdot f(t_n) + \frac{1}{2} \cdot h \cdot \left[ f(t_{n+1}) - f(t_n) \right] = \frac{1}{2} \cdot h \cdot \left[ f(t_n) + f(t_{n+1}) \right] \end{equation}

Ficamos assim com uma solução mais precisa para a iteracção da solução, \begin{equation} y_{n+1} = y_n + \frac{h}{2} \cdot \left[ f(t_n, y(t_n)) + f(t_{n+1}, y(t_{n+1}) \right] \end{equation} Há aqui no entanto o problema da solução da iteracção envolver o valor a calcular, $y(t_{n+1})$, que não conhecemos. Podemos no entanto fazer a aproximação de Euler para a determinação intercalar de $y(t_{n+1}) \simeq y(t_n) + h f(t_n)$. Ficamos então com a solução, \begin{equation} y_{n+1} = y_n + \frac{h}{2} \cdot \left[ f(t_n, y(t_n)) + f(t_n + h, y_n + h f(t_n, y_n)) \right] \end{equation}

Este é o chamado método de Runge-Kutta de 2a ordem, que possui uma precisão superior ao de Euler e que é suficiente para resolver o pêndulo gravítico.

Solução da equação do pêndulo com algoritmo Runge-Kutta de 2a ordem

Ir para o Índice

Retomemos a equação do movimento do pêndulo, mas desta vez vamos optimizar numericamente a equação. A solução da equação do movimento do pêndulo é uma função periódica, que possui uma escala de tempo associada ao seu período e que depende do problema. Podemos assim trabalhar a equação do movimento de forma a expressá-la em termos de uma variável adimensional $\tau$, tal que,
$ t = T_0 \; \tau$
onde $T_0$ é a escala de tempo do problema. Determinemos esta escala de tempo.

\begin{equation} \frac{1}{T_0^2} \; \frac{d^2\theta}{d\tau^2} \; + \; \frac{g}{L} \; \sin\theta = 0 \quad \Rightarrow \quad \frac{d^2\theta}{d\tau^2} \; + \; T_0^2 \; \frac{g}{L} \; \sin\theta = 0 \end{equation} Para a precisão numérica, é importante manter tanto quanto possível os valores próximos da unidade. Deriva-se assim a escala de tempo, \begin{equation} T_0^2 \; \frac{g}{L} = 1 \quad \Rightarrow \quad T_0 = \sqrt{\frac{g}{L}} \end{equation}

A equção do movimento a resolver fica então muito mais simples, \begin{equation} \frac{d^2\theta}{d\tau^2} \; + \; \sin\theta = 0 \end{equation} Podemos agora escrever o sistema de equações diferenciais de primeira ordem de forma vectorial, escolhendo como variáveis do problema o ângulo $\theta$ e a velocidade angular $\omega$. Teremos assim o vector de variáveis a resolver, \begin{equation} y = \begin{pmatrix} \theta \\ \omega \end{pmatrix} \end{equation} e o vector de funções, \begin{equation} f = \begin{pmatrix} \omega \\ -\sin\theta \end{pmatrix} \end{equation}

e portanto a equação diferencial escreve-se como, \begin{equation} \frac{d}{d\tau} \; \begin{pmatrix} \theta \\ \omega \end{pmatrix} = \begin{pmatrix} \omega \\ -\sin\theta \end{pmatrix} \end{equation} Vindo para a solução de Runge-Kutta 2, \begin{equation} \begin{pmatrix} \theta \\ \omega \end{pmatrix}_{n+1} = \begin{pmatrix} \theta \\ \omega \end{pmatrix}_n + \frac{h}{2} \; \left[ \begin{pmatrix} \omega \\ -\sin\theta \end{pmatrix}_n + \begin{pmatrix} \omega \\ -\sin\theta \end{pmatrix}_{n+1} \right] \end{equation} sendo o termo $f(n+1)$ calculado em $y(t_{n+1}) \simeq y(t_n) + h f(t_n)$.

Algoritmo

Vamos utilizar a linguagem python e em particular as suas bibliotecas:

Adicionalmente existe a biblioteca:

Os passos do algoritmo a implementar para resolver a equação pelo método RK2 resumem-se ao seguinte:

Implementação em python

ficheiro: pendulum_numerical_solution.py

# -*- coding: utf-8 -*-
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pendulum_tools as tools

"""
Determina o movimento (theta, velocidade, tempo) de um pendulo, descrito pela
ODE:
theta''(t) + alpha*sin(theta(t)) = 0

Esta equacao e convertida num sistema de duas equacoes diferenciais de 1a ordem:
theta'(t) = omega
omega'(t) = -alpha*sin(theta(t))

definindo a variavel y = (theta, omega), podemos escrever o sistema como:
y'(t) = (omega, -alpha*sin(theta(t)))

variaveis:
T : tempo total de 0 a T
n : numero de intervalos de tempo
dt : intervalo de tempo
theta0 : angulo inicial
v0: velocidade inicial
alpha : sqrt(g/L)
"""

# --- definicao dos valores iniciais

T = 10. # segundos
n = 10000
dt = T/float(n)
theta0 = np.deg2rad(60.) # 10 graus convertidos para radianos
v0 = 0.

print('dt (segs)={}'.format(dt))

g = 9.806 # m/s^2
L = 0.50  # m
alpha = g/L


# --- definir funcoes ODE

def ODEfuncs(y, time, alpha):
    """
    definicao das funcoes f(y,t): dy/dt = f(y,t)
    """
    f0 = y[1]
    f1 = -alpha*np.sin(y[0])
    return np.array([f0,f1])

# --- metodo de Euler: obter a proxima iteraccao

def euler(y, time, dt, ODEfuncs, alpha):
    """
    metodo de euler
    """
    y_next = y + ODEfuncs( y, time, alpha ) * dt
    return y_next

# --- metodo de RK2: obter a proxima iteraccao

def RK2(y, time, dt, ODEfuncs, alpha):
    """
    metodo de runge-kutta ordem 2
    """
    K1 = dt * ODEfuncs( y, time, alpha)
    K2 = dt * ODEfuncs( y+K1, time+dt, alpha)
    y_next = y + 0.5 * (K1 + K2)
    return y_next

# --- SOLUCAO: ciclo nos intervalos de tempo

t = np.linspace(0,T,n+1)
yeuler = np.zeros((n+1,2)) #array com n+1 linhas e 2 colunas preenchido com zero's
yrk2 = np.zeros((n+1,2))
yeuler[0] = theta0, v0
yrk2[0] = theta0, v0

for k in range(n):
    yeuler[k+1] = euler(yeuler[k], t[k], dt, ODEfuncs, alpha)
    yrk2[k+1]   = RK2(  yrk2[k],   t[k], dt, ODEfuncs, alpha)

# --- find period from zero crossings

period1 = tools.freq_zero_crossing(yeuler[:,0], dt) # recuperar array com valores da 1a coluna (theta)
period2 = tools.freq_zero_crossing(yrk2[:,0], dt)
print('period1={} segs period2={}'.format(period1, period2))

# --- plot resultados

fig1 = plt.figure(1, figsize=(16,12))

ax11 = fig1.add_subplot(221)
ax11.plot(t,yeuler[:,0])
ax11.set_title('euler: angle as function of time')

ax12 = fig1.add_subplot(222)
ax12.plot(t,yeuler[:,1],label='velocity(t)')
ax12.set_title('euler: velocity as function of time')

ax13 = fig1.add_subplot(223)
ax13.plot(t, yrk2[:,0], color='red')
ax13.set_title('rk2: angle as function of time', color='red')

ax14 = fig1.add_subplot(224)
ax14.plot(t, yrk2[:,1], label='velocity(t)', color='red')
ax14.set_title('rk2: velocity as function of time', color='red')

fig2 = plt.figure(2, figsize=(10,6))

ax21 = fig2.add_subplot(111)
ax21.plot(t, yeuler[:,0], color='blue', label='euler')
ax21.plot(t, yrk2[:,0], color='red', label='rk2')
ax21.set_title('angle as function of time')
ax21.legend(loc='upper left')

fig2.savefig('FIG_pendulum_numerical_solution_angle.png')

plt.show()

O ficheiro pode ser obtido aqui.

Mostra-se de seguida a solução obtida usando os métodos de Euler e RK2, para um ângulo inicial $\theta=60$ graus. O período de oscilação é de $1.533$ segundos.

Utilização de bibliotecas numéricas python

Ir para o Índice

A biblioteca scipy pode ser utilizada para a resolução de equações diferenciais.
A função disponível é: scipy.integrate.odeint()

Lancemos uma sessão de python e façamos:

 >>> # importamos o modulo matematico onde esta a funcao odeint()
 >>> from scipy.integrate.odepack import odeint
 >>> # vejamos o help da função
 >>> help(odeint)
Help on function odeint in module scipy.integrate.odepack:

odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
    Integrate a system of ordinary differential equations.

    Solve a system of ordinary differential equations using lsoda from the
    FORTRAN library odepack.

    Solves the initial value problem for stiff or non-stiff systems
    of first order ode-s::

        dy/dt = func(y, t0, ...)

    where y can be a vector.

    *Note*: The first two arguments of ``func(y, t0, ...)`` are in the
    opposite order of the arguments in the system definition function used
    by the `scipy.integrate.ode` class.

    Parameters
    ----------
    func : callable(y, t0, ...)
        Computes the derivative of y at t0.
    y0 : array
        Initial condition on y (can be a vector).
    t : array
        A sequence of time points for which to solve for y.  The initial
        value point should be the first element of this sequence.
    args : tuple, optional
        Extra arguments to pass to function.


Apêndice: o uso da linguagem Python

Ir para o Índice

Definir variáveis em python e utilizar o módulo matemático (math)

No início da primeira lição aprendemos a usar o interpretador de python e concluímos que o mais confortável é criarmos um ficheiro com o auxílio do editor de texto, nano, e introduzirmos lá o código que pretendemos correr. É isso que vamos fazer agora criando o ficheiro t1.py, onde vamos calcular o coseno de um dado valor em radianos. O valor é fornecido ao programa a partir do teclado.

Para editarmos o ficheiro t1.py escrevemos na linha de comando linux,

 $ nano t1.py

e de seguida introduzimos as seguintes linhas de programação python no ficheiro:

  # importamos o modulo matematico onde esta a funcao cos()
  import math
  # associamos a variável x ao valor de entrada / nota: conversão para real feita pela função float()
  x = float(input("angulo: "))
  # calculamos o coseno
  cx = math.cos(x)
  # imprimimos uma mensagem que seja clara
  print("O coseno do angulo ", x, " em radianos é = ", cx )

E finalmente para corrermos o ficheiro e vermos o resultado,

 $ python3 t1.py
 angulo: 2.3
 O coseno do angulo  2.3  em radianos é =  -0.6662760212798241

Definir ciclos e testar condições

Vamos agora criar um programa python que seja capaz de somar os números inteiros de 1 a 100. No final, queremos testar se o resultado é superior ao número 1000 ou não.

  # inicilaizar a variavel que contem a soma
  soma = 0
  for x in range(1,101):
     soma += x
  if soma>1000:
     print(soma)

Definir funções e módulos do utilizador

Em python as funções que podem retornar algo ou não, criam-se com def. Vamos criar um ficheiro com o nome pyfun.py, com o auxílio do editor nano, que contenha todas as funções que vamos desenvolver neste estágio. Poderemos assim criar o nosso próprio módulo pyfun que podemos importar (import) em programas python.

Vamos então editar o pyfun.py integrando lá as seguintes funções:

 import math

 def cos(x):
   x = float(input("angulo: "))
   cx = math.cos(x)
   print("O coseno do angulo ", x, " em radianos é = ", cx )
   return cx

 def soma(x,y):
   c = x + y
   print("a soma de ", x, " e de ", y, " dá a soma ", c)
   return c

Para utilizarmos as funções definidas no módulo pyfun basta fazermos o seguinte programa teste, t2.py:

 import math
 import pyfun

 x = float(input("angulo: "))
 pyfun.cos(x)
 x,y = float(input("soma de a,b :"))
 pyfun.soma(x,y)


Fernando Barão
Instituto Superior Técnico - Dep. de Física
Agosto 2018