Código em Python para simulações de osciladores lineares e FPUT quadrático: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
Sem resumo de edição
 
(8 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 61: Linha 61:
def modo_de_osc(poseq, l, m, am, N):
def modo_de_osc(poseq, l, m, am, N):
   return poseq + am*np.sin(np.pi*(l+1)*m/(N+1))
   return poseq + am*np.sin(np.pi*(l+1)*m/(N+1))
#Só para corda do FPUT, não precisa de flag
def condicao_inicial_1d_vertical(L, N, m):         
  pos = np.zeros((N,1), dtype="float")
  poseq = np.zeros((N,1), dtype="float")
  for i in range(0, N):
    pos[i] = modo_de_osc(poseq[i], i, m, 0.5, N)         
 
  return pos


def condicao_inicial_2d(N):
def condicao_inicial_2d(N):
Linha 178: Linha 188:


   return figname
   return figname
def energia_do_modo(x, x_old, N, m, dt):
  #Energia do modo m
  omega = 2*np.sin(m*np.pi/(2*(N+1)))
  omega2 = omega**2
  tmp = np.arange(1, N+1, 1)
  tmp[-1:] = 0
  pos = np.reshape(x, (1,N))
  pos_old = np.reshape(x_old, (1,N))
   
  produto_escalar = np.sqrt(2.0 / (N+1)) * np.dot(pos, np.sin( tmp * m * np.pi / (N+1) ))
  produto_escalar_ant = np.sqrt(2.0 / (N+1)) * np.dot(pos_old, np.sin( ( tmp * m * np.pi / (N+1) ) ))
   
  return 0.5 * (produto_escalar - produto_escalar_ant)**2 / dt**2 + (0.5 * omega2 * produto_escalar**2)
</source>
Imprimir os cinco primeiros modos de oscilação para N=32 partículas (Figura 12 da Wiki).
<source lang="python">
N=32
modo = np.zeros((5,N), dtype="float")
nodes = [i for i in range(1,33)]
for i in range(0,5):
  for j in range(0, N):
    modo[i,j] = modo[i,j] + np.sin((j+1)*np.pi*(i+1)/(N+1))
 
plt.plot(nodes, modo[0,:], "k.", label="Modo 1")
plt.plot(nodes, modo[1,:], "r-.", label="Modo 2")
plt.plot(nodes, modo[2,:], "b-", label="Modo 3")
plt.plot(nodes, modo[3,:], "g--", label="Modo 4")
plt.plot(nodes, modo[4,:], "y:", label="Modo 5")
plt.xlabel("Número da partícula")
plt.ylabel("Deslocamento em relação à posição de equilíbrio")
plt.legend()
plt.show()
</source>
</source>
<gallery class="center" mode=packed heights=250px>
modos.png
</gallery>


Exemplo de como podem ser feitos os gráficos dos modos e como comparar com as frequências esperadas. Também estou calculando a energia total do sistema. O exemplo é do modo 1, do caso N=3. Note que a energia não é exatamente conservada, é uma dificuldade do método de Verlet. Para N=2 a conservação é bem obedecida, uma reta próxima a 0.2. Seria possível melhorar diminuindo o passo de integração, mas preferi deixar assim para ilustrar.
Exemplo de como podem ser feitos os gráficos dos modos e como comparar com as frequências esperadas (gráficos 2 a 11 da Wiki). Também estou calculando a energia total do sistema. O exemplo é do modo 1, do caso N=3. Note que a energia não é exatamente conservada, é uma dificuldade do método velocity Verlet. Para N=2 a conservação é bem obedecida, uma reta próxima a 0.2. Seria possível melhorar diminuindo o passo de integração, mas preferi deixar assim para ilustrar.


<source lang="python">
<source lang="python">
Linha 268: Linha 323:
compara_simetricoN3.png ‎
compara_simetricoN3.png ‎
</gallery>
</gallery>
'''FPUT unidimensional, deslocamento horizontal:'''
Setar alpha diferente de zero e flag=4 (último parâmetro das funções condicao_inicial(), plota_figura_1d() e velocity_verlet()).
Salva animação no arquivo fput.gif. Bem grande o gif gerado.
As animações não têm linhas representando as molas e o tamanho das partículas é bem pequeno. Fiz isso de propósito pra visualizar melhor o caso com N=32 partículas, do artigo da RBF.
'''Lembrar de alterar as variáveis globais de acordo. Diferentemente do caso bidimensional, o valor de L tem q ser fornecido e deve ser obviamente diferente de zero.'''
<source lang="python">
alpha = 0.25
k = 1.0
massa = 1.0
L = 33
N = 32
images=[]
t=0.0
dt = 0.1
tmax = 30
posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")
#Por enquanto não estou usando posições em Y
posY = np.zeros((N,1), dtype="float")
velX = np.zeros((N,1), dtype="float")
acelX = np.zeros((N,1), dtype="float")
posX, posXeq = condicao_inicial_1d(L, N, 1, 4)
while(t < tmax):
  figname = plota_figura_1d(L, posX, posY, t, 4)
  images.append(imageio.imread(figname))
  os.remove(figname)
 
  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 4)
  t += dt
imageio.mimsave('fput.gif', images, fps=30)
</source>
FPUT unidimensional, deslocamento vertical (como no programa em C):
Setar alpha diferente de zero e flag=4 (último parâmetro das funções condicao_inicial(), plota_figura_1d() e velocity_verlet()).
<source lang="python">
alpha = 0.25
k = 1.0
massa = 1.0
L = 33
N = 32
images=[]
t=0.0
dt = 0.1
tmax = 100
posX = np.zeros((N,1), dtype="float")
posY = np.zeros((N,1), dtype="float")
posYeq = np.zeros((N,1), dtype="float")
velY = np.zeros((N,1), dtype="float")
acelY = np.zeros((N,1), dtype="float")
step = L/(N+1)
for i in range(0, N):
  posX[i] = (i+1)*step
posY = condicao_inicial_1d_vertical(L, N, 1)
while(t < tmax):
  figname = plota_figura_1d(L, posX, posY, t, 4)
  images.append(imageio.imread(figname))
  os.remove(figname)
 
  posY, velY, acelY = velocity_verlet(posY, posYeq, velY, acelY, k, massa, alpha, 4)
  t += dt
imageio.mimsave('fput.gif', images, fps=30)
</source>
Cálculo e gráficos das energias por modo do FPUT quadrático:
<source lang=python>
alpha = 0.25
k = 1.0
massa = 1.0
L = 33
N = 32
images=[]
t=0.0
dt = 0.1
tmax = 14000
energia1 = []
energia2 = []
energia3 = []
energia4 = []
energia5 = []
posX = np.zeros((N,1), dtype="float")
posY = np.zeros((N,1), dtype="float")
posYeq = np.zeros((N,1), dtype="float")
velY = np.zeros((N,1), dtype="float")
acelY = np.zeros((N,1), dtype="float")
step = L/(N+1)
for i in range(0, N):
  posX[i] = (i+1)*step
posY = condicao_inicial_1d_vertical(L, N, 1)
while(t < tmax):
  pos_old = posY
 
  posY, velY, acelY = velocity_verlet(posY, posYeq, velY, acelY, k, massa, alpha, 4)
  energia1.append(energia_do_modo(posY, pos_old, N, 1, dt)) 
  energia2.append(energia_do_modo(posY, pos_old, N, 2, dt))
  energia3.append(energia_do_modo(posY, pos_old, N, 3, dt))
  energia4.append(energia_do_modo(posY, pos_old, N, 4, dt))
  energia5.append(energia_do_modo(posY, pos_old, N, 5, dt))
  t += dt
</source>
<source lang=python>
time = np.linspace(0, tmax, int(tmax/dt))
plt.figure(figsize=(15,7))
plt.title('Energias por modo de vibração')
plt.xlabel("t")
plt.ylabel("Energia (u.a.)")
plt.plot(time, energia1,label="Modo 1")
plt.plot(time, energia2,label="Modo 2")
plt.plot(time, energia3,label="Modo 3")
plt.plot(time, energia4,label="Modo 4")
plt.plot(time, energia5,label="Modo 5")
plt.legend()
</source>
<gallery class="center" mode=packed heights=250px>
energias_por_modo.png ‎
</gallery>
Osciladores acoplados lineares unidimensionais:
Salva animação no arquivo OA1D-N#.gif, com # igual ao número de partículas da simulação. Bem grande o gif gerado. Fica ainda maior quanto mais partículas tiver a simulação. '''Lembrar de alterar as variáveis globais de acordo. Diferentemente do caso bidimensional, o valor de L tem q ser fornecido e deve ser obviamente diferente de zero.'''
<source lang=python>
alpha = 0.0
k = 1.0
massa = 1.0
L = 12
N = 4
images=[]
t=0.0
dt = 0.1
tmax = 50
posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")
posY = np.zeros((N,1), dtype="float")
velX = np.zeros((N,1), dtype="float")
acelX = np.zeros((N,1), dtype="float")
 
posX, posXeq = condicao_inicial_1d(L, N, 1, 3)
while(t < tmax):
  figname = plota_figura_1d(L, posX, posY, t, 3)
  images.append(imageio.imread(figname))
  os.remove(figname)
 
  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 3)
 
  t += dt
imageio.mimsave('OA1D-N{:.0f}.gif'.format(N), images, fps=60)
</source>
Osciladores acoplados em uma rede bidimensional L x L (Aqui o L é calculado de maneira a ter um grid quadrado com espaçamento 1 entre as partículas, só precisa ser fornecida a quantidade de partículas N). '''N TEM QUE SER UM QUADRADO PERFEITO.''' LEMBRAR DE ALTERAR AS VARIÁVEIS GLOBAIS, FAZER L=0.0 E N = 4, 9, 16, 25,..., ETC.
<source lang=python>
alpha = 0.0
k = 1.0
massa = 6.0
L = 0.0
N = 49
images=[]
t=0.0
dt = 0.1
tmax = 30
posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")
posY = np.zeros((N,1), dtype="float")
posYeq = np.zeros((N,1), dtype="float")
velX = np.zeros((N,1), dtype="float")
velY = np.zeros((N,1), dtype="float")
acelX = np.zeros((N,1), dtype="float")
acelY = np.zeros((N,1), dtype="float")
 
L, posX, posY, posXeq, posYeq = condicao_inicial_2d(N)
while(t < tmax):
  figname = plota_figura_2d(L, posX, posY, t, 3)
  images.append(imageio.imread(figname))
  os.remove(figname)
  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 3)
  posY, velY, acelY = velocity_verlet(posY, posYeq, velY, acelY, k, massa, alpha, 3)
  t += dt
imageio.mimsave('OA2D-N{:.0f}.gif'.format(N), images, fps=60)
</source>

Edição atual tal como às 17h14min de 4 de maio de 2022

Bibliotecas usadas:

import numpy as np
import random
import matplotlib.pyplot as plt
import os
import imageio

Funções:

#flag == 1 -> Duas massas acopladas por mola central (livres nos lados) -> Não está sendo usado neste trabalho
#flag == 2 -> Oscilador acoplado com duas massas e três molas (linear) -> Não precisaria mais existir... o flag 3 já #             resolve o caso N=2 também
#flag == 3 -> Oscilador acoplado com N>=2 massas e N+1 molas (linear)
#flag == 4 -> FPUT com termo quadrático

# L é o tamanho unidimensional de todo o sistema, de parede a parede;
# N é o número de partículas
# m é o número do modo que se quer inicializar: m=1,...,N (pode ser maior q N, mas vai repetir o primeiro modo)

def condicao_inicial_1d(L, N, m, flag):
  #Espaçamento das posições de equilíbrio das partículas
  step = L/(N+1)

  posx = np.zeros((N,1), dtype="float")
  posxeq = np.zeros((N,1), dtype="float")

  for i in range(0, N):
    posxeq[i] = (i+1)*step

  l = posxeq[1] - posxeq[0]

  #Duas massas acopladas por uma mola, ligeiramente comprimidas. Deixei por razões sentimentais, foi o primeiro que 
  #fiz.
  if(flag==1):                   
    posx[0] = posxeq[0]+0.5
    posx[1] = posxeq[1]-0.5
  #Duas massas, três molas, sistema acoplado:
  elif(flag==2):
    for i in range(0, N):
      posx[i] = modo_de_osc(posxeq[i], i, m, 0.5, N) 
  #N>=2 massas, N+1 molas, sistema acoplado, condição randômica de deslocamento do equilíbrio comentada:
  elif(flag==3):
    for i in range(N):
      #posx[i] = posxeq[i] + random.uniform(-l/5, +l/5) -> Se quiser condição inicial randômica, descomentar essa 
      #linha e comentar a linha seguinte
      posx[i] = modo_de_osc(posxeq[i], i, m, 0.5, N)
  #FPUT:
  elif(flag==4):
    for i in range(0, N):
      posx[i] = modo_de_osc(posxeq[i], i, m, 0.5, N)          

  return posx, posxeq

# Para inicializar o sistema em um determinado modo de oscilação no tempo t=0
# l é a posição da partícula na cadeia (de 0 a N-1), m é o número do modo que se quer: 1, 2, 3,...,N
# am é a amplitude da oscilação (tecnicamente, a amplitude do modo vai ser am*np.sin(...))
def modo_de_osc(poseq, l, m, am, N):
  return poseq + am*np.sin(np.pi*(l+1)*m/(N+1))

#Só para corda do FPUT, não precisa de flag
def condicao_inicial_1d_vertical(L, N, m):           
  pos = np.zeros((N,1), dtype="float")
  poseq = np.zeros((N,1), dtype="float")

  for i in range(0, N):
    pos[i] = modo_de_osc(poseq[i], i, m, 0.5, N)          
  
  return pos

def condicao_inicial_2d(N):
  N_sqrt = int(np.sqrt(N))          # N tem q ser um quadrado perfeito!

  L = N_sqrt + 1
  l = step = 1

  posx = np.zeros((N,1), dtype="float")
  posxeq = np.zeros((N,1), dtype="float")

  posy = np.zeros((N,1), dtype="float")
  posyeq = np.zeros((N,1), dtype="float")

  for i in range(N):
    posxeq[i] = (i%N_sqrt)*step + 1
    posyeq[i] = int(i/N_sqrt)
    posx[i] = posxeq[i] + random.uniform(-l/5, +l/5)
    posy[i] = posyeq[i] + random.uniform(-l/5, +l/5)
  
  #Oscilador 2D também poderia ter inicialização por modo... mas não implementei.

  return L, posx, posy, posxeq, posyeq

def acel_quad(pos, poseq, k, m, alpha, flag):
  size = len(pos)
  l = poseq[1] - poseq[0]
  a = np.zeros((size,1), dtype="float")
  
  if(flag==1):
    a[0] = -k/m*(pos[0]-pos[1]+l)   #l é a distância entre as massas quando a mola está relaxada
    a[1] = k/m*(pos[0]-pos[1]+l)
  elif(flag==2):
    a[0] = k/m*((pos[1]-poseq[1])-2*(pos[0]-poseq[0]))     
    a[1] = k/m*((pos[0]-poseq[0])-2*(pos[1]-poseq[1])) 
  elif(flag==3 or flag==4):         #Se flag==3, cuidar para fazer alpha=0!
    x = pos-poseq                                      
    for i in range(size):                                 
      #Cuidar as fronteiras
      if(i==0):
        a[i] = k/m * (x[i+1] - 2*x[i] - alpha*(x[i]**2-(x[i+1]-x[i])**2))
      elif(i==(size-1)):
        a[i] = k/m * (x[i-1] - 2*x[i] - alpha*((x[i]-x[i-1])**2-x[i]**2))
      else:
        a[i] = k/m * (x[i+1] + x[i-1] - 2*x[i]) * (1.0 + alpha * (x[i+1] - x[i-1]));

  return a

def velocity_verlet(pos, poseq, vel, acel, k, m, alpha, flag):
  size = len(pos)
  new_acel = np.zeros((size,1), dtype="float")
  new_pos = np.zeros((size,1), dtype="float")
  new_vel = np.zeros((size,1), dtype="float")
  
  new_pos = pos + vel*dt + 0.5*acel*dt**2
  new_acel = acel_quad(new_pos, poseq, k, m, alpha, flag)
  new_vel = vel + 0.5*(acel + new_acel)*dt

  return new_pos, new_vel, new_acel

def plota_figura_1d(L, posX, posY, t, flag):
  figname = 'f{:.1f}.jpeg'.format(t)

  plt.figure(figsize=(5,3))
  plt.xlim(0, L)
  plt.ylim(-1, 1)
  #plt.yticks([])
  plt.axis("off")
  #plt.grid()

  if(flag==1):
    plt.title('Duas massas acopladas por uma mola - Tempo {:.1f}'.format(t))
    plt.scatter(posX, posY, s=300, c="black")
  elif(flag==2 or flag == 3):
    plt.title('Oscilador Acoplado Linear - Tempo {:.1f}'.format(t))
    plt.scatter(posX, posY, s=300, c="black")
  elif(flag==4):
    plt.title('FPUT quadrático - Tempo {:.1f}'.format(t))
    plt.scatter(posX, posY, s=5, c="black")
  
  plt.axvline(x=0)
  plt.axvline(x=L)

  if(flag != 4):
    size = len(posX)
    aj = 0.4            #Ajuste para plotar as linhas horizontais, varia com o tamanho do scatter, aqui, s=300

    for i in range(size):
      if(i==0):
        plt.axhline(y=0, xmin=0, xmax=(posX[i]-aj)/L, color='r', linestyle=':')
        plt.axhline(y=0, xmin=(posX[i]+aj)/L, xmax=(posX[i+1]-aj)/L, color='k', linestyle=':')
      elif(i==(size-1)):
        plt.axhline(y=0, xmin=(posX[i]+aj)/L, xmax=1, color='r', linestyle=':')
      else:
        plt.axhline(y=0, xmin=(posX[i]+aj)/L, xmax=(posX[i+1]-aj)/L, color='k', linestyle=':')

  plt.savefig(figname, dpi=100)           
  plt.close()
  return figname

def plota_figura_2d(L, posX, posY, t, flag):
  figname = 'f{:.1f}.jpeg'.format(t)

  plt.figure(figsize=(5,3))
  plt.xlim(0, L)
  plt.ylim(-1, L)
  #plt.yticks([])
  plt.axis("off")
  #plt.grid()

  plt.title('Oscilador Acoplado Linear 2D - Tempo {:.1f}'.format(t))
  
  plt.scatter(posX, posY, s=300, c="black")

  plt.savefig(figname, dpi=100)           
  plt.close()

  return figname

def energia_do_modo(x, x_old, N, m, dt):
  #Energia do modo m
  omega = 2*np.sin(m*np.pi/(2*(N+1)))
  omega2 = omega**2

  tmp = np.arange(1, N+1, 1)
  tmp[-1:] = 0

  pos = np.reshape(x, (1,N))
  pos_old = np.reshape(x_old, (1,N))
    
  produto_escalar = np.sqrt(2.0 / (N+1)) * np.dot(pos, np.sin( tmp * m * np.pi / (N+1) ))
  produto_escalar_ant = np.sqrt(2.0 / (N+1)) * np.dot(pos_old, np.sin( ( tmp * m * np.pi / (N+1) ) ))
    
  return 0.5 * (produto_escalar - produto_escalar_ant)**2 / dt**2 + (0.5 * omega2 * produto_escalar**2)

Imprimir os cinco primeiros modos de oscilação para N=32 partículas (Figura 12 da Wiki).

N=32

modo = np.zeros((5,N), dtype="float")

nodes = [i for i in range(1,33)]

for i in range(0,5):
  for j in range(0, N):
    modo[i,j] = modo[i,j] + np.sin((j+1)*np.pi*(i+1)/(N+1))
  
plt.plot(nodes, modo[0,:], "k.", label="Modo 1")
plt.plot(nodes, modo[1,:], "r-.", label="Modo 2")
plt.plot(nodes, modo[2,:], "b-", label="Modo 3")
plt.plot(nodes, modo[3,:], "g--", label="Modo 4")
plt.plot(nodes, modo[4,:], "y:", label="Modo 5")

plt.xlabel("Número da partícula")
plt.ylabel("Deslocamento em relação à posição de equilíbrio")

plt.legend()
plt.show()

Exemplo de como podem ser feitos os gráficos dos modos e como comparar com as frequências esperadas (gráficos 2 a 11 da Wiki). Também estou calculando a energia total do sistema. O exemplo é do modo 1, do caso N=3. Note que a energia não é exatamente conservada, é uma dificuldade do método velocity Verlet. Para N=2 a conservação é bem obedecida, uma reta próxima a 0.2. Seria possível melhorar diminuindo o passo de integração, mas preferi deixar assim para ilustrar.

alpha = 0.0
k = 1.0
massa = 1.0
L = 12
N = 3

t=0.0
dt = 0.1
tmax = 30

posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")

posY = np.zeros((N,1), dtype="float")

velX = np.zeros((N,1), dtype="float")

acelX = np.zeros((N,1), dtype="float")
  
posX, posXeq = condicao_inicial_1d(L, N, 1, 3)

x1 = []
x2 = []
x3 = []
dotx1 = []
dotx2 = []
dotx3 = []
energia = []

x1.append(posX[0]-posXeq[0])
x2.append(posX[1]-posXeq[1])
x3.append(posX[2]-posXeq[2])
dotx1.append(velX[0])
dotx2.append(velX[1])
dotx3.append(velX[2])

energia.append(0.5*((posX[0]-posXeq[0])**2+(posX[1]-posXeq[1])**2+(posX[2]-posXeq[2])**2))

while(t < tmax):
  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 3)

  x1.append(posX[0]-posXeq[0])
  x2.append(posX[1]-posXeq[1])
  x3.append(posX[2]-posXeq[2])
  dotx1.append(velX[0])
  dotx2.append(velX[1])
  dotx3.append(velX[2])

  energia.append(0.5*((posX[0]-posXeq[0])**2+(posX[1]-posXeq[1])**2+(posX[2]-posXeq[2])**2) + 0.5*(velX[0]**2+velX[1]**2+velX[2]**2))
  
  t += dt
time = np.linspace(0, tmax, len(energia))
plt.plot(time, energia)
plt.ylim(0, 1)
plt.show()
time = np.linspace(0, 30, len(x1))
plt.plot(time, 0.35*np.cos(np.sqrt(2-np.sqrt(2))*time), "g-", label=r"$cos(\sqrt{2-\sqrt{2}}t)$")
plt.plot(time, x1, "k.", label=r"$x_{1}$")
#plt.plot(time, x2, "r-", label=r"$x_{2}$")
#plt.plot(time, x3, "b.", label=r"$x_{2}$")


plt.xlim(0,30)

plt.xlabel("t")
plt.ylabel("Deslocamento em relação à posição de equilíbrio")

plt.title("Oscilador Linear (Modo 1) - N=3")
plt.legend(loc="lower right")
plt.show()

FPUT unidimensional, deslocamento horizontal:

Setar alpha diferente de zero e flag=4 (último parâmetro das funções condicao_inicial(), plota_figura_1d() e velocity_verlet()).

Salva animação no arquivo fput.gif. Bem grande o gif gerado.

As animações não têm linhas representando as molas e o tamanho das partículas é bem pequeno. Fiz isso de propósito pra visualizar melhor o caso com N=32 partículas, do artigo da RBF.

Lembrar de alterar as variáveis globais de acordo. Diferentemente do caso bidimensional, o valor de L tem q ser fornecido e deve ser obviamente diferente de zero.

alpha = 0.25
k = 1.0
massa = 1.0
L = 33
N = 32
images=[]

t=0.0
dt = 0.1
tmax = 30

posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")

#Por enquanto não estou usando posições em Y
posY = np.zeros((N,1), dtype="float")

velX = np.zeros((N,1), dtype="float")

acelX = np.zeros((N,1), dtype="float")

posX, posXeq = condicao_inicial_1d(L, N, 1, 4)

while(t < tmax):
  figname = plota_figura_1d(L, posX, posY, t, 4)
  images.append(imageio.imread(figname))
  os.remove(figname)
  
  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 4)

  t += dt

imageio.mimsave('fput.gif', images, fps=30)

FPUT unidimensional, deslocamento vertical (como no programa em C):

Setar alpha diferente de zero e flag=4 (último parâmetro das funções condicao_inicial(), plota_figura_1d() e velocity_verlet()).

alpha = 0.25
k = 1.0
massa = 1.0
L = 33
N = 32
images=[]

t=0.0
dt = 0.1
tmax = 100

posX = np.zeros((N,1), dtype="float")

posY = np.zeros((N,1), dtype="float")
posYeq = np.zeros((N,1), dtype="float")

velY = np.zeros((N,1), dtype="float")

acelY = np.zeros((N,1), dtype="float")

step = L/(N+1)

for i in range(0, N):
  posX[i] = (i+1)*step

posY = condicao_inicial_1d_vertical(L, N, 1)

while(t < tmax):
  figname = plota_figura_1d(L, posX, posY, t, 4)
  images.append(imageio.imread(figname))
  os.remove(figname)
  
  posY, velY, acelY = velocity_verlet(posY, posYeq, velY, acelY, k, massa, alpha, 4)
  t += dt

imageio.mimsave('fput.gif', images, fps=30)

Cálculo e gráficos das energias por modo do FPUT quadrático:

alpha = 0.25
k = 1.0
massa = 1.0
L = 33
N = 32
images=[]

t=0.0
dt = 0.1
tmax = 14000

energia1 = []
energia2 = []
energia3 = []
energia4 = []
energia5 = []

posX = np.zeros((N,1), dtype="float")

posY = np.zeros((N,1), dtype="float")
posYeq = np.zeros((N,1), dtype="float")

velY = np.zeros((N,1), dtype="float")

acelY = np.zeros((N,1), dtype="float")

step = L/(N+1)

for i in range(0, N):
  posX[i] = (i+1)*step

posY = condicao_inicial_1d_vertical(L, N, 1)

while(t < tmax):
  pos_old = posY
  
  posY, velY, acelY = velocity_verlet(posY, posYeq, velY, acelY, k, massa, alpha, 4)

  energia1.append(energia_do_modo(posY, pos_old, N, 1, dt))  
  energia2.append(energia_do_modo(posY, pos_old, N, 2, dt))
  energia3.append(energia_do_modo(posY, pos_old, N, 3, dt))
  energia4.append(energia_do_modo(posY, pos_old, N, 4, dt))
  energia5.append(energia_do_modo(posY, pos_old, N, 5, dt))

  t += dt
time = np.linspace(0, tmax, int(tmax/dt))

plt.figure(figsize=(15,7))
plt.title('Energias por modo de vibração')
plt.xlabel("t")
plt.ylabel("Energia (u.a.)")
plt.plot(time, energia1,label="Modo 1")
plt.plot(time, energia2,label="Modo 2")
plt.plot(time, energia3,label="Modo 3")
plt.plot(time, energia4,label="Modo 4")
plt.plot(time, energia5,label="Modo 5")
plt.legend()

Osciladores acoplados lineares unidimensionais:

Salva animação no arquivo OA1D-N#.gif, com # igual ao número de partículas da simulação. Bem grande o gif gerado. Fica ainda maior quanto mais partículas tiver a simulação. Lembrar de alterar as variáveis globais de acordo. Diferentemente do caso bidimensional, o valor de L tem q ser fornecido e deve ser obviamente diferente de zero.

alpha = 0.0
k = 1.0
massa = 1.0
L = 12
N = 4
images=[]

t=0.0
dt = 0.1
tmax = 50

posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")

posY = np.zeros((N,1), dtype="float")

velX = np.zeros((N,1), dtype="float")

acelX = np.zeros((N,1), dtype="float")
  
posX, posXeq = condicao_inicial_1d(L, N, 1, 3)

while(t < tmax):
  figname = plota_figura_1d(L, posX, posY, t, 3)
  images.append(imageio.imread(figname))
  os.remove(figname)
  
  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 3)
  
  t += dt

imageio.mimsave('OA1D-N{:.0f}.gif'.format(N), images, fps=60)

Osciladores acoplados em uma rede bidimensional L x L (Aqui o L é calculado de maneira a ter um grid quadrado com espaçamento 1 entre as partículas, só precisa ser fornecida a quantidade de partículas N). N TEM QUE SER UM QUADRADO PERFEITO. LEMBRAR DE ALTERAR AS VARIÁVEIS GLOBAIS, FAZER L=0.0 E N = 4, 9, 16, 25,..., ETC.

alpha = 0.0
k = 1.0
massa = 6.0
L = 0.0
N = 49
images=[]

t=0.0
dt = 0.1
tmax = 30

posX = np.zeros((N,1), dtype="float")
posXeq = np.zeros((N,1), dtype="float")

posY = np.zeros((N,1), dtype="float")
posYeq = np.zeros((N,1), dtype="float")

velX = np.zeros((N,1), dtype="float")
velY = np.zeros((N,1), dtype="float")

acelX = np.zeros((N,1), dtype="float")
acelY = np.zeros((N,1), dtype="float")
  
L, posX, posY, posXeq, posYeq = condicao_inicial_2d(N)

while(t < tmax):
  figname = plota_figura_2d(L, posX, posY, t, 3)
  images.append(imageio.imread(figname))
  os.remove(figname)

  posX, velX, acelX = velocity_verlet(posX, posXeq, velX, acelX, k, massa, alpha, 3)
  posY, velY, acelY = velocity_verlet(posY, posYeq, velY, acelY, k, massa, alpha, 3)
  t += dt

imageio.mimsave('OA2D-N{:.0f}.gif'.format(N), images, fps=60)