Численное интегрирование. Сеточные методы. Метод Монте-Карло

TODO: определение интеграла через предел суммы

Определённый интеграл

\[\int_a^b f(x) dx\]

Можно понимать как площадь под графиком функции f(x) от на отрезке [a, b]. Такая площадь понимается с положительным знаком, на участках, где f(x) > 0, и с отрицательным на участках f(x) < 0. При этом площадь можно посчитать приближённо, если достаточно хорошо подогнать под её форму некоторое множество геометрических фигур, площадь которых мы можем вычислить точно, и просуммировать их площади.

Для упрощения последующих вычислений разобъём отрезок [a, b] на N подотрезков длиной \(h = \frac{b - a}{N}\). Тогда начало i-го подотрезка находится в точке \(a + ih\), а конец в точке \(a + (i + 1)h\)

TODO: картинка

Метод левых прямоугольников состоит в том, что в качестве геометрической фигуры, дающей нам приближённую площадь под функцией, мы берём прямоугольник с высотой, равной значению f(x) в левой точке отрезка.

\[S_{left\ rect} = \sum_{i = 0}^{N - 1} h f(a + ih) = h \sum_{i = 0}^{N - 1} f(a + ih)\]

Ошибка метода (отличие приближённого значения интеграла от точного аналитического) в этом случае будет линейно зависеть от h.

\[| S_{left\ rect} - S_{exact} | = O(h) = O(\frac{1}{N})\]

При увеличении количества отрезков в 10 раз количество вычислений также увеличится в 10 раз, ошибка аналогично упадёт в 10 раз.

Метод правых прямоугольников обладает такой же ошибкой и отличается лишь тем, что в качестве высоты прямоугольника берётся значение функции на правом краю отрезка.

\[S_{right\ rect} = h \sum_{i = 1}^{N} f(a + ih)\]

Метод средних прямоугольников берёт высоту прямоугольника в средней точке отрезка и за счёт большей симметрии получает в два раза меньшую ошибку, однако она по-прежнему так же зависит от N:

\[S_{middle\ rect} = \sum_{i = 0}^{N - 1} h f(a + (i + 0.5) h) = h \sum_{i = 0}^{N - 1} f(a + (i + 0.5)h)\]\[| S_{middle\ rect} - S_{exact} | \approx 0.5 | S_{left\ rect} - S_{exact} |\]\[| S_{middle\ rect} - S_{exact} | = O(h) = O(\frac{1}{N})\]

Если подумать о том, как можно улучшить саму форму геометрической фигуры, которую мы подгоняем мод форму функции на каждом отрезке, можно придти к выводу, что следующим шагом на пути уменьшения ошибки будет использование трапеций. Это действительно так:

\[S_{trapezium} = \sum_{i = 0}^{N - 1} \frac{h}{2} \left(f(a + ih) + f(a + (i + 1)h)\right) = h \left(\frac{f(a)}{2} + \left( \sum_{i = 1}^{N - 1} f(a + ih) \right) + \frac{f(b)}{2}\right)\]\[| S_{trapezium} - S_{exact} | = O(h^2) = O(\frac{1}{N^2})\]

Как видим, тут мы получаем достаточно большое увеличение точности: количество подотрезков и вычислений увеличивается в N раз, ошибки уменьшается в \(N^2\).

Метод Симпсона

Следующим шагом для улучшения точности будет использовать в качестве фигуры фрагмент параболы. Для того, чтобы построить параболу, нам нужно три коэффициента, и три точки на каждом отрезке.

На отрезке длиной h выберем три точки: две на краях отрезка и одну по центру. Проведём через значения функции в этих точках параболу. Определив коэффициенты параболы и проинтегрировал её, получим формулу Симпсона:

\[S = \frac{h}{6} \left( f(-h/2) + 4 f(0) + f(h/2) \right)\]

Вернувшись к задаче интегрирования на отрезке из многих шагов и применив формулу Симпсона, получим:

\[S_{Simpson} = \sum_{i = 0}^{N - 1} \frac{h}{6} \left( f(a + ih) + 4 f(a + (i + 0.5) h)) + f(a + (i + 1) h) \right)\]

После раскрытия суммы и приведения однородных слагаемых получим:

\[S_{Simpson} = \frac{h}{6} \left( f(a) + 2 \sum_{i=1}^{N-1}f(a+ih) +4 \sum_{i=0}^{N-1}f(a+(i + \frac{1}{2})h) + f(b) \right)\]

Ошибка в методе Симпсона ведёт себя очень хорошо: падает в \(N^4\) при увеличении количества вычислений в N раз.

\[| S_{Simpson} - S_{exact} | = O(h^4) = O(\frac{1}{N^4})\]

Генераторы случайный чисел

Несмотря на то, что в разговоре часто упоминаются генераторы случайных чисел, важно понимать, что случайных чисел на самом деле в компьютере нет, просто потому, что процессор при одних и тех же входных данных и алгоритмах выдаёт всегда один и тот же ответ, а источников случайности в компьютере либо нет вообще, либо их недостаточно для получения большого количества случайных чисел. Например, температура на процессоре меняется с точностью до одной десятой градуса несколько раз в секунду, очевидно этого физического источника случайности мало для создания тысяч случайных чисел.

Большинство случайных чисел являются на самом деле псевдослучайными и получаются по жёстко заданным математическим формулам. Однако, эти формулы специально составлены математиками так, чтобы получаемая последовательность чисел была максимально похожа на случайную.

  1. Вероятность выпадения каждого числа из некоторого диапазона, в котором выбираются числа, должна быть одинакова. Это достаточно легко проверить, если запустить ГПСЧ столько раз, чтобы каждое число из диапазона было выброшено по несколько раз, и измерить вероятность выпадения каждого. Она должна быть приблизительно одинакова.

  2. Вероятность выпадения каждого числа должна слабо зависеть от предыдущих чисел.

  3. TODO рисунок с состоянием

    Если знать начальное состояние ГПСЧ, последовательность чисел можно восстановить.

  4. Так как состояние имеет конечный размер, то последовательность случайных числе периодически повторяется. Это следует из того, что при обходе конечного множества состояний одно из них неминуемо повторится, и после этого начнётся периодическое повторение последовательности, так как математическая формула, изменяющее состояние, неизменна.

ГПСЧ языка Си

TODO период

TODO независимых измерений

ГПСЧ MT19937 (Вихрь Мерсена 19937)

Период генератора составляет \(2^{19937}\) чисел.

Гарантируется независимость 623 сгенерированных случайных чисел.

MT19937 доступен отдельной библиотекой для языка Си, а так же встроен в некоторые другие языки как генератор по умолчанию.

Работа с ГПСЧ в Си

Напечатать 10 случайных чисел можно так:

#include <stdio.h>

int main(void)
{
    int i;
    for (i = 0; i < 10; ++i) {
        int x = rand();
        printf("%d\n");
    }
    return 0;
}

Прототип функции rand находится в stdlib.h и выглядит так:

int rand(void);

То есть, rand ничего не принимает и возвращает int. В документации к стандартной библиотеке можно прочесть, что rand возвращает случайное целое положительно число от нуля до константы RAND_MAX включительно. Константу RAND_MAX поменять нельзя, потому что она уже вкомпилирована в машинной код функции rand. Для того, чтобы наша программа печатала числа от 0 до 10 включительно, можно проделать трюк с остатком от деления:

#include <stdio.h>

int main(void)
{
    int i;
    for (i = 0; i < 10; ++i) {
        int x = rand() % 11; // Остаток от деления на 11 -- число от 0 до 10
        printf("%d\n");
    }
    return 0;
}

Этот метод хорошо работает только в том случае, если верхняя граница N требуемых случайных чисел достаточно мала по сравнению с RAND_MAX. Если же это не так, то числа будут распределены в диапазоне от 0 до N неравномерно.

Например, если RAND_MAX == 150 и N == 100, то очевидно, что вероятность выпадения чисел от 0 до 50 в два раза выше, чем чисел от 51 до 99.

Значение константы RAND_MAX на разных компиляторах различно, но оно не меньше 32767.

Если запустить один из тех двух примеров, которые мы только что обсудили, несколько раз, можно увидеть, что последовательности случайных чисел получаются одни и те же. Это происходит из-за того, что при старте программы на Си начальное состояние ГПСЧ всегда одно и то же. Если сделать это состояние зависимым от времени, мы сможем получать различные последовательности случайных чисел. Задать начальное состояние мы можем с помощью функции srand:

#include <stdio.h>

int main(void)
{
    srand(time(NULL)); // Вызов time(NULL) возвращает время в секундах,
            // прошедшее с 1 января 1970 года, которое является началом
            // отсчёта времени в большинстве современных компьютеров.
            // Его мы используем как начальное состояние ГПСЧ.
    int i;
    for (i = 0; i < 10; ++i) {
        int x = rand() % 11; // Остаток от деления на 11 -- число от 0 до 10
        printf("%d\n");
    }
    return 0;
}

Конечно, такая программа всё ещё будет выводить одинаковые последовательности, если её дважды запустить за одну секунду, но всё же это намного лучше, чем одна фиксированная последовательность.