Brent的方法被包装在while循环中

Question

我试图用python2来解决一个非常复杂的函数，它依赖于两个变量ne_b和T，使用Brent的方法和scipy。为了保证输入函数的准确性，我把它分解成几个函数，所有的变量都是ne_b，对于某个温度T，我成功地用Brent方法求主函数f的最小值

现在，我想找出，在f（ne_b）=0时，不同温度值下，ne_b的解。因此，我创建了一个超过1500<；t峎3<；25000的交互循环，并列出了完全相同的函数，只是为了避免冗余而重新命名，因此对于每个值t峎3，计算一个新的ne_b函数。在

整个代码：

import scipy.optimize as op
from scipy.optimize import fsolve 
from scipy.optimize import minimize_scalar
from scipy.optimize import bisect

#Defining constants for computation.
kB = 8.6173303e-5 # eV/K
k = 1.38065e-16 #dyne * cm / K
a = 7.566e-15 #erg cm^-3 K^-4
Rn = 13.598 #eV
X = 0.71
Y = 0.27
Z = 0.02
me = 9.110e-28 #grams
ma = 1.661e-24 #grams
C = 4.83e15 #Saha constant in cm^-3 K^-3/2
He_1 = 24.5873876 #Ionization potential of helium in eV.
He_2 = 54.4177630 #Ionization potential of ionized helium in eV.
Ca_1 = 6.113158 #ionization potential of neutral calcium in eV.
Ca_2 = 11.87172 #Ionization potential of ionized calcium.
T=[] #Temperature in kelvin. 
#U_11=0 #Partition function for neutral hydrogen.
U_21 = 1. #Partition function for ionized hydrogen.
U_12 = 1. #Partition function for neutral helium.
U_22 = 1.9999 #Partition function for ionized helium.
U_32 = 1. #Partition function of doubly ionized helium.
#U_120=0 #Partition function for neutral calcium.
#U_220=0 #Partition function for ionized calcium.
U_320 = 1. #Partition function for doubly ionized calcium.
P_gas = 100. #dyne cm^-2
A_1 = 1.00790 #Atomic mass of hydrogen.
A_2 = 4.00260 #Atomic mass of helium.
A_20 = 40.08 #Atomic mass of calcium. 
a_1 = float((X/A_1) / ((X/A_1) + (Y/A_2) + (Z/A_20))) 
a_2 = float((Y/A_2) / ((X/A_1) + (Y/A_2) + (Z/A_20)))
a_20 = float((Z/A_20) / ((X/A_1) + (Y/A_2) + (Z/A_20)))
mu_n = float(1. / ((X/A_1) + (Y/A_2) + (Z/A_20)))

#Defining arrays for problem 3 plots. 
N = [] #total number density of gas.
T_3 = [] #Temperature in kelvin.
ne_3 = [] #Brent-solved electron density.


t_3 = 1500 #Kelvin.
while t_3 <= 25000:
    xp20 = [100, 2520, 2800, 3150, 3600, 4200, 5040, 6300, 8400, 12600, 25200]
    fp11 = [1.9999, 1.9999, 1.9999, 1.9999, 1.9999, 1.9999, 1.9999, 1.9999, 1.9999, 2.0091, 2.3334]
    U_11b = np.interp(t_3, xp20, fp11)

    fp120 = [1, 1, 1.0023, 1.0069, 1.02329, 1.06660, 1.183, 1.5171, 2.9174, 21.4783, 172981.636]
    U_120b = np.interp(t_3, xp20, fp120)

    fp220 = [1.9999, 1.9999, 2.0045, 2.0137, 2.037, 2.0893, 2.208, 2.4604, 3.0409, 4.5499, 6.6834]
    U_220b = np.interp(t_3, xp20, fp220) #Partition function of ionized calcium.

    n_3 = (P_gas) / (k*t_3) #Finding total number density.
    N.append(n_3)

    #Neutral hydrogen.
    prod_11b = Rn / (kB * t_3)
    exp_11b = np.exp(-prod_11b)
    def Y_11b(ne_b):
            return (1/ne_b) * (C*t_3**1.5*(U_21 / U_11b)*exp_11b)
    def f_11b(ne_b):
            return (1 + Y_11b(ne_b))**-1

    #Ionized hydrogen.
    def f_21b(ne_b):
            return Y_11b(ne_b) / (1 + Y_11b(ne_b))

    #Neutral helium.
    prod_12b = He_1 / (kB * t_3)
    exp_12b = float(np.exp(-prod_12b))
    def Y_12b(ne_b):
            return (1/ne_b) * (C * t_3**1.5 * (U_22 / U_12) * exp_12b)

    #Ionized helium.
    prod_22b = He_2 / (kB * t_3)
    exp_22b = float(np.exp(-prod_22b))
    def Y_22b(ne_b):
            return (1/ne_b) * (C * t_3**1.5 * (U_32 / U_22) * exp_22b)

    #For neutral helium.
    def f_12b(ne_b):
            return (1 + Y_12b(ne_b) + Y_22b(ne_b)*Y_12b(ne_b))**-1

    #Final ionized helium expression.
    def f_22b(ne_b):
            return Y_12b(ne_b) * (1 + Y_12b(ne_b) + (Y_12b(ne_b) * Y_22b(ne_b)))**-1

    #Doubly ionized helium.
    def f_32b(ne_b):
            return (Y_22b(ne_b) * Y_12b(ne_b)) / (1 + Y_12b(ne_b) + (Y_22b(ne_b) * Y_12b(ne_b)))

    #Neutral calcium.
    prod_120b = Ca_1 / (kB * t_3)
    exp_120b = float(np.exp(-prod_120b))
    def Y_120b(ne_b):
            return (1/ne_b) * (C * t_3**1.5 * (U_220b / U_120b) * exp_120b)

    #Singly ionized calcium.
    prod_220b = Ca_2 / (kB * t_3)

    #For neutral calcium:
    def f_120b(ne_b):
            return (1 + Y_120b(ne_b) + (Y_120b(ne_b) * Y_220b(ne_b)))**-1

    #For ionized calcium.
    def f_220b(ne_b):
            return f_120b(ne_b) * Y_120b(ne_b)

    #Doubly ionized calcium.
    def f_320b(ne_b):
            return f_220b(ne_b) * Y_220b(ne_b)

    #Creating the master sum function.
    def sum(ne_b):
            return (a_1 * f_21b(ne_b)) + (a_2 * (f_22b(ne_b) + 2.*f_32b(ne_b))) + (a_20 * (f_220b(ne_b) + 2.*f_320b(ne_b)))

    #Creating the final function g(ne_b).
    def g(ne_b):
            return ((n_3 - ne_b) * sum(ne_b)) - ne_b

    #Attempting to minimize complicated function using Brent's method.
    #ne_bf = op.minimize_scalar(g, bounds=(1, 1e15), method='brent', tol=1e-08)
    ne_bf = op.brentq(g, -1e15, 1e15) #Just guessing at values, one near total suspected density, one opposite for bounds.
    ne_3.append(ne_bf)

    t_3 = t_3 + 50 #Steps of 50 kelvin.
    T_3.append(t_3)

print "End of code so far."

当我运行循环时，错误是：

Traceback (most recent call last):

File "hw6.py", line 444, in ne_bf = float(op.brentq(g, -1e15, 1e15)) #Just guessing at values, one near total suspected density, one opposite for bounds.

File"/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/zeros.py", line 415, in brentq r = _zeros._brentq(f,a,b,xtol,rtol,maxiter,args,full_output,disp)

File "hw6.py", line 440, in g return ((n_3 - ne_b) * sum(ne_b)) - ne_b

File "hw6.py", line 436, in sum return (a_1 * f_21b(ne_b)) + (a_2 * (f_22b(ne_b) + 2.*f_32b(ne_b))) + (a_20 * (f_220b(ne_b) + 2.*f_320b(ne_b)))

File "hw6.py", line 384, in f_21b return Y_11b(ne_b) / (1 + Y_11b(ne_b))

File "hw6.py", line 378, in Y_11b return (1/ne_b) * (C*t_3**1.5*(U_21 / U_11b)*exp_11b)

ZeroDivisionError: float division by zero

我需要为每个函数追加一个新的值。在

我已经测试了这些函数中的每一个随机浮点值，比如g（12.78432），它们分别计算。只有当我试着把每个函数归零时，我才得到这个错误。如有任何建议，我们将不胜感激。如果需要，我还可以提供更多信息，如果需要，我可以提供成功的代码（为t_3离散值运行的代码）。谢谢您！在

Answer 1

问题是你试图除以零。找出分母为零时每个函数返回的值，可能是None，float("Inf")，或{}。然后让每个函数在除法前检查分母是否为零。在