Foucaultovo kyvadlo

Výpočet trajektorie Foucaultova kyvadla v neinerciální soustavě otáčející se spolu se Zemí

import numpy as np
import matplotlib.pyplot as plt

T=1 #perioda otoceni Zeme
omega=2*np.pi/T #uhlova frekvence otaceni Zeme
dt=T/100000 #časový krok pro numerické řešení pohybové rovnice
m=1 #hmotnost zavazi
Tk=T/4#perioda kmitu 
k=(2*np.pi)**2*m/Tk**2 #tuhost pruziny

pole časů, souřadnic a rychlostí

t=np.arange(0,T,dt)
x=np.zeros(np.size(t))
y=np.zeros(np.size(t))
vx=np.zeros(np.size(t))
vy=np.zeros(np.size(t))

počáteční podmínky

x[0]=0
y[0]=0
vx[0]=0
vy[0]=1

numerické řešení pohybové rovnice

for i in range(np.size(t)-1):
    x[i+1]=x[i]+vx[i]*dt
    y[i+1]=y[i]+vy[i]*dt
    vx[i+1]=vx[i]+(omega**2*x[i]+2*omega*vy[i]-k/m*x[i])*dt
    vy[i+1]=vy[i]+(omega**2*y[i]-2*omega*vx[i]-k/m*y[i])*dt

vykreslení výsledku

fig,ax=plt.subplots(figsize=(5,5))
ax.plot(x,y)    
ax.set_xlabel('x',fontsize=12)
ax.set_ylabel('y',fontsize=12)
plt.show()