使用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中的每个元素都是强度值，则为q = B_1/B_2。 brentq将为q中的每个元素解决T（温度）的分析上不可解的方程，并制作与q相同尺寸的新温度阵列。换句话说，我试图使用普朗克方程来解决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返回一个值，一组args。有一些求解器（如ode）可让您定义一个函数，该函数使用矢量变量并返回匹配的矢量导数。它看起来不像这个根发现者允许的。所以你坚持迭代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标志，因为这增加了复杂性。