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]:
In [ ]:
In [48]:
integrate(func,-5,5,5000)
Out[48]:
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))
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]:
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]:
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]:
In [ ]: