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 = PI/2.0, 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

1. ustaw $c = 0$
2. wylosuj punkt $(x_i , y_i )$, gdzie wartości losowe $ x_p < x_i < x_k$ oraz $ y_p < y_i < y_k$
3. 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,
4. 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,
5. jeżeli wylosowany punkt $(x_i , y_i)$ nie spełnia żadnego z powyższych warunków wówczas, pozostawiamy zmienną $c$ bez zmian
6. powtórz $n$ krotnie kroki od 2 do 5
7. 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 6

Treść zadania znajduje się w Moodle pod adresem Zadanie 7

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$