# Partikkel skytes i tyngdefeltet

En partikkel skytes ut i tyngdefeltet med fart $v_0=|\vec{v}(0)|$ i en retning som danner en vinkel $\alpha$ med horisontalen (z-aksen). Vi antar at det kun er tyngdekraften som virker, så Newton's annen lov gir oss at 

$$
\begin{equation}
m \vec{a} = -mg \mathbf{k},
\end{equation}
$$

der $m$ er massen til partikkelen. Vi kan se bort ifra massen siden vi kan dele begge sider på $m$.



Vi kan integrere Newton's annen low to ganger og finne posisjonen $\vec{r}(t)=x \mathbf{i} + z \mathbf{k}$ for partikkelen. Bruker først at $\vec{a}=d\vec{v}/dt$ slik at vi får

$$
\begin{equation}
\frac{d \vec{v}}{dt} = -g\mathbf{k}.
\end{equation}
$$

Integrerer opp og får

$$
\begin{equation}
\vec{v}(t) - \vec{v}(0) = -g t \mathbf{k}.
\end{equation}
$$

Bruker så at $\vec{v} = d \vec{r} / dt$ og integrerer på samme måte fra

$$
\begin{equation}
\frac{d \vec{r}}{dt} = \vec{v}(0) - g t \mathbf{k}.
\end{equation}
$$

Ender opp med at

$$
\begin{equation}
\vec{r}(t) - \vec{r}(0) = \vec{v}(0) t - \frac{1}{2}g t^2 \mathbf{k},
\end{equation}
$$

der vi kan ta bort $\vec{r}(0)$ siden vi regner at utskytningen starter i origo. 

Vektoren $\vec{v}(0)$ kan dekomponeres langs $x$ og $y$-akser ved å bruke vinkelen $\alpha$. Får at 

$$
\begin{equation}
\vec{v}(0) = v_0 \cos(\alpha) \mathbf{i} + v_0 \sin(\alpha) \mathbf{k},
\end{equation}
$$

og dermed har vi funnet posisjonsvektoren $\vec{r}$ med komponenter

$$
\begin{align}
  x &= v_0 \cos(\alpha) t,\\
  z &= v_0 \sin(\alpha) t - \frac{1}{2}gt^2.
\end{align}
$$

slik at posisjonsvektoren blir

$$
\begin{align}
\vec{r}(t) &= x(t) \mathbf{i} + z(t) \mathbf{k}, \notag\\
           &= v_0 \cos(\alpha) t \mathbf{i} + (v_0 \sin(\alpha) t - \frac{1}{2}gt^2) \mathbf{k}.
\end{align}
$$

Hastighetsvektoren kan også samles sammen til

$$
\begin{equation}
\label{eq:v}
\vec{v}(t) = v_0 \cos(\alpha) \mathbf{i} + (v_0 \sin(\alpha) - gt) \mathbf{k}.
\end{equation}
$$

Vi skal nå modellere dette med et interaktivt plot som viser partikkelbanen og vektoren. Til dette trenger vi å gjøre et par observasjoner. 

 1. Hva er domenet til partikkelen?
   1. Hvor høyt skytes partikkelen?
   2. Hvor lang tid tar det før den kommer ned igjen?

Finner svar på A ved å se på når hastigheten i z-retningen er lik null (fra $z$-komponenten til $\vec{v}$ over)

$$
\begin{align*}
v_0 \sin(\alpha) - g t_0 &= 0 \\
t_0 &= \frac{v_0 \sin(\alpha)}{g}.
\end{align*}
$$

Høyden ved $t_0$ er gitt ved $z(t_0)$, som vi kaller $H$

$$
z(t_0) = H = \frac{(v_0 \sin(\alpha))^2}{2 g}
$$

Så domenet for $z$ må være minst $[0, H]$ for å få med hele partikkelbanen. 

Tiden det tar for partikkelen å komme ned igjen finner vi ved å løse $z(t) = 0$

$$
\begin{align*}
v_0 \sin(\alpha) t - \frac{1}{2}gt^2 &= 0 \\
t\left(v_0 \sin(\alpha) - \frac{1}{2}g t \right) &= 0
\end{align*}
$$

som løses enkelt til

$$
t = 0 \lor t = \frac{2 v_0 \sin(\alpha)}{g}
$$

Vi kaller tiden da partikkelen igjen treffer bakken ($z=0$) for $T=2 v_0 \sin(\alpha)/g$. Vi innser at simuleringen må gå minst i $t \in [0, T]$.

Vi finner domenenet for $x$-aksen ved å sette sluttiden inn i uttrykket for $x$

$$
x\left(T\right) = v_0 \cos(\alpha) T.
$$

Så $x$-domenet blir $[0, v_0 \cos(\alpha) T]$.

Da er vi klare til å starte plottingen!

Vi starter med å importere verktøyene vi trenger. Vi vil lage et interaktivt plot som viser partikkelen langs en bane. Da kan vi bruke [ipywidgets](https://ipywidgets.readthedocs.io/en/latest/).


```python
from ipywidgets import interact
import matplotlib.pyplot as plt
import numpy as np
```


```python
%matplotlib inline

g = 9.81
v0 = 10
a = np.pi/4
T = 2*v0*np.sin(a)/g
N = 100 # Bruker interval med N punkter
t0 = np.linspace(0, T*1.1, N) 
```


```python
# Beregner hele partikkelbanen først
x0 = v0*np.cos(a)*t0
z0 = v0*np.sin(a)*t0 - 0.5*g*t0**2
    
def partikkelbane(t):    
    # plot hele partikkelbanen
    plt.figure()
    plt.plot(x0, z0, 'b')
    
    # plot partikkelen
    x = v0*np.cos(a)*t
    z = v0*np.sin(a)*t - 0.5*g*t**2
    plt.plot(x, z, 'ok')
    
    # plot en hastighetsvektor ved t0 og t
    plt.arrow(0, 0, v0*np.cos(a)*T/4, v0*np.sin(a)*T/4, head_width=v0*T/100)
    plt.arrow(x, z, v0*np.cos(a)*T/4, (v0*np.sin(a)-g*t)*T/4, head_width=v0*T/100)
    
    # posisjonsvektor
    plt.arrow(0, 0, x, z, length_includes_head=True)
    
    plt.ylim(z0.min(), 2*z0.max())
    
    plt.text(0.25*v0*np.cos(a), 0.3*v0*np.sin(a), r'$\vec{v}(0)$')
    plt.text(x+0.25*v0*np.cos(a), z+0.2*(v0*np.sin(a)-g*t), r'$\vec{v}(t)$')
    plt.text(0.5*x, 0.4*z, r'$\vec{r}(t)$')
    
    plt.show()
    
interact(partikkelbane, t=(0, T, T/20))
```


    interactive(children=(FloatSlider(value=0.7208020195581524, description='t', max=1.4416040391163047, step=0.07…





    <function __main__.partikkelbane(t)>



Vi kan nå beregne buelengden, eller strekningen partikkelen tilbakelegger, ved å integrere numerisk

$$
\begin{equation}
L = \int_{0}^{t_m}|\vec{v}(t)|dt
\end{equation}
$$

Numerisk integrasjon kan gjøres i `numpy` ved å bruke trapes-metoden [trapz](https://docs.scipy.org/doc/numpy/reference/generated/numpy.trapz.html). Den krever som input funksjonsverdiene $|\vec{v}(t_i)|$ og punktene der funksjonsverdiene er samplet, her $t_i=i T/(N-1)$ for $i=0, 1, \ldots, N-1$ og vi har brukt $N=100$ punkter i koden over.  


```python
t0 = np.linspace(0, T, N) 
x0 = v0*np.cos(a)*t0
z0 = v0*np.sin(a)*t0 - 0.5*g*t0**2
L = np.trapz(np.sqrt((v0*np.cos(a))**2 + (v0*np.sin(a) - g*t0)**2), t0)     
print('Buelengden er %2.3f'%(L))
```

    Buelengden er 11.700


# Sympy

Vi skal nå gjøre samme implemetering ved å benytte [sympy](https://sympy.org) for eksakt symbolsk manipulering av vektorer. Til det skal vi bruke `sympy` sin `vector`-modul som inneholder en klasse for et koordinatsystem `CoordSys3D`.


```python
import sympy as sp
from sympy.vector import CoordSys3D
```

Vi antar nå først at både `v0,a` og `t` er variable


```python
g = 9.81
v0, a, t = sp.symbols('v0,a,t')
N = CoordSys3D('N')

x = v0*sp.cos(a)*t 
z = v0*sp.sin(a)*t - 0.5*g*t**2

# Definer posisjonsvektoren
r = x*N.i + z*N.k

# Deriver mhp t for å finne hastigheten
v = sp.diff(r, t)

print(r)
print(v)
```

    (t*v0*cos(a))*N.i + (-4.905*t**2 + t*v0*sin(a))*N.k
    (v0*cos(a))*N.i + (-9.81*t + v0*sin(a))*N.k



```python
# Legg parametre i dictionary
v0 = 10
a = np.pi/4
d = {'v0': v0, 'a': a}
T = 2*v0*np.sin(a)/g

# Redefiner variable med disse parameterne
x = x.subs(d)
z = z.subs(d)
r = r.subs(d)
v = v.subs(d)

# Gjør x og z kallbare funksjoner for numpy-arrays
x = sp.lambdify(t, x)
z = sp.lambdify(t, z)

# Evaluer x og z for hele tidsdomenet
t0 = np.linspace(0, T*1.1, 100)

x0 = x(t0)
z0 = z(t0)
vx0 = np.float(v.subs(t, 0).dot(N.i))
vz0 = np.float(v.subs(t, 0).dot(N.k))
```


```python
%matplotlib inline
    
def partikkelbane(t):    
    # plot hele partikkelbanen
    plt.figure()
    plt.plot(x0, z0, 'b')
    
    # plot partikkelen
    plt.plot(x(t), z(t), 'ok')
    
    # plot en hastighetsvektor ved t0 og t
    vx = np.float(v.subs('t', t).dot(N.i))
    vz = np.float(v.subs('t', t).dot(N.k))
    scale = T/4
    plt.arrow(0, 0, vx0*scale, vz0*scale, head_width=v0*T/100)
    plt.arrow(x(t), z(t), vx*scale, vz*scale, head_width=v0*T/100)
    
    # plot posisjonsvektor
    plt.arrow(0, 0, x(t), z(t), head_width=v0*T/100, length_includes_head=True)
    
    plt.ylim(z0.min(), 2*z0.max())
    
    plt.text(0.25*vx0, 0.3*vz0, r'$\vec{v}(0)$')
    plt.text(x(t)+0.25*vx0, z(t)+0.2*vz, r'$\vec{v}(t)$')
    plt.text(0.5*x(t), 0.4*z(t), r'$\vec{r}(t)$')
    
    plt.show()
    
interact(partikkelbane, t=(0, T*1.1, T*1.1/20))
```


    interactive(children=(FloatSlider(value=0.7928822215139677, description='t', max=1.5857644430279354, step=0.07…





    <function __main__.partikkelbane(t)>




```python
Ls = sp.Integral(sp.sqrt(v.dot(v)), (t, 0, T)).evalf()
print('Buelengde = %2.4f'%(Ls))
```

    Buelengde = 11.7002


0c12c2e11368ac4be0def9bc0774aeb894d21af3
ipynb
Jupyter Notebook
notebooks/partikkelbane.ipynb
mikaem/MEK1100-21
b2dfb4dc3598f57989dbf1f397179ced9a8e39b9
BSD-2-Clause
notebooks/partikkelbane.ipynb
mikaem/MEK1100-21
b2dfb4dc3598f57989dbf1f397179ced9a8e39b9
BSD-2-Clause
notebooks/partikkelbane.ipynb
mikaem/MEK1100-21
b2dfb4dc3598f57989dbf1f397179ced9a8e39b9
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES