求解由 Python 中的函数定义的数组中每个元素的定积分

Question

我正在尝试整合数组的值。函数中使用的数组是“H”，积分的积分限制是“surf”和“base”。这些值是从 .txt 文件中读取的。 'def flow_pramB(z)' 定义要积分的函数，'def integrand(a,b)' 执行积分。然后我想调用函数“integrand”来生成一个新的已集成的单个值数组。

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)

此代码只是使其工作的许多尝试中的一个示例。非常感谢您提供解决方案或只是让我知道我需要尝试其他模块才能使其正常工作！

Answer 1

函数flow_pramB 定义不正确。它应该接受一个浮点数并返回一个浮点数，而不是像现在这样的列表。例如，这是一种可能的方法（注意我们使用 numpy 数组而不是 for 循环），

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)

这会运行，但会产生一些 NaNs 和警告，因为 flow_param 不会产生正确的结果（因为我不确定它应该产生什么）。这仍然让您了解如何使这项工作发挥作用，

1. 验证flow_paramB 采用浮点数并返回正确的浮点数结果。
2. 检查integrate.quad(flow_paramB, a_float, b_float) 是否运行并产生正确的结果。
3. 然后像你做的那样用np.vecorize 装饰器包装它。