Física computacional: Números aleatorios.

Curso 2016-17 1

Números aleatorios y método de Monte Carlo.

Vamos a analizar ahora un método para evaluar integrales definidas utilizando números aleatorios, mé-
todo que es muy poderoso sobre todo cuando se trata de resolver integrales en múltiples dimensiones.

El teorema básico de Monte Carlo permite estimar integrales definidas multidimensionales de la forma
Z b1 Z b2 Z bn Z
I= dx1 dx2 · · · f (x1 , x2 , . . . , xn )dxn = f (~x)dn x (1)
a1 a2 an V
mediante la expresión
r
hf 2 i − hf i2
I ≈ V hf i ± V (2)
N
con
N N
1 X 2 1 X 2
hf i ≡ f (xi ) hf i ≡ f (xi ) (3)
N N
i=1 i=1
siendo V el volumen multidimensional, f (xi ) el valor de la función en cada uno de los N puntos xi
aleatoriamente distribuidos en V.

El término que aparece después del ± en la ecuación (2) representa una estimación del error de la
integral.

Cálculo de integrales de linea

EJERCICIO 1:
Implementar un programa que calcule la integral
Z 1
I= (1 − x2 )1.5 dx
0

mediante el método de Monte Carlo.1

RESULTADO:
Obtener un conjunto N de números aleatorios xi en el intervalo [0,1] y, de acuerdo con la ecuación
(2), aproxime la integral mediante la expresión

N
1 X
I ≈ V hf i ≡ f (xi ) (V = 1 − 0) (4)
N
i=1

En la tabla siguiente se da el valor de la integral para distintos números de puntos aleatorios.

N I Error
100 0.542766011609 3.15e-02
1,000 0.596106490640 1.04e-02
10,000 0.589884472906 3.29e-03
100,000 0.589156327124 1.05e-03
1,000,000 0.588452115452 3.32e-04
10,000,000 0.588974322757 1.05e-04
100,000,000 0.589047857125 3.32e-05

1 3π
La solución exacta es 16
≈ 0.58904862260 . . .
Física computacional: Números aleatorios. Curso 2016-17 2

EJERCICIO 2:
Implementar un programa que calcule la integral
Z ∞
I= e−x dx
0

mediante el método de Monte Carlo. 2

RESULTADO:
Hagamos el siguiente cambio de variable

1
y=
x+1
con lo que
1 dy
x= − 1, dx = −
y y2
y como cuando x = 0 ⇒ y = 1 y para x = ∞ ⇒ y = 0, nos queda finalmente
Z ∞ Z 1
−x 1 −( y1 −1)
I= e dx = e dy
0 0 y2

Obteniendo un conjunto N de números aleatorios yi en el intervalo [0,1] podemos aproximar la

integral por la ecuación (4). En la siguiente tabla se da la integral para distintos valores de puntos
aleatorios.

N I Error
100 0.908366745442 5.57e-02
1,000 0.991364933854 1.59e-02
10,000 1.009085861080 4.93e-03
100,000 0.999059706743 1.58e-03
1,000,000 1.000329142065 5.00e-04
10,000,000 1.000022080063 1.58e-04
100,000,000 1.000005711056 5.00e-05

Cálculo de integrales de superficie

EJERCICIO 3: y
Implementar un programa que calcule, mediante el
método de Monte Carlo, el área de un circulo de
radio R=1.5 cuyo centro se encuentra en el origen
de coordenadas. x

RESULTADO:
Inscribamos el círculo en un cuadrado de lado 2R
como muestra la figura 1. La ecuación (1) adopta
la forma
Fig. 1

2 1
La solución exacta es e1− x
Física computacional: Números aleatorios. Curso 2016-17 3

Z R Z R
I= dx f (x, y)dy ≈ V hf i (5)
−R −R

donde V = (R + R)(R + R) = 4R2 representa el área del cuadrado. Definamos f(x,y) de tal forma
que
1 si x2 + y 2 ≤ R2
(
f (x, y) =
0 si x2 + y 2 > R2
Se trazan N pares de puntos aleatorios de coordenadas xi , yi tales que −R < xi < R e −R <
yi < R y se obtiene el área del circulo por
N
4R2 X
Scir = 4R2 hf i = f (xi , yi )
N
i=1

En la siguiente tabla se da el valor del área del circulo para diferentes valores de N

N S Error
100 7.200000000000 3.60e-01
1,000 7.110000000000 1.16e-01
10,000 7.049700000000 3.71e-02
100,000 7.061940000000 1.17e-02
1,000,000 7.074495000000 3.69e-03
10,000,000 7.066345500000 1.17e-03
100,000,000 7.068223620000 3.70e-04
valor real 7.068583470577. . .

EJERCICIO 4: y

Implementar un programa que calcule, mediante el

método de Monte Carlo, el centro de masas de la
figura 2 supuesta homogénea. 6

RESULTADO: r=
2
x
Las coordenadas del centro de masas viene dado por
8
R R
xdm ydm Fig. 2
xG = R , yG = R
dm dm

siendo dm = ρdxdy

En la siguiente tabla se dan los resultados obtenidos para diferentes valores de N

N xG Error yG Error
100 3.322788417 2.75e-01 3.390963979 2.28e-01
1,000 3.436734471 1.01e-01 3.238444091 8.25e-02
10,000 3.759905273 3.15e-02 3.337816926 2.47e-02
100,000 3.715825414 1.01e-02 3.325536089 7.95e-03
1,000,000 3.694471614 3.19e-03 3.323785383 2.51e-03
10,000,000 3.698428933 1.01e-03 3.323793502 7.95e-04
100,000,000 3.698447315 3.19e-04 3.324003108 2.52e-04
Real: xG = 3.698769 . . . yG = 3.323999 . . .
Física computacional: Números aleatorios. Curso 2016-17 4

Cálculo de integrales de volumen

EJERCICIO 5: z

Implementar un programa que calcule, mediante el

método de Monte Carlo, el volumen de una esfera
de radio R=1.5 cuyo centro se encuentra en el
origen de coordenadas.
y

RESULTADO: x
Inscribamos la esfera en un cubo de lado 2R como
muestra la figura 3. La ecuación (1) adopta la forma
Z R Z R Z R Fig. 3
I= dx dy f (x, y, z)dz ≈ Vcub hf i (6)
−R −R −R

donde Vcub representa el volumen del cubo. Definamos f(x,y,z) de tal forma que

1 si x2 + y 2 + z 2 ≤ R2
(
f (x, y, z) =
0 si x2 + y 2 + z 2 > R2

Si trazamos N pares de puntos aleatorios de coordenadas xi , yi , zi tal que −R < xi < R, −R <
yi < R y −R < zi < R, obteniendo cuantos de ellos, Nint , caen dentro de la esfera, resulta que
Nint
Vesf = Vcub hf i = Vcub
N
En la siguiente tabla se da el valor del área del circulo para diferentes valores de N

N S Error
100 16.470000000000 1.32e+00
1,000 13.878000000000 4.27e-01
10,000 14.253300000000 1.35e-01
100,000 14.135310000000 4.26e-02
1,000,000 14.124915000000 1.35e-02
10,000,000 14.136473700000 4.26e-03
100,000,000 14.137902270000 1.35e-03
valor real 14.137166941154. . .
y
EJERCICIO 6: 4

Implementar un programa que calcule, mediante el 2

método de Monte Carlo, el volumen del trozo de toro 1

circular representado en la figura 4. 0 1 2 3 4 x

La región de interés se encuentra contenida entre los

planos x = 1 y x = 4, y = −3 e y = 4, z = −1 y z = 1.

Fig. 4

RESULTADO:
La ecuación de un toro circular centrado en el origen de coordenadas y de eje de revolución el
eje z es
Física computacional: Números aleatorios. Curso 2016-17 5

p 2
x2 + y 2 − R + z 2 = r 2

En la siguiente tabla se da el valor del volumen para diferentes valores de N

N V Error
100 24.360000000000 2.07e+00
1,000 21.042000000000 6.64e-01
10,000 22.016400000000 2.10e-01
100,000 22.027740000000 6.63e-02
1,000,000 22.062852000000 2.10e-02
10,000,000 22.096582200000 6.63e-03
100,000,000 22.098762420000 2.10e-03