Simulace matematického kyvadla

import numpy as np
import matplotlib.pyplot as plt

teoretická hustota pravděpodobnosti je $f(x)=\frac{1}{\pi\sqrt{x_0^2-x^2}}$

#teoretická hustota pravděpodobnosti výskytu kyvadla
def f(x,x0):
    if abs(x)<x0:
        return 1.0/(np.pi*np.sqrt(x0**2-x**2))
    else:
        return 0.0

#vektorizovaná funkce f(x) aby byla schopna pracovat najednou s celým polem x
vf=np.vectorize(f)        

Protože ještě neznáme metodu inverzní funkce pro generování náhodné proměnné se zadaného rozdělení, budeme simulovat časy kdy mačkáme spoušť fotoaparátu.

ndata=100000 #počet simulovaných dat
nbins=100 #počet bunů histogramu
l=1 #délka závěsu 1 m
g=9.81 #tíhové zrychlení
alpha_0=5*np.pi/180 #=uhel max. výchylky 5 stupňů (převedeno na radiány) 
x0=alpha_0*l #max. výchylka
omega=np.sqrt(l/g)
T=2*np.pi/omega #perioda kmitu
t=T*np.random.random(ndata)

časová závislost $x$-ové souřadnice kyvadla je $x(t) = x_0\cos\omega t$

x=x0*np.cos(omega*t)

Nakreslíme normovaný histogram simulovaných dat a pro srovnání také teoretickou hustotu pravděpodobnosti

plt.hist(x,bins=nbins,density=True,label="simulace")
#do grafu nakreslíme pro srovnání teoretickou hustotu pravděpodobnosti f(x)
x_teor=np.linspace(-x0,x0,1000)
plt.plot(x_teor,vf(x_teor,x0),c='red',label="teoretická f(x)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.show()