Lab. 13 Równania różniczkowe - zagadnienie początkowe

Równania różniczkowe

Przykład:

$$\frac{dy}{dx} = 3 xy $$ $$x^2 y'' + axy' + y = 0$$

Równanie różniczkowe zwyczajne - jedna zmienna niezależna
Rząd równania - rząd najwyższej pochodnej

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

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 $$

Metoda Eulera

Wartości funkcji $y(x_i)$ w kolejnych punktach $x_i$ ze stałym krokiem $h = x_{i+1} -x_{i}$

$$ y_{i+1} = y_{i} + h \cdot f(x_i, y_i) $$

Interpretacja graficzna

$$ y(x+h)-y(h)=\int_x^{x+h} y^{\prime} d x=\int_x^{x+h} f(x, y) d x $$

Błędy

Błąd lokalny (pojedynczego kroku) - proporcjonalny do $h^2$

$$ E_a = \frac{f'(x_i,y_i)\cdot h^2}{2!}$$

Błąd globalny - akumuluje błędy z wszystkich kroków, można pokazać, że jest proporcjonalny do $h$

Im mniejszy krok tym lepsza dokładność metody ale dla bardzo małych $h$ zaczynają dominować błędy zaokrągleń

Przykład

Rozwiążmy problem początkowy $$\frac{dy}{dt} = y, \qquad y(0)=1 $$

Rozwiązanie analityczne jest postaci $y = e^t$

euler.c

#include <stdio.h>
#include <math.h>
 
double f(double t, double y)
{
   return y;
}
 
int main()
{
   double h=0.1, y, t;
   int i;
 
   t=0.0; y = 1.0;
   printf("t        Euler    Exact\n");
   for (i=0; i<= 10; i++)
   {
      printf("%lf %lf %lf\n", t, y, exp(t));
      y=y+h*f(t, y);
      t=t+h;
   }
   return 0;
}

Metoda midpoint

$$ y_{i+1} = y_i + h\cdot f\left(x_i + \frac{1}{2}h, \, y_i + \frac{1}{2} h f(x_i, y_i)\right)$$

Metoda Heuna

$$ y_{i+1} = y_i + h\cdot \left( \frac{1}{2} f(x_i, y_i) + \frac{1}{2} f\left(x_i + h, y_i + h f(x_i, y_i)\right) \right)$$

Metoda Ralstona

$$ y_{i+1} = y_i + h\cdot \left( \frac{1}{3} f(x_i, y_i) + \frac{2}{3} f\left(x_i + \frac{3}{4} h, y_i + \frac{3}{4}h f(x_i, y_i)\right) \right)$$

Zmodyfikowane metody Eulera

Metoda midpoint $$ y_{i+1} = y_i + h\cdot f\left(x_i + \frac{1}{2}h, \, y_i + \frac{1}{2} h f(x_i, y_i)\right)$$

Metoda Heuna

$$ y_{i+1} = y_i + h\cdot \left( \frac{1}{2} f(x_i, y_i) + \frac{1}{2} f\left(x_i + h, y_i + h f(x_i, y_i)\right) \right)$$

Metoda Ralstona

$$ y_{i+1} = y_i + h\cdot \left( \frac{1}{3} f(x_i, y_i) + \frac{2}{3} f\left(x_i + \frac{3}{4} h, y_i + \frac{3}{4}h f(x_i, y_i)\right) \right)$$

Zadanie 11

Treść zadania znajduje się w Moodle pod adresem Zadanie 11

Dodatkowe ćwiczenia

Ćw. 1 Rozwiąż dowolną metodą numeryczną równanie różniczkowe na odcinku $[0,1]$ z krokiem $h=0.1$ $$ y^{\prime}=-x^2 y, \quad y(0)=1, $$

Ćw. 2 Napisz program wyznaczający trajektorię układu Lorentza danego układem równań różniczkowych $$ \begin{aligned} \frac{d u}{d t} & =p(v-u), \\ \frac{d v}{d t} & =-u w+R u-v \\ \frac{d w}{d t} & =u v-b w . \end{aligned} $$

Przyjmij wartości współczynników $p=16$, $b=4$ i $R=35$. Warunki początkowe wynoszą $u=5$, $v=5$, $w=5$. Program prosi o podanie wielkości kroku $h$ oraz liczby iteracji $N$ po czym wypisuje kolejne współrzędne trajektorii uzyskane dowolną metodą rozwiązywania równań różniczkowych zwyczajnych.