Cvičení z Numerické matematiky 2

Metoda konečných diferencí

Oberhuberova prezentace

Cauchyho úloha a Rungovy-Kuttovy metody

Oberhuberova prezentace

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

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

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

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

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:

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

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
                         );

}

Řešení okrajových úloh

Oberhuberova prezentace

Mame-li úlohu se dvěma stupni volnosti, potřebujeme dvě počáteční podmínky. Místo toho ale často máme zadané okrajové podmínky. Z přednášky známe metodu střelby, ale ta nejde zobecnit do více dimenzí (parciální diferenciální rovnice). Místo toho použijeme metodu sítí: provedeme diskretizaci a vyřešíme soustavu lineárních rovnic. Detaily viz prezentace a přednáška.