Python/Scipy Integration array

Question

I am trying to write a program that does the following:

Takes values of V from an array
Passes the V values into an integral that is with respect to E
Output integral results into an array I
Plot I against V

The equation looks quite nasty, but everything is a constant other than V. Here is the equation. The equation isn't very important.

How should I go about this problem? My attempt (as shown below) does not calculate the integral for each value of V read from the file.

from scipy import integrate #integrate.quad
from numpy import *
import pylab
import datetime
import time
import os
import math

# import V
fn = 'cooltemp.dat'
V = loadtxt(fn,unpack=True,usecols=[1])

# variables
del1, del2, R, E, fE, fEeV = 1,2,1,2,1,1
e = 1.602176565*10**-19

# eqn = dint(abc)
a = E/( math.sqrt( E**2 - del1**2 ) )
b = ( E+ e*V )/( math.sqrt( ( E + e*V )**2) - del2**2)
c = fE-fEeV
d = 1/(e*R) # integration constant
eqn = a*b*c

# integrate 
result = quad(lambda E: eqn,-inf,inf)

# current
I = result*d

# plot IV curve
pylab.plot(V,I,'-r')

## customise graph
pylab.legend(['degree '+str(n),'degree '+str(q),'data'])
pylab.axis([0,max(x),0,max(y)])
pylab.xlabel('voltage (V)')
pylab.ylabel('current (A)')
tc = datetime.datetime.fromtimestamp(os.path.getmtime(fn))
pylab.title('IV curve\n'+fn+'\n'+str(tc)+'\n'+str(datetime.datetime.now()))
pylab.grid(True)
pylab.show()

* Updated attempt:

from scipy import integrate
from numpy import *
import pylab
import datetime
import time
import os
import math

# import V
fn = 'cooltemp.dat'
V = loadtxt(fn,unpack=True,usecols=[1])
# print V

# variables
del1, del2, R, E, fE, fEeV = 1.0,2.0,1.0,2.0,1.0,1.0
e = 1.602176565*10**-19

I=[]
for n in range(len(V)):

    constant = 1/(e*R) # integration constant
    eqn = (E/( math.sqrt( E**2 - del1**2 ) ))*(( E + e*V[n] )/( math.sqrt( ( E + e*V[n] )**2) - del2**2))*(fE-fEeV)

    # integrate 
    result,error = integrate.quad(lambda E: eqn,-inf,inf)
    print result
    # current
    I.append(result*constant)

I = array(I)

# plot IV curve
pylab.plot(V,I,'-b')

Looks like the error is in result = quad(lambda E: eqn, -inf, inf). lambda should contain the same variable on both sides of the colon, and that is not happening. By the way, that expression is not very readable. — heltonbiker
– heltonbiker, Commented Aug 15, 2012 at 14:50
Thanks heltonbiker. Does that mean I Should def the eqn and then try quad(lambda E: eqn, inf, -inf) ? I have since updated my code, I understand that the eqn is hard to follow, so it is no longer in chunks. — 8765674
– 8765674, Commented Aug 15, 2012 at 16:27
If eqn is a function, then you could use quad(lambda E: eqn(E), -inf, inf) but that would mean you don't need a lambda in the first place, using then directly quad(eqn, -inf, inf), supposing E varies continuously between -inf and inf. — heltonbiker
– heltonbiker, Commented Aug 15, 2012 at 16:37

ChrisB · Accepted Answer · 2012-08-15 18:27:38Z

5

You have a few problems:

The "function" you pass to quad always returns eqn, which is just a precalculated number. You need to define a proper function that takes a given value for E as input and returns the integrand. This function will also need to assume a fixed value for V. Assuming the code you provided calculates the proper quantity for a given value of V and E (I haven't checked, just copy-pasting):

# import V
fn = 'cooltemp.dat'
V = loadtxt(fn,unpack=True,usecols=[1])
# print V

@np.vectorize
def result(x):
    def integrand(E):
        del1, del2, R, fE, fEeV = 1.0,2.0,1.0,1.0,1.0
        e = 1.602176565*10**-19
        a = E/( math.sqrt( E**2 - del1**2 ) )
        b = ( E+ e*x )/( math.sqrt( ( E + e*x )**2) - del2**2)
        c = fE-fEeV
        d = 1/(e*R) # integration constant
        return a * b * c
    return quad(integrand, -inf, inf)

I = result(V)

To Summarize:

result(v) evaluates the full integral (over E) for a fixed value of v
integrand(E) evaluates the integrand at fixed E (the integration variable) and, incidentally, fixed V (it grabs the value from outside the function, which is why the definition for integrand is nested inside the definition for result)
the @np.vectorize trick is just a nice convenience function that allows you to pass arrays for V into result. Numpy will loop over the values for you, and return an array instead of a scalar

edited Aug 15, 2012 at 18:27

answered Aug 15, 2012 at 15:00

ChrisB

4,7587 gold badges32 silver badges43 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

8765674 Over a year ago

Thank you. I added a loop to my program as shown in the update. Do I still need to def the function prior to integrating so that E is a variable and remove E=1?

ChrisB Over a year ago

I've edited my answer to be a bit more explicit. I'm also not sure whether quad will know what to do with (-inf, inf) -- you may need to help it out and provide finite input (big negative and positive numbers)

8765674 Over a year ago

Thanks. Just tried that method and got the following errors: UserWarning: The integral is probably divergent, or slowly convergent. and ValueError: x and y must have same first dimension

BenMorel · Accepted Answer · 2014-06-24 17:04:35Z

You should use np.vectorize to pass arrays into the equations and get arrays back. For example, this calculates the following equation (Comoving distance if you are curious...):

comoving distance


import numpy as np
from scipy.integrate import quad

spl=299792458.0    #speed of light in m/s
Mpc=3.0856E22      # Mpc in m
pc=3.0856E16       # pc in m

def HubbleTime(H0):
    return 3.0856e17/(H0/100.0)

def HubbleDist(H0):
        """returns the Hubble Distance (in Mpc) for given H_0"""
    return spl*HubbleTime(H0)/Mpc

def Integrand(z, Om, OLam):
        """ This is the E(z) function from Hogg (2000)
        Integrand(z, Om, OLam)
    """
    return ( Om*(1+z)3 + OLam)(-0.5)

def CosmComDist(z, H0=70, Om=0.30, OLam=0.70):
    """Gives the comoving distance at redshift z
    CosmComDist(z, H0=70, Om=0.30, OLam=0.70)
    """
    CMD=HubbleDist(H0)*quad(Integrand, 0, z, args=(Om, OLam))[0]
    return CMD

CosmComDist=np.vectorize(CosmComDist)
redshifts = np.linspace(0,1,100)
distances = CosmComDist(redshifts)

Collectives™ on Stack Overflow

Python/Scipy Integration array

2 Answers 2

3 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related