I know how to use the Monte-Carlo to calculate integrals, but I was wondering if it is possible to use the trapezium rule combined with numpy for efficiency to get the same integral, I am not sure which one is the fastest or is the latter possible?
for e.g. to integrate e**-x**2 > y I could use the monte-carlo method like this:
import numpy as np
import matplotlib.pyplot as plt
X = np.random.rand(500000,2)
X[:,0] = X[:,0]*4-2
J = np.where(X[:,1] < np.exp(-X[:,0]**2))[0]
Xp = X[:2000]
Ip = [i for i in range(len(Xp)) if i in J]
Inp = [i for i in range(len(Xp)) if i not in J]
plt.plot(Xp[Ip,0],Xp[Ip,1], 'bd', Xp[Inp,0],Xp[Inp,1], 'rd')
plt.show()
And this can calculate very easily :
print len(J) / 500000.0 * 4
Which gives:
1.767784
In this case it was easy but if the intervals are not specified given in like [a,b] , n, and I want to make a function then I think the method above is not really efficient, at least I think so.
So , my question is can I integrate a continuous function like e.g.cos(x)/x for a determined interval like [a,b] in a function with trapezium rule?
And is it better than the method i used here?
Every advice is welcome.
Just use scipy.integrate.quad:
from scipy import integrate
from np import inf
from math import exp, sqrt, pi
res, errEstimate = integrate.quad(lambda x: exp(-x**2), -inf, +inf)
print(res) #output: 1.7724538509055159
print(sqrt(pi)) #output: 1.7724538509055159
The last line simply checks that the evaluated integral is indeed square root of Pi (it's the Gaussian integral).