使用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標誌，因爲這增加了複雜性。