Sistemi NDE višjega reda

Diferencialne enačbe višjega reda lahko rešujemo z istimi metodami kot sisteme enačb 1. reda.

Prevod enačbe višjega reda na sistem enačb 1. reda

Navadno diferencialno enačbe (NDE) reda $n$

\[\begin{equation}\label{eq:nde} y^{(n)}(t) = f(y^{(n-1)}, \ldots y'(t), y(t), t) \end{equation}\]

lahko s preprosto vpeljavo novih spremenlivk prevedemo v sistem enačb 1. reda. Za vsak odvod $y^{(k)}(t)$ reda $k < n$ vpeljemo novo spremenljivko

\[\begin{equation}\label{eq:spremenljivke} \begin{aligned} z_1(t) &= y(t)\cr z_2(t) &= y'(t)\cr \ldots \cr z_{n}(t) & = y^{(n-1)}(t) \end{aligned} \end{equation}\]

Če upoštevamo, da je $y^{(n)}(t) = z_n'(t)$, lahko enačbo \eqref{eq:nde} lahko v novih koordinatah zapišemo

\[ z_n'(t) = f(z_n(t), \ldots z_2(t),z_1(t), t)\]

enačbe za $z_k$ dobimo tako, da v odvajamo enačbe \eqref{eq:spremenljivke} in dobimo sistem $n$ enačb prvega reda

\[\begin{aligned}\label{eq:sistem} z_1'(t) &= z_2(t)\cr z_2'(t) &= z_3(t)\cr \ldots \cr z_{n}'(t) & = f(z_n(t), \ldots z_2(t),z_1(t), t) \end{aligned}\]

Poševni met

Teniško žogico vržemo v zrak pod odločenim kotom. Kako se žogica giblje? Enačbe gibanja lahko zapišemo z 2. Newtonovim zakonom

\[\begin{equation}\label{eq:newton2} \ddot{\vec{x}}(t) = \vec{a} = \frac{1}{m}\vec{F}, \end{equation}\]

kjer je

vektor $\vec{x}(t)=(x(t),y(t),z(t))$ vektor, ki označuje položaj podan s koordinatami,
vektor $\vec{F}$ vsota vseh sil, ki delujejo na žogico,
parameter $m$ pa masa žogice.

Sistem 3 enačb 2. reda \eqref{eq:newton2} lahko spremnimo v sistem 6 enač 1. reda, če vpeljemo nove spremnljivke za odvode koordinat. V tem primeru so nove spremenljivke kar komponente vektorja hitrosti

\[\vec{v(t)} = \dot{\vec{x}}(t) = (v_x(t), v_y(t), v_z(t))\]

Sistem \eqref{eq:newton2} lahko prevedemo v sistem

\[\begin{equation}\label{eq:enacbe_gibanja} \dot{\vec{z}}(t) = \begin{bmatrix}\dot{\vec{x}}\cr \dot{\vec{v}}\end{bmatrix} = \begin{bmatrix}\vec{v}\cr \frac{1}{m}\vec{F}\end{bmatrix} \end{equation}\]

Na kroglico delujeta dve sili:

gravitacijska sila $\vec{F}_g = -mg\vec{e}_z$, kjer $g$ gravitacijski pospešek na površju Zemlje in $\vec{e}_z$ enotski vektor v vertikalni smeri.
sila upora $\vec{F}_u$, ki jo opiše kvadratni zakon upora.

Če določimo začetni položaj $\vec{x}(t_0) = \vec{x}_0$ in začetno hitrost $\vec{v}(t_0)=\vec{v}_0$ v določenem trenutku $t_0$, sistem NDE \eqref{eq:enacbegibanja} natanko določa gibanje žogice. Trajetorjija žogice je rešitev [začetnega problema](https://en.wikipedia.org/wiki/Initialvalueproblem) za enačbe \eqref{eq:enacbegibanja} in začetni pogoj

\[\vec{z}(t_0) = \begin{bmatrix}\vec{x}_0\cr \vec{v}_0\end{bmatrix}.\]

Kvadratni zakon upora

Sila upora, je sila s katero sredstvo(zrak, voda) zavira pot telesa, ko se le to premika po sredstvu. Način, kako je upor odvisen od hitrosti in dimenzij telesa in lastnosti sredstva je zelo odvisen od Raynoldsovega števila. Za naš primer, ko žogica potuje skozi zrak velja kvadratni zakon upora, po katerem je velikost sile upora sorazmerna kvadratu velikosti hitrosti in kaže v nasprotno smer hitrosti

\[\vec{F}_u(\vec{v}) = -\frac{1}{2}c_v\rho S\|\vec{v}\|\vec{v}.\]

Parametri, ki nastopajo v kvadratnem zakonu so

brezdimenzijski koeficient upora $c_v$, ki je odvisen od oblike (za kroglo je $c_v=0.47$)
gostota sredstva $qrho$ (za zrak je $\rho=1.2\mathrm{kgm}^{-3}$)
ploščina $S$ prečnega prereza telesa v smeri, ki je pravokotna na gibanje telesa

Naloga

Napiši funkcije

upor(v; c=0.47, ρ=1.2, S=0.01), ki izračuna zračni upor v odvisnoti od vektorja hitrosti v
gravitacija(m; g=9.81), ki izračuna vektor sile teže na telo z maso m in
posevni_met((x,v), t, parametri), ki vrne desne strani sistema enačb \eqref{eq:enacbe_gibanja}.

Funkcijo resi_nde predelaj tako, da računa zaporedje približkov, dokler žogica ne pade na tla.

Napiši funkcijo nicla_nde(f, y), ki poišče rešitev enačbe f(y(t))=0, pri čemer je y(t) rešitev sistema NDE. Uporabiš lahko Newtonovo metodo ali kubično interpolacijo.

Vprašanja

Kako daleč leti teniška žogica, če jo bržemo z začetne višine 2m s hitrostjo $30\mathrm{kmh}^{-1}$ pod kotom 45°? Kolikšna je dolžina trajektorije, ki jo je žogica prepotovala?

Pod kolikšnim kotom naj vržemo žogico, da bo letela najdlje. Ali je kot odvisen od začetne hitrosti?

Perioda geosinhrone orbite

Zaradi vrtenja okrog svoje osi in razgibanega reliefa, Zemlje ne moremo obravnavati kot popolno kroglo. Sploščenost na polih in relief se pozna tudi na Zemljinem gravitacijskem polju. Gravitacijsko polje lahko opišemo kot gradient potenciala.

Gravitacijski potencial Zemlje lahko razvijemo z multipolnim razvojem

\[u=-{\frac {\mu }{r}}+\sum _{n=2}^{N_{z}}{\frac {J_{n}P_{n}^{0}(\sin \theta )}{r^{n+1}}}+\sum _{n=2}^{N_{t}}\sum _{m=1}^{n}{\frac {P_{n}^{m}(\sin \theta )(C_{n}^{m}\cos m\varphi +S_{n}^{m}\sin m\varphi )}{r^{n+1}}}\]

kjer je $\mu = GM$ produkt mase zemlje in gravitacijske konstante. Glavnemu delu potenciala $-\frac{\mu}{r}$ ustreza sila

\[\vec{F} = \frac{\mu}{r^2}\vec{x}.\]

Upoštevali bomo le en dodatni člen (imenovan tudi kvadrupolni moment)

\[{\frac {J_{2}\ P_{2}^{0}(\sin \theta )}{r^{3}}}=J_{2}{\frac {1}{r^{3}}}{\frac {1}{2}}(3\sin ^{2}\theta -1)=J_{2}{\frac {1}{r^{5}}}{\frac {1}{2}}(3z^{2}-r^{2})\]

višje člene v razvoju pa zanemarimo. Sila, ki je posledica kvadrupolnega momenta, je enaka

\[\begin{aligned} &F_{x}=-{\frac {\partial u}{\partial x}}=J_{2}\ {\frac {x}{r^{7}}}\left(6z^{2}-{\frac {3}{2}}(x^{2}+y^{2})\right)\\ &F_{y}=-{\frac {\partial u}{\partial y}}=J_{2}\ {\frac {y}{r^{7}}}\left(6z^{2}-{\frac {3}{2}}(x^{2}+y^{2})\right)\\ &F_{z}=-{\frac {\partial u}{\partial z}}=J_{2}\ {\frac {z}{r^{7}}}\left(3z^{2}-{\frac {9}{2}}(x^{2}+y^{2})\right)\end{aligned}\]

Natančne vrednosti koeficientov $\mu$ in $J_2$ lahko razberemo iz geopotencialnega modela:

\[\mu = 398600,440 km^3⋅s^{−2}, J_2 = 1,75553 × 10^{10} km^5⋅s^{−2}\]