CODEWITHSUNDEEP
pythonnumerical-integration
user3346077
user3346077

Reputation: 13

Integration on python

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.

Upvotes: 1

Views: 3493

Answers (3)

gg349
gg349

Reputation: 22691

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)

Upvotes: 1

wflynny
wflynny

Reputation: 18531

I don't see it recommended much on this site, but you could try sympy:

In [1]: import sympy as sp
In [2]: x, a = sp.symbols(('x', 'a'))
In [3]: f = sp.cos(a*x) * sp.exp(-x**2)
In [4]: res = sp.integrate(f, (x, -sp.oo, sp.oo))
In [5]: res
Out[5]: sqrt(pi) * exp
In [6]: sp.pprint(res)
Out[6]:
         2 
       -a  
       ────
  ___   4  
╲╱ π ⋅ℯ

For numerical integration, try the scipy package.

Upvotes: 2

duffymo
duffymo

Reputation: 308898

Your computer will have a very hard time representing those boundary limits.

Start by plotting your function.

It also helps to know the answer before you start.

I'd recommend breaking it into two integrals: one from minus-infinity to zero and another from zero to plus-infinity. As noted by flebool below, it's an even function. Make sure you know what that means and the implications for your solution.

Next you'll need an integration scheme that can deal with boundary conditions at infinity. Look for a log quadrature scheme.

A naive Euler integration would not be my first thought.

Upvotes: 0

Related Questions