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

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

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

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

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

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

Выполнил:

Студент 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

Приложения

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

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

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

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

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

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

 

 

Введение

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

 

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

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

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

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

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

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

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

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

Распределение	Формула плотности	Моделирующая формула
Равномерное непрерывное распределение на
Стандартное нормальное распределение
Экспоненциальное (показательное) распределение
Гамма- распределение
Бета-распределение

 

 

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

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

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

 

 

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

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

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

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

 

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

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

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

 

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

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

 

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

Интегральное уравнение Вольтерра  имеет вид  , где и - известные функции, . Для нахождения приближённого значения функции используется схема Неймана-Улама. Также, как и при вычислении интегралов методом Монте-Карло, проводятся 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; {случайный вектор}