\documentclass[10pt,a4paper]{article}

%% Polonizacja
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}
\usepackage{polski}
\usepackage{indentfirst}

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts}
\usepackage{latexsym}

%% Różne
%\usepackage{listings}
%\lstloadlanguages{C}
%\usepackage{comment}
%\usepackage{enumerate}
%\usepackage[dvips]{hyperref}


%% Numerowanie równań
\numberwithin{equation}{section}

%% Twierdzenia
\newtheorem{thm}{Twierdzenie}[section]
\newtheorem{df}{Definicja}[section]

%% Własne polecenia
\DeclareMathOperator{\sign}{sign}

\newcommand*{\en}{\enspace}
\newcommand*{\cc}[1]{{#1}^{\ast}}
\newcommand*{\n}[1][n]{$#1$\nobreakdash-\hspace{0pt}}
\newcommand*{\abs}[1]{\lvert #1 \rvert}
\newcommand*{\norm}[1]{\ensuremath{\lVert #1 \rVert}}
\newcommand*{\ABS}[1]{\ensuremath{\left\lvert #1 \right\rvert}}
\newcommand*{\NORM}[1]{\ensuremath{\left\lVert #1 \right\rVert}}
\newcommand*{\field}[1]{\mathbb{#1}}
\newcommand*{\vect}[1]{\ensuremath{\mathbf{#1}}}
\newcommand*{\matr}[1]{\ensuremath{\mathbf{#1}}}

\begin{document}

\title{Przykłady równań różniczkowych zwyczajnych}
\author{Jakub Narębski \and Tomasz Pawłowski \and Tomasz Werner}
\maketitle



\section{Oscylator Duffinga z~wymuszeniem}
\label{sec:duffing}

\begin{equation}
  \label{eq:duffing}
  \ddot{x} + k \dot{x} + x^3 = A \cos \omega t
\end{equation}

Start z~$x(0) = \dot{x}(0) = 0$, dla parametrów $\omega = 1$, $k =
0.05$.  Zbadać zachowanie rozwiązania w~zależności od wartości
parametru $A$ dla $A \in [1, 12]$, wykreślając wykres w~przestrzeni
fazowej dla $t \in [0, 100 \pi]$.  Zobaczyć dla jakiego $A$ występują
gwałtowne zmiany wyglądu trajektorii.

Można także wziąć ogólniejszą postać równania Duffinga
\begin{equation}
  \label{eq:duffing-more}
  \ddot{x} + \gamma \dot{x} + \omega_0^2 x + \varepsilon x^3 
  = A \cos \omega t
\end{equation}



\section{Oscylator van der Pola z~wymuszeniem}
\label{sec:van-der-pol}

\begin{equation}
  \label{eq:van-der-pol}
  \ddot{x} - \varepsilon (1 - \alpha x^2) \dot{x} + \omega_0^2 x
  = \omega^2 A \cos \omega t
\end{equation}

Warto zbadać zarówno zachowanie z, jak i~bez siły wymuszającej.  Przy
braku siły wymuszającej tzn. dla $A = 0$ warto wziąć $\alpha =
1/\varepsilon$ i~$\omega_0 = 1$ i~zbadać zachowanie układu dla
wartości parametru $\varepsilon$ zmieniających się od 0.01 do 4.
Zobaczyć, że niezależnie od warunków początkowych dochodzimy do ruchu
quasiperiodycznego.

Przy występowaniu wymuszenia, dla małej amplitudy siły wymuszającej
\begin{math}
  \dfrac{A \omega}{\varepsilon (2/\sqrt{\alpha})} \ll 1
\end{math}
dostajemy ruch quasiperiodyczny, zmieniający się łagodnie ze zmianą
parametrów.  Dla dużej amplitudy siły wymuszającej 
\begin{math}
  \dfrac{A \omega}{\varepsilon (2/\sqrt{\alpha})} \gg 1
\end{math}
dostajemy ruch okresowy.  Pojawia się on w~sposób nagły przy wzroście
$A$.  

Jeśli dodamy do równania~\eqref{eq:van-der-pol} nieliniową sztywność
typu Duffinga to otrzymamy równanie van der Pola -- Duffinga
o~postaci
\begin{equation}
  \label{eq:van-der-pol--duffing}
  \ddot{x} 
  - \varepsilon (1 - \alpha x^2) \dot{x} 
  + \omega_0^2 (1 + \gamma x^2) x
  = \omega^2 A \cos \omega t
\end{equation}
Równanie to przechodzi od zachowania quasiharmonicznego, przez orbity
okresowe po zachowanie chaotyczne przy zmianie parametrów.



\section{Wahadło z~napędem}
\label{sec:wahadlo}

Przez wybór odpowiednich jednostek $t$ oraz $A$ i~$B$ można sprowadzić
równanie wahadła matematycznego z~wymuszeniem do postaci
\begin{equation}
  \label{eq:wahadlo}
  \ddot{\theta} + \gamma \dot{\theta} + \sin \theta
  = A \cos \omega t + B 
\end{equation}
Wprowadźmy oznaczenie
\begin{equation*}
  \Omega = \frac{2 \pi}{\omega} \frac{A + B}{\gamma}
\end{equation*}
Dla wartości parametru $A = 1.0$, $\omega = 1.76$ mamy: (a) dla
$\gamma = 1.576$, $\Omega = 1.4$ ruch okresowy, (b) dla $\gamma =
1.253$, $\Omega = 1.2$ przejście do chaosu, (c) dla $\gamma = 1.081$,
$\Omega = 1.094$ zachowanie chaotyczne.  Jeśli położymy $B = 0$, to
zachowanie wahadła w~funkcji parametrów $A$ i~$\omega$ miota się
miedzy ewolucja regularną i~chaotyczną.  Wrażliwośc na warunki
początkowe pojawia się, gdy wahadło przechodzi przez punkt szczytowy
z~bardzo małą prędkością.



\section{Równanie drapieżnik -- ofiara}
\label{sec:prey--predator}

Równanie
\begin{equation}
  \label{eq:prey--predator}
  \left\{
  \begin{aligned}
    \dot{X} &= \alpha X  - \beta X Y \\
    \dot{Y} &= \beta X Y - \gamma Y
  \end{aligned}
  \right.
\end{equation}
nazywamy równaniem drapieżnik--ofiara (ang: \emph{prey--predator});
$X$~oznacza (odpowiednią unormowaną) liczbę zajęcy, zaś $Y$~oznacza
liczbę lisów.  Dla parametrów $\alpha = \beta = \gamma = 1$ równanie
to ma punkt stały $X = Y = 1$.  Można zobaczyć, że jest to punkt
stabilny dla cyklu; dla dowolnych parametrów początkowych (np. $X(0) =
Y(0) = \frac{1}{4}$) dostajemy w~przestrzeni fazowej krzywą
zamkniętą. 



\section{Równanie Lorentza}
\label{sec:lorentz}

Równanie Lorentza powstało jako obcięcie (uproszczenie) równania ruchu
dla rolki konwekcyjnej (w~badaniach dynamiki atmosfery).  Równanie to
ma następująca postać:
\begin{equation}
  \label{eq:lorentz}
  \left\{
  \begin{aligned}
    \dot{x}_1 &= \sigma (x_3 - x_1)\\
    \dot{x}_2 &= r x_1 - x_2 - x_1 x_2 \\
    \dot{x}_3 &= - b x_3 + x_1 x_2
  \end{aligned}
  \right.
\end{equation}
gdzie $\sigma$ jest liczbą Prandtla i~jest proporcjonalne do lepkości,
$b$ jest proporcjonalne do rozmiarów komórki konwekcyjej zaś $r$ jest
parametrem kontrolnym (proporcjonalnym do przewodnictwa cieplnego).
Równanie to ma trzy punkty stałe:
\begin{subequations}
  \label{eq:lorentz:punkty-stale}
  \begin{align}
    \label{eq:lorentz:punkty-stale:1}
    x_1 &= x_2 = x_3 = 0 \\
    \label{eq:lorentz:punkty-stale:2}
    x_1 &= x_2 = \pm \sqrt{b(r - 1)},\quad x_3 = r - 1
    \tag{\ref{eq:lorentz:punkty-stale}bc}
  \end{align}
\end{subequations}
Punkt stały $x_1 = x_2 = x_3 = 0$ odpowiadający przewodnictwu cieplemy
bez konwekcji jest punktem stabilnym dla~$0 < r < 1$; dla $r > 1$
punkt ten staje się punktem siodłowym.  Drugi i~trzeci punkt krytyczny
(które istnieją dla $r \ge 1$) jest stabilny dla $1 < r < r_1$, dla
$r_1 < r < r_c$ mamy dwa cykle graniczne, zaś dla $r > r_c$ pojawia
się chaos.  Krytyczna wartośc parametru kontrolnego dana jest wzorem
\begin{equation}
  \label{eq:lorentz:r_c}
  r_c = \frac{\sigma(\sigma + b + 3)}{\sigma - b - 1}
\end{equation}
Dla najczęsciej stosowanych wartości parametrów $\sigma = 10, b =
\frac{8}{3}$ dostajemy $r_1 \simeq 13.54$, zaś $r_c = 24.7368$.
Trajektorie chaotyczne istnieją do $r \simeq 27.74$; najczęściej do
badania chaotycznego zachowania równania Lorentza bierze się $r =
24.74$.



\section{Zredukowane zagadnienie trzech ciał}
\label{sec:satelita}

Zredukowane zagadnienie trzech ciał to równania ruchu ciała próbnego
(satelity) w~polu grawitacyjnym dwu dużo cięższych od niego obiektów;
na tyle ciężkich że ruch ciało próbne nie wpływa na ich ruch
i~te dwa ciała poruszają się po elipsach (w~szczególnym przypadku: po
okręgach).  Jeżeli dwa cięższe ciała poruszają sie po okręgach wokół
wspólnego środka masy to ruch satelity można rozważać w~nieinercjalnym
układzie w~którym dwa cięższe ciała spoczywają.  Oznaczmy cięższe
ciała indeksami 1 i~2.  Niech $M_i$ oznacza masę \n[i]tego ciała,
a~$\mu_i$ jego masę zredukowaną
\begin{align}
  \label{eq:masa-zredukowana}
  \mu_i &= \frac{M_i}{M_1 + M_2}
\end{align}
Wówczas w~układzie w~którym ciała 1 i~2 spoczywają mają one
współrzędne $(-\mu_2, 0)$ i~$(\mu_1, 0)$.

Dobierając odpowiednio jednostki dostajemy następujące równanie ruchu
dla satelity:
\begin{equation}
  \label{eq:satelita}
  \left\{
  \begin{aligned}
    \ddot{x} &= x + 2\dot{y} 
    - \frac{\mu_2}{r_1^2} \frac{(x + \mu_1)}{r_1^{\vphantom{2}}}
    - \frac{\mu_1}{r_2^2} \frac{(x - \mu_2)}{r_2^{\vphantom{2}}} \\
    \ddot{y} &= y - 2\dot{x} 
    - \frac{\mu_2}{r_1^2} \frac{y}{r_1^{\vphantom{2}}}
    - \frac{\mu_1}{r_2^2} \frac{y}{r_2^{\vphantom{2}}} 
  \end{aligned}
  \right.
\end{equation}
gdzie $r_i$ jest odległością ciała $i$ od położenia satelity.

Dla układu Ziemia--Księżyc 
\begin{math}
  \mu_{\text{Z}} 
  = \dfrac{M_{\text{K}}}{M_{\text{K}} + M_{\text{K}}} 
  \approx \dfrac{1}{82.45}
\end{math}.
Jeśli się zacznie z~warunku początkowego $x(0) = 1.2, y(0) = 0$,
$\dot{x}(0) = 0, \dot{y}(0) = -1.0493570983$ to wychodzi krzywa
zamknięta z~okresem $T = 6.192169331396$.

Do sprawdzenia dokładności rachunków można zastosować badanie stałości
całki Jacobiego
\begin{equation}
  \label{eq:satelita:jacobi}
  J = \frac{1}{2} [ \dot{x}^2 + \dot{y}^2 - x^2 - y^2 ] 
  - \frac{\mu_1}{r_1} - \frac{\mu_2}{r_2}
\end{equation}
Dla powyższych danych początkowych (i~wartości parametrów) jest ona
równa $J = -1.041589$.

\emph{Uwaga:} do analizy tego zagadnienia należy zastosowac metode ze
zmiennym krokiem, gdyż wynik składa się z~częsci wolnozmiennej (dla
której potrzebna byłaby bardzo duża liczba kroków gdyby krok był mały)
i~szybkozmiennej (dla której potrzebny jest bardzo mały krok aby nie
utracić dokładności).


\end{document}
