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:

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:
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.
  1. Napisz funkcję:
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.
  1. 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]\).

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

  3. 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.

  4. 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.

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