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

clear, clc
% определяем кривые
f1 = @(x) x.^2 + 2;
f2 = @(x) -exp(abs(x)); 
% пределы интегрирования 
a = -3;
b = 2;
% определяем диапазон значений y
x = linspace(a,b,500);
y1 = f1(x);
y2 = f2(x);
maxy = max( [max(y1), max(y2)] ); % максимальное значение ф-ций
miny = min( [min(y1), min(y2)] ); % минимальное значение ф-ций
S = (b-a)*(maxy-miny); % площадь прямойгольной области
 
% начинаем Монте-Кралить
N = 5000; % число точек
x = a + (b-a)*rand(1,N);
y = miny + (maxy-miny)*rand(1,N);
L1 = f1(x)<=y; % логический вектор: если точка ниже перв. кривой, то 1
n1 = sum( L1 ); % кол-во точек под первой кривой
L2 = f2(x)<=y; % логический вектор: если точка ниже второй. кривой, то 1
n2 = sum( L2 ); % кол-во точек под второй кривой
m = abs(n1-n2); % кол-во точек между кривыми
I = m/N*S % значение интеграла
 
% рисуем облака! ^_^
x1 = x(L1);
y1 = y(L1);
x2 = x(L2);
y2 = y(L2);
plot(x1,y1,'.r',x2,y2,'.g','EraseMode','xor');

clear, clc
% определяем кривые
f1 = @(x) cos(x./2) + 5;
f2 = @(x) 0; 
% пределы интегрирования 
a = 0; b = 4;
% определяем диапазон значений y
x = linspace(a,b,500);
y1 = f1(x); y2 = f2(x);
maxy = max( [max(y1), max(y2)] ); % максимальное значение ф-ций
miny = min( [min(y1), min(y2)] ); % минимальное значение ф-ций
S = (b-a)*(maxy-miny); % площадь прямойгольной области
 
N = 2000; % число точек
x = a + (b-a)*rand(1,N);
y = miny + (maxy-miny)*rand(1,N);
L1 = f1(x)<=y; % логический вектор: если точка ниже перв. кривой, то 1
n1 = sum( L1 ); % кол-во точек под первой кривой
L2 = f2(x)<=y; % логический вектор: если точка ниже второй. кривой, то 1
n2 = sum( L2 ); % кол-во точек под второй кривой
m = abs(n1-n2); % кол-во точек между кривыми
I = m/N*S % значение интеграла
 
% Построение точек
x1 = x(L1); y1 = y(L1);
x2 = x(L2); y2 = y(L2);
plot(x1,y1,'.r',x2,y2,'.g','EraseMode','xor'); 
 
% Построение графика функции
x3 = 0:0.01:4;
y3 = cos(x./2)+5;
plot(x3,y3)

% Построение графика функции
hold on;
f = inline('cos(x./2)+5');
x3 = 0:0.01:4;
y3 = f(x3);
plot(x3,y3)

clear, clc
% определяем кривые
f1 = @(x) 2*(2*x-1);
f2 = @(x) -2*(2*(x-1)-1);
f3 = @(x) -2*(sin(2*pi*(x-2))+1);
% пределы интегрирования 
a = 0;
b = 3;
% определяем диапазон значений y
x = linspace(a,b,500);
y1 = f1(x);
y2 = f2(x);
y3 = f3(x);
maxy = max( [max(y1), max(y2), max(y3)] ); % максимальное значение ф-ций
miny = min( [min(y1), min(y2), min(y3)] ); % минимальное значение ф-ций
S = (b-a)*(maxy-miny); % площадь прямойгольной области
 
% Метод Монте-Карло
N = 5000; % число точек
x = a + (b-a)*rand(1,N);
y = miny + (maxy-miny)*rand(1,N);
L1 = f1(x)<= y; % логический вектор: если точка ниже перв. кривой, то 1
n1 = sum( L1 ); % кол-во точек под первой кривой
L2 = f2(x)<= y; % логический вектор: если точка ниже второй. кривой, то 1
n2 = sum( L2 );% кол-во точек под второй кривой
L3 = f3(x)<= y;
n3 = sum( L3 );
m = abs(n1-n2-n3); % кол-во точек между кривыми
I = m/N*S % значение интеграла
 
% заполняем точками
x1 = x(L1);
y1 = y(L1);
x2 = x(L2);
y2 = y(L2);
x3 = x(L3);
y3 = y(L3);
plot(x1,y1,'.r',x2,y2,'.g',x3,y3,'.r','EraseMode','xor');

clear, clc
% определяем функцию
f1 = @(x) 2*(2*x-1).*((0<=x)&(x<1)) + ...
    -2*(2*(x-1)-1).*((1<=x)&(x<2)) + ...
    -2*(sin(2*pi*(x-2))+1).*((2<=x)&(x<3));
% пределы интегрирования 
a = 0;
b = 3;
% определяем диапазон значений y
x0 = linspace(a,b,500);
y0 = f1(x0);
maxy =max(y0); % максимальное значение ф-ций
miny = min(y0); % минимальное значение ф-ций
S = (b-a)*(maxy-miny); % площадь прямойгольной области
 
% Метод Монте-Карло
N = 5000; % число точек
x = a + (b-a)*rand(1,N);
y = miny + (maxy-miny)*rand(1,N);
L1 = (y<=f1(x)); % логический вектор: если точка ниже перв. кривой, то 1
m = sum( L1 ); % кол-во точек под кривой
I = m/N*S; % значение интеграла
% заполняем точками
x1 = x(L1);
y1 = y(L1);
plot(x,y,'.y',x1,y1,'.r',x0,y0,'b','linewidth',2); 
 
% А теперь расчет по длинам
for i = 1:14
    N = 2^i;
    x = a + (b-a)*rand(1,N);
    y = miny + (maxy-miny)*rand(1,N);
    L1 = (y<=f1(x)); % логический вектор: если точка ниже перв. кривой, то 1
    m = sum( L1 ); % кол-во точек под кривой
    F(i) = m/N*S; % значение интеграла
end
F'