Введение

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

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

Метод Монте-Карло используется в реальной жизни, например, в задачах, связанных с физикой, создании искусственного интеллекта, прогнозировании погоды и так далее, а также имеет огромное применение в финансах, где числовой метод Монте-Карло используется для расчёта стоимости акций, прогнозировании продаж, управления проектами и многого другого.^[1]

Основное преимущество использования Монте-Карло заключается в том, что этот метод обеспечивает множество возможных результатов и вероятность каждого из большого пула случайных выборок данных, однако, метод зависит от предположений, и это иногда может быть сложной задачей. Некоторые другие преимущества Монте-Карло: он изучает поведение системы без её построения, обеспечивает в целом точные результаты, по сравнению с аналитическими моделями, помогает обнаружить неожиданное явление и поведение системы, а также выполнить анализ "что, если". ^[2]

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

Описание Метода

Как известно, чтобы вычислить. интеграл по графику, нам нужно вычислить площадь под кривой. Давайте рассмотрим график функции, который выглядит как 1/4 круга с радиусом, равным единице. Его площадь равна π * радиус². Теперь, эту часть круга, мы поместим в квадрат со стороной, равной единице, и, соответственно, площадью равной единице. Используя эту информацию, мы можем вычислить площадь круга, которая будет равна π * радиус²/4 = π/4 = 0,785.

Теперь давайте выберем несколько точек на нашем графике, которые находятся внутри квадрата. Примечание: они не должны быть внутри круга, просто, внутри квадрата. Допустим, мы использовали 18 разных точек, где 14 из них находятся внутри круга, а другие 4 - только внутри квадрата. Примечание: поскольку круг занимает значительную часть квадрата, логично, что количество точек в круге будет больше, чем количество точек вне его. Теперь, чтобы вычислить площадь круг, нам нужно разделить количество точек внутри круга (14) на общее количество точек (18). Площадь круга = 14/18 = 0,77777778, что примерно равно 0,785. Этот метод аппроксимации интеграла называется Монте-Карло, и его формула^{[3], [4]}:

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

Но что, если нам нужно аппроксимировать многомерный интеграл? Что, если форма нашей функции не такая простая? Допустим, нам нужно вычислить объём сферы. Мы можем сделать это, вычислив тройной интеграл:

Мы собираемся интегрировать, используя метод Монте-Карло.

Первым шагом будет выбор N - количества случайных точек для генерации и r - радиуса сферы. Следующим шагом будет генерация точек, значений для r_i, θ_i, φ_i относительно пределов интеграла, поэтому 0≤r_i≤r, 0≤θ_i≤π, 0≤ф_i≤2π.

Для каждой случайной точки r_i, θ_i, ф_i, мы оценим подынтегральное выражение f (r_i,θ_i, ф_i) = r_i² * sin(θ_i) и после этого вычислим среднее значение подынтегрального выражения, суммируя все значения, которые мы получили на предыдущем шаге, деленные на N.

Теперь мы готовы к нашему последнему шагу, где мы вычисляем объём сферы. Умножим окончательный результат из предыдущей части на трехмерный объём, который можно вычислить как Объём = (b_x - a_x) * (b_y - a_y) * (b_z - a_z) = (r - 0) * (π - 0) * (2π - 0). Объём приближения =

является окончательным приближением тройного интеграла для объёма сферы с использованием численного метода Монте-Карло. ^[5]

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

Интеграция Монте-Карло для пространства D:

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

Как мы можем заметить, выражение под корнем погрешности метода - это стандартное отклонение, то есть в целом, ошибка метода будет зависеть от выражения

Когда N - количество точек - увеличивается, погрешность метода уменьшается.

Оригинальный Код

import numpy as np

N=1000 #Количество точек для случайной выборки
approximations = np.zeros(N)

def f(x, y, z): #определить функцию для интегрирования (здесь функция в интеграле для расчета объёма сферы)
    return 1
  
def random_points(a_x, b_x, a_y, b_y, a_z, b_z): #Определить функцию для генерации случайных точек
    r=np.random.uniform(a_x, b_x) #Сгенерировать случайные значения радиуса
    theta=np.random.uniform(a_y, b_y) #Сгенерировать случайные значения угла тэта
    phi=np.random.uniform(a_z, b_z) #Сгенерировать случайные значения угла фи
    return r, theta, phi 
  
def spherical_coordinates(r, theta, phi): #Определить функцию для конвертации в сферические координаты
    x=r*np.sin(theta)*np.cos(phi)
    y=r*np.sin(theta)*np.sin(phi) 
    z=r*np.cos(theta)
    return x, y, z  
  
for i in range (N):
    r, theta, phi = random_points(0, 5, 0, np.pi, 0, 2*np.pi)
    x, y, z = spherical_coordinates(r, theta, phi)
    dV=r**2*np.sin(theta) #Оставшийся компонент от dV=r**2*sin(theta) dr dtheta dphi"
    approximations[i]=f(x, y, z)*dV #Вычислить подынтегральную функцию
    
Monte_Carlo_Approximation=(np.sum(approximations)/N)*(5-0)*(np.pi-0)*(2*np.pi-0) 
print ("Объём сферы подсчитанный с помощью метода Монте-Карло =", Monte_Carlo_Approximation)  

Применение

Мы будем вычислять тройной интеграл:

Сначала мы определим нашу функцию f (x, y, z) = x² + y².

Из описания задачи мы видим, что радиус сферы равен 2, так как в сферических координатах x² + y² + z² = r², значит, r = 2. Поэтому, мы собираемся изменить наш оригинальный код и сгенерировать случайные значения для r от 0 до 2. Затем, поскольку мы работаем только с половиной сферы, которая находится там, где y ≥ 0, наше значение φ будет 0 ≤ φ ≤ π. Для θ значение по-прежнему 0 ≤ θ ≤ π.

Затем мы преобразуем наши координаты x, y, z в сферические координаты, где
x = r ∗ sin θ cos φ
y = r ∗ sin θ sin φ
z = r ∗ cos θ
dV = r² ∗ sin θ dr dθ dφ

Мы оценим подынтегральные функции f (r_i, θ_i, φ_i) = ((r_i ∗ sin θ_i ∗ cos φ_i)² + (r_i ∗ sin θ_i ∗ sin φ_i)²) ∗ r²_i ∗ sin(θ_i), а затем вычислим приближение интеграла с помощью численного метода Монте-Карло.
Приближение =

Мы также вычислим интеграл вручную и получим точный результат, который равен 128π/15 ≈ 26,8082573106.

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

import numpy as np

N=1000

approximations=np.zeros(N)
Monte_Carlo_Approximation=[]

def f(x, y, z): 
    return x**2+y**2
for i in range (N):
    r=np.random.uniform(0, 2) 
    theta=np.random.uniform(0, np.pi) 
    phi=np.random.uniform(0, np.pi)
    x=r*np.sin(theta)*np.cos(phi) 
    y=r*np.sin(theta)*np.sin(phi) 
    z=r*np.cos(theta)
    dV=r**2*np.sin(theta) 
    approximations[i]=f(x, y, z)*dV
    Monte_Carlo_Approximation.append(np.sum(approximations)/(i+1)*(2-0)*(np.pi-0)*(np.pi-0)) 
    
absolute_error=abs(Monte_Carlo_Approximation[N-1]-(128*np.pi)/15)

print("Интеграл, вычисленный с помощью метода Монте-Карло =", Monte_Carlo_Approximation[N-1])
print("Абсолютная погрешность =", absolute_error)

Заключение

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

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