布伦特的方法包裹在一个while循环中

Question

我正在尝试使用 python 2 来解决一个非常复杂的函数，该函数依赖于两个变量 ne_b 和 T，使用带有 scipy 的 Brent 方法。为了保证输入函数的准确性，我把它分成了几个函数，都是变量ne_b。对于某个温度 T，我已经成功地使用 Brent 方法找到了主函数 f 的最小值。

现在，我想在 f(ne_b) = 0 处找到 ne_b 的解，但是对于不同的温度值。所以我创建了一个超过 1500

整个代码：

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（最近一次通话最后一次）：

文件“hw6.py”，第 444 行，在 ne_bf = float(op.brentq(g, -1e15, 1e15)) #只​​是猜测值，一个接近总怀疑密度，一个相反的边界。

文件“/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/zeros.py”，第 415 行，在 brentq r = _zeros._brentq(f,a,b,xtol,rtol,maxiter,args,full_output,disp)

文件“hw6.py”，第 440 行，以 g 为单位 返回 ((n_3 - ne_b) * sum(ne_b)) - ne_b

文件“hw6.py”，第 436 行，总和 返回 (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)))

文件“hw6.py”，第 384 行，在 f_21b 中 返回 Y_11b(ne_b) / (1 + Y_11b(ne_b))

文件“hw6.py”，第 378 行，在 Y_11b 中 返回 (1/ne_b) * (C*t_3**1.5*(U_21 / U_11b)*exp_11b)

ZeroDivisionError：浮点除以零

我只需要将每个新函数归零，然后将该值附加到数组中以便稍后绘制。

我已经针对一些随机浮点值（例如 g(12.78432)）测试了这些函数中的每一个，它们分别计算。只有当我尝试将每个函数归零时才会出现此错误。任何建议表示赞赏。如果需要，我还可以提供更多信息，或者如果需要，我可以提供成功的代码（运行离散值 t_3 的代码）。谢谢！

Answer 1

问题是您试图除以零。弄清楚如果分母为零，您希望每个函数返回什么，可能是None、float("Inf") 或numpy.nan。然后让每个函数在尝试除法之前检查分母是否为零。。