In [1]:
import numpy as np
from matplotlib import pylab as plt
import numpy.random as rand
from tqdm import tqdm
from mpl_toolkits.mplot3d.axes3d import Axes3D
%matplotlib inline

In [2]:
def integrate(func, a, b, n):
    x = np.random.uniform(a, b, (n, 2))
    #print(x)
    #print(x.shape)
    res = np.sum(func(x))
    res = res*((np.abs(b-a))**2)/n
    return res

def func(x):
    res = np.square(x[:,0]) + np.square(x[:,1])
    return res

def grid_vect(a,b,N):
    x = np.linspace(a, b, N, endpoint=True)
    h = x[1] - x[0]
    return np.mgrid[a:b+h:h, a:b+h:h].reshape(2,-1).T

In [3]:
a = -5
b = 5
n = 50

In [4]:
xy = grid_vect(a,b,n)
Z = func(xy)
Z = Z.reshape(1,Z.size)
fig = plt.figure(figsize=(14,6))
ax = fig.add_subplot(1, 2, 1, projection='3d')
x = np.linspace(a, b, n, endpoint=True)
ax.plot_wireframe(xy[:,0], xy[:,1], Z)


Out[4]:
<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x232297bd4a8>

In [ ]:


In [48]:
integrate(func,-5,5,5000)


Out[48]:
1668.5064927483595

In [14]:
true_val = 1666.67
n_max = 1500000
err_l = []
err_l_theor = []
y_l = []
i_l = []

D = np.square(b-a)/12

for i in tqdm(range(1,50000)):
    i_l.append(i)
    y = integrate(func, a, b, i)
    y_l.append(y)
    err = np.abs(true_val - y)
    err_l.append(err)
    err_l_theor.append(np.sqrt(D/i))


100%|██████████████████████████████████████████████████████████████████████████| 49999/49999 [00:33<00:00, 1491.72it/s]

In [15]:
fig = plt.figure(figsize=(14,7))
axes = fig.add_subplot (1, 1, 1)
plt.title("Error vs N, Monte-Carlo method test")
plt.xlabel("N, points")
plt.ylabel("Abs. error")
axes.set_yscale ('log')
axes.set_xscale ('log')
plt.plot(i_l, err_l, 'rx')
plt.plot(i_l, err_l_theor)


Out[15]:
[<matplotlib.lines.Line2D at 0x2322a103b00>]

In [8]:
xxx = np.linspace(1,500000000,30000)
yyy = np.sqrt(1/xxx)
fig = plt.figure(figsize=(14,7))
axes = fig.add_subplot (1, 1, 1)
axes.set_yscale ('log')
axes.set_xscale ('log')
plt.plot(xxx,yyy)


Out[8]:
[<matplotlib.lines.Line2D at 0x112805377b8>]

In [11]:
fig = plt.figure(figsize=(14,7))
axes = fig.add_subplot (1, 1, 1)
plt.title("Error vs N, Monte-Carlo method test")
plt.xlabel("N, points")
plt.ylabel("Abs. error")
axes.set_yscale ('log')
axes.set_xscale ('log')
plt.plot(i_l, np.array(err_l)*(1/true), 'rx')
plt.plot(i_l, np.array(err_l_theor)*(1/true))


Out[11]:
[<matplotlib.lines.Line2D at 0x112809d4f98>]

In [ ]: