Метод Монте Карло
.pdfЛабораторная работа №5 Метод Монте-Карло
Найдем методом Монте-Карло площадь одного сектора круга:
Впишем сектор круга в квадрат :
Сформируем |
набор случайных точек |
внутри квадрата, то есть |
и |
. Для этого нам необходимо получить определенное |
количество, обозначить через NN, случайных чисел. В maple это можно сделать с помощью команды rand(K..M). Эта команда формирует целые случайные числа от K до M. Нам же нужны дробные числа от нуля до единицы. Мы получим целые числа от 0 до 1000, а затем делением на 1000
получим необходимые ( |
) нам числа. Те |
точки, |
для которых |
||
справедливо: |
, очевидно, лежат внутри |
сектора. |
Подсчитаем |
||
количество таких точек и с помощью формулы |
|
|
найдем площадь |
||
|
|
||||
сектора. |
|
|
|
|
|
> restart; |
|
|
|
|
|
> NN:=1000; |
|
|
|
|
|
|
NN := 1000 |
|
|
|
|
Формируем массивы X и Y случайных чисел 0≤X[i] |
|
и |
. |
>r:=rand(0.. NN);
>X:=[seq(evalf(r()/NN),i=1..NN)]:
>Y:=[seq(evalf(r()/NN),i=1..NN)]:
Включим графические пакеты:
>with(plots);
>with(plottools);
N1 и N2 переменные для подсчета точек попавших внутрь и вне сектора. Напомним, что по условиям формирования случайных точек все точки попадают внутрь квадрата.
1
Введем переменные для подсчета точек внутри сектора и вне его. Очевидно, что N1+N2=NN. Первоначально N1 и N2 равны нулю.
> N1:=0;
N1 := 0
> N2:=0;
N2 := 0
Зададим графические объекты дугу (S) и линии (l1 и l2), ограничивающие квадрат :
>S:=arc([0,0],1,0..Pi/2,thickness=3,color=green):
>l1:=line([0,1],[1,1],thickness=2,color=red):
>l2:=line([1,0],[1,1],thickness=2,color=red):
Зададим цикл, формирующий графические объекты SS[i]- диски радиусом 0.005 и центром в точках X[i], Y[i] и имеющие черный цвет, если их центр находится внутри сектора и черный вне его. С помощью команд N1:=N1+1; и N2:=N2+1;будем считать точки внутри и вне сектора.
>for i from 1 to NN do if X[i]^2+Y[i]^2<1 then
SS[i]:=disk([X[i],Y[i]],0.005,color=black); N1:=N1+1;
>else SS[i]:=disk([X[i],Y[i]],0.005,color=red); N2:=N2+1;
>fi;
>od:
Изобразим точки-диски:
> display(S,l1,l2,seq(SS[i],i=1..NN));
2
Возьмем полученные результаты:
> N1:=N1;
N1 := 792
> N2:=N2;
N2 := 208
Возьмем найденную нашим методом площадь:
> Sq:=evalf(N1/1000);
Sq := 0.7920000000
Сравним с известным результатом:
> evalf(Pi/4);
0.7853981635
Применение метода Монте-Карло для решения задачи Дирихле
Задачей Дирихле называется задача:
, (1)
с краевыми условиями:
. (2)
То есть необходимо найти значения функции u(x,y),удовлетворяющей уравнению (1) внутри некоторой области D, значения которой на границе этой области задаются (2).
Пример: Проверим как методом Монте-Карло решается задача Дирихле, если
и известны значения функции на круге радиуса 1, с центром в начале координат.
3
>restart;
>f:=(x,y)->x^2-y^2;
f:= ( x, y ) x2 y2
>diff(diff(f(x,y),x),x)+diff(diff(f(x,y),y),y);
0
> plot3d(f(x,y),x=-10..10,y=-10..10,axes=normal);
> r:=rand(1..4); r := pro c()
lo calt;
glo bal_seed;
_seed := irem( a _seed, p );
t := _seed;
to concats do _seed := irem( a _seed, p ); t := s t _seed end do; irem( t, divisor) offset
end proc
> h:=0.01;
h := 0.01
> N:=500;
4
N := 500
> S:=array(1..N);
S := array( 1 .. 500, [ ] )
>for i from 1 to N do
>X[1]:=0; Y[1]:=0.5; k:=1;
>while X[k]^2+Y[k]^2<1 do
>if r()=1 then X[k+1]:=X[k]+h; Y[k+1]:=Y[k]; k:=k+1;
>else
>if r()=2 then X[k+1]:=X[k]; Y[k+1]:=Y[k]+h; k:=k+1;
>else
>if r()=3 then X[k+1]:=X[k]-h; Y[k+1]:=Y[k]; k:=k+1;
>else X[k+1]:=X[k]; Y[k+1]:=Y[k]-h; k:=k+1;
>fi; fi; fi;
>od;
>S[i]:=f(X[k],Y[k]);
>od:
>SS:=add(S[n],n=1..N)/N;
SS:= -0.2668754000
>f(X[1],Y[1]);
-0.25
>
> SS := -.2588660000;
|
|
SS := -0.2588660000 |
|
|
|
|
|
|
|
Задания |
|
|
|
|
|
1. |
Методом Монте-Карло найти площадь эллипса: |
|
|
|
. |
||
|
|
||||||
2. |
Методом Монте-Карло найти площадь фигуры ограниченной |
||||||
|
параболами |
и y 2 x2 . |
|
|
|
|
|
3. |
Получить у преподавателя функцию |
и область и границу ее, на |
|||||
|
которой известны значения функции, и найти значение (проверить) в |
||||||
|
точке внутри области методом Монте-Карло. |
|
|
|
|
5