Całkowanie numeryczne: Metoda Monte Carlo

Aby zrozumieć zasadę metody Monte Carlo wyobraźcie sobie, iż chcecie wyznaczyć pole koła wpisanego w kwadrat o boku równym 2 (pole to co do wartości jest równe liczbie pi, ale na razie udajmy, że o tym nie wiemy). W tym celu wyznaczamy wewnątrz kwadratu dużo losowych punktów. Następnie zliczamy te punkty, które wpadają do wnętrza koła. Pole koła jest w przybliżeniu równe:

	- pole koła
	- pole kwadratu
	- liczba punktów w kole
	- liczba wszystkich punktów

Przykład:

Oto odpowiedni skrypt JavaScript symulujący obliczanie powierzchni koła dla podanego przykładu:

JavaScript

<html>
<head>
  <title>Wyznaczanie liczby PI
  metodą Monte Carlo</title>
</head>
<body>

<div style="overflow-x: auto;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frm1"
        style="text-align:
               center;
               background-color:
               #E7E7DA">
          <b>
          Obliczanie pola
          koła wpisanego
          <br/>
          w kwadrat o boku 2
          <br/>&nbsp;&nbsp;
          za pomocą metody
          Monte Carlo
          &nbsp;&nbsp;</b><br/>
          (C)2026 mgr
          Jerzy Wałaszek<hr>
          Podaj liczbę
          punktów<br/>
          do wygenerowania:<br/>
          <input
          type="text"
          name="inp_n"
          size="16"
          value="10000"
          style="text-align:
                 right">
          <hr>
          <input
          onclick="js_p();"
          value="Oblicz pole koła"
          name="B1"
          type="button">
          <hr>
          Pole koła wynosi:
          <div id="out1">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script language="javascript">

// Przykładowa aplikacja
// obliczająca pole koła
// wpisanego w kwadrat
// za pomocą metody
// Monte Carlo
//----------------------
// (C)2004 mgr Jerzy Wałaszek

function js_p()
{
  var n = parseInt(document
          .frm1.inp_n.value);
  var nk,s,x,y,i;

  if(isNaN(n) || (n < 2))
    s = "<font color=red><b>" +
        "Popraw dane" +
        "</b></font>";
  else
  {
    nk = 0
    for(i = 0; i < n; i++)
    {
      x = Math.random() * 2;
      y = Math.random() * 2;
      if(Math.sqrt((x-1) *
        (x-1) + (y-1) *
        (y-1)) <= 1) nk++;
    }
    s = 4 * nk / n;
    s = Math.round(s * 100000) /
        100000;
  }
  document.getElementById("out1")
  .innerHTML = s
}

</script>

</body>
</html>

Oczywiście wynik jest bliski liczbie π = 3,1415926535... dopiero dla dużych wartości n. Zaczynają wtedy obowiązywać prawa dużych liczb i pomimo przypadkowości wyboru punktów, pojawia się prawidłowość. Ponieważ punkty rozkładają się równomiernie w obrębie pola kwadratu, to stosunek liczby punktów wewnątrz koła do liczby wszystkich punktów w kwadracie jest równy stosunkowi pola koła do pola kwadratu. Stąd właśnie pochodzi nasz wzór:

Wzór ten jest podstawą wyznaczania wartości całki oznaczonej za pomocą metody Monte Carlo, czyli losowania punktów. Zasada jest następująca:

Dla danej funkcji f(x), której całkę oznaczoną chcemy obliczyć w przedziale całkowania <x_p, x_k>, wyznaczamy prostokąt obejmujący pole pod wykresem tej funkcji o wysokości h i długości podstawy (x_k - x_p). W dalszej kolejności losujemy n punktów i zliczamy te punkty n_w, które wpadają w pole pod wykresem funkcji. Wartość całki wyraża się wzorem przybliżonym:

Otrzymany wzór ma kilka wad. Na przykład w ogólnym przypadku trudno wyznaczyć wysokość h. Również kłopoty pojawiają się, gdy funkcja zmienia znak w przedziale całkowania. Dlatego częściej jako metodę Monte Carlo przyjmuje się metodę, która wyznacza średnią z wartości funkcji w przedziale całkowania na podstawie serii losowo wybranych współrzędnych x. Następnie średnia ta jest mnożona przez długość przedziału całkowania i otrzymujemy przybliżoną wartość całki oznaczonej. Wzór ma następującą postać:

x_los jest wartością pseudolosową zmiennoprzecinkową z przedziału <x_p, x_k> Wartość tę otrzymujemy w przykładowych językach programowania wg poniższych wzorów:

Język	Instrukcja
Pascal	xlos := xp + random * (xk - xp);
C++	xlos = xp + (double)rand()/(double)(RAND_MAX+1) * (xk - xp);
Basic	xlos = xp + rnd(1) * (xk - xp)
JavaScript	xlos = xp + Math.random() * (xk - xp)
Python (dodatek)	import random ... xlos = random.uniform(xp, xk)

do podrozdziału do strony

Specyfikacja problemu

Dane wejściowe

x_p	- początek przedziału całkowania, x_p ∈ R
x_k	- koniec przedziału całkowania, x_k ∈ R
n	- liczba punktów podziałowych, n ∈ N
f(x)	- funkcja rzeczywista, której całkę liczymy

Dane wyjściowe

s	- przybliżona wartość całki oznaczonej funkcji f(x) w przedziale <x_p,x_k>, s ∈ R

Zmienne pomocnicze

dx	- odległość między dwoma sąsiednimi punktami podziałowymi, dx ∈ R
i	- licznik punktów podziałowych, i ∈ N
x_los	- punkt wybrany losowo z przedziału całkowania, x_los ∈ R

Lista kroków

K01:	ii
K02:	Dla i = 1,2,...,n : wykonuj kroki K03...K04
K03:	x_los ← liczba losowa z <x_p,x_k>
K04:	s ← s + f(x_los)
K05:
K08:	Zakończ

Schemat blokowy

Obliczenia rozpoczynamy od pobrania informacji o przedziale całkowania oraz o ilości punktów losowych, które należy wygenerować w celu obliczenia wartości średniej funkcji w tym przedziale. Dokładność metody rośnie wraz ze wzrostem n.

W zmiennej s będziemy obliczać sumy wartości funkcji. Zmienna ta posłuży później do wyliczenia średniej oraz samej całki oznaczonej. Na początku obliczeń ustawiamy ją na 0. W zmiennej dx zapamiętujemy szerokość przedziału całkowania. Wartość ta jest wykorzystywana zwykle przy generacji liczby losowej oraz na końcu przy obliczeniu wartości całki.

Rozpoczynamy pętlę iteracyjną kontrolowaną przez zmienną i. Pętla ta wykona się n-razy. Wewnątrz pętli generujemy liczbę pseudolosową x_los leżącą w przedziale <x_p.x_k>. Metoda generacji zależy od wybranego języka programowania, który udostępnia odpowiednie funkcje pseudolosowe. Właściwe procedury generacji tej liczby podaliśmy wyżej w tabelce. Po wyznaczeniu liczby pseudolosowej x_los obliczamy wartość funkcji w tym punkcie i dodajemy ją do sumy s. Gdy pętla się zakończy, wyliczamy średnią wartość funkcji w przedziale całkowania, mnożymy ją przez długość tego przedziału i otrzymujemy przybliżoną wartość całki oznaczonej. Wypisujemy wyniki i kończymy algorytm.

do podrozdziału do strony

W przedziale <0, 1> całka funkcji f(x) = x² + 2x ma dokładną wartość 1.333... W naszym programie liczymy n = 10000 losowych punktów.

C++

// Obliczanie całki oznaczonej
// metodą Monte Carlo
// ---------------------------
// (C)2004 mgr Jerzy Wałaszek

#include <iomanip>
#include <iostream>
#include <cstdlib>
#include <time.h>

using namespace std;

// Tutaj definiujemy funkcję
double f(double x)
{
  return(x * x + 2 * x);
}

// Program główny

int main()
{
  // liczba punktów losowych
  const int N = 10000;
  double xp,xk,s,dx;
  int i;

  // 3 cyfry po przecinku
  cout << setprecision(3)
  // format stałoprzecinkowy
       << fixed;

  cout << "Obliczanie calki oznaczonej\n"
          "    Metoda Monte Carlo\n"
          "---------------------------\n"
          "(C)2004 mgr J.Walaszek  I LO\n\n"
          "f(x) = x * x + 2 * x\n\n"
          "Poczatek przedzialu calkowania\n\n"
          "xp = ";
  cin >> xp;
  cout << "\nKoniec przedzialu calkowania\n\n"
          "xk = ";
  cin >> xk;
  cout << endl;
  srand(time(NULL));
  s  = 0;
  dx = xk - xp;
  for(i = 1; i <= N; i++)
    s += f(xp+((double)rand()/(double)(RAND_MAX+1)*dx));
  s = dx * s / N;
  cout << "Wartosc calki wynosi : " << s
       << endl << endl;
  system("pause");
  return 0;
}

Pascal

// Obliczanie całki oznaczonej
// metodą Monte Carlo
// ---------------------------
// (C)2004 mgr Jerzy Wałaszek

program int_montecarlo;

// Tutaj definiujemy funkcję
function f(x : double) : double;
begin
  f := x * x + 2 * x;
end;

// Program główny

// Liczba punktów losowych
const N = 10000;
var
  xp,xk,s,dx : double;
  i : integer;

begin
  writeln('Obliczanie calki oznaczonej');
  writeln('    Metoda Monte Carlo');
  writeln('---------------------------');
  writeln('(C)2004 mgr J.Walaszek I LO');
  writeln;
  writeln('f(x) = x * x + 2 * x');
  writeln;
  writeln('Poczatek przedzialu calkowania');
  writeln;
  write('xp = '); readln(xp);
  writeln;
  writeln('Koniec przedzialu calkowania');
  writeln;
  write('xk = '); readln(xk);
  writeln;

  // Inicjujemy generator liczb
  // pseudolosowych
  randomize;

  s  := 0;
  dx := xk - xp;
  for i := 1 to N do
    s := s + f(xp + random * dx);
  s := dx * s / N;
  writeln('Wartosc calki wynosi : ',s:0:3);
  writeln;
  writeln('Nacisnij Enter...');
  readln;
end.

Basic

' Obliczanie całki oznaczonej
' metodą Monte Carlo
' ---------------------------
' (C)2004 mgr Jerzy Wałaszek

Declare Function f(x As Double) _
                     As Double

' Program główny

' Liczba punktów losowych
Const N = 10000

Dim As Double xp,xk,s,dx
Dim As Integer i

Print "Obliczanie  calki oznaczonej"
Print "za pomoca metody Monte Carlo"
Print "----------------------------"
Print "(C)2004 mgr J.Walaszek  I LO"
Print
Print "f(x) = x * x + 2 * x"
Print
Print "Poczatek przedzialu calkowania"
Print
input "xp = ", xp
Print
Print "Koniec przedzialu calkowania"
Print
Input "xk = ", xk
Print

Randomize

s = 0
dx = xk - xp
For i = 1 To N
  s += f(xp + Rnd() * dx)
Next
s *= dx / N
Print Using _
      "Wartosc calki wynosi : #####.###";s
Print
Print "Nacisnij klawisz Enter..."
Sleep
End

' Tutaj definiujemy funkcję **
function f(x As Double) _
             As Double
  f = x * x + 2 * x
End Function

Python (dodatek)

# Obliczanie całki oznaczonej
# metodą Monte Carlo
# ---------------------------
# (C)2004 mgr Jerzy Wałaszek

import random

def f(x):
    return x * x + 2 * x

# Program główny

# Liczba punktów losowych
n = 10000

print("Obliczanie  całki oznaczonej")
print("za pomocą metody Monte Carlo")
print("----------------------------")
print("(C)2004 mgr J.Wałaszek  I LO")
print()
print("f(x) = x * x + 2 * x")
print()
print("Początek przedziału całkowania")
print()
xp = float(input("xp = "))
print()
print("Koniec przedziału całkowania")
print()
xk = float(input("xk = "))
print()

s = 0
dx = xk - xp
for i in range(1, n + 1):
    s += f(random.uniform(xp, xk))
s *= dx / n
print("Wartość całki wynosi : ",end="")
print("%.3f" % s)
print()
input("Naciśnij Enter...")

Wynik:

Obliczanie  całki oznaczonej
za pomocą metody Monte Carlo
----------------------------
(C)2004 mgr J.Wałaszek  I LO

f(x) = x * x + 2 * x

Początek przedziału całkowania

xp = 0

Koniec przedziału całkowania

xk = 1

Wartość całki wynosi : 1.331

Naciśnij Enter...

JavaScript

<html>
    <head>
        <title>Całkowanie
        numeryczne metodą
        Monte Carlo</title>
    </head>
    <body>

<div style="overflow-x: auto;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frm2"
        style="text-align:
               center;
               background-color:
               #E7E7DA">
          <b>Obliczanie całki
          oznaczonej<br>
          &nbsp;&nbsp;
          za pomocą metody
          Monte Carlo
          &nbsp;&nbsp;
          </b><br>
          (C)2026 mgr
          Jerzy Wałaszek
          <hr>
          Całkowana funkcja:
          <br>
          f(x) = x<sup>2</sup> + 2x
          <hr>
          Przedział całkowania<br>
          x<sub>p</sub> =
          <input
          name="xp_inp"
          size="15"
          value="0"
          type="text"
          style="text-align:
                 right">
          <br>
          x<sub>k</sub> =
          <input
          name="xk_inp"
          size="15"
          value="1"
          type="text"
          style="text-align:
                 right">
          <hr>
          <input
          onclick="js_montecarlo();"
          value="Oblicz całkę"
          name="B2"
          type="button">
          <hr>
          Wartość całki wynosi:
          <div id="out2">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script language="JavaScript">

// Obliczanie całki oznaczonej
// metodą Monte Carlo
// ---------------------------
// (C)2004 mgr Jerzy Wałaszek

// Tutaj definiujemy funkcję
function f(x)
{
  return(x * x + 2 * x);
}

function js_montecarlo()
{
  //liczba punktów losowych
  var N = 10000;
  var xp,xk,s,dx,i,t;

  xp = parseFloat(document
       .frm2.xp_inp.value);
  xk = parseFloat(document
       .frm2.xk_inp.value);
  if(isNaN(xp) || isNaN(xk))
    t = "<font color=red><b>" +
        "Popraw dane wejściowe!" +
        "</b></font>";
  else
  {
    s  = 0;
    dx = xk - xp;
    for(i = 1; i <= N; i++)
      s += f(xp + Math.random() *
           dx);
    s = dx * s / N;
    t = Math.floor(s * 1000) /
        1000;
  };
  document.getElementById("out2")
  .innerHTML = t;
}

</script>

    </body>
</html>

do podrozdziału do strony

Metoda Monte Carlo

Algorytm całkowania metodą Monte Carlo

Opis algorytmu

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

Schemat blokowy

Przykładowe programy