Using the Monte Carlo numerical method to calculate multidimensional integrals
Introduction
Back in the 1940s, John von Neumann and Stanislaw Ulam invented Monte Carlo simulation or Monte Carlo numerical method. They named it after the famous gambling place in Monaco, because the method has the same random characteristics as the game of roulette.
Monte Carlo methods are a broad class of computational algorithms that rely on repeated random sampling to produce numerical results. The basic concept is to use randomness to solve problems that could in principle be deterministic. Numerical Monte Carlo methods use three classes of problems, such as optimization, numerical integration, and generating results based on probability distributions.
The Monte Carlo method is used in real life, such as in physics problems, artificial intelligence, weather forecasting, and so on, and also has huge applications in finance, where the Monte Carlo numerical method is used to calculate stock prices, forecast sales, manage projects, and much more.[1]
The main advantage of using Monte Carlo is that it provides many possible outcomes and the probability of each from a large pool of random data samples, however, the method depends on assumptions and this can sometimes be challenging. Some other advantages of Monte Carlo are: it studies the behavior of a system without constructing it, provides generally accurate results compared to analytical models, helps to detect unexpected phenomena and system behavior, and perform “what if” analysis. [2]
In this paper, the Monte Carlo method will be used to approximate both one-dimensional and multidimensional integrals.
Description of the Method
As you know, to calculate the integral over a graph, we need to calculate the area under the curve. Let's consider the graph of a function that looks like 1/4 of a circle with a radius of one. Its area is π * radius2. Now, we will place this part of the circle in a square with a side equal to one, and, accordingly, an area equal to one. Using this information, we can calculate the area of the circle, which will be equal to π * radius2/4 = π/4 = 0.785.
Now let's pick a few points on our graph that are inside the square. Note: They don't have to be inside the circle, just inside the square. Let's say we used 18 different points, where 14 of them are inside the circle and the other 4 are just inside the square. Note: Since the circle takes up a significant portion of the square, it makes sense that the number of points inside the circle would be greater than the number of points outside it. Now, to calculate the area of the circle, we need to divide the number of points inside the circle (14) by the total number of points (18). Area of the circle = 14/18 = 0.77777778, which is approximately 0.785. This method of approximating the integral is called Monte Carlo, and its formula is[3], [4]:
Let's mention that the points we take are random, so if we use the same curve to calculate its area, we may get a slightly different result. And as we can summarize, more points will give us a more accurate result.
But what if we want to approximate a multidimensional integral? What if the form of our function is not so simple? Let's say we want to calculate the volume of a sphere. We can do this by calculating a triple integral:
We are going to integrate using the Monte Carlo method.
The first step is to select N – the number of random points to generate and r – the radius of the sphere. The next step is to generate points, values for riθiφi relatively limits of the integral, therefore 0≤ri≤r, 0≤θi≤π, 0≤фi≤2π.
For each random point ri, θifiwe will estimate the integrand f (ri,θifi) = ri2 * sin(θi) and after that we calculate the average value of the integrand by summing up all the values that we obtained in the previous step, divided by N.
Now we are ready for our last step where we calculate the volume of the sphere. We multiply the final result from the previous part by the three-dimensional volume, which can be calculated as Volume = (bx – ax) * (by – ay) * (bz – az) = (r – 0) * (π – 0) * (2π – 0). Approximation volume =
is the final approximation of the triple integral for the volume of a sphere using the Monte Carlo numerical method. [5]
Method Error
Monte Carlo integration for D space:
Method Error =
As we can see, the expression under the root of the method error is the standard deviation, that is, in general, the method error will depend on the expression
As N – the number of points – increases, the error of the method decreases.
Original Code
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) Application
We will calculate the triple integral:
First we define our function f(x, y, z) = x2 + y2.
From the description of the problem we see that the radius of the sphere is 2, since in spherical coordinates x2 + y2 + z2 = r2which means r = 2. So, we're going to modify our original code and generate random values for r between 0 and 2. Then, since we're only working with the half of the sphere that's where y ≥ 0, our value of φ will be 0 ≤ φ ≤ π. For θ, the value is still 0 ≤ θ ≤ π.
We then transform our x,y,z coordinates into spherical coordinates, where
x = r ∗ sin θ cos φ
y = r ∗ sin θ sin φ
z = r ∗ cos θ
dV = r2 ∗ sin θ dr dθ dφ
We will estimate the integrands f (riθiφi) = ((ri ∗ sin θi ∗ cos φi)2 + (ri ∗ sin θi ∗ sin φi)2) ∗ r2i ∗ sin(θi), and then we calculate the approximation of the integral using the numerical Monte Carlo method.
Approximation =
We will also calculate the integral manually and get the exact result, which is 128π/15 ≈ 26.8082573106.
In the code, we will output the approximation value as well as the absolute error.
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)Conclusion
The Monte Carlo method is an effective method for calculating multidimensional integrals. It is easy to use, and its error can be reduced by increasing the number of samples.
Like any method, Monte Carlo has its drawbacks. For example, when working with a small amount of information, the dispersion is too high, and in such cases, it is better to turn to other methods.
