Ирина Эланс

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

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

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ВОЛОГОДСКИЙ ГОСУДАРСТВЕННЫЙ ПЕДАГОГИЧЕСКИЙ УНИВЕРСИТЕТ

ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ И КОМПЬЮТЕРНЫХ ТЕХНОЛОГИЙ

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

Выполнил:

Студент 3-го курса ФПМиКТ

Гудков Т.Н.

Научный руководитель:

Кандидат физико-математических наук, доцент А.С. Сипин

2012 г.

г. Вологда

Содержание

Введение…………………………………………………………………..3

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

1.1. Распределение случайной величины…………………………4

1.2. Моделирующие формулы……………………………………..5

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

2.1. Общая схема метода Монте-Карло…………………………...6

2.2. Оценка погрешности метода…………………………………..7

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

3.1. Вычисление стандартных интегралов………………………...8

3.2. Вычисление кратных интегралов...…....……………………...9

3.3. Решение интегрального уравнения Вольтерра………………10

Приложения

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

Монте-Карло (Delphi)…………………………………………..11

Программа для решения интегрального уравнения Вольтерра

методом Монте-Карло (Delphi)…………….………………….20

Заключение……………………………………………………………....22

Литература…………………………………………………………….....23

Введение

Метод Монте-Карло – это численный метод, известный с 50-х годов нашего века и имеющий широкое применение в различных областях математики, физики, химии и многих других. Изначально метод использовался при разработке водородной бомбы. Идея была развита польским математиком Станиславом Уламом, именно он предложил использовать ЭВМ для расчётов методом Монте-Карло. Метод носит название в честь района в стране Монако, который известен своими казино. Суть метода заключается в использовании датчика псевдослучайных чисел, аналогичного рулетке. Одно из важнейших применений метода – вычисление интегралов и решение интегральных уравнений. В настоящее время метод Монте-Карло продолжает развиваться и становится более эффективным за счёт появления параллельного программирования.

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

Распределение случайной величины

Пусть – случайная величина, тогда .

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

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

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

Одно из самых часто используемых распределений случайных величин – равномерное распределение на отрезке . Для моделирования случайной величины, распределённой по такому закону, не нужен специальный алгоритм моделирования, поскольку датчик псевдослучайных чисел изначально выдаёт случайные величины из отрезка . Но если нас интересует случайная величина , равномерно распределённая на отрезке , где и - любые действительные числа ( < ), то потребуется простой алгоритм моделирования, который сгенерирует случайную величину при помощи случайной величины , равномерно распределённой на . Моделирующая формула выглядит следующим образом: .

1.2 Моделирующие формулы

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

Распределение

Формула плотности

Моделирующая формула

Равномерное непрерывное распределение на

Стандартное нормальное распределение

Экспоненциальное (показательное) распределение

Гамма- распределение

Бета-распределение

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

Общая схема метода Монте-Карло

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

Результат вычислений напрямую зависит от объёма выборки, то есть от количества проведённых испытаний. Чем больше проводится испытаний, тем точнее результат. В теории метода можно узнать о том, как более разумно выбрать случайные величины, как находить её значения. Одна из важнейших задач данного метода – нахождение способов уменьшения дисперсии, которые позволяют уменьшить получаемую в расчётах ошибку.

Оценка погрешности метода

Пусть было проведено n независимых испытаний, по результатам которых было найдено выборочное среднее . Оно было принято в качестве оценки математического ожидания: . Если повторить эти же n испытаний, то получится уже другая оценка . Из этого можно сделать вывод, что точную оценку мат. ожидания получить просто невозможно, поэтому в процессе расчётов методом Монте-Карло оценивают допущенную ошибку.

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

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

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

Вычисление стандартных интегралов

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

После выбора распределения, необходимо смоделировать случайную величину по этому распределению. Приближённое значение определённого интеграла вычисляется по формуле , где - сгенерированное значение случайной величины, смоделированной по выбранному распределению, - плотность выбранного распределения, n – объём выборки.

Вычисление кратных интегралов

Для получения приближённого значения кратного интеграла вида моделируется s-мерный случайный вектор. На каждом шаге значения компонент , смоделированного по выбранному распределению вектора, подставляются в подынтегральную функцию. Поскольку смоделированные случайные величины независимы, то совместная плотность выбранного распределения от случайного вектора равна произведению плотностей от значений каждой компоненты: . Итоговая формула для вычисления кратного определённого интеграла методом Монте-Карло выглядит следующим образом: , где - значения компонент -мерного случайного вектора, смоделированного по выбранному распределению при испытании, – кратность вычисляемого интеграла, n – объём выборки.

Решение интегрального уравнения Вольтерра

Интегральное уравнение Вольтерра имеет вид , где и - известные функции, . Для нахождения приближённого значения функции используется схема Неймана-Улама. Также, как и при вычислении интегралов методом Монте-Карло, проводятся n испытаний, в каждом из которых вычисляется значение функции и приближённо находится значение интеграла по формуле . Здесь - равномерно распределённая на случайная величина, где , – первый номер, для которого выполнено неравенство . Общая формула для решения интегрального уравнения Вольтерра методом Монте-Карло имеет вид , где – это я траектория, на которой оценивается интеграл.

Приложение 1

Программа для вычисления кратных интегралов методом Монте-Карло (Delphi)

program CalcMultipleIntegrals;

{$APPTYPE CONSOLE}

uses

SysUtils;

const

nmax = 10; {максимально возможная кратность интеграла}

type

tBorder = record

bottom: LongInt; {нижняя граница интегрирования}

top: LongInt; {верхняя граница интегрирования}

end;

tParameters = array [1..nmax] of double; {параметры подынтегральной функции}

tBorders = array [1..nmax] of tBorder; {массив границ кратного интеграла}

var

alpha, alpha1, alpha2: double;

vars: tParameters;

{функция для возведения числа в степень}

function Pow(x, y: double): double;

begin

if (x < 0) then

Pow := (-1) * exp(y * ln(abs(x)))

else if (x > 0) then

Pow := exp(y * ln(abs(x)))

else

Pow := 0;

end;

{функция для вычисления факториала}

function Factorial(n: longint): longint;

var

i, Factor: longint;

begin

Factor := 1;

if n > 0 then

begin

for i := 1 to n do

Factor := Factor * i;

Factorial := Factor;

end

else

Factorial := 1;

end;

{плотность равномерного распределения}

function pUniformDist(k: LongInt; bord: tBorders): double;

begin

pUniformDist := 1 / (bord[k].right - bord[k].left)

end;

{функция для моделирования нормально распределённой с.в.}

function randomValueNormDist: double;

begin

alpha1 := Random;

alpha2 := Random;

randomValueNormDist := Sqrt(-2*ln(alpha1)) * cos(2*Pi*alpha2);

end;

{плотность нормального распределения}

function pNormDist(x: double): double;

begin

pNormDist := 1 / sqrt(2 * pi) * exp(-x * x / 2);

end;

{функция для моделирования с.в., распределённой экспоненциально}

function randomValueExpDist(L: double): double;

begin

randomValueExpDist := ((-1/L) * ln(alpha));

end;

{плотность экспоненциального распределения}

function pExpDist(x: double; L: integer): double;

begin

pExpDist := L * exp(-L * x);

end;

{гамма-функция}

function Gamma(n: LongInt): LongInt;

begin

Gamma := Factorial(n-1);

end;

{функция для моделирования с.в., распределённой по гамма-распределению}

function randomValueGammaDist(A: longint; B: double): double;

var

i: longint;

sum: double;

begin

sum := 0;

for i := 1 to A do

begin

alpha := random;

sum := sum + randomValueExpDist(1/B);

end;

randomValueGammaDist := sum;

end;

{плотность гамма-распределения}

function pGammaDist(x: double; A: longint; B: double): double;

begin

if x >= 0 then

pGammaDist := (Pow(x, A - 1)*exp(-x/B))/(Gamma(A)*Pow(B, A))

else

pGammaDist := 0;

end;

{функция для моделирования с.в., распределённой по бета-распределению}

function randomValueBetaDist(m: longint; p: double): double;

var

i: longint;

ksi: double;

begin

ksi := 1;

for i := 1 to m do

begin

alpha := random;

ksi := ksi * pow(alpha, 1 / (p + i - 1));

end;

randomValueBetaDist := ksi;

end;

{бета-функция}

function Beta(x, y: longint): double;

begin

Beta := Gamma(x) * Gamma(y) / Gamma(x + y);

end;

{плотность бета-распределения}

function pBetaDist(x: double; p, m: longint): double;

begin

pBetaDist := pow(x, p - 1) * pow(1 - x, m - 1) / Beta(p, m);

end;

{подынтегральная функция}

function f(p: tParameters; bord: tBorders; num: LongInt):double;

var

i: LongInt;

rez: Double;

begin

rez := p[1]*p[1]+p[2]*p[2]*p[3]; {f(x, y) = x*x + y*y*y}

for i := 1 to num do begin

if not((bord[i].bottom<=p[i]) and (p[i]<=bord[i].top)) then

{если значение одной из переменных функции выходит за границы интегрирования, то функция обнуляется}

rez := 0;

end;

f := rez;

end;

{начало основной программы}

var

i, j:LongInt; {счётчики}

r: LongInt; {номер выбранного распределения}

n: LongInt; {объём выборки}

mul: LongInt; {кратность интеграла}

L: LongInt; {параметр экспоненциального распределения}

p, m: LongInt; {параметры бета-распределения}

W, B: longint; {параметры гамма-распределения}

sum, sum2: double; {суммы вычисленных значений и квадратов этих значений}

disp: Double; {дисперсия}

pMult: double; {произведение плотностей распределения компонент вектора}

val: Double; {вычисленное значение интеграла за одно испытание}

ksi: tParameters; {случайный вектор}

bord: tBorders; {массив границ интегрирования}

begin

randomize;

writeln('Выберите распределение случайной величины:');

writeln('1 - Равномерное распределение');

writeln('2 - Нормальное распределение');

writeln('3 - Экспоненциальное распределение');

writeln('4 - Бета распределение');

writeln('5 - Гамма-распределение');

readln(r);

Readln(mul);

for i := mul downto 1 do begin

write('нижняя граница: ', mul-i+1, ' = ');

Readln(bord[i].bottom); {i-я нижняя граница интегрирования}

write('верхняя граница ', mul-i+1, ' = ');

Readln(bord[i].top); {i-я верхняя граница интегрирования}

end;

write('Количество испытаний = ');

readln(n);

case r of

1: {вычисление интеграла по равномерному распределению}

begin

sum := 0;

sum2 := 0;

for i := 1 to n do

begin

pMult := 1;

for j := 1 to mul do begin

ksi[j]:= bord[j].bottom + (bord[j].top - bord[j].bottom) * random;

pMult := pMult * pMultipleUniformDist(j, bord);

end;

val := f(ksi, bord, mul) / pMult ;

sum := sum + val;

sum2 := sum2 + val * val;

end;

writeln('Результат: ', sum / n: 0: 8);

disp := (sum2 - sum * sum / n) / n;

writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

readln;

end;

2: {вычисление интеграла по нормальному распределению}

begin

sum := 0;

sum2 := 0;

for i := 1 to n do

begin

pMult := 1;

for j := 1 to mul do begin

ksi[j] := randomValueNormDist;

pMult := pMult * pNormDist(ksi[j]);

end;

val := f(ksi, bord, mul) / pMult;

sum := sum + val;

sum2 := sum2 + val * val;

end;

writeln('Результат: ', sum / n: 0: 8);

disp := (sum2 - sum * sum / n) / n;

writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

readln;

end;

3: {вычисление интеграла по экспоненциальному распределению}

begin

write('L = ');

readln(L);

sum := 0;

sum2 := 0;

for i := 1 to n do

begin

pMult := 1;

for j := 1 to mul do begin

alpha := random;

ksi[j] := randomValueExpDist(L);

pMult := pMult * pExpDist(ksi[j], L);

end;

val := f(ksi, bord, mul) / pMult;

sum := sum + val;

sum2 := sum2 + val * val;

end;

writeln('Результат: ', sum / n: 0: 8);

disp := (sum2 - sum * sum / n) / n;

writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

readln;

end;

4: {вычисление интеграла по бета-распределению}

begin

write('p = ');

readln(p);

write('m = ');

readln(m);

sum := 0;

sum2 := 0;

for i := 1 to n do

begin

pMult := 1;

for j := 1 to mul do begin

alpha := random;

ksi[j] := randomValueBetaDist(m, p);

pMult := pMult * pBetaDist(ksi[j], p, m);

end;

val := f(ksi, bord, mul) / pMult;

sum := sum + val;

sum2 := sum2 + val * val;

end;

writeln('Результат: ', sum / n: 0: 8);

disp := (sum2 - sum * sum / n) / n;

writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

readln;

end;

5: {вычисление интеграла по гамма-распределению}

begin

write('W= ');

readln(W);

write('B= ');

readln(B);

sum := 0;

sum2 := 0;

for i := 1 to n do

begin

pMult := 1;

for j := 1 to mul do begin

ksi[j] := randomValueGammaDist(W, B);

pMult := pMult * pGammaDist(ksi[j], W, B);

end;

val := f(ksi, bord, mul) / pMult;

sum := sum + val;

sum2 := sum2 + val * val;

end;

writeln('Результат: ', sum / n: 0: 8);

disp := (sum2 - sum * sum / n) / n;

writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

readln;

end;

Пример работы программы:

Вычислим интеграл

Точное значение интеграла: 0,25

Результаты расчётов программы:

	Приближённое значение интеграла
	Равномерное распределение		Нормальное распределение		Экспоненциальное распределение		Гамма-распределение		Бета-Распределение
Количество испытаний	Результат	Ошибка	Результат	Ошибка	Результат	Ошибка	Результат	Ошибка	Результат	Ошибка
n = 10000	0,24817	0,00661	0,24441	0,03325	0,24338	0,01992	0,24791	0,00655	0,25613	0,02060
n = 100000	0,25088	0,00209	0,25316	0,01058	0,25070	0,00638	0,25013	0,00208	0,25177	0,00641
n = 1000000	0,25029	0,00066	0,25021	0,00332	0,25041	0,00202	0,25025	0,00066	0,25126	0,00202

Приложение 2

Программа для решения интегрального уравнения Вольтерра методом Монте-Карло (Delphi)

var

n: longint; {объём выборки}

t: double; {верхняя граница интегрирования}

sum, sum2: double; {суммы вычисленных значений и квадратов этих значений}

disp: double; {дисперсия}

t0, t1, q: double;

begin

write('top border(t) = ');

readln(t);

sum := 0;

sum2 := 0;

for i := 1 to n do begin

sum := sum + y(t);

sum2 := sum2 + y(t)*y(t);

t0 := t;

t1 := Random;

q := 1;

while (t1 < t0) do begin

q := q * a(t1);

val := y(t1)*q;

sum := sum + val;

sum2 := sum2 + val * val;

t0 := t1;

t1 := Random;

end;

writeln('Результат: ', sum/n : 0: 8);

disp := (sum2 - sum * sum / n) / n;

writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

readln;

end;

Пример работы программы:

Вычислим интегральное уравнение Вольтерра при

Точное значение уравнения:

Объём выборки	Приближённое решение уравнения	Ошибка
n = 100000	0,35856	0,00227
n = 1000000	0,35995	0,00082
n = 10000000	0,36006	0,00008

Результаты расчётов программы:

Заключение

В процессе исследований было выяснено, что метод Монте-Карло имеет несколько важных преимуществ:

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

Также были замечены очевидные недостатки метода:

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

Литература

С. М. Ермаков, С. А. Михайлов «Статистическое моделирование», 1982г.
С. М. Ермаков «Метод Монте-Карло в вычислительной математике», 2009г.

Вычисление интегралов методом Монте-Карло 📙 Курсовая → 🆔 59112