Cvičení z Numerické matematiky 2

Metoda konečných diferencí

Definice. Nechť

x_{0} \in ℝ

H_{x_{0}}

je jeho okolí a

f, g : H_{x_{0}} \to ℝ

f

aproximuje

g

s přesností řádu

r \in ℕ

, pokud existuje

C \in ℝ^{+}

takové, že

\lim_{x \to x_{0}} \frac{| f (x) - g (x) |}{{| x - x_{0} |}^{r}} = C

Definice. Nechť

x_{0} \in ℝ

H_{x_{0}}

je jeho okolí,

f : H_{x_{0}} \to ℝ

h \in ℝ^{+}

l_{1}, l_{2} \in ℕ

. Potom budeme značit

$n ≔ l_{1} + l_{2}$
$x_{i} ≔ x_{0} + i h, i \in {- l_{1}, \dots, l_{2}}$
diference $f_{i} ≔ f (x_{i})$

Věta. Nechť

f \in 𝒞^{2} ((x_{0}, x_{1}))

. Potom

\frac{f_{1} - f_{0}}{h} = f^{'} (x_{0}) + 𝒪 (h)

Věta. Nechť

f \in 𝒞^{2} ((x_{- 1}, x_{0}))

. Potom

\frac{f_{0} - f_{- 1}}{h} = f^{'} (x_{0}) + 𝒪 (h)

Věta. Nechť

f \in 𝒞^{2} ((x_{- 1}, x_{1}))

. Potom

\frac{f_{1} - f_{- 1}}{2 h} = f^{'} (x_{0}) + 𝒪 (h)

Věta. Nechť

f \in 𝒞^{3} ((x_{- 1}, x_{1}))

. Potom

\frac{f_{1} - f_{- 1}}{2 h} = f^{'} (x_{0}) + 𝒪 (h^{2})

Věta. Nechť

f \in 𝒞^{4} ((x_{- 1}, x_{1}))

. Potom

\frac{f_{1} - 2 f_{0} + f_{- 1}}{h^{2}} = f^{''} (x_{0}) + 𝒪 (h^{2})

Věta. Nechť

f \in 𝒞^{(n + 1)} ((x_{- l_{1}}, x_{l_{2}}))

. Potom pro každé

k \in \hat{n}

lze

k

-tou derivaci aproximovat pomocí konečných diferencí s přesností řádu

n + 1 - k

Cauchyho úloha a Rungovy-Kuttovy metody

Oberhuberova prezentace

Mějme obyčejnou diferenciální rovnici s počáteční podmínkou:

x^{'} (t) = f (t, x (t)), t \in (0, T)

x (0) = x_{0}

Zkusíme ji vyřešit „simulací“ s diskrétními časovými intervaly. Zvolíme si časový krok $τ$ , tedy budeme mít $N ≔ \frac{T}{τ}$ kroků. Budeme počítat:

x_{i + 1} ≔ x_{i} + Δ x (i τ, τ)

Otázkou je, jak zvolit vhodné $Δ x$ . Chtěli bychom, aby platilo:

Δ x (t, τ) = x (t + τ) - x (t) \overset{\overset{}{Taylor}}{=} τ x^{'} (t) + \frac{τ^{2}}{2} x^{''} (t) + \dots

První derivaci si můžeme spočítat přímo ze zadání. S druhou už to bude horší:

x^{''} (t) = \partial_{1} f (t, x (t)) + \partial_{2} f (t, x (t)) \cdot x^{'} (t) = \partial_{1} f (t, x (t)) + \partial_{2} f (t, x (t)) \cdot f (t, x (t))

K tomu už potřebujeme znát derivace $f$ . A s dalšími derivacemi by to bylo ještě horší. Pokud se vykašleme na všechny derivace kromě první, dostaneme známou Eulerovu metodu. Ale tento postup je obecně dost numericky nestabilní. Radši použijeme Rungovu-Kuttovu metodu. Ta funguje tak, že přírůstky budeme hledat ve tvaru:

Δ x (t, τ) = \sum_{i = 1}^{n} p_{i} k_{i} (t, τ)

k_{1} (t, τ) ≔ τ f (t, x (t))

k_{i} (t, τ) ≔ τ f (t + α_{i} τ, x (t) + \sum_{j = 1}^{i - 1} β_{i, j} k_{j} (t, τ)), i > 1

Teď je samozřejmě otázkou, jak vhodně zvolit koeficienty $p_{i}, a_{i}, β_{i, j}$ . Obecně je chceme zvolit tak, aby se přírůstek shodoval s Taylorem do co nejvyššího stupně. Detaily viz prezentace.

Rungova-Kuttova metoda pro $n = 1$

Δ x (t, τ) = p_{1} k_{1} (t, τ) = p_{1} τ f (t, x (t))

\begin{aligned} φ_{1} (τ) & = x (t + τ) - x (t) - p_{1} τ f (t, x (t)) \\ = τ x^{'} (t) + \frac{τ^{2}}{2} x^{''} (s) - p_{1} τ f (t, x (t)) \\ = (1 - p_{1}) τ f (t, x (t)) + \frac{τ^{2}}{2} x^{''} (s) \end{aligned}

φ_{1}^{'} (τ) = (1 - p_{1}) f (t, x (t)) + τ x^{''} (s)

φ_{1}^{'} (0) = (1 - p_{1}) f (t, x (t)) = 0 ⟹ p_{1} = 1

Δ x (t, τ) = τ f (t, x (t))

Tím jsme opět dostali známou Eulerovu metodu.

Rungova-Kuttova metoda pro $n = 2$

Δ x (t, τ) = p_{1} k_{1} (t, τ) + p_{2} k_{2} (t, τ) = p_{1} τ f (t, x (t)) + p_{2} τ f (t + α_{2} τ, x (t) + β_{2,1} k_{1} (t, τ))

\begin{aligned} φ_{2} (τ) & = τ x^{'} (τ) + \frac{τ^{2}}{2} x^{''} (t) + \underset{z}{\underset{⏟}{\frac{τ^{3}}{6} x^{'''} (s)}} - p_{1} τ f (t, x (t)) - p_{2} τ f (t, x (t) + β_{2,1} k_{1} (t, τ)) \\ = τ f (t, x (t)) + \frac{τ^{2}}{2} (\partial_{t} f (t, x (t)) + \partial_{x} f (t, x (t)) f (t, x (t))) + z - p_{1} τ f (r, x (t)) - p_{2} τ f (t, x (t) + β_{2,1} τ f (t, x (t))) \end{aligned}

\begin{aligned} φ_{2}^{'} (τ) = f (t, x (t)) & + τ (\partial_{t} f (t, x (t)) + \partial_{x} f (t, x (t)) f (t, x (t))) \\ - p_{1} f (t, x (t)) \\ - p_{2} f (t, x (t) + β_{2,1} τ f (t, x (t))) \\ - p_{2} τ (\partial_{t} f (\dots) α_{2} + \partial_{x} f (\dots) β_{2,1} f (t, x (t))) \end{aligned}

φ_{2}^{''} (τ) = \partial_{τ} f (t, x (t)) + \partial_{x} f (t, x (t)) f (t, x (t)) - 2 p_{2} (\partial_{t} f (\dots) α_{2} + \partial_{x} f (\dots) β_{2,1} f (t, x (t))) - p_{2} τ \partial_{t} (\dots)

φ_{2}^{'} (0) = \dots \overset{♡}{=} 0 ∴ 1 - p_{1} - p_{2} = 0

φ_{2}^{''} (0) = \dots \overset{♡}{=} 0 ∴ 1 - 2 α_{2} p_{2} = 1 - 2 β_{2,1} p_{2} = 0

První možnost:

p_{1} = p_{2} = \frac{1}{2}, α_{2} = β_{2,1} = 1

Druhá možnost:

p_{1} = 0, p_{2} = 1, α_{2} = β_{2,1} = \frac{1}{2}

Vidíme, že pro každý řád se to počítá výrazně pracněji.

Cauchyho úloha vyššího řádu

x^{(n)} (t) = f (t, x (t), \dots, x^{(n - 1)} t)

Použijeme substituci, čímž to převedeme na soustavu:

u_{i} (t) ≔ x^{(i)} (t)

u_{n}^{'} (t) = f (t, u_{0} (t), \dots, u_{n - 1} (t))

u_{i - 1}^{'} (t) = u_{i} (t)

Musíme znát počáteční podmínky pro každé $u$ , tedy počáteční derivace až do řádu $n - 1$ .

Runge-Kuttovy metody můžeme přepsat do vektorového tvaru:

Δ \vec{u} (t, τ) = \sum_{i = 1}^{n} p_{i} {\vec{k}}_{i} (t, τ)

{\vec{k}}_{i} (t, τ) = τ \vec{F} (t + \sum_{j = 2}^{i} α_{j} τ, \vec{u} (t) + \sum_{j = 1}^{i - 1} β_{i, j} {\vec{k}}_{j} (t, τ))

Objektově orientované programování

Oberhuberova prezentace

Když chceme v C++ vykreslit graf funkce, normální člověk by napsal něco takového:

void drawFunction(std::function<double(double)> f, double left, double right, double step) {
  for (double x = left; x <= right; x += step) {
    std::cout << x << " " << f(x) << "\n";
  }
}

int main() {
  drawFunction([](double x) { return (2 * x + 3) * x + 1; }, 0.0, 1.0, 0.1);
  drawFunction([](double x) { return 2 * sin(5 * x + 1); }, 0.0, 1.0, 0.1);
}

Ale Oberhuber radši napíše tohle:

struct Function
{
    virtual double getFx( double x );
}

struct Quadratic : public Function    // !!!!
{
    double a, b, c;

    double getFx ( double x )
    {
        return ( a * x + b ) * x + c;
    }
}

struct Sinus : public Function    // !!!!
{
    double A, f, p;

    double getFx ( double x )
    {
        fx = A * sin ( f * x + p );
    }
}

void drawFunction( Function * function ,
                   double left , double right , double step )
{
    for ( double x = left ; x <= right ; x += step )
    {
        printf( " %f %f \n", x, function−>getFx( x ) );
    }
}

int main ( int argc , char * argv[] )
{
        /****
         * Evaluate f (x) = 2 x^2 + 3 x + 1 on <0,1>
         */
         Quadratic f1;
         f1.a = 2.0;
         f1.b = 3.0;
         f1.c = 1.0;
         drawFunction ( &f1, // WE DO NOT NEED TO CARE ABOUT FUNCTION TYPE
                        0.0, 1.0,     // on interval <0, 1>
                        0.1,          // we draw it with resolution 0.1
                                      //    −> 11 points on the interval
                         );

        /****
         * Evaluate f(x) = 2 sin( 5 *x + 1 ) on <0,1>
         */
         Sinus f2;
         f2.A = 2.0;
         f2.f = 5.0;
         f2.p = 1.0;
         drawFunction ( &f2, // WE DO NOT NEED TO CARE ABOUT FUNCTION TYPE
                        0.0, 1.0,     // on interval <0, 1>
                        0.1,          // we draw it with resolution 0.1
                                      //    −> 11 points on the interval
                         );

}

Cvičení z Numerické matematiky 2

Metoda konečných diferencí

Cauchyho úloha a Rungovy-Kuttovy metody

Rungova-Kuttova metoda pro n=1

Rungova-Kuttova metoda pro n=2

Cauchyho úloha vyššího řádu

Objektově orientované programování

Řešení okrajových úloh

Rungova-Kuttova metoda pro $n = 1$

Rungova-Kuttova metoda pro $n = 2$