hi i have been given a question by my lecturer to integrate a function through python and he gave us very little information. the boundaries are +infinity and -infinity and the function is
(cos a*x) * (e**-x**2)
so far I have
def gauss_cosine(a, n):
sum=0.0
dx = ((math.cosine(a*x)*math.exp(-x**2)))
return
for k in range (0,n):
x=a+k*dx
sum=sum+f(x)
return dx*sum
not sure if this is right at all.
Well, your integral has an analytical solution, and you can calculate it with sympy, as @Bill pointed out, +1.
However, what I think was the point of the question is how to numerically calculate this integral, and this is what I discuss here.
The integrand is even. We reduce the domain to [0,+inf], and will multiply by 2 the result.
We still have an oscillatory integral on an unbounded domain. This is often a nasty beast, but we know that it is convergent, and well behaved at +- inf. In other words, the exp(-x**2) decays to zero fast enough.
The trick is to change variable of integration, x=tan(t), so that dx=(1+x**2)dt. The domain becomes [0,pi/2], it is bounded and the numerical integration is then a piece of cake.
Example with the simpson's rule from scipy, with a=2. With just 100 discretization points we have a 5 digits precision!
from scipy.integrate import simps
from numpy import seterr, pi, sqrt, linspace, tan, cos, exp
N = 100
a = 2.
t = linspace(0, pi / 2, N)
x = tan(t)
f = cos(a * x) * exp(-x ** 2) * (1 + x ** 2)
print "numerical solution = ", 2 * simps(f, t)
print "analytical solution = ",sqrt(pi) * exp(-a ** 2 / 4)