Testing nquad

We just try to see if python’s scipy.integrate library works in the correct way. For example, we try to test nquad on numerical integrations, especially multiple integrations.

import numpy as np
from scipy import integrate

Maybe an easy example is to calculate the area of a circle. Here we have the formula:

Specifically, if the radius of the circle is 2, we calculate the blue area in the graph.

We try to calculate it, using nquad of scipy.integrate library:

def integrandA(x0, y0, r0):
 return 1;
def domyA(r0):
 return [-r0, r0]
def domxA(y0, r0):
 return [- np.sqrt(r0**2 - y0**2 ), np.sqrt(r0**2 - y0**2)]
r0v=2
vA=integrate.nquad(integrandA, [domxA, domyA], args=(r0v,))

Then we get:

(area, error)  = (12.566370614359187, 8.001883600172732e-09) 
expected value =  12.566370614359172

So we see that the area is calculated accurately up to 13 decimal places.

The next example is to calculate the volume of a 3-dimensional sphere. The formula is presented as follows:

And we consider the case where the radius of the sphere is 2. The corresponding code is:

def integrandB(x0, y0, z0, r0):
 return 1;
def domzB(r0):
 return [-r0, r0]
def domyB(z0, r0):
 return [- np.sqrt(r0**2 - z0**2 ), np.sqrt(r0**2 - z0**2)]
def domxB(y0, z0, r0):
 return [- np.sqrt(r0**2 - z0**2 - y0**2 ), np.sqrt(r0**2 - z0**2 - y0**2 ) ]
r0v=2
vB=integrate.nquad(integrandB, [domxB, domyB, domzB], args=(r0v,))

Then we get:

volume, error  = (33.510321638291174, 8.001883600172732e-09) 
expected value =  33.510321638291124

where we can see that the calculation is accurate to 13 decimal places.