info2_lab05.Rmd

---
number: 5
course: Informatyka 2
material: Instrukcja 5
---

# Całkowanie numeryczne układów równań różniczkowych zwyczajnych

Pliki do wykorzystania w poniższym ćwiczeniu można pobrać za pomocą poniższych linków:

 - [Plik nagłówkowy rk4.h](http://ccfd.github.io/courses/code/info2/rk4.h)
 - [Plik źródłowy rk4.cpp](http://ccfd.github.io/courses/code/info2/rk4.cpp)

## Wstęp

Celem ćwiczenia jest zastosowanie metody **Eulera** oraz metody **Rugego-Kutty 4 rzędu** do numerycznego rozwiązania równań ruchu dynamiki Newtona.
Jako przykład takiego zagadnienia posłuży nam wahadło matematyczne.

## Równania ruchu

Ruch wahadła matematycznego najwygodniej jest opisać w układzie współrzędnych biegunowych związanych z jego punktem zaczepienia.
Otrzymamy wtedy równanie różniczkowe wraz z warunkami początkowymi:
$$
\left\{
\begin{array}{l}
  \dfrac{d^2\alpha}{dt^2} = -\dfrac{g}{l}\sin(\alpha) \\
  \alpha(t_0) = \alpha_0 \\
  \dfrac{d\alpha}{dt}(t_0) = \omega_0
\end{array}
\right.
$$
gdzie:

- $\alpha$ - kąt wychylenia wahadła z położenia równowagi,
- $g$ - przyśpieszenie ziemskie,
- $l$ - długość wahadła,
- $m$ - masa kulki zaczepionej na końcu wahadła,
- $\alpha_0$ - początkowe wychylenie wahadła,
- $\omega_0$ - początkowa prędkość wahadła.

Równanie to możemy sprowadzić do układu równań różniczkowych zwyczajnych pierwszego rzędu za pomocą podstawienia:
$$
\dfrac{d\alpha}{dt} = \omega
$$
Układ równań ma teraz postać:
$$
\tag{*}
\left\{
\begin{array}{l}
  \dfrac{d\omega}{dt} = -\dfrac{g}{l}\sin(\alpha) \\
  \dfrac{d\alpha}{dt} = \omega \\
  \alpha(t_0) = \alpha_0 \\
  \omega(t_0) = \omega_0
\end{array}
\right.
$$

## Rozwiązanie układu równań różniczkowych metodą Eulera

Mamy układ dwóch równań różniczkowych pierwszego rzędu:
$$
\left\{
\begin{array}{l}
  \dfrac{d\omega}{dt} = F_1(\alpha, \omega, t) \\
  \dfrac{d\alpha}{dt} = F_2(\alpha, \omega, t) \\
\end{array}
\right.
$$
Szukanymi funkcjami są $\omega = \omega(t)$ oraz $\alpha = \alpha(t)$.
Układ ten można rozwiązać metodą Eulera.
Jedna iteracja całkowania z krokiem $h$ będzie miała postać:
\begin{equation}
\left\{
\begin{array}{l}
  \omega(t_{i+1}) = \omega(t_i) + h \cdot F_1(\alpha(t_i), \omega(t_i), t_i) \\
  \alpha(t_{i+1}) = \alpha(t_i) + h \cdot F_2(\alpha(t_i), \omega(t_i), t_i) \\
\end{array}
\right.
\end{equation}

## Ćwiczenia

Dla wahadła opisanego układem równań (*):

1. Napisz funkcję o nagłówku:
```c++
void rhs_fun(double t, double *X, double *F);
```
która oblicza wartości prawych stron równań różniczkowych.
Argumenty funkcji to:

- `t` - zmienna niezależna (czas),
- `X` - tablica zmiennych zależnych ($\alpha$ i $\omega$),
- `F` - tablica do której zapisane zostaną obliczone prawe strony równań różniczkowych.

2. Napisz funkcję:
```c++
void veuler(double t, double *X, double h, int n,
            void (* fun)(double, double *,double *), double *X1);
```
która wykonuje jeden krok całkowania układu równań różniczkowych zwyczajnych pierwszego rzędu metodą Eulera.
Argumenty funkcji to:

- `t` - zmienna niezależna,
- `X` - tablica wartości zmiennych zależnych w kroku `t`,
- `h` - krok całkowania,
- `n` - rozmiar tablicy,
- `fun` - wskaźnik do funkcji obliczającej prawe strony równań,
- `X1` - tablica do której zapisane zostaną wartości zmiennych zależnych w kroku `t + h`.

2. Napisz program, który używając metody Eulera wyznacza zależności kąta wychylenia wahadła $\alpha$ oraz prędkości kątowej $\omega$ od czasu dla $t \in \left[0, \ldots 10\right]$.

3. Narysuj wykres trajektorii układu w przestrzeni fazowej ($\alpha-\omega$).

4. Powtórz obliczenia korzystając z metody Rungego-Kutty 4 rzędu, która jest zaimplementowana w bibliotece `rk4.cpp`.
Odpowiednia funkcja nazywa się `vrk4` a jej nagłówek jest analogiczny do nagłówka funkcji `veuler`.

5. Wyznacz zależność energii całkowitej wahadła od czasu $E(t)$.
Energia całkowita wahadła wyraża się wzorem:
$$
E = \frac{ml^2}{2}\left(\frac{d\alpha}{dt}\right)^2 + mgl(1-\cos(\alpha))
$$
**Uwaga**: Przy braku dyssypacji, energia mechaniczna powinna być stała.

6. Powtórz obliczenia dla różnych kroków czasowych.