How to deal with numerical integration in python with small results?-CodePudding

I am running into an issue with integration in Python returning incorrect values for an integral with a known analytical solution. The integral in question is

LaTex expression for the integral (can't post photos yet)

For the value of sigma I am using (1e-15),the solution to this integral has a value of ~ 1.25e-45. However when I use the scipy integrate package to calculate this I get zero, which I believe has to do with the precision required from the calculation.

#scipy method
import numpy as np
from scipy.integrate import quad

sigma = 1e-15
f = lambda x: (x**2) * np.exp(-x**2/(2*sigma**2))

#perform the integral and print the result 
solution = quad(f,0,np.inf)[0]
print(solution)
0.0

And since precision was an issue I tried to also use another recommended package mpmath, which did not return 0, but was off by ~7 orders of magnitude from the correct answer. Testing larger values of sigma result in the solution being very close to the corresponding exact solution, but it seems to get increasingly incorrect as sigma gets smaller.

#mpmath method
import mpmath as mp

sigma = 1e-15
f = lambda x: (x**2) * mp.exp(-x**2/(2*sigma**2))

#perform the integral and print the result 
solution = mp.quad(f,[0,np.inf])
print(solution)
2.01359486678988e-52

From here I could use some advice on getting a more accurate answer, as I would like to have some confidence applying python integration methods to integrals that cannot be solved analytically.

CodePudding user response：

you should add extra points for the function as 'mid points', i added 100 points from 1e-100 to 1 to increase accuracy.

#mpmath method
import numpy as np
import mpmath as mp

sigma = 1e-15
f = lambda x: (x**2) * mp.exp(-x**2/(2*sigma**2))

#perform the integral and print the result
solution = mp.quad(f,[0,*np.logspace(-100,0,100),np.inf])
print(solution)

1.25286197427129e-45

Edit: turns out you need 10000 points instead of 100 points to get a more accurate result, of 1.25331413731554e-45, but it takes a few seconds to calculate.

CodePudding user response：

Most numerical integrators will run into issues with numbers that small due to floating point precision. One solution is to scale the integral before calculating. Letting q -> x/sigma, the integral becomes:

f = lambda q: sigma**3*(q**2) * np.exp(-q**2/2)
solution = quad(f, 0, np.inf)[0]

# solution: 1.2533156529417088e-45