Brent-Methode in eine While-Schleife gewickelt

Question

Ich versuche, Python 2 zu verwenden, um eine sehr komplizierte Funktion zu lösen, die auf zwei Variablen, ne_b und T, mit Brent-Methode mit scipy beruht. Um die Genauigkeit der Eingabe der Funktion zu gewährleisten, habe ich sie in mehrere Funktionen aufgeteilt, alle Variablen ne_b. Für eine bestimmte Temperatur T habe ich erfolgreich die Methode von Brent benutzt, um das Minimum der Hauptfunktion f zu finden.

Nun möchte ich die Lösungen von ne_b bei f (ne_b) = 0 finden, aber für unterschiedliche Temperaturwerte. Also habe ich eine interative Schleife über 1500 < t_3 < 25000 erstellt, und genau die gleichen Funktionen aufgeführt, nur umbenannt, um Redundanz zu vermeiden, so dass für jeden Wert t_3 eine neue Funktion von ne_b berechnet wird.

Der gesamte Code:

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." 

Wenn ich die Schleife ausgeführt, der Fehler ist:

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

Ich brauche nur jede neue Funktion auf Null und diesen Wert in einem Array append zum Plotten später.

Ich habe jede dieser Funktionen für einige zufällige Float-Wert wie g (12.78432) getestet, und sie berechnen separat. Nur wenn ich versuche, jede Funktion auf Null zu setzen, erhalte ich diesen Fehler. Jeder Rat wird geschätzt. Ich kann auch weitere Informationen bereitstellen, wenn angefordert, oder den erfolgreichen Code (der, der für einen diskreten Wert von t_3 ausgeführt wurde), wenn angefordert. Vielen Dank!

Answer 1

Das Problem ist, dass Sie versuchen, durch Null zu teilen. Finde heraus, was jede Funktion zurückgeben soll, wenn der Nenner Null ist, wahrscheinlich None, float("Inf"), oder numpy.nan. Dann muss jede Funktion prüfen, ob der Nenner 0 ist vor es Division versucht.