In [2]: 
from pylab import *
%matplotlib inline
from scipy.integrate import odeint

# nastavitve za izris grafov (http://matplotlib.org/1.3.1/users/customizing.html)
rc('text', usetex=True)
rc('font', size=12, family='serif', serif=['Computer Modern'])
rc('xtick', labelsize='small')
rc('ytick', labelsize='small')
rc('legend', fontsize='medium')
rc('figure', figsize=(5, 3))
rc('lines', linewidth=2.0)
rc('axes', color_cycle=['k'])
rc('contour', negative_linestyle='solid')

Navadne diferencialne enačbe

Pri reševanju diferencialne enačbe gre za določitev obnašanja neke količine, odvisne od ene ali več spremenljivk, pri čemer v povezavi med iskano količino \(y\) in neodvisnimi spremenljivkami \(x_i\) nastopajo tudi odvodi iskane količine po neodvisnih spremenljivkah, npr. \(\frac{\partial y}{\partial x_i}\), \(\frac{\partial^2 y}{\partial x_i\partial x_j},\ldots\) Če je neodvisna spremenljivka ena sama, govorimo o navadnih diferencialnih enačbah. Če v diferencialni enačbi nastopajo samo prvi odvodi, govorimo o navadni diferencialni enačbi prveg reda. Če poleg njih nastopijo tudi drugi odvodi, je to navadna diferencialna enačba drugega reda.

Navadno diferencialno enačbo prvega reda lahko numerično rešimo tako, da odvod nadomestimo z diferenčnim približkom:

\[\frac{\mathrm{d}y(x)}{\mathrm{d}x}=f\left(y(x),x\right)\]

\[\frac{y(x+\Delta x)-y(x)}{\Delta x}=f\left(y(x),x\right)\]

\[y(x+\Delta x)=y(x)+f\left(y(x),x\right)\Delta x\]

Na ta način lahko s poznavanjem začetnega pogoja \(y(0)\) z zaporedno uporabo zgornje zveze poiščemo vrednosti količine \(y\) pri poznejših časih:

\[y(0)\rightarrow y(\Delta x)\rightarrow y(2\Delta x)\ldots\]

Navadno diferencialno enačbo višjega reda lahko z vpeljavo novih spremenljivk spremenimo v sistem navadnih diferencialnih enačb prvega reda:

\[\frac{\mathrm{d}^2y(x)}{\mathrm{d}x^2}=f\left(y(x),\frac{\mathrm{d}y(x)}{\mathrm{d}x}, x\right),\]

\[\frac{\mathrm{d}y(x)}{\mathrm{d}x}=u(x),\]

\[\frac{\mathrm{d}u(x)}{\mathrm{d}x}=f\left(y(x),u(x),x\right),\]

ki ga nato numerično rešimo tako, kot smo se naučili pri navadni diferencialni enačbi prvega reda:

\[y(x+\Delta x)=y(x)+u(x)\Delta x\]

\[u(x+\Delta x)=u(x)+f(y(x), u(x), x)\Delta x\]

Tu potrebujemo dva začetna pogoja, vrednost iskane funkcije \(y(0)\) in njenega odvoda \(u(0)=\left.\frac{\mathrm{d}y}{\mathrm{d}x}\right|_{x=0}\):

\[\left(y(0), u(0)\right)\rightarrow\left(y(\Delta x), u(\Delta x)\right)\rightarrow\left(y(2\Delta x), u(2\Delta x)\right)\ldots\]

Algoritmi, ki jih za reševanje navadnih diferencialnih enačb ponujajo profesionalna računalniška orodja, so veliko bolj sofisticirani od zgoraj opisanjega. Z uporabo diferenčnih približkov višjega reda, spremenljivega koraka \(\Delta x\), ..., dosežejo večjo stabilnost in hitrost iskanja rešitve.

Kot primer si oglejmo izračun poti \(s\) in hitrosti \(\frac{\mathrm{d}s}{\mathrm{d}t}\) telesa pri padanju le-tega v zemeljskem gravitacijskem polju. Na telo poleg gravitacijske sile deluje zaviralna sila, za katero velja kvadratni zakon upora. Zapišimo drugi Newtonov zakon:

\[m\frac{\mathrm{d}^2s}{\mathrm{d}t^2}=mg-c\left(\frac{\mathrm{d}s}{\mathrm{d}t}\right)^2\]

Začetna pogoja sta \(s(0)=0\) in \(\left.\frac{\mathrm{d}s}{\mathrm{d}t}\right|_{t=0}=0\).

In [9]: 
t=linspace(0, 10, 1000)
m=1.0
g=9.81
c=0.01
def deriv(y, t): 
    return array([y[1], g - c / m * y[1] ** 2])

y0=array([0, 0])
y=odeint(deriv, y0, t)
s=y[:, 0]
v=y[:, 1]

fig, (axa, axb) = subplots(2, 1, sharex='col', sharey='row', figsize=(5, 6))

axa.plot(t, s);
axa.set_ylabel(r'Pot (m)')
axa.grid(alpha=0.5)
axa.annotate('(a)', xy=(0.03, 0.92), xycoords='axes fraction')

axb.plot(t, v);
axb.set_ylabel(r'Hitrost (m/s)')
axb.set_xlabel(r'\v Cas (s)')
axb.grid(alpha=0.5)
axb.annotate('(b)', xy=(0.03, 0.92), xycoords='axes fraction')

subplots_adjust(hspace=0.1)

Naloge

  1. Izračunaj in nariši, kako sta odklon iz ravnovesne lege in kotna hitrost matematičnega nihala odvisna od časa, če ob času \(t=0\) nihalo miruje, odklon od ravnovesne lege pa je a) \(30^\circ\) in b) \(150^\circ\) (nihalo naj ima togo vrvico). Primerjaj z analitičnim rezultatom za majhne odklone iz ravnovesne lege. Krožna frekvenca nihala naj bo \(1\,\)s\(^{-1}\).

  2. Izračunaj in nariši tira Zemlje in Lune pri gibanju okoli Sonca. Podatki: \(M_S=2\times 10^{30}\,\)kg, \(M_Z=6\times 10^{24}\,\)kg, \(M_L=7.35\times 10^{22}\,\)kg, \(d_{ZS}=1.5\times 10^{8}\,\)km, \(d_{LZ}=3.84\times 10^5\,\)km, \(T_Z=365.25\,\)dni in \(T_L=29.53\,\)dni. Začetni položaj Zemlje naj bo \(\mathbf{r}_Z(0)=(d_{ZS}, 0)\), Lune pa \(\mathbf{r}_L(0)=(d_{ZS}, d_{LZ})\). Začetna hitrost Zemlje naj bo \(\mathbf{v}_Z(0)=(0, 2\pi d_{ZS}/T_Z)\), Lune pa \(\mathbf{v}_L(0)=(-2\pi d_{LZ}/T_L, 2\pi d_{ZS}/T_Z)\). Ker ima Sonce bistveno večjo maso od Zemlje in Lune, lahko predpostaviš, da miruje v izhodišču koordinatnega sistema.