Solve for a definite integral for each element within an array that is defined by a function in Python

Question

I am attempting to integrate the values of an array. The Array used within the function is "H" and limits of integration for the integral are "surf" and "base". The values are read in from a .txt file. 'def flow_pramB(z)' defines the function to be integrated and 'def integrand(a,b)' performs the integral. Then I would like to call the function 'integrand' to produce a new array of individual values that have been integrated.

import numpy as np
from scipy import integrate

datafile = np.genfromtxt(" filename/flow_line_num9.txt",delimiter='\t',skiprows=1, dtype=float)

#The first four parameters are all arrays of numbers 
ptID = datafile[:,0] 
surf = datafile[:,9] #m
base = datafile[:,10] #m
H = surf - base #m

Bo = 2.207 #Pa yr^1/3
Ct = .16612 #K^k
Tr = 273.39 #K
k = 1.17
Ts = -19.0 #celsius
Tb = -2.0 #celsius

@np.vectorize
def integrand(a,b):
    def flow_pramB(z):
        temp = []
        for i in range(0,len(ptID)):
            tempA = ((Ts-Tb)*pow((z/H[i]),.333333))+Tb
            temp.append(tempA)
        B = []
        for i in range(0,len(ptID)):
            Bpram = (Bo*np.exp((3155/temp[i]) - (Ct/(pow((Tr-temp[i]),k)))))
            B.append(Bpram)
        return B
    return integrate.quad(flow_pramB, a, b)

B = integrand(surf, base)

This code is simply one example of many many tries to get this to work. A solution or just letting me know that I need to try other modules to get this to work would be greatly appreciated!

Answer 1

The function flow_pramB is incorrectly defined. It should take a float and return a float, not a list as it does now. For instance, here is a possible approach (note that we use numpy arrays instead of for loops),

def flow_pramB(z):
    temp = ((Ts-Tb)*pow((z/H[:]),.333333))+Tb
    B = (Bo*np.exp((3155/temp[:]) - (Ct/(pow((Tr-temp[:]),k)))))
    return B.sum() # suming all the elements of B (probably not the right thing to do)

@np.vectorize
def integrand(a,b):
    return integrate.quad(flow_pramB, a, b)

B = integrand(surf, base)

this runs, but produces some NaNs and warning because flow_param does not produce the correct results (as I'm not sure what it should produce). Still this gives you an idea, of how to make this work,

1. Verify that flow_paramB takes a float and returns the correct float result.
2. Check that integrate.quad(flow_paramB, a_float, b_float) runs and produces the correct results.
3. Then wrap it with the np.vecorize decorator as you did.