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 0x112c7a0d7b8>






    

In [ ]:

    




    

In [5]:

    

integrate(func,-5,5,5000)


    

    
Out[5]:





1666.3885134827892

In [6]:

    

true = 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,5000)):
    i_l.append(i)
    y = integrate(func, a, b, i)
    y_l.append(y)
    err = true - y
    err_l.append(err)
    err_l_theor.append(np.sqrt(D/i))

for i in tqdm(range(5000,n_max,150)):
    i_l.append(i)
    y = integrate(func, a, b, i)
    y_l.append(y)
    err = true - y
    err_l.append(err)
    err_l_theor.append(np.sqrt(D/i))

for i in tqdm(range(n_max+1, 100*(n_max+1), n_max)):
    i_l.append(i)
    y = integrate(func, a, b, i)
    y_l.append(y)
    err = true - y
    err_l.append(err)
    err_l_theor.append(np.sqrt(D/i))


    

    




100%|███████████████████████████████████████████████████████████████████████████| 4999/4999 [00:00<00:00, 10940.94it/s]
100%|██████████████████████████████████████████████████████████████████████████████| 9967/9967 [07:35<00:00, 10.83it/s]
100%|████████████████████████████████████████████████████████████████████████████████| 100/100 [07:54<00:00, 12.84s/it]

In [7]:

    

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[7]:





[<matplotlib.lines.Line2D at 0x112c7fb7860>]






    

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 [ ]: