Efficient way for calculating integrals?

Question

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.

Uhm, no, integrating 1-dimensional, nicely bounded, infinitely differentiable functions with rapidly decaying derivatives with Monte-Carlo is indeed a *very* bad idea. But where did you get the idea that trapezium rule isn't complete garbage either? There are much better methods for numerical integration. What's wrong with [`scipy.integrate`](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html)? — Andrey Tyukin, May 17 '18 at 15:44
@AndreyTyukin I am still very beginner in programming, the other methods I found were using too many new commands to comprehend, so I am starting from the simplest atm, but if you recommend any certain method I am listening . — JKM, May 17 '18 at 15:47
"very beginner" is relative... It's not at all obvious that an average "beginner" would be able to implement a Monte Carlo method for computing integrals. It's not really programming-specific either. — Andrey Tyukin, May 17 '18 at 16:01

score 1 · Accepted Answer · answered May 17 '18 at 15:57

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).

score 0 · Answer 2 · answered Oct 12 '18 at 21:44

You can also use the Riemann approximation. Below code is in Java

package math;

import java.util.Optional;
import java.util.function.*;
import java.util.stream.IntStream;
import static java.lang.Math.*;

public class IntegralJava8
{
    public interface Riemann extends BiFunction<Function<Double, Double>, Integer, 
    BinaryOperator<Double>> { }

    public static void main(String args[])
    {
        int N=100000;
        Riemann s = (f, n) -> (a, b) -> 
        IntStream.range(0, n)
        .mapToDouble(i -> f.apply(a + i * ((b - a) / n)) * ((b - a) / n)).sum();
        Optional<Double> gaussIntegral = 
            Optional.of(s.apply(x -> exp(-pow(x, 2)), N).apply(-1000.0, 1000.0));
        gaussIntegral.ifPresent(System.out::println);
    }
}

In above class, it will calculate the Gaussian Integral from -infinity to infinity which is equal to square root of PI (1.772)

Efficient way for calculating integrals?

2 Answers2