Автор работы: Пользователь скрыл имя, 08 Апреля 2013 в 23:07, курсовая работа
Методами Монте-Карло называется семейство методов решения численных задач, основанных на моделирование различных случайных величин. Главная идея методов Монте-Карло заключается в следующем: пусть поставлена некая задача и скалярная величина является ее решением (это может быть решение СЛАУ или же значение интеграла , и т.д.). Постараемся придумать такую случайную величину , чтобы . Если это удастся сделать, то естественно ожидать, что решение поставленной задачи можно получить как при достаточно большом .
Введение 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
Замечание: не смотря на свою простоту, простейший метод требует слишком большие мощности выборки , что с одной стороны одновременно создает необходимость иметь генератор случайных чисел с большим отрезком апериодичности, а также требует внушительных временных затрат на построение оценки . Картину дополняет вероятностный характер оценки точности: существует не нулевая вероятность, что заданная рассчетная точность не будет достигнута даже на самых качественных случайных числах и нужном объеме выборки .
Как было показано ранее, простейший
метод в своей наивной
имеющем место с вероятностью . Однако, как уже отмечалось, точность зависит от двух критериев величин; увеличение хоть и увеличило точность решения, но сделало расчет чрезвычайно долгим. Постараемся, теперь, уменьшить дисперсию .
Как и ранее, будем полагать, что - равномерно распределенная в случайная величина. Ее дисперсия , откуда легко видно, что при . На этом простом факте основан метод существенной выборки. Опишем его в общем случае.
Вновь рассмотрим задачу приближенного вычисления интеграла
Разобьем произвольно
Очевидно, что и . В каждой из областей будем рассматривать случайную точку с плотностью . Для оценки воспользуемся простейшим методом Монте-Карло: так как
то, выбрав независимых реализаций точки можем запистаь оценку:
Складывая такие оценки, получим новую оценку для отыскиваемого интеграла:
Для построения которой необходимо использовать количество случайных чисел.
Численный пример 2 Пересчитаем численный пример 1 с использованием ускоренного метода График подынтегральной функции
При и имеем При и имеем При и имеем |
|
Численный пример 3 Пусть требуется взять интеграл
При и имеем |
|
Численный пример 4 Пусть требуется взять интеграл
При и имеем При и имеем При и имеем |
|
Численный пример 5 Пусть требуется взять интеграл
При и имеем При и имеем При и имеем |
|
Численный пример 6 Пусть требуется взять интеграл
При и имеем При и имеем При и имеем |
Замечание: как видно из примеров выше, метод существенной выборки заметно ускоряет сходимость простейшего метода Монте-Карло для вычисления одномерных определенных интегралов и позволяет уменьшить размеры выборки случайных чисел. Такая экономизация делает ускоренный метод крайне привлекательным для использования.
Построим расчетную схему ускоренного метода Монте-Карло с использованием равномерно распределенных случайных величин для приближенного вычисления интеграла
Для этого, прямоугольную область разобьем на частей по с шагом , и на частей по с шагом . Возьмем случайных величин , равномерно распределенных на отрезке ( ), а также случайных величин , равномерно распределенных на ( ). Случайная величина же будет иметь математическое ожидание:
при
|
Численный пример 7 Пусть требуется взять интеграл
При и имеем При и имеем |
|
Численный пример 8 Пусть требуется взять интеграл
При и имеем |
|
Численный пример 9 Пусть требуется взять интеграл
При и имеем |
Замечание: Ускоренный метод Монте-Карло замечательно сходится даже на двумерных интегралах!
В данной работе мы рассмотрели методы
генерации псевдослучайных
Был рассмотрен метод Монте-Карло в приложении к вычислению одномерных и двумерных интегралов. Построен общий ход решения методов Монте-Карло, сделана вероятностная оценка погрешности методов.
Построен и непосредственно апробован на одномерном интеграле простейший метод Монте-Карло, после чего была выяснена его несостоятельность. Простейший метод был ускорен с помощью метода существенной выборки, что дало приемлемые вычислительные затраты и хорошую точность на тестовых примерах вычисления одномерных и двумерных интегралов.
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; } } |
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; } } } } |
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(Function { // шаг разбиения 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(Function 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( } } } |
Информация о работе Основы численных методов курсовая работа