Monte Carlo method¶
Shared by Antoine Rideau
On this page you will find how to calculate using Castor the value of \(\pi\) with Monte Carlo method.
Monte Carlo methods are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle. They are mainly used in physical and mathematical problems related to optimization, numerical integration, and generating draws from a probability distribution.
First, draw 1 by 1 square and inscribe a quadrant within it.
Then n points are scattered randomly with \((x,y)\) coordinate over the square.
// Parameters
double n = 5000.; //Number of points
auto MX = rand(1, n);
auto MY = rand(1, n);
As the quadrant’s surface is equal to
and the square’s surface is
the probability to be in the quadrant is
So, points inside the quadrant i.e points with \(x^2 + y^2 \leq 1\) are counted.
// Pi computation
matrix<std::size_t> Incircle;
Incircle = find(pow(MX, 2) + pow(MY, 2) <= 1);
double Pi = 4. * (numel(Incircle) / n);
Another way of doing it is by using sum
size_t Nin = sum<size_t>(pow(MX, 2) + pow(MY, 2) <= 1);
double Pi = 4. * (Nin / n);
See sum
Code¶
#include "castor/matrix.hpp"
#include "castor/graphics.hpp"
using namespace castor;
int main(int argc, char const *argv[])
{
// Parameters
double n = 5000.; //Number of points
auto MX = rand(1, n);
auto MY = rand(1, n);
// Pi computation
matrix<std::size_t> Incircle;
Incircle = find(pow(MX, 2) + pow(MY, 2) <= 1);
double Pi = 4. * (numel(Incircle) / n);
std::cout << "Calculated value of pi: " << Pi << endl;
// Visu
auto X = linspace(0, 1, 1000);
auto Y = sqrt(1 - pow(X, 2));
figure fig;
plot(fig, MX, MY, {"b"});
plot(fig, eval(MX(Incircle)), eval(MY(Incircle)), {"r"});
plot(fig, X, Y, {"r-"});
drawnow(fig);
return 0;
}
With this code you should get these outputs :
Calculated value of pi: 3.148
