Lab. 14 Równania różniczkowe c.d.

  • Metody Rungego-Kutty
  • układy równań różniczkowych
  • równania różniczkowe wyższego rzędu

Zagadnienie początkowe

Równanie różniczkowe zwyczajne rzędu pierwszego $$ \frac{dy}{dx} = f(x, y)$$ warunek początkowy $$ y(x_0) = y_0 $$

Szukamy funkcji $y(x)$. Rozwiązanie numeryczne: wartości funkcji $y(x_i)$ w węzłach $x_i$

Metody Rungego-Kutty

Ogólna postać $$y_{n+1}=y_{n}+\sum_{i=1}^{s} w_{i} k_{i}$$ $$\begin{array}{l} k_{1}=h f\left(x_{n}, y_{n}\right) \\ k_{i}=h f\left(x_{n}+a_{i} h, y_{n}+\sum_{j=1}^{i-1} b_{i j} k_{j}\right), \, i>1 \\ \end{array}$$ gdzie $ w_{i}, a_{i}, b_{ij}$ stałe

  • Metoda Eulera, pierwszego rzędu, błąd obcięcia $O(h)$
  • Metoda Heuna, midpoint, Ralstona, drugiego rzędu, błąd obcięcia $O(h^2)$
  • Metada Rungego-Kutty 4 rzędu, błąd obcięcia $O(h^4)$, najlepszy stosunek między dokładnością a liczbą obliczeń

Układy równań różniczkowych

Szukamy $y_{i}=y_{i}(x), \quad i=1,2, \ldots, m$

spełniających

$$\left\{\begin{array}{l} \frac{\mathrm{d} y_{1}}{\mathrm{~d} x}=f_{1}\left(x, y_{1}, \ldots y_{m}\right) \\ \frac{\mathrm{d} y_{2}}{\mathrm{~d} x}=f_{2}\left(x, y_{1}, \ldots y_{m}\right) \\ \ldots \\ \frac{\mathrm{d} y_{m}}{\mathrm{~d} x}=f_{m}\left(x, y_{1}, \ldots, y_{m}\right) \end{array}\right.$$

z warunkiem początkowym

$$y_{i}\left(x_{0}\right)=y_{i, 0}, \qquad i=1,2, \ldots, m$$

Rozwiązanie układu równań różniczkowych za pomocą metod R-K

  1. wyznaczamy $k_1$ dla wszystkich $m$ równań
  2. wyznaczamy $k_2$ dla wszystkich $m$ równań
  3. …
  4. mając wszystkie $k_i$ wyznaczamy przybliżenie $y_i$ dla każdego z równań

Przykład: model epidemii

Model SIR (S - susceptible, I - infected, R - recovered)

$$\begin{aligned} \frac{dS}{dt} & = -\alpha \cdot S \cdot I \\ \frac{dR}{dt} & = I(t) \big/ T \\ \frac{dI}{dt} & = \alpha \cdot S \cdot I - I \big/T \\ \end{aligned}$$

  • $\alpha$ prawd. zarażenia przy spotkaniu dwóch osobników
  • $T$ średni czas zdrowienia

Załóżmy, że $\alpha=0.0002, \, T=12$ dni

Przyjmijmy warunki początkowe $ I(0) = 1, R(0)=0, S(0)=1000$.

Użyjmy metody Heuna $$\begin{array}{c} y_{i+1}=y_{i}+\left(\frac{1}{2} k_{1}+\frac{1}{2} k_{2}\right) h \\ k_{1}=f\left(x_{i}, y_{i}\right) \\ k_{2}=f\left(x_{i}+h, y_{i}+k_{1} h\right) \end{array}$$

do wyznaczenia $I$ oraz $S$ w czasie $$\begin{aligned} \frac{dS}{dt} & = f_S(t, S, I)\\ \frac{dI}{dt} & = f_I(t, S, I)\\ \end{aligned}$$

$$\begin{aligned} k_{S1} & = f_S(t_i, S_i, I_i) \\ k_{I1} & = f_I(t_i, S_i, I_i) \\ k_{I2} & = f_I(t_i + h, S_i + k_{S1}h, I_i + k_{I1}h) \\ k_{S2} & = f_S(t_i + h, S_i + k_{S1}h, I_i + k_{I1}h) \\ S_{i+1} & =S_{i} + \left(\frac{1}{2} k_{S1}+\frac{1}{2} k_{S2}\right) h \\ I_{i+1} & =I_{i} + \left(\frac{1}{2} k_{I1}+\frac{1}{2} k_{I2}\right) h \\ \end{aligned}$$

Przykład w C

sir.c
#include <stdio.h>
#include <math.h>

double a = 0.0002;
double T = 12;

double f_s(double t, double s, double i)
{
   return - a*s*i;
}

double f_i(double t, double s, double i)
{
   return  a*s*i - i/T;
}

int main()
{
   double h=0.05, t0=0.0, tk=200, t;
   double n = 1000;    /* populacja */
   double i0 = 1;      /* jeden zainfekowany */ 
   double s0 = n - i0; /* liczba podatnych */ 
   double i = i0, s = s0;
   double k_s1, k_s2, k_i1, k_i2;
   
   printf("# t      S         I       R\n");
   for (t=t0; t<= tk; t=t+h)
   {
      printf("%6.2f  %6.1f  %6.1f  %6.1f\n", t, s, i, n-s-i);
      
      k_s1 = f_s(t, s, i);
      k_i1 = f_i(t, s, i);
      
      k_s2 = f_s(t+h, s + k_s1 * h, i + k_i1 * h);
      k_i2 = f_i(t+h, s + k_s1 * h, i + k_i1 * h);
      
      s = s + 0.5 * (k_s1 + k_s2 ) * h;
      i = i + 0.5 * (k_i1 + k_i2 ) * h;      
   }

   return 0;
}

Równania wyższych rzędów

Równanie różniczkowe rzędu $n$

$$y^{(n)}=f\left(x, y, y^{\prime}, \ldots, y^{(n-1)}\right)$$

można zamienić na układ $n$ równań pierwszego rzędu podstawiając

$$\begin{array}{lllll} y_{0}=y & y_{1}=y^{\prime} & y_{2}=y^{\prime \prime} & \ldots & y_{n-1}=y^{(n-1)} \end{array}$$

otrzymujemy

$$y_{0}^{\prime}=y_{1} \\ y_{1}^{\prime}=y_{2} \\ \ldots \\ y_{n}^{\prime}=f\left(x, y_{0}, y_{1}, \ldots, y_{n-1}\right)$$

Zadanie

Napisz program rozwiązujący równanie różniczkowe opisujące drgania tłumione sprężyny $$ \frac{d^{2} x}{d t^{2}}+2 \beta \frac{d x}{d t}+\omega_{0}^{2} x=0$$ z warunkiem początkowym $x(t=0)=1, \quad x^{\prime}(t=0)=0 $.
Zaimplementuj rozwiązanie z użyciem metody Rungego-Kutty czwartego rzędu z krokiem całkowania $h=0.01$.
Przyjmij wartość współczynnika tłumienia $\beta = 1.8$, częstotliwość drgań własnych $\omega_0 = 10 $ Hz.
Dla podanej przez użytkownika wartości czasu końcowego $T_k$ program wypisuje w kolejnych liniach wartości czasu $t_i$, amplitudy $x(t_i)$ oraz prędkości $x^{\prime}(t_i)$ w kolejnych węzłach całkowania poczynając od $t=0$ kończąc na $t=T_k$.

Metody Rungego-Kutty rzędu IV

$$\begin{array}{l} y_{i+1}=y_{i}+\frac{1}{6}\left(K_{1}+2 K_{2}+2 K_{3}+K_{4}\right) h \\ K_{1}=f\left(x_{i}, y_{i}\right) \\ K_{2}=f\left(x_{i}+\frac{1}{2} h, y_{i}+\frac{1}{2} K_{1} h\right) \\ K_{3}=f\left(x_{i}+\frac{1}{2} h, y_{i}+\frac{1}{2} K_{2} h\right) \\ K_{4}=f\left(x_{i}+h, y_{i}+K_{3} h\right) \end{array}$$

Dodatkowe ćwiczenia

Ćw. 1