# evaluvating variable inside a function while integrating

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 close merge delete

Sort by » oldest newest most voted

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))
more

@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?

( 2012-12-13 07:26:37 +0100 )edit

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

( 2012-12-13 07:59:31 +0100 )edit

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))
more

@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.

( 2012-12-14 02:29:32 +0100 )edit

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)
more