Vizsgáljuk meg, hogy a $v_0=30~\mathrm{m}/\mathrm{s}$-os kezdősebességű, $\alpha=50^\circ$-os szögű ferde hajítás differenciálegyenlet-rendszerének megoldása ugyanazt a pályát adja-e, mint az analitikus megoldás! Ábrázoljuk a numerikus és az egzakt megoldást ugyanazon koordinátarendszerben!
Segítségül:
\begin{align} \ddot{x}&=0&\dot{x}_0&=v_0\cdot \cos\alpha &x_0&=0\\ \ddot{y}&=-g& \dot{y}_0&=v_0\cdot \sin\alpha & y_0&=0 \end{align}Először töltsük be a szükséges csomagokat
%pylab inline
from scipy.integrate import *
Ezután definiáljuk a növekményt megadó függvényt, majd oldjuk meg a differenciálegyenletet. Mivel a probléma két darab másodrendű differenciálegyenletet tartalmaz, négy darab elsőrendű egyenletre lehet visszavezetni, amit egy négydimenziós vektorra vonatkozó differenciálegyenletként írhatunk:
$$ \frac{\textrm{d}}{\textrm{d}t} \left(\begin{array}{c} x\\ v_x\\ y\\ v_y \end{array}\right) = \left(\begin{array}{c} v_x\\ 0\\ v_y\\ -g \end{array}\right) $$Erre már definiálhatunk függvényt, és megoldhatjuk az egyenletet:
g = 10 # Megadjuk a kezdeti feltételeket
v_0 = 30
alpha = 50/180 * pi
r_0 = [0, v_0 * cos(alpha), 0, v_0 * sin(alpha)]
def f(r, t): # Majd a vektor novekményét visszaadó függvényt is definiáljuk
return [r[1], 0, r[3], -g]
t = linspace(0, 4.6, 100)
u = odeint(f, r_0, t) # végül megoldjuk az egyenletet
Ezzel, az u
listában tárolt vektorok 0. és 2. komponenseik a a test $x$ és $y$ koordinátái vannak eltárolva a mozgás során. Hozzunk létre hasonló x
és y
tömböket, de az ezekben tárolt értékeket az egzakt képletekkel számoljuk ki:
x = v_0 * t * cos(alpha)
y = v_0 * t * sin(alpha) - g/2 * t**2
Végül ábrázoljuk mindkét eredményt:
plot(x, y, "b-", lw=2) # Az exakt eredményt folytonos vonallal
plot(u[:, 0], u[:, 2], "rx", markersize=5) # A numerikus eredmént piros x-ekkel
xlabel("x-tengely", size=12)
ylabel("y-tengely", size=12)
title("50°-os ferdehajítás pályája", size=18, y=1.05)
Mint látszik, az eredmény nagyon közel azonos.