Integrazione numerica con metodo Monte Carlo

import random

def f(x): return 3*x**2+2*x+1

a = 1
b = 5
fmax = 100
AREA = (b-a)*fmax

punti = 0
PUNTI = 5000
for p in range(PUNTI):
    x = random.uniform(a, b)
    y = random.uniform(0, fmax)
    if(y <= f(x)):
        punti += 1

rapporto = punti/PUNTI
area     = rapporto*AREA

print("PUNTI = %d, punti = %d, rapporto = %.4f, area = %.4f" %(PUNTI, punti, rapporto, area))

import matplotlib.pyplot as plt
import numpy as np
import random

def f(x): return 3*x*x+2*x+1
a=1
b=5
fmax=100
AREA=(b-a)*fmax
punti=0
PUNTI=5000

xp=[] # ASCISSE CASUALI
yp=[] # ORDINATE CASUALI --- SOTTO IL GRAFICO
xP=[] # ASCISSE CASUALI
yP=[] # ORDINATE CASUALI --- SOPRA IL GRAFICO

for p in range(PUNTI):
    x=random.uniform(a,b)
    y=random.uniform(0,fmax)
    if(y <= f(x)):
        punti += 1
        xp.append(x) # punto sotto
        yp.append(y) 
    else:
        xP.append(x) # punto sopra
        yP.append(y)

rapporto=punti/PUNTI 
area    =rapporto*AREA 

print("PUNTI=%d, punti=%d, rapporto=%.4f, area=%.4f" %(PUNTI,punti,rapporto,area))

X=np.linspace(a-1,b+1)
yf=[]
for x in X:
    yf.append(f(x))

plt.grid(which="major")
plt.plot(X , yf, color="blue" , linewidth="2")
plt.scatter(xp, yp, color="green", marker =".")
plt.scatter(xP, yP, color="red", marker ="."  )
plt.title("Integrazione numerica")

plt.show()

import matplotlib.pyplot as plt
import numpy as np
import random

def f(x): return 3*x*x+2*x+1
a=1
b=5
fmax=100
AREA=(b-a)*fmax
punti=0
PUNTI=1000

xp=[]  # ASCISSE CASUALI
yp=[]  # ORDINATE CASUALI --- SOTTO IL GRAFICO
xP=[]  # ASCISSE CASUALI
yP=[]  # ORDINATE CASUALI --- SOPRA IL GRAFICO
ya=[]  # Approssimazioni area

for p in range(1,PUNTI+1):
    x=random.uniform(a,b)
    y=random.uniform(0,fmax)
    if(y <= f(x)):
        punti += 1
        xp.append(x)  # punto sotto
        yp.append(y)        
    else:
        xP.append(x)  # punto sopra
        yP.append(y)
    area=punti/p*AREA # approssimazione
    ya.append(area)

rapporto=punti/PUNTI 
area    =rapporto*AREA 

print("PUNTI=%d, punti=%d, rapporto=%.4f, area=%.4f" %(PUNTI,punti,rapporto,area))

X=np.linspace(a-1,b+1)
yf=[]
for x in X:
    yf.append(f(x))

plt.subplot(2,1,1)
plt.grid()
plt.plot(ya, linewidth="2")
plt.title("Approssimazione area")
plt.subplot(2,1,2)
plt.grid()
plt.plot(X , yf, color="blue", linewidth="2")
plt.scatter(xp, yp, color="green", marker=".")
plt.scatter(xP, yP, color="red", marker=".")

plt.show()