\documentclass[10pt,a4paper]{article}

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

%% Wzorki matematyczne
\usepackage{amsmath}

\newcommand{\vect}[1]{\ensuremath{\mathbf{#1}}}

\begin{document}

\title{
  Zadania do Metod Numerycznych II A\\
  Grupy 2, 3. Seria 2.\\[1 em]
  \large
  Równania różniczkowe zwyczajne}
\author{Jakub Narębski}
\date{03.04.2001}
\maketitle

We wszystkich programach powinna być możliwość zadania przez
użytkownika kroku całkowania, czy to przez jawne jego zadanie, czy
przez podanie zakresu i~liczby kroków.

Poniżej przypomnienie niektórych wzorów na metody jednokrokowe:
\begin{itemize}
\item Metoda Eulera --- pierwszego rzędu
  \begin{equation}
    \label{eq:Euler}
    u_{n+1} = u_n + \Delta t\: f(u_n, t_n)
  \end{equation}

\item Ulepszona metoda Eulera (metoda Heuna) --- drugiego rzędu
  \begin{equation}
    \label{eq:Heun}
    u_{n+1} = u_n + \tfrac{1}{2} \Delta t\, 
    \bigl( f(u_n, t_n) + f(u_n + \Delta t\: f(u_n, t_n), t_n + \Delta t) \bigr)
  \end{equation}

\item Zmodyfikowana metoda Eulera (metoda Runge--Kutty 2. rzędu)
  \begin{equation}
    \label{eq:RK2}
    u_{n+1} = u_n + \Delta t \,
    f\bigl(u_n + \tfrac{1}{2} \Delta t\: f(u_n, t_n), t_n + \tfrac{1}{2} \Delta t \bigr)
  \end{equation}
\end{itemize}

Wykresy w~programie \texttt{gnuplot} robimy za pomocą polecenia
\begin{flushleft}
  \hspace{2 em}
  \texttt{plot} \textit{"nazwapliku"} \texttt{using} 1:2 
\end{flushleft}
gdzie 1:2 oznacza, że wartości $x$ mają być brane z~pierwszej kolumny
w~pliku o~nazwie \texttt{nazwapliku}, zaś $y$ z~drugiej kolumny (są to
wartości domyślne).  Warto przejrzeć skrypt \texttt{wykresik.sh},
dostępny na mojej stronie domowej.  W~programie tym dostępna jest
wbudowany system pomocy, wywoływany za pomoca polecenia \texttt{help},
np. \texttt{help plot}.



\section{Ruch w polu siły lepkości}
\label{sec:lepki}

Należy rozwiązać zagadnienie \textbf{ruchu w~polu siły oporu lepkiego}, t.j.
\begin{equation}
  \label{eq:lepki}
  \left\{
    \begin{aligned}
      \frac{dv(t)}{dt} &= - \alpha v(t) \\
      v(0) &= v_0
    \end{aligned}
  \right.
\end{equation}
Podać po jakim czasie cząstka osiągnie ustaloną prędkość (tzn. jej
prędkość $v$ będzie równa $v_k$). Zamiast tego można podać jaką
prędkość będzie miała cząstka po zadanym czasie. Prędkość początkowa
oraz współczynnik oporu lepkiego powinny być zadawane przez
użytkownika bez potrzeby rekompilowania programu.

Zadanie należy rozwiązać za pomocą metody Eulera, ulepszonej metody
Eulera (metody Heuna), zmodyfikowanej metody Eulera (metody
Runge-Kutty 2. rzędu) oraz klasycznej metody Runge-Kutty 4. rzędu.
Należy zrobić wykres rozwiązania od czasu oraz błędu rozwiązania
numerycznego (różnicy między rozwiązaniem numerycznym a~analitycznym)
od czasu. Pod Linuksem można to zrobić zapisując dane do pliku,
a~następnie robiąc wykres za pomocą programu \texttt{gnuplot}.

Uwaga: ponieważ równanie to jest równaniem liniowym można łatwo
zastosować metody niejawne, np. metodę Eulera wstecz 
\(\vect{u}_{n+1} = \vect{u}_n + \Delta t \: \vect{u}_{n+1}\)


\section{Wahadło matematyczne i~oscylator harmoniczny}
\label{sec:wahadlo-matematyczne}

Należy metodą Runge--Kutty rozwiązać równanie \textbf{wahadła
  matematycznego} (zrobić symulację).  Wahadłem matematycznym nazywamy
punkt materialny umiesczony na końcu sztywnego, nieważkiego pręta
w~polu stałej siły (siły grawitacji).  Siła dana jest wzorem $F(x) = -
k' sin(x)$.  Zagadnieniu temu odpowiada następujący układ równań
różniczkowych zwyczajnych
\begin{equation}
  \label{eq:wahadlo}
  \left\{
    \begin{aligned}
      \frac{d v}{d t} &= -k \sin(x) \\
      \frac{d x}{d t} &= v
    \end{aligned}
  \right.
  \quad\text{gdzie }
  x \in [-\pi,\pi]
\end{equation}
Należy zbadać zachowanie układu dla \(x_0 = 0\), \(k = 1\) i~różnych
wartości początkowych prędkości \(v_0\).  Maksymalny krok to \(0.1\)
okresu drgań równoważnego oscylatora harmonicznego (czyli \(0.1\cdot 2
\pi\) dla \(k = 1\)), minimalny czas symulacji to pięć okresów (czyli
należy podać wyniki od \(t = 0\) do \(t = 5\cdot 2\pi\), przy \(k =
1\)).  

Wyniki powinny być przedstawione na wykresie położenia i~prędkości od
czasu ($x(t)$ vs. $t$, $v(t)$ vs. $t$) oraz tzw. wykresie fazowym:
zależność prędkości od położenia, czyli wykresie parametrycznym
($v(t)$ vs $x(t)$).

Wymagane jest \emph{przynajmniej} wypisywanie na ekran (bądź do pliku,
do dalszej obróbki, np. przez \texttt{gnuplota}) wartości $t, x(t),
v(t)$.  Wypisywanie powinno następować co ileś kroków, bądź co zadany
okres czasu. Można również zadać ile wartości ma wypisać (zapisać).

\medskip

Dla małych wielkości prędkości początkowej układ powinien się
zachowywać podobnie jak oscylator harmoniczny: wykresy $x(t)$ i~$v(t)$
powinny być odpowiednio sinusoidą i~cosinusoidą, zaś na wykresie
fazowym powinno być widać okrąg. Przy zbyt dużej wartości kroku
(a~zatem dużym błędzie numerycznym) okrąg będzie się rozszerzał.
Metoda Eulera jest niestabilna już dla oscylatora harmonicznego, więc
nie należy jej stosować.

Dla dużych prędkości początkowych prędkość będzie oscylowała, ale
nigdy nie będzie zmieniać znaku, zaś kąt będzie rósł; stąd potrzebna
jest normalizacja kąta (zmiennej $x$) do przedziału $[-\pi,\pi]$,
np. wzorem \(x = \bigl((x + \pi) \bmod 2 \pi\bigr) - \pi\)).

\bigskip

Nalezy także rozwiązać za pomocą tego samego (bądź zmodyfikowanego)
programu równanie \textbf{oscylatora harmonicznego}.  Użyć metody
Eulera oraz metody Runge--Kutty 4. rzędu.  Porównać obie metody,
korzystając z~analitycznego (ścisłego) rozwiązanie.  W~tym zagadnieniu
nie jest potrzebna normalizacja położenia do przedziału $[-\pi, \pi]$
i~nie należy jej wykonywać.  Wyniki przedstawić jak dla wahadła
matematycznego.



\section{Metody wielokrokowe jawne}
\label{sec:wielokrokowe}

Zbadać zbieżność i koszt (ilość obliczeń funkcji prawej strony
równania) dla zagadnienia o znanym rozwiązaniu, np. oscylatora
harmonicznego bądź równania ruchu wzmocnionego $u' = \alpha u$, gdzie
$\alpha > 0$.  Uwaga: w~zagadnieniu ruchu wymuszonego dobrać
odpowiednio mały współczynnik $\alpha$ oraz $t_{\text{kon}}$.  

Proszę także rozwiązać zagadnienie:
\begin{equation}
  \label{eq:jakies}
  \left\{
    \begin{aligned}
      \frac{d u}{d t} &= \left(\frac{2 \log(u) + 8}{t} - 5\right)u \\
      u(1) &= 1
    \end{aligned}
  \right.
\end{equation}
dla $t \in [1, 6]$, z~warunkiem początkowym $u(1) = 1$.  Ścisłe
rozwiązanie tego równania na zadanym przedziale dane jest wzorem
\begin{equation*}
  u(t) = \exp\left( -t^2 + 5 t - 4 \right)
\end{equation*}
Należy zrobić wykres rozwiązania od czasu (czy to bezpośrednio w
programie, czy za pomocą programu \texttt{gnuplot}).  Porównać
dokładność i~ilość obliczeń funkcji z~klasyczną metodą Runge--Kutty.
Porównać błąd dla różnych wielkości kroku i~ewentualnie wyliczyć
(wyznaczyć) rząd metody (podpowiedź: $\log(x^n) = n \log(x)$).
\emph{Uwaga:} metoda Milne'a nie jest stabilna, i~jej rząd wyliczany
na podstawie równania oscylatora harmonicznego może być niższy niż 4.
Może ona nawet dawać błędne wyniki dla tego równania.


\bigskip

\textbf{\textit{Dla chętnych:}} metody typu
\textbf{predyktor--korektor}. Porównać właściwości z~metodami
wielokrokowymi jawnymi, stosując je do tych samych zagadnień.

\vfill

\begin{flushright}
  Jakub Narębski
\end{flushright}

\end{document}


