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

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 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 <= a: continue
if k <= a and b <= k:
return numerical_integral(D[k],a,b,**kwds)
elif k <= a:
sum += numerical_integral(D[k],a,k,**kwds)
elif a <= k and k <= b:
sum += numerical_integral(D[k],k,k,**kwds)
elif k <= b <= k:
sum += numerical_integral(D[k],k,b,**kwds)
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 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 <= a1: continue
if k <= a1 and a2 <= k:
print 'one: key=',k, 'interval from', a1, 'to', a2
return numerical_integral(2*pi*x*B( Lambda,T[k]),a1,a2)
if k <= a1:
print 'two: key=',k, 'interval from', a1, 'to', k
sum += numerical_integral(2*pi*x*B( Lambda,T[k]),a1,k)
continue
if a1 <= k and k <= a2:
print 'three: key=',k, 'interval from', k, 'to', k
sum += numerical_integral(2*pi*x*B( Lambda,T[k]),k,k)
continue
if k <= a2 <= k:
print 'four: key=',k, 'interval from', k, 'to', a2
sum += numerical_integral(2*pi*x*B( Lambda,T[k]),k,a2)
if a1 >= k:
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. 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)

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?

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