[Pascal] Liczba Pi metodą Monte Carlo

[WhiteB]

Czy mógłby mi ktoś napisać prosty programik szukający przybliżonej wartości liczby Pi metodą Monte Carlo w Pascalu? Jestem kompletna noga w tym i proszę Was o pomoc

[szw1710]

Koło wpisane w kwadrat jednostkowy ma pole \(\displaystyle{ \pi}\). Należy więc w kwadracie \(\displaystyle{ 0\le x\le 1,\;0\le y\le 1}\) wylosować dużo punktów \(\displaystyle{ (x,y),}\) a następnie sprawdzić, ile z nich leży w kole, tzn. ile punktów \(\displaystyle{ (x,y)}\) spełnia nierówność

\(\displaystyle{ \left(x-\frac{1}{2}\right)^2+\left(y-\frac{1}{2}\right)^2\le 1}\)

Środkiem tego koła jest punkt \(\displaystyle{ \left(\frac{1}{2},\frac{1}{2}\right),}\) stąd ta nierówność - odległość punktu od środka nie przekracza promienia.

Stosunek liczby punktów leżących w kole do liczby wylosowanych punktów daje przybliżenie liczby \(\displaystyle{ \pi.}\)

Program robiłbym tak:

1. Określenie liczby punktów do wylosowania.
2. Losowanie punktu \(\displaystyle{ (x,y)}\)
3. Sprawdzenie powyższej nierówności; jeśli zachodzi, zwiększenie licznika punktów leżących w kole o jeden.
4. Powrót do punktu 2, pętla aż do liczby wylosowanych punktów.

W pseudokodzie:

Kod: Zaznacz cały

licznik:=0;
podaj liczbę n punktów do wylosowania;
dla k:=1 do n powtarzaj
{
  losuj punkt (x,y)
  jeśli zachodzi nierówność (x-0.5)^2+(y-0.5)^2<=1, to licznik := licznik+1
};
pi:=licznik/n

[WhiteB]

Dziękuje Ci wielce. Mam nadzieję, że dalej sobie poradzę już

[WhiteB]

Kod: Zaznacz cały

PROGRAM MONTE;
USES GRAPH;
VAR

ILPROB,ILSUKCESOW,X,Y,KARTA,TRYB:LONGINT;
POLE_C:REAL;
BEGIN
KARTA:=DETECT;INITGRAPH(KARTA,TRYB, 'kkk');
RANDOMIZE;
FILLELLIPSE(320,240,200,200);
FOR ILPROB:=1 TO 200000 DO
BEGIN
X:=RANDOM(640);
Y:=RANDOM(480);
IF (GETPIXEL(X,Y)=15)OR(GETPIXEL(X,Y)=4)THEN
BEGIN
PUTPIXEL(X,Y,4);
ILSUKCESOW:=ILSUKCESOW+1;
END
ELSE PUTPIXEL(X,Y,YELLOW);
END;
POLE_C:=640*480/200.5/200.5;
WRITE('POLE=',POLE_C/ILPROB*ILSUKCESOW:1:2);
READLN;
CLOSEGRAPH;
END.

Czy mógłby ktoś wskazać mi błąd i jak poprawić?

[szw1710]

Tam pewien błąd zrobiłem. Moja metoda nie liczy przybliżenia liczby \(\displaystyle{ \pi,}\) ale \(\displaystyle{ \frac{\pi}{4}.}\) Do \(\displaystyle{ \pi}\) potrzeba koła o promieniu \(\displaystyle{ 2}\), np. wpisanego w kwadrat \(\displaystyle{ -1\le x\le 1,\;-1\le y\le 1.}\) Stosowna nierówność przybierze postać \(\displaystyle{ x^2+y^2\le 1.}\)

[WhiteB]

Kod: Zaznacz cały

program Losowanie;


var
licznik, n, k: integer;
iks, ygrek : integer;
pi, q : double;

begin

licznik:=0;

WriteLN('Witaj w losowaniu! Podaj liczbe punktow do wylosowania');
ReadLN(n);

for k:=2 to n do
begin
iks:=random(2);
ygrek:=random(2);

	if ((iks)*(iks))+((ygrek)*(ygrek)) <= 1 then
    begin
	licznik:=licznik+1;
    WriteLn('x: ', iks , ' - y: ', ygrek, ' - licznik: ', licznik);
    end
end;

pi:=(4*licznik)/n;

WriteLN('Przyblizona wartosc liczby pi to: ', pi:0:4);


Read(q);
 end.

może być coś takiego?

[szw1710]

Uruchom program i zobacz jakie generuje wyniki.

[WhiteB]

4,000

[szw1710]

No właśnie. Więc coś musisz poprawić. Losujesz punkty w kwadracie \(\displaystyle{ [0,2]\times [0,2].}\) A zatem koło ma środek \(\displaystyle{ (1,1)}\) i promień \(\displaystyle{ 1,}\) a jego pole to \(\displaystyle{ \pi.}\) Więc sprawdzamy warunek

\(\displaystyle{ (x-1)^2+(y-1)^2\le 1.}\)

Wartość licznik/n podać powinna przybliżenie \(\displaystyle{ \pi.}\)

Nie widzę w Twoim kodzie zainicjowania generatora liczb pseudolosowych.

[zeus_156]

Kod: Zaznacz cały

...
var
licznik, n, k: integer;
iks, ygrek : integer;
pi, q : double;
...

uważam, że wszystkie zmienne powinny być typu double, inaczej zawsze będziesz sprawdzał tylko 3 wartości: 0 i 1 oraz 2 bo integer nie może mieć wartości pośrednich.


Na pascal-u się nie znam ale napisałem program w C#:

Kod: Zaznacz cały

using System;
using System.Text;

namespace CalcPiByMonteCarlo
{
 class Program
 {
  static void Main(string[] args)
  {
   int licznik = 0;
   int raportowanie = 10;
   double x, y, tmp;
   var rnd = new Random();
   for (int i = 1; i <= 10000000; i++)
   {
    x = 2 * rnd.NextDouble();  /* NextDouble() zwraca wartosc z przedzialu <0,0;1,0) */
    y = 2 * rnd.NextDouble();
    if ((x - 1) * (x - 1) + (y - 1) * (y - 1) <= 1.0)
    {
     licznik += 1;
    }
    if (i == raportowanie)
    {
     raportowanie *= 10;
     tmp = 4 * licznik;
     Console.WriteLine("i={0} pi={1}", i, tmp / i);
    }
   }
   Console.Read();
  }
 }
}

i otrzymałem takie wyniki:
i=10 pi=2.8
i=100 pi=3.32
i=1000 pi=3.168
i=10000 pi=3.1812
i=100000 pi=3.14312
i=1000000 pi=3.142008
i=10000000 pi=3.1415628

[szw1710]

Witam archeologa Temat ma półtora roku.

Ładne wyniki otrzymałeś. Na tym metoda Monte Carlo polega. A może teraz spróbuj napisać programik wyliczający przybliżoną wartość całki \(\displaystyle{ \int_0^1 e^{x^2}\dd x}\). Też metodą Monte Carlo. Skorzystaj z tego, że całka oznaczona funkcji nieujemnej jest polem obszaru, jaki ogranicza wykres tej funkcji i oś pozioma (potocznie nazywana \(\displaystyle{ x}\)).

[zeus_156]

\(\displaystyle{ \int_0^1 e^{x^2}\dd x}\)

Po lekkiej modyfikacji programu udało mi się obliczyć wartość całki. Kod programu:

Kod: Zaznacz cały

using System;
using System.Text;

namespace MonteCarlo
{
 class Program
 {
  static void Main(string[] args)
  {
   int licznik = 0;
   int raportowanie = 10;
   double x, y, tmp;
   var rnd = new Random();
   for (int i = 1; i <= 100000000; i++)
   {
    x = rnd.NextDouble();  /* NextDouble() zwraca wartosc z przedzialu <0,0;1,0) */
    y = 3 * rnd.NextDouble();
    if (Math.Exp(x*x) >= y)
    {
     licznik += 1;
    }
    if (i == raportowanie)
    {
     raportowanie *= 10;
     tmp = 3 * licznik;  /* Liczba 3 jest tu dlatego poniewaz obszar z jakiego losuje to 3 (1*3) */
     Console.WriteLine("i={0} calka={1}", i, tmp / i);
    }
   }
   Console.WriteLine("Done");
   Console.Read();
  }
 }
}

A wyniki tego programu to:

i=10 calka=1.5
i=100 calka=1.38
i=1000 calka=1.425
i=10000 calka=1.4328
i=100000 calka=1.46313
i=1000000 calka=1.464888
i=10000000 calka=1.4631342
i=100000000 calka=1.46265048
Done

[szw1710]

Świetnie. Myślę, że rozumiesz, na czym rzecz polega