Ask Your Question
0

evaluvating variable inside a function while integrating

asked 2012-12-12 18:54:10 -0500

indiajoe gravatar image

updated 2012-12-13 00:58:34 -0500

I got a strange problem. The code to reproduce the problem is given below

from scipy.constants import h, c, k
def T2(x):
    a=11717  
    if x < 21500 : return a*(x**-0.53)
    else :
        print(x)  # Just for debugging
        return a*(x**-0.75)
# Blackbody Planky function
def B(Lambda,Temp):
    return 2*h*c**2/(Lambda**5 *(exp(h*c/(Lambda*k*Temp))-1))

def flux2(Lambda):
    return numerical_integral(2*pi*x*B(Lambda,T2(x)),7,2150)[0]

print 2.5*log(flux2(9000*10**-10))

The problem is inside the T2() function. Since the integral in x is from 7 to 2150, The if condition should get satisfied. and return the a(x*-0.53) . But instead it is evaluating the else condition. print x is printing the alphabet 'x' instead of the value of variable x it is supposed to take during each point in integral. I guess i have understood how these functions work inside an integral wrongly. What is it that I am doing wrong here?


Update: I instead tried the Piecewise() function to define T2() as follows,

a=11717
T2= Piecewise([[(0,215),a*(x**-0.53)],[(215,21500),a*(x**-0.75)]])

but inside the integral function I am getting ValueError

ValueError: Value not defined outside of domain.

edit retag flag offensive close merge delete

3 answers

Sort by ยป oldest newest most voted
0

answered 2012-12-12 21:21:17 -0500

ndomes gravatar image

When calling T2(x) inside the function flux2 x is an element of Symbolic Ring and bool(x < 21500) gives False.

The range 7 to 2150 doesn't matter because T2(x) is evaluated before any substitution of x happens.

You can get the desired behavior with an assumption ( try with and without forget() )

assume(x<=2150)
#print assumptions()
#forget()
from scipy.constants import h, c, k
def T2(x):
    a=11717  
    if x < 21500 : 
        return a*(x**-0.53)
    else :
        print(x.parent())  # Just for debugging        
        return a*(x**-0.75)
# Blackbody Planky function
def B(Lambda,Temp):
    return 2*h*c**2/(Lambda**5 *(exp(h*c/(Lambda*k*Temp))-1))

def flux2(Lambda):
    return numerical_integral(2*pi*x*B(Lambda,T2(x)),7,2150)[0]

print 2.5*log(flux2(9000*10**-10))
edit flag offensive delete link more

Comments

@ndomes : Thanks for the help. But what i really want is the function T2(x) to give different function output corresponding to the value of x. It is a piecewise defined function. How can i do that?

indiajoe gravatar imageindiajoe ( 2012-12-13 00:26:37 -0500 )edit

I tried the Piecewise() function but now getting a Value error (Details I have updated in question)

indiajoe gravatar imageindiajoe ( 2012-12-13 00:59:31 -0500 )edit
0

answered 2012-12-13 06:06:56 -0500

ndomes gravatar image

another try:

put the functions in a python dictionary with the respective intervals as keys.

run through the sorted keys and add up the integrals for each fitting interval

from scipy.constants import h, c, k

a=11717

f1(x) = a*(x**-0.53)
f2(x) = a*(x**-0.75)

T = {(1,10000):f1,(10000,20000):f2, (20000,50000):f1}


# Blackbody Planky function
def B(Lambda,Temp):
    return 2*h*c**2/(Lambda**5 *(exp(h*c/(Lambda*k*Temp))-1))

a1 = 300; a2 = 22000  # interval

def flux2(Lambda):
    sum = 0   
    L = copy(T.keys())
    L.sort()
    for k in L:
        if k[1] <= a1: continue
        if k[0] <= a1 and a2 <= k[1]:
            print 'one: key=',k, 'interval from', a1, 'to', a2
            return numerical_integral(2*pi*x*B( Lambda,T[k]),a1,a2)[0]            
        if k[0] <= a1: 
            print 'two: key=',k, 'interval from', a1, 'to', k[1]
            sum += numerical_integral(2*pi*x*B( Lambda,T[k]),a1,k[1])[0]
            continue
        if a1 <= k[0] and k[1] <= a2:
            print 'three: key=',k, 'interval from', k[0], 'to', k[1]
            sum += numerical_integral(2*pi*x*B( Lambda,T[k]),k[0],k[1])[0]
            continue
        if k[0] <= a2 <= k[1]:
            print 'four: key=',k, 'interval from', k[0], 'to', a2
            sum += numerical_integral(2*pi*x*B( Lambda,T[k]),k[0],a2)[0]
        if a1 >= k[1]:
            pass            
    return sum

print 2.5*log(flux2(9000*10**-10))
edit flag offensive delete link more

Comments

@ndomes: Thankyou, This is fine. but i would have preferred something more simply to integrate a piecewise defined function. So that all the limits of function changes are defined inside the function and the outside integrations commands are independent of what was inside the function.

indiajoe gravatar imageindiajoe ( 2012-12-13 19:29:32 -0500 )edit
0

answered 2012-12-15 00:21:21 -0500

ndomes gravatar image

To get it more convenient I wrapped numerical_integral for use with piecewise functions:

def numerical_integral_piecewise(f,a,b,**kwds):

    """
    Returns a numerical integral from a piecewise function 

    ``f``  - piecewise function

    ``a``  - lower bound of the interval

    ``b``  - upper bound of the interval

    ``**kwds`` -  options passed to numerical_integral


    EXAMPLES::

        sage: f1(x) = 1
        sage: f2(x) = 2
        sage: f3(x) = 3
        sage: f = piecewise([[(0,4),f1],[(4,6),f2],[(6,10),f3]])
        sage: numerical_integral_piecewise(f,2,9)
        15.0

        if a (or b ) is outside the domain of f, the interval (a,b) is cut 
        to the lower bound (or upper bound) of the domain of f

        sage: numerical_integral_piecewise(f,-2,20)
        20.0

    TODO
        return an error estimate as numerical_integral does

    """    

    D = dict(f.list())
    sum = 0   
    for k in D.keys():
        if k[1] <= a: continue
        if k[0] <= a and b <= k[1]:
            return numerical_integral(D[k],a,b,**kwds)[0]            
        elif k[0] <= a: 
            sum += numerical_integral(D[k],a,k[1],**kwds)[0]
        elif a <= k[0] and k[1] <= b:
            sum += numerical_integral(D[k],k[0],k[1],**kwds)[0]
        elif k[0] <= b <= k[1]:
            sum += numerical_integral(D[k],k[0],b,**kwds)[0]        
    return sum


f1(x) = 1
f2(x) = 2
f3(x) = 3
f = piecewise([[(0,4),f1],[(4,6),f2],[(6,10),f3]])
numerical_integral_piecewise(f,-2,20)
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-12-12 18:54:10 -0500

Seen: 213 times

Last updated: Dec 15 '12