Lab. 7 Całkowanie Monte Carlo

Całka oznaczona

$$F = \int_{a}^{b} f(x) d x$$

pole powierzchni pod funkcją na odcinku $[a,b]$

Zwykły estymator Monte Carlo

$$\left\langle F^{N}\right\rangle=(b-a) \frac{1}{N} \sum_{i=1}^{N} f\left(X_{i}\right)$$

gdzie $X_{i}$ jest zmienną losową z rozkładu jednostajnego na odcinku $[a,b]$.

Dla próby $N=4$

Warto zauważyć, że $$\operatorname{Pr}\left(\lim _{N \rightarrow \infty}\left\langle F^{N}\right\rangle=F\right)=1$$

Generator liczb losowych w C

random.c

#include<stdio.h>
#include<stdlib.h>
#include<time.h>

int main()
{
   double a=0.0, b=1.0, x;
   int i, n = 10;
   
   srand(time(0));

   for(i=0; i<n; i++)
   {
       x = a + (b - a) * rand()/RAND_MAX;
       printf("%f\n", x);
   }
}

Przykład

Wyznaczenie całki $F=\int_{0}^{\frac{\pi}{2}} \sin (x) d x$

Wartość analityczna $F=[-\cos (x)]_{0}^{\frac{\pi}{2}}=-\cos \left(\frac{\pi}{2}\right)--\cos (0)=1$

mc1.c

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>

#define PI  3.14159265359

double f(double x)
{
   return sin(x);
}

double mc1(double a, double b, int n, double (*func)(double))
{
   double s = 0.0, x, fx;
   int j;

   for (j=0; j < n; j++){
      x = a + (b-a) * ((double)rand() / RAND_MAX);
      s += (*func)(x);
    }
   s = (b - a) * s / n;
   return s;
}

int main()
{
   double a = 0.0, b = 2*PI, s;
   int n;

   printf("Ile losowan n=");
   scanf("%d", &n);
   
   srand(time(0));

   s = mc1(a, b, n, f);
   printf("n=%6d s=%10.6f \n", n, s );
}

Metoda MC akceptacji-odrzuceń

Pole prostokąta wyznaczonego przez zakres całkowania

$$P=\left|x_{k}-x_{p}\right|\left|y_{k}-y_{p}\right|$$

$$\frac{P_{\text {prostokata }}}{\text { Calka }}=\frac{n}{c}$$

$$\text { calka }=P_{\text {prostokata }} \cdot \frac{c}{n}=\left|x_{k}-x_{p}\right| \cdot\left|y_{k}-y_{p}\right| \cdot \frac{c}{n}$$

Algorytm

ustaw $c = 0$
wylosuj punkt $(x_i , y_i )$, gdzie wartości losowe $ x_p < x_i < x_k$ oraz $ y_p < y_i < y_k$
jeżeli wylosowany punkt $(x_i , y_i )$ leży nad osią OY i jednocześnie pod wykresem funkcji całkowanej, czyli spełnia nierówność $0 < y_i \leq f(x_i)$, wówczas zwiększamy zmienną $c$ o jeden,
jeżeli wylosowany punkt $(x_i , y_i )$ leży pod osią OY i jednocześnie nad wykresem funkcji całkowanej, czyli spełnia nierówność $0 > y_i \geq f(x_i )$, wówczas zmniejszamy zmienną $c$ o jeden,
jeżeli wylosowany punkt $(x_i , y_i)$ nie spełnia żadnego z powyższych warunków wówczas, pozostawiamy zmienną $c$ bez zmian
powtórz $n$ krotnie kroki od 2 do 5
zwróć wynik $$\left|x_{k}-x_{p}\right| \cdot\left|y_{k}-y_{p}\right| \cdot \frac{c}{n}$$

Własności metody MC

wartość całki jest przybliżeniem (estymatorem) o rozkładzie normalnym
dokładność rośnie (tj. wariancja maleje) ze wzrostem rozmiaru próbki losowej $\sigma\left[\left\langle F^{N}\right\rangle\right] \propto \frac{1}{\sqrt{N}}$
bardzo prosta metoda, która z powodzeniem może być stosowana do obliczania całek wielowymiarowych
procedura prosta do zrównoleglenia
powolna zbieżność
oszacowanie błędu nie zawsze proste
wyniki zależą od jakości generatora liczb
w metodzie odrzuceń musimy znać zakres wartości maksymalnych i minimalnych funkcji całkowanej na odcinku $[x_p, x_k]$

Zadanie

Napisz program, który za pomocą metody Monte Carlo wyznaczy przybliżoną wartość liczby $\pi$ poprzez wyznaczenie pola powierzchni ćwiartki koła o promieniu $r=1$ oraz oszacuje odchylenie standardowe wyniku. Pole ćwiartki koła równe jest wartości całki

$$P = \int_0^1{\sqrt{1-x^2}} dx$$

Z wzoru na pole koła wynika też, że $P = \frac{1}{4}\pi r^2 = \frac{\pi}{4}$

Procedurę Monte Carlo dla danej wartości losowań $k$ powtórz $n=100$ i wyznacz wartość średnią $\mu$ oraz odchylenie standardowe $\sigma$ z uzyskanej serii przybliżeń $x_1, x_2, \ldots, x_n$ zgodnie ze wzorami $$\mu = \frac{1}{n} \sum_{i=1}^n x_i \qquad \sigma_n=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(x_{i}-\mu\right)^{2}}=\sqrt{\left(\frac{1}{n}\sum_{i=1}^{n} x_{i}^{2} \right)-\mu^{2}}$$

Wartość $k$ określającą ilość losowań pojedynczej procedury MC podaje użytkownik na początku działania programu. Program wypisz wyniki dla dwóch poznanych metod: Monte Carlo estymującej wartość oczekiwaną oraz za pomocą metody akceptacji-odrzuceń.

Dodatkowe ćwiczenia

Ćw. 1

Używając metody Monte Carlo wyznacz objętość sfery o promieniu 1

$$x^{2}+y^{2}+z^{2} \leq 1$$

Ćw. 2

Za pomocą metody Monte Carlo rozwiąż zadanie z poprzednich zajęć aby wyznaczyć przybliżoną wartość liczby $\pi$