# Implementing numerical integration in Python

I have this Python code that implements a rectangular numerical integration. It evaluates the (K-1)-dimensional integral for arbitrary integer $K \geq 1$

$$\int_{u_K = 0}^{\gamma}\int_{u_{K-1} = 0}^{\gamma-u_K}\cdots\int_{u_2}^{\gamma-u_K-\cdots-u_3}F_U(\gamma-\sum_{k=2}^Ku_k)f_U(u_2)\cdots f_U(u_K)\,du_2\cdots du_K$$

where $F_U$ corresponds to the cdf function in the code, and $f_U$ to the pdf. I implemented it using recursion as follows:

    #************************** Import necessary libraries******************************************
import numpy as np
import matplotlib.pyplot as plt
import time
#******************************Set the constant scalars and vectors***************************
start_time = time.time()
KU = 3
eta_dB = 3
eta = 10**(eta_dB/10)
ExpanF = 50

tdB = np.arange(-5,11,4)
tVec = 10**(tdB/10)
thVec = (tVec/eta)*(ExpanF-eta*(KU-1))

N = 10000 # For simulation

du = 0.01

#******************************Define functions to be used***************************************
#Define the CDF of U
def CDF(u):
return 1-1/(u+1)

#Define the PDF of U
def pdf(u):
return 1/((1+u))**2

def FK(h, k):
#print(f'h is {h}, and k is {k}')
if k == 1:
res = CDF(h)
else:
#n_iter = int(h/du)
res =  0
u = 0
while u < h:
res += FK(h-u, k-1)*pdf(u)*du
u += du
return res

#*******************Find the numerical and simulation values of the integral******************
ResultNum = []
ResultSim = []

if (ExpanF-eta*(KU-1)) > 0:
for t in thVec:
# Numerical
ResultNum.append(1-FK(t, KU))
# Simulation
count = 0
for n in range(0, N):
if np.sum(np.random.exponential(1, (1, KU))/np.random.exponential(1, (1, KU))) <= t:
count += 1
ResultSim.append(1-count/N)


The parameter du is the resolution of the numerical integral, and the smaller, the more accurate the numerical approximation is. However, decreasing du comes at the expense of more computational time using the above code. For example I tried to run the above code on my machine Intel i5 @2.5 GH with 6 GB of RAM for du=0.0001, and it took 15 hours and didn't finish, at which case I had to abort it.

Is there anyway the above code can be optimized to run significantly faster? I tried C++, and I got a 30x speed factor using the same code and algorithm, but I was wondering if I can get a similar or better speed up factor in Python.

• Please fix your indentation. The easiest way to post code is to paste it into the question editor, highlight it, and press Ctrl-K to mark it as a code block. – 200_success Jul 29 '18 at 17:56
• Did you try scipy.integrate.nquad? – Gareth Rees Jul 30 '18 at 8:05
• @GarethRees Not really. Can it be used for an arbitrary number of nested integrals? – BlackMath Jul 30 '18 at 14:37
• See the documentation: "Function signature should be func(x0, x1, ..., xn, t0, t1, ..., tm). Integration is carried out in order. That is, integration over x0 is the innermost integral, and xn is the outermost." – Gareth Rees Jul 30 '18 at 15:15
• @GarethRees Doesn't sound like what I am looking for. I don't want to specify the functions and limits each time I change K. I need a dynamic program where all I change is K. That is why I am trying to program the numerical integration from scratch. Otherwise, I would have used the built-in numerical integrations in MATLAB or any other platform. But these functions are not flexible for an arbitrary number of nested integrals. – BlackMath Jul 30 '18 at 16:00