使用 brentq 函数迭代 Python 中的数组

Question

我在使用 brentq 函数迭代数组的每个元素时遇到问题。下面定义的函数中的 q 是一个 FITS 文件数组，我们通过 brentq 函数使用此数组中的每个元素作为运行的输入，以便求解 T.

本质上，我的问题在于不是特别了解在哪里或如何实现适当的 for 循环来迭代 q.

的每个元素的函数

关于如何解决这个问题有什么建议吗？

def f(T,q,coeff1,coeff2,coeff3):
    return q*const3 - ((exp(const2/T)-1)/(exp(const/T)-1))

a = brentq(f, 10, 435.1, args=(q,4351.041,4262.570,0.206))
print a

newhdu = fits.PrimaryHDU(a)
newhdulist = fits.HDUList([newhdu])
newhdulist.writeto('Temp21DCOT.fits')

进一步解释：我正在尝试做的基础是使用brentq使用我们初始数组的强度值求解温度值（我们的 FITS 文件）。

这个方程是从普朗克方程的两个波长的比率推导出来的，所以q = B_1/B_2如果我们要符合物理学，q中的每个元素都是强度值。 brentq 将为 q 中的每个元素求解 T（温度）的解析不可解方程，并生成与 q 大小相同的新温度数组。换句话说，我正在尝试使用 Plank 方程求解 FITS 文件中每个像素的温度。

注意：我重新发布了这个以更有效地澄清问题。

Answer 1

您是否遇到迭代问题或效率问题？

这个迭代对我有用：

In [485]: from scipy import optimize

In [486]: def f(T,q,coeff1,coeff2,coeff3):
        return q*coeff3 - ((np.exp(coeff2/T)-1)/(np.exp(coeff1/T)-1))
        # corrected the coeff use

In [487]: q=np.linspace(1,3,10)
# q range chosen to avoid the different signs ValueError

In [488]: A=[optimize.brentq(f, 10, 435.1, args=(i,4351.041,4262.570,0.206),full_output=True) for i in q]

In [489]: A
Out[489]: 
[(55.99858839149886, <scipy.optimize.zeros.RootResults at 0xa861284c>),
 (64.14621536172528, <scipy.optimize.zeros.RootResults at 0xa861286c>),
 (72.98658083834341, <scipy.optimize.zeros.RootResults at 0xa861288c>),
 (82.75638321495505, <scipy.optimize.zeros.RootResults at 0xa86128ac>),
 (93.73016750496367, <scipy.optimize.zeros.RootResults at 0xa86128cc>),
 (106.25045004489637, <scipy.optimize.zeros.RootResults at 0xa86128ec>),
 (120.76612665626851, <scipy.optimize.zeros.RootResults at 0xa861290c>),
 (137.88917389176325, <scipy.optimize.zeros.RootResults at 0xa861292c>),
 (158.4854607193551, <scipy.optimize.zeros.RootResults at 0xa861294c>),
 (183.82941862839408, <scipy.optimize.zeros.RootResults at 0xa861296c>)]

In [490]: [a[1].iterations for a in A]
Out[490]: [8, 9, 10, 10, 10, 10, 10, 9, 8, 10]

在 brentq 文档中 f returns 一个值，一组 args。有一些求解器，例如 ode 求解器，可以让您定义一个函数，该函数采用矢量变量和 returns 匹配矢量导数。看起来这个根查找器不允许这样做。所以你不得不迭代 args 值，并解决每种情况。我将迭代写成列表理解。其他迭代格式也是可能的（for 循环等）。我们甚至可以将此 brentq 调用包装在一个可以通过 np.vectorize 传递的函数中。但这仍然是一次迭代，节省的时间很少。

有多种处理多维数组的方法。一个简单的方法是 flatten 输入，进行 1d 迭代，然后重塑结果。例如：

In [517]: q1=q.reshape(2,5)

In [518]: q1
Out[518]: 
array([[ 1.        ,  1.22222222,  1.44444444,  1.66666667,  1.88888889],
       [ 2.11111111,  2.33333333,  2.55555556,  2.77777778,  3.        ]])

In [519]: np.array([optimize.brentq(f, 10, 435.1, args=(i,4351.041,4262.570,0.206)) for i in q1.flat]).reshape(q1.shape)
Out[519]: 
array([[  55.99858839,   64.14621536,   72.98658084,   82.75638321,
          93.7301675 ],
       [ 106.25045004,  120.76612666,  137.88917389,  158.48546072,
         183.82941863]])

我放弃了 full_output 标志，因为这会增加复杂性。

使用 brentq 函数迭代 Python 中的数组

Iterating an Array in Python using the brentq Function

python

arrays

numpy

mathematical-optimization

scipy