我有一段代码正试图从IDL复制到Python中。我已经复制了它的大部分,并且它可以工作,除了当我试图通过broyden1非线性解算器运行我写的函数(称为broyfunc)时,我得到了一个错误消息:我不能有一个包含infs或nan的数组。我也不知道怎么检查这个。我已经包括了我的错误信息和代码。谢谢您!在
import numpy as np
import scipy
import sympy
import sympy.mpmath as mpmath
# Inputs needed for AN_RC_MOD
p0 = 1.5 # reference pressure (surface pressure) [bar]
T0 = 94 # temperature at reference pressure [K]
Ttoa = 175 # temperature at top of atmosphere (assumed as stratopause) [K]
F1 = 1.5 # flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2]
F2 = 1.1 # flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2]
F20 = 0.75 # flux absorbed in channel '2' down to reference level p0 [W/m**2]
Fi = 0 # internal heat flux [W/m**2]
gam = 1.4 # ratio of specific heats, cp/cv (see Eq. 8 of RC12)
alpha = 0.77 # empirical adjustment to dry adiabatic lapse rate (see Eq. 10 of RC12)
n = 4/3 # scaling parameter that relates tau ~ p^n, where tau is gray thermal optical depth (see Eq.6 of RC12)
beta = alpha*((gam-1)/gam)
# Outputs:
# taurc - gray thermal optical depth of radiative-convective boundary
# tau0 - gray thermal optical depth down to reference pressure p0
# k1 - attenuation parameter for channel '1' [dimensionless]
# k2 - attenuation parameter for channel '2' [dimensionless]
# Important constants
sigma = 5.67e-8 # Stefan-Boltzmann constant [J/s/m**2/K**4]
D = 1.66 # diffusivity factory, (see, e.g., Rodgers & Walshaw 1966)
# Need inputs for: (numbers are for Titan, from Example_w_fluxes.pro)
#F1C = F1 ;flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2]
#F20C = F20 ;flux absorbed in channel '2' down to reference level p0 [W/m**2]
#F2C = F2 ;flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2]
#FiC = Fi ;internal heat flux [W/m**2]
#T0C = T0 ;temperature at reference pressure [K]
#p0C = p0 ;reference pressure (e.g., surface pressure) [bar]
#TtoaC = Ttoa ;temperature at top of atmosphere [K]
#nC = n ;scaling power between gray optical depth and pressure
#betaC = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12
# Need inputs for: (numbers are for Titan, from Example_w_fluxes.pro)
#F1C = F1 ;flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2]
#F20C = F20 ;flux absorbed in channel '2' down to reference level p0 [W/m**2]
#F2C = F2 ;flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2]
#FiC = Fi ;internal heat flux [W/m**2]
#T0C = T0 ;temperature at reference pressure [K]
#p0C = p0 ;reference pressure (e.g., surface pressure) [bar]
#TtoaC = Ttoa ;temperature at top of atmosphere [K]
#nC = n ;scaling power between gray optical depth and pressure
#betaC = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12
# x[0] = optical depth of Radiative-Convective boundary, taurc
# x[1] = total thermal optical depth down to reference pressure p0 (tau0)
# initial guess that R-C boundary is at D*taurc~1, and tau0~10*p0^n
x = [1/D, 10*p0**n] # initial guesses of taurc and tau0, respectively
print x
# Define function F for Broyden optimization
def broyfunc(x):
# Variables
taurc = x[0] # current value of optical depth of r-c boundary
tau0 = x[1] # current value of total thermal optical depth down to reference pressure p0
# Compute fluxes and temperatures at r-c boundary ;;;
# first check if solver is attempting to use negative values for the
# optical depths (which is unphysical), and, if so, return large numbers
# to the solver to prevent it from trying such values as solutions
if x[0] < 0 or x[1] < 0:
x = [1.e3,1.e3]
print x
else:
x = x
# Shortwave optical depth in channel "2" is given by
# LOG( F2/(F2-F20) ), and k2 is defined as the ratio of
# this to the thermal optical depth at p0 (i.e., tau0)
if F2 > 0:
k2 = math.log(F2/(F2-F20))/tau0
else:
# if the flux in channel "2" is zero,
# then no attenuation in this channel
k2 = 0
# k1 is computed by requiring the radiative temperature profile
# to equal the top-of-atmosphere temperature (Ttoa, set by user),
# which is accomplished by evaluating Eq. 18 in RC12 at tau=0
if F1 > 0:
k1 = D*( 2/F1*( sigma*Ttoa**4 - Fi/2 - F2/2*(1 + k2/D) ) - 1 )
else:
# if the flux in channel "1" is zero,
# then no attenuation in this channel
k1 = 0
# return values of k1, k2 to user
k1C = k1 # attenuation parameter for channel '1' (assumed as "stratosphere")
k2C = k2 # attenuation parameter for channel '2' (assumed as "troposphere")
print 'k1C', k1C
print 'k2C', k2C
# four cases: (1) both k1, k2 are greater than zero
# (2) k1 greater than zero, k2 equal to (or less than) zero
# (3) k2 greater than zero, k1 equal to (or less than) zero
# (4) both k1, k2, are equal to (or less than) zero
# case 1: both k1, k2 are greater than zero #
if k1 > 0 and k2 > 0:
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = (F1/2)*(1 + (D/k1) + (1 - D/k1)*math.exp(-k1*taurc)) + (F2/2)*(1 + (D/k2) + (1 - D/k2)*math.exp(-k2*taurc)) + (Fi/2)*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D/k1 + (k1/D - D/k1)*math.exp(-k1*taurc)) + F2/2*(1 + D/k2 + (k2/D - D/k2)*math.exp(-k2*taurc)) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# case 2: k1 greater than zero, k2 equal to (or less than) zero #
if k1 > 0 and k2 < 0:
k2C = 0
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = F1/2*(1 + D/k1 + (1 - D/k1)*math.exp(-k1*taurc)) + F2/2*(2 + D*taurc) + Fi/2*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D/k1 + (k1/D - D/k1)*math.exp(-k1*taurc)) + F2/2*(1 + D*taurc) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# case 3: k2 greater than zero, k1 equal to (or less than) zero #
if k1 < 0 and k2 > 0:
k1C = 0
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = F1/2*(2 + D*taurc) + F2/2*(1 + D/k2 + (1 - D/k2)*math.exp(-k2*taurc)) + Fi/2*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D*taurc) + F2/2*(1 + D/k2 + (k2/D - D/k2)*math.exp(-k2*taurc)) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# case 4: both k1, k2, are equal to (or less than) zero
if k1 < 0 and k2 < 0:
k1C = 0
k2C = 0
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = F1/2*(2 + D*taurc) + F2/2*(2 + D*taurc) + Fi/2*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D*taurc) + F2/2*(1 + D*taurc) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# Upwelling thermal flux at R-C boundary from convective region (RC12 Eq. 13)
#Fup_c = sigma*T0**4*math.exp(D*taurc)*((math.exp(-D*tau0)) + (1/((D*tau0)**(4*beta/n)))*(mpmath.gammainc(1+(4*beta/n),D*tau0) - mpmath.gammainc(1+(4*beta/n),D*taurc)))
#Fup_c = (sigma*T0**4)*math.exp(D*taurc)*(math.exp((-D*tau0) + 1/(D*tau0)**(4*beta/n)*mpmath.gammainc(1+4*beta/n)*(mpmath.gammainc(1+4*beta/n, D*tau0) - mpmath.gammainc(1+4*beta/n,D*taurc))))
a1 = D*taurc
a2 = -D*tau0
a3 = D*taurc
a4 = D*tau0
b1 = (4*beta)/n
b2 = b1 + 1
Fup_c1 = (sigma*T0**4)*(math.exp(a1))
Fup_c2 = math.exp(a2)
Fup_c3 = 1/(a4**b1)
Fup_c4 = mpmath.gammainc(b2, a3)
Fup_c5 = mpmath.gammainc(b2, a4)
Fup_c = Fup_c1*(Fup_c2 + (Fup_c3*(Fup_c4 - Fup_c5)))
# temperature at R-C boundary from convective region (RC12 Eq. 11)
T_c = T0*(taurc/tau0)**(beta/n)
print 'D*taurc', D*taurc
print 'Fup_r', Fup_r
print 'Fup_c', Fup_c
print 'T_r', T_r
print 'T_c', T_c
# system is solved when the temperatures and upwelling fluxes are continuous
# at the R-C boundary, so the solver attempts to minimize the difference
# between these quantities as computed from the radiative expressions and the
# convective expressions
#result = np.dtype('f8')
result = [T_r - T_c, Fup_r - Fup_c]
print result
# essentially we solve simultaneous equations:
# T_r - T_c = 0 and Fup_r - Fup_c = 0
broyfunc(x)
#y = np.dtype('f8')
y = scipy.optimize.broyden1(broyfunc, x)
ValueError回溯(最近一次调用) 在() 183 y=np.D类型(“f8”) 184打印x -->;185 y=scipy.optimize.broyden1(布罗伊芬克,x,iter=10) 186个
//Python/lib/python2.7/site-packages/scipy/optimize/非林.pyc在broyden1中(F,xin,iter,alpha,reduction_method,max_rank,verbose,maxiter,F_tol,F_rtol,x_rtol,tol_norm,line_search,callback,**kw)
//Python/lib/python2.7/site-packages/scipy/optimize/非林.pyc在nonlin_solve中(F,x0,jacobian,iter,verbose,maxiter,F_tol,F_rtol,x_rtol,tol_norm,line_search,回调,完整输出,引发异常) 279 dx=np.inf公司 280 Fx=功能(x) -->;281 Fx_norm=标准值(Fx) 282 283雅可比=雅可比(雅可比)
//Python/lib/python2.7/site-packages/scipy/linalg/杂项pyc正常(a,ord) 115#与numpy的区别仅在于非有限处理和 116#布拉斯 --&117加仑np.asarray chkfinite公司(一) 118如果ord在(None,2)和(a.ndim==1)和(a.ndim==1)中。数据类型.char在“fdFD”)中: 119#使用blas实现快速稳定的欧几里德范数
//Python/lib/python2.7/site-packages/numpy/lib/function_基.pyc在asarray chkfinite中(a,dtype,order) 588如果a。数据类型.char在typecodes['AllFloat']和notnp.是有限的(a) .all(): 589提升值错误( -->;590“数组不能包含inf或nan”) 591回路a 592个
ValueError:数组不能包含inf或nan
最直接的问题是您没有从
broyfunc
返回任何内容:所以布罗登积分器所看到的是
^{pr2}$(不可否认,这是一条可怕的错误信息。)
在添加
return result
之后,我得到了类似于但这是一个单独的问题,需要真正思考代码应该做什么。;—)
还有一个值得警告的问题:你写作
但您使用的是Python 2,因此有整数除法:
要么升级到python3(更好),通过添加一个小数点(
4./3
)来确保所有内容都是一个浮点数,或者,如果您坚持使用python2并且正在编写数字代码,我建议您将from __future__ import division
添加到程序的开头,以便相关问题 更多 >
编程相关推荐