CODEWITHSUNDEEP
pythonscipy
coalescedance
coalescedance

Reputation: 25

Using scipy to integrate a function?

I am trying to use SciPy to integrate this function:

y(x) = (e^-ax)*cos(x) between 0 and 4pi.

Here is the code I have so far:

from numpy import *
from scipy.integrate import simps
a = 0 
x = linspace(0 , 4*pi, 100)
y = (e^-(a*x))*cos(x)
integral_value = simps(y,x)
print integral_value

However it doesn't seem to be working. Any help would be appreciated, thanks!

Upvotes: 2

Views: 2914

Answers (3)

Bill Bell
Bill Bell

Reputation: 21663

Cracking a nut with a sledgehammer ...

>>> from sympy import *
>>> var('a x')
(a, x)
>>> y = exp(-a*x)*cos(x)
>>> y.subs(a,0)
cos(x)

Should it surprise any of us to find that this integrates to 0 (zero) over the given interval?

Upvotes: 1

willeM_ Van Onsem
willeM_ Van Onsem

Reputation: 477883

Well if you run the program you obtain the following error:

TypeError: ufunc 'bitwise_xor' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

So you know the problem is with the ^ (bitwise xor) in your function. In Python one uses ** to take the exponent.

If one writes:

y = (e**-(a*x))*cos(x)

instead, one gets:

>>> print integral_value
-0.000170200006112

The full program:

from numpy import *
from scipy.integrate import simps
a = 0 
x = linspace(0 , 4*pi, 100)
y = (e**-(a*x))*cos(x)
integral_value = simps(y,x)
print integral_value

You can also make explicit use of numpy functions with:

from numpy import *
from scipy.integrate import simps
a = 0 
x = linspace(0 , 4*pi, 100)
y = exp(-a*x)*cos(x)
integral_value = simps(y,x)
print integral_value

In order to increase the precision, you can increase the number of points (100 is not that much).

Upvotes: 5

ewcz
ewcz

Reputation: 13097

import numpy as np
import math
from scipy.integrate import simps
a = 0
x = np.linspace(0 , 4*math.pi, 100)

#create a vectorized function which can be applied directly to an array    
fn = np.vectorize(lambda x: math.exp(-a*x)*math.cos(x))
y = fn(x)

integral_value = simps(y, x)
print integral_value

this indeed yields value of -0.000170200006112. Note that for a=0, the integral is equal to zero, so in order to get "closer" you would need to refine the grid...

Upvotes: 1

Related Questions