Основы численных методов курсовая работа

Федеральное агентство по образованию

Национальный Исследовательский  Технологический Университет «МИСиС»


 

 

Институт информационных технологий и автоматизированных систем управления


 

 

Кафедра инженерной кибернетики


 

 

Курс:

Основы численных методов

 

Курсовая работа:

Метод Монте-Карло.

Применение метода для решения задач интегрирования

 

 

Исполнитель:

Сенченко Роман Владимирович,

студент группы ММ-09-01

Руководитель:

Гопенгауз Владимир Израилевич

 

 

 

 

 

 

 

 

 

Москва, 2012 г. 

Оглавление

Введение 3

§1. Моделирование случайных  величин 4

1.1. Генерация стандартной  случайной величины 4

1.2. Статистическая  проверка генератора 5

Длина апериодичности 5

Критерий Пирсона 6

Частотный тест 7

Тест серий 8

§2. Применение метода Монте-Карло для вычисления интегралов 9

2.1. Общий метод 9

Оценка математического  ожидания 9

Погрешность метода 9

2.2. Простейший метод  Монте-Карло 10

Одномерный определенный интеграл 11

2.3. Простейший метод  с повышенной скоростью сходимости 12

Двумерный определенный интеграл 17

Выводы 20

Список литературы 21

Приложения и листинги 22

Приложение 1. Листинг  генератора Д.Лемера 22

Приложение 2. Листинг  статистических функций 24

Приложение 3. Листинг  основной части программы 27

 

Введение

Методами Монте-Карло называется семейство методов решения численных  задач, основанных на моделирование  различных случайных величин. Главная  идея методов Монте-Карло заключается  в следующем: пусть поставлена некая  задача и скалярная величина является ее решением (это может быть решение СЛАУ или же значение интеграла , и т.д.). Постараемся придумать такую случайную величину , чтобы . Если это удастся сделать, то естественно ожидать, что решение поставленной задачи можно получить как при достаточно большом .

 

 

§1. Моделирование  случайных величин

1.1. Генерация стандартной случайной величины

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

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

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

Не вдаваясь в дальнейшие подробности, укажем сразу простую и достаточно хорошую функцию, предложенную Д.Лемером:

,

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

 

 

А именно:

Тройка чисел  являются параметрами функции. Генератор стандартной случайной величины, основанный на этой функции носит название генератор Д.Лемера.

1.2. Статистическая проверка  генератора

Разумеется, не для всех параметров генератор Д.Лемера, также называемый линейным конгруэнтным генератором, генерирует последовательности, обладающие хорошими свойствами. Обстоятельные исследования по этому поводу приводятся в [2]. Мы же лишь воспользуемся рекомендациями, данными в [1] для проверки качества генерируемой последовательности при уже выбранных параметрах генератора.

Длина апериодичности

Прежде, чем приступать к каким-либо проверкам получаемой последовательности, следует определить длину отрезка  апериодичности. Как легко видеть из самой формулировки алгоритма  Д.Лемера:

Числа берутся по модулю . Это означает, что гарантированно наступит момент времени, когда на очередной итерации (обозначим ее ) получится число равное числу ( ). Число называют длиной отрезка апериодичности.

Ясно, что использовать числа, полученные алгоритмом Д.Лемера при  , не имеет смысла.

Критерий Пирсона

Рассмотрим произвольную случайную  величину , которая может быть одномерной или многомерной, дискретной или непрерывной. Обозначим через множество всевозможных значений . Фиксируем какое-нибудь разбиение множества на попарно непересекающихся множеств таких, что при .

Выберем независимых значений величины и обозначим через количество значений, принадлежащих . Легко видеть, что .

В качестве меры отклонения «истинных» значений от «теоретических» удобно выбрать величину

Согласно теореме Пирсона, каковы бы ни были исходная величина и (такое, что все ), при каждом

,

Где плотность  , называется плотностью распределения с степенями свободы (аналитическую формулу для можно найти в [3]).

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

Строим величину для конкретной гипотезы (меняется разбиение и соответствующие вероятности ), вычисляем значение теоретической величины . Если , то принимаем гипотезу о законе распределения; иначе – отвергаем.

 

Для удобства проверки гипотезы о  стандартности чисел, генерируемых алгоритмом Д.Лемера, возьмем аппроксимацию  Голдштейна из [3]:

Частотный тест

Пусть первые десятичные разряды всех чисел, генерируемых генератором Д.Лемера на всем промежутке апериодичности .

Во-первых, если распределены равномерно, то частоты встречи цирф будут примерно равны. Аналогично, примерно равны будут и частоты независмых пар

Сосчитаем матрицу частот при : - частота встречи пары ; а также сосчитаем вектора , и . По этим данным посчитаем величины

  1. с 9-ю степенями свободы – тест проверки частот (frequency test)
  2. с 99-ю степенями свободы – тест проверки пар (serail test)

И проверим гипотезу о равномерности  распределения с помощью критерия .

Тест серий

Числа образуют серию длины , если . Пусть - это количество серий длины , а - количество всех остальных серий выборки . Общее количество серий обозначим через . Величина с степенями свободы вычисляется по формуле

,

 и 

 

Замечание 1: в конечном счете, эти тесты не гарантируют хорошего распределения вероятностей. Посему, стоит также построить гистограмму всех чисел, полученных на отрезке апериодичности и сравнить ее с ожидаемой.

 

Замечание 2: эмпирически удалось подобрать параметры генератора Д.Лемера, при которых длина отрезка апериодичности составляет и чиста обладают хорошими случайными свойствами.

 

§2. Применение метода Монте-Карло для вычисления интегралов

2.1. Общий метод

Оценка математического ожидания

Рассмотрим произвольную случайную  величину , у которой существует конечное математическое ожидание . Чтобы оценить величину , выберем независимых реализаций случайной величины и вычислим математическое ожидание этой выборки:

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

Более того, при достаточно больших  можно положить приближенное равенство

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

Погрешность метода

Предположим, дополнительно к 2.1., что  у случайной величины существует еще и конечная дисперсия. Тогда

И по центральной предельной теореме (а величина подчиняется ей, коль скоро ) имеем:

,

Положим, что  . Тогда из последнего соотношения получим, что

Следовательно, при достаточно больших  значениях 

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

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

2.2. Простейший метод Монте-Карло

Обозначим через  произвольную область (ограниченную или неограниченную, связную или несвязную), в которой определены функции и , причем,

.

Рассмотрим задачу о приближенном вычислении интеграла

.

Чтобы построить метод Монте-Карло  для рассчета интеграла  , рассмотрим случайную точку с плотностью вероятностей , а также, на ряду с ней, рассмотрим действительную случайную величину . Легко видеть, что

.

Но, согласно оценки из пункта 2.1, и действительна оценка точности вычисления из пункта 2.2.

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

Построим метод Монте-Карло на основе простейшего, позволяющий брать  интегралы вида

,

который всегда можно представить  в виде

,

где - плотность равномерно распределенной случайной величины на отрезке .

Согласно простейшему методу, построим случайную величину с плотностью , и сконструируем еще одну случайную величину . Тогда, с учетом оценок математического ожидания имеем

при достаточно больших  .

 

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

 

Численный пример 1

Пусть требуется взять интеграл

График подынтегральной функции

 

При , имеем

При , имеем

При , и выбранный генератор Д.Лемера не способен сгенерироваьт такую выборку.


 

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

2.3. Простейший метод с повышенной скоростью сходимости

Как было показано ранее, простейший метод в своей наивной реализации совершенно не пригоден для решения  численных задач. Напомним, что это  было связано с неравенством

,

имеющем место с вероятностью . Однако, как уже отмечалось, точность зависит от двух критериев величин; увеличение хоть и увеличило точность решения, но сделало расчет чрезвычайно долгим. Постараемся, теперь, уменьшить дисперсию .

Как и ранее, будем полагать, что  - равномерно распределенная в случайная величина. Ее дисперсия , откуда легко видно, что при . На этом простом факте основан метод существенной выборки. Опишем его в общем случае.

 

Вновь рассмотрим задачу приближенного  вычисления интеграла

,

Разобьем произвольно множество  ; далее, пусть

 и 

Очевидно, что  и . В каждой из областей будем рассматривать случайную точку с плотностью . Для оценки воспользуемся простейшим методом Монте-Карло: так как

,

то, выбрав независимых реализаций точки можем запистаь оценку:

 

 

 

Складывая такие оценки, получим  новую оценку для отыскиваемого  интеграла:

Для построения которой необходимо использовать количество случайных чисел.

 

Численный пример 2

Пересчитаем численный пример 1 с  использованием ускоренного метода

График подынтегральной функции

 

При и имеем

При и имеем

При и имеем


 

 

 

 

 

 

 

 

 

 

Численный пример 3

Пусть требуется взять интеграл

 

При и имеем


 

Численный пример 4

Пусть требуется взять интеграл

 

При и имеем

При и имеем

При и имеем


 

 

 

Численный пример 5

Пусть требуется взять интеграл

 

При и имеем

При и имеем

При и имеем


 

Численный пример 6

Пусть требуется взять интеграл

 

При и имеем

При и имеем

При и имеем


 

Замечание: как видно из примеров выше, метод существенной выборки заметно ускоряет сходимость простейшего метода Монте-Карло для вычисления одномерных определенных интегралов и позволяет уменьшить размеры выборки случайных чисел. Такая экономизация делает ускоренный метод крайне привлекательным для использования.

Двумерный определенный интеграл

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

Для этого, прямоугольную область  разобьем на частей по с шагом , и на частей по с шагом . Возьмем случайных величин , равномерно распределенных на отрезке ( ), а также случайных величин , равномерно распределенных на ( ). Случайная величина же будет иметь математическое ожидание:

при

 

 

 

 

 

 

 

 

Численный пример 7

Пусть требуется взять интеграл

 

При и имеем

При и имеем


 

Численный пример 8

Пусть требуется взять интеграл

 

При и имеем


 

 

 

 

 

Численный пример 9

Пусть требуется взять интеграл

 

При и имеем


 

Замечание: Ускоренный метод Монте-Карло замечательно сходится даже на двумерных интегралах!

 

Выводы

В данной работе мы рассмотрели методы генерации псевдослучайных чисел  при помощи генератора Д.Лемера, а  также статистические тесты, позволяющие  проводить контроль качества получаемых последовательностей.

Был рассмотрен метод Монте-Карло  в приложении к вычислению одномерных и двумерных интегралов. Построен общий ход решения методов  Монте-Карло, сделана вероятностная  оценка погрешности методов.

Построен и непосредственно  апробован на одномерном интеграле  простейший метод Монте-Карло, после  чего была выяснена его несостоятельность. Простейший метод был ускорен  с помощью метода существенной выборки, что дало приемлемые вычислительные затраты и хорошую точность на тестовых примерах вычисления одномерных и двумерных интегралов.

 

Список литературы

  1. Численные методы Монте-Карло, И.М. Соболь. Главная редакция физико-математической литературы изд-ва «Наука», 1973
  2. Дональд Кнут Искусство программирования, том 2. Получисленные методы = The Art of Computer Programming, vol.2. Seminumerical Algorithms. — 3-е изд. — М.: «Вильямс», 2007. — С. 832. — ISBN 0-201-89684-2
  3. Кобзарь А.И. Прикладная математическая статистика. Для инженеров и научных работников. – М.: ФИЗМАТЛИТ, 2006. – 816 с. – ISBN 5-9221-0707-0

 

 

Приложения и  листинги

Приложение 1. Листинг генератора Д.Лемера

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

 

namespace MonteCarlo

{

    /* Генератор реализации базовой случайной величины

     * (равномерно  распределенной в [0, 1])

     * на  базе алгоритма Д.Лемера

     */

    class LehmerGenerator

    {

        /* Конструктор генератора

         * (m0, M, g) - параметры генератора

         * capacity - объем генератора

         */

        public LehmerGenerator(ulong m0, ulong M, ulong g, ulong capacity)

        {

            m_m0 = m_m = m0;

            m_M = M;

            m_g = g;

 

            m_totalCapacity = m_capacity = capacity;

        }

 

        /* Установка генератора в начальное состояние

         */

        public void SetInitialState()

        {

            m_m = m_m0;

            m_capacity = m_totalCapacity;

        }

 

 

        /* Взятие очередной реализации базовой случайной величины

         */

        public double GetValue()

        {

            m_capacity--;

            if (m_capacity < 0)

            {

                throw new CapacityOverException();

            }

 

            m_m = (m_m * m_g) % m_M;

 

            return ((double)m_m) / ((double)m_M);

        }

 

        /* Получение оставшегося запаса генератора

         */

        public ulong GetCapacity()

        {

            return m_capacity;

        }

 

        /* Исключени: генератор исчерпан

         */

        public class CapacityOverException

            : Exception

        {

        };

 

        /* Текущее значение параметра m тройки (m, M, g)

         */

        private ulong m_m;

 

        /* Оставшийся объем генератора

         */

        private ulong m_capacity;

 

        /* Начальное значение параметра m0

         */

        private ulong m_m0;

 

        /* Параметр M тройки (m, M, g)

         */

        private ulong m_M;

 

        /* Параметр g тройки (m, M, g)

         */

        private ulong m_g;

 

        /* Начальный объем генератора

         */

        private ulong m_totalCapacity;

    }

}


 

 

 

Приложение 2. Листинг статистических функций

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

 

namespace MonteCarlo

{

    /* Класс-обертка статистических функций

     */

    static class Statistics

    {

        /* Распределение Пирсона

         * (аппроксимация Голдштейна квантилей хи-квадрат)

         * n - количество степеней свободы

         * alpha - коэффиицент значимости из [0.001, 0.999]

         */

        static public double PearsonDistribution(int n, double alpha)

        {

            double[] a = {1.0000886, 0.4713941, 0.0001348028, -0.008553069, 0.00312558, -0.0008426812, 0.00009780499};

            double[] b = {-0.2237368, 0.02607083, 0.01128186, -0.01156761, 0.005169654, 0.00253001, -0.001450117};

            double[] c = {-0.01513904, -0.008986007, 0.02277679, -0.01323293, -0.006950356, 0.001060438, 0.001565326};

           

            double d = 0;

            if (0.001 <= alpha && alpha < 0.5)

            {

                d = -2.0637 * Math.Pow(Math.Log(1 / alpha) - 0.16, 0.4274) + 1.5774;

            }

            else if (0.5 <= alpha && alpha <= 0.999)

            {

                d = 2.0637 * Math.Pow(Math.Log(1 / (1 - alpha)) - 0.16, 0.4274) - 1.5774;

            }

 

            double Chi = 0;

            for (int i = 0; i <= 6; i++)

            {

                Chi += Math.Pow(n, -i / 2.0) * Math.Pow(d, (double)i) *

                    (a[i] + b[i] / n + c[i] / Math.Pow((double)n, 2.0));

            }

            Chi = n * Math.Pow(Chi, 3.0);

 

            return Chi;

        }

 

        /* Частотный тест генератора Лемера:

         * проверка частот и проверка  частот пар

         * generator - генератор Лемера

         * alpha - коэффициент значимости

         */

        static public bool FrequencyTest(LehmerGenerator generator, double alpha)

        {

            generator.SetInitialState();

 

            // мощность генератора

            ulong N = generator.GetCapacity();

 

            // полная табличка частот

            ulong[,] V = new ulong[10, 10];

            for (ulong i = 0; i < N / 2; i++)

            {

                V[(ulong)(generator.GetValue() * 10), (ulong)(generator.GetValue() * 10)]++;

            }

 

            // вектора частот пар

            ulong[] a = new ulong[10];

            ulong[] b = new ulong[10];

            ulong[] c = new ulong[10];

            for (int i = 0; i < 10; i++)

            {

                for (int j = 0; j < 10; j++)

                {

                    a[i] += V[i, j];

                    b[i] += V[j, i];

                }

                c[i] = a[i] + b[i];

            }

 

            // проверка пар

            double empiricalChi1 = 0;

            for (int i = 0; i < 10; i++)

            {

                empiricalChi1 += Math.Pow(c[i] - 0.1 * N, 2);

            }

            empiricalChi1 *= 10.0 / (double)N;

 

            // проверка частот

            double empiricalChi2 = 0;

            for (int i = 0; i < 10; i++)

            {

                for (int j = 0; j < 10; j++)

                {

                    empiricalChi2 += Math.Pow(V[i, j] - 0.005 * N, 2.0);

                }

            }

            empiricalChi2 *= 50.0 / (double)N;

 

            generator.SetInitialState();

 

            // вывод ответа

            if (empiricalChi1 < PearsonDistribution(9, alpha) && empiricalChi2 < PearsonDistribution(99, alpha))

            {

                return true;

            }

            else

            {

                return false;

            }

        }

 

        /* Тест серий генератора Лемера

         * generator - генератор Лемера

         * m - количество степеней свободы

         * alpha - коэффициент значимости

         */

        static public bool RunTest(LehmerGenerator generator, int m, double alpha)

        {

            generator.SetInitialState();

 

            // полная матрица частот

            ulong[] V = new ulong[m + 1];

            for (int i = 0; i <= m; i++)

            {

                V[i] = 0;

            }

 

            ulong l = 1;

            int currentValue = 0,

                previousValue = (int)(generator.GetValue() * 10);

            bool endFlag = false;

            while (endFlag == false)

            {

                if (generator.GetCapacity() > 0)

                {

                    currentValue = (int)(generator.GetValue() * 10);

                }

                else

                {

                    endFlag = true;

                }

 

                if (currentValue == previousValue && endFlag == false)

                {

                    l++;

                }

                else

                {

                    if (l > (ulong)m)

                    {

                        V[m]++;

                        l = 1;

                    }

                    else

                    {

                        V[l - 1]++;

                        l = 1;

                    }

                    previousValue = currentValue;

                }

            }

 

            // общее количество серий

            ulong n = 0;

            for (int i = 0; i <= m; i++)

            {

                n += V[i];

            }

 

            // проверка серий

            double empiricalChi = 0;

            for (int i = 0; i < m; i++)

            {

                double p = 9 * Math.Pow(10.0, -(double)(i + 1));

                empiricalChi += Math.Pow(V[i] - n * p, 2.0) / (n * p);

            }

            empiricalChi += Math.Pow(V[m] - n * Math.Pow(10.0, (double)(-m)), 2.0) /

                (n * Math.Pow(10.0, (double)(-m)));

 

            generator.SetInitialState();

 

            // вывод ответа

            if (empiricalChi < PearsonDistribution(m, alpha))

            {

                return true;

            }

            else

            {

                return false;

            }

        }

    }

}


 

 

Приложение 3. Листинг основной части программы

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.IO;

 

namespace MonteCarlo

{

    class Program

    {

        /* Делегаты одномерной и двумерной функций

         */

        public delegate double Function1(double x);

        public delegate double Function2(double x, double y);

 

        /* Метод простейшего интегрирования по методу Монте-Карло

         * f - одномерная функция, интегрируемая на отрезке [a, b]

         * generator - генератор Д.Лемера, источник стандартных случайных чисел

         * n - мощность выборки

         */

        public static double SimpleIntegrate(Function1 f, double a, double b, LehmerGenerator generator, ulong n)

        {

            // случайная величина, равномерно распределенная в [a, b]

            UniformVariates Q = new UniformVariates(generator, a, b);

 

            // вычисление интеграла

            double I = 0;

            for (ulong i = 0; i < n; i++)

            {

                I += f(Q.GetValue());

            }

            I *= (b - a) / n;

 

            return I;

        }

 

        /* Метод простейшего интегрирования по методу Монте-Карло

         * для одномерной функции

         * (ускоренный методом существенной выборки)

         * f - одномерная функция, интегрируемая на отрезке [a, b]

         * generator - генератор Д.Лемера, источник стандартных случайных чисел

         * m - мощность разбиения

         * n - мощность выборки

         */

        public static double AcceleratedIntegrate1(Function1 f, double a, double b, LehmerGenerator generator, int m, ulong n)

        {

            // шаг разбиения

            double h = (b - a) / m;

 

            // вычисление интеграла

            double I = 0;

            for (int i = 0; i < m; i++)

            {

                // случайная величина, равномерно распределенная в m-ом отрезке

                double alpha = a + i * h,

                       beta = a + (i + 1) * h;

                UniformVariates Q = new UniformVariates(generator, alpha, beta);

 

                // вычисление частичного интеграла

                double S = 0;

                for (ulong j = 0; j < n; j++)

                {

                    S += f(Q.GetValue());

                }

                S *= (beta - alpha) / n;

                I += S;

            }

 

            return I;

        }

 

        /* Метод простейшего интегрирования по методу Монте-Карло

         * для двумерной функции

         * (ускоренный методом существенной выборки)

         * f - двумерная функция, интегрируемая на отрезке [a, b]

         * generator - генератор Д.Лемера, источник стандартных случайных чисел

         * m - мощность разбиения

         * n - мощность выборки

         */

        public static double AcceleratedIntegrate2(Function2 f, double a, double b,

            double c, double d, LehmerGenerator generator, int mx, int my, ulong n)

        {

            // шаг разбиения

            double hx = (b - a) / mx,

                   hy = (d - c) / my;

 

            // вычисление интеграла

            double I = 0;

            for (int i = 0; i < mx; i++)

            {

                for (int j = 0; j < my; j++)

                {

                    // случайная величина, равномерно распределенная в m-ом квадратике

                    double alpha = a + i * hx,

                           beta = a + (i + 1) * hx,

                           gamma = c + j * hy,

                           delta = c + (j + 1) * hy;

 

                    UniformVariates Qx = new UniformVariates(generator, alpha, beta);

                    UniformVariates Qy = new UniformVariates(generator, gamma, delta);

 

                    // вычисление частичного интеграла

                    double S = 0;

                    for (ulong k = 0; k < n; k++)

                    {

                        S += f(Qx.GetValue(), Qy.GetValue());

                    }

                    S *= (hx * hy) / n;

                    I += S;

                }

            }

 

            return I;

        }

 

        /* Интегрируемые функции одной переменной

         */

        public static double I1(double x)

        {

            return Math.Pow(x, 3) * Math.Sin(x + 1);

        }

 

        public static double I2(double x)

        {

            return Math.Cos(x) / Math.Pow(x, 2) + Math.Pow(x, 3) * Math.Sin(x);

        }

 

        public static double I3(double x)

        {

            return 8 * Math.Pow(x, 5) - 10 * Math.Pow(x, 4) + 5 * Math.Pow(x, 3) + x - 18;

        }

 

        public static double I4(double x)

        {

            return x / Math.Sqrt(x + 1);

        }

 

        public static double I5(double x)

        {

            return Math.Exp(x) * Math.Sin(10 * x);

        }

 

        /* Интегрируемые функции двух переменных

         */

        public static double J1(double x, double y)

        {

            return Math.Pow(x, 2) + Math.Pow(y, 2) * Math.Sin(x);

        }

 

        public static double J2(double x, double y)

        {

            return Math.Pow(x, 2) - Math.Pow(y, 2) - Math.Cos(3 * x) + y * Math.Sin(4 * x);

        }

 

        public static double J3(double x, double y)

        {

            return Math.Pow(x, 3) - y + Math.Cos(x) * Math.Exp(y) + Math.Sqrt(x * y);

        }

 

        static void Main(string[] args)

        {

            ulong m = 1,

                  M = (ulong)Math.Pow(2, 26),

                  g = (ulong)Math.Pow(5, 17),

                  L = (ulong)Math.Pow(2, 24);

 

 

            LehmerGenerator G = new LehmerGenerator(m, M, g, L);

 

            Console.WriteLine(AcceleratedIntegrate2(J3, 0, 7, 0, 7, G, 120, 120, 5000));

        }   

    }

}