被零除？（以牛顿迭代法）

Question

我正在进行牛顿迭代以找到 T_a。除了第一个定义之外，代码中的一切似乎都很好。我的 rho(T_a) returns 除以零（它假定 T_a 为零，而它只是一个变量。如果我将等式中的 T_a 更改为类似100，一切顺利。知道为什么它返回除以零吗？

from numpy import *
import numpy as np
import pylab
import scipy
from scipy.optimize import leastsq
from math import *
import matplotlib.pyplot as plt
from scipy import integrate

#     THETA NOTATION:
#pi/2: substellar point
#-pi/2: antistellar point
#0: terminators


#define constants used in equations:
alb = 0.2 #constant albedo
F = 866 #J/s*m**2
R = 287.0 #J/K*kg
U = 5.0 #m/s
C_p = 1000 #J/K*kg
C_d = 0.0015
p1 = 10**4
p2 = 10**5.0
p3 = 10**6.0 #Pa
sig = 5.67*10**-8.0 #J/s*m**2*K**4 #Stefan-Boltzmann cst


def rho(T_a):
    p1=10000.0
    R=287.0 #J/K*kg
    return (p1/(T_a*R))

def a(T_a):
    U = 5 #m/s
    C_p = 1000 #J/K*kg
    C_d = 0.0015
    return rho(T_a)*C_p*C_d*U




#################################################
##### PART 2 : check integrals equality
#################################################


#define the RHS and LHS of integral equality

def LHS(theta):
    return (1-alb)*F*np.sin(theta)*np.cos(theta)


#define the result of each integral

Left = integrate.quad(lambda theta: LHS(theta), 0, pi/2)[0]


#define a function 1-(result LHS/result RHS)  >>> We look for the zero of this


x0=130.0 #guess a value for T_a
#T_a = 131.0

#Python way of solving for the zero of the function

#Define T_g in function of T_a, have RHS(T_a) return T_g**4 etc, have             result_RHS(T_a) return int.RHS(T_a),
#have func(T_a) return result_LHS/result_RHS


def T_g(T_a,theta):
    return np.roots(array([(sig),0,0,(a(T_a)),((-a(T_a)*T_a)-LHS(theta))]))[3]

def RHS(theta,T_a):
    return sig*T_g(T_a,theta)**4*np.cos(theta)

def result_RHS(T_a,theta):
    return integrate.quad(lambda theta: RHS(T_a,theta), -pi/2, pi/2)[0]

def function(T_a,theta):
    return 1-((Left/result_RHS(T_a,theta)))

theta = np.arange(-pi/2, pi/2, pi/20)

T_a_0 = scipy.optimize.newton(function,x0,fprime=None,args=(theta,),tol=    (10**-5),maxiter=50000)

输出：

Traceback (most recent call last):
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 85, in <module>
    T_a_0 = scipy.optimize.newton(function,x0,fprime=None,args=(theta,),tol=(10**-5),maxiter=50000)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/optimize/zeros.py", line 120, in newton
    q0 = func(*((p0,) + args))
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 81, in function
    return 1-((Left/result_RHS(T_a,theta)))
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 78, in result_RHS
    return integrate.quad(lambda theta: RHS(T_a,theta), -pi/2, pi/2)[0]
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 312, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 78, in <lambda>
    return integrate.quad(lambda theta: RHS(T_a,theta), -pi/2, pi/2)[0]
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 75, in RHS
    return sig*T_g(T_a,theta)**4*np.cos(theta)
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 72, in T_g
    return np.roots(array([(sig),0,0,(a(T_a)),((-a(T_a)*T_a)-LHS(theta))]))[3]
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 38, in a
    return rho(T_a)*C_p*C_d*U
  File "/Users/jadecheclair/Documents/PHY479Y/FindT_a.py", line 32, in rho
    return (p1/(T_a*R))
ZeroDivisionError: float division by zero

Answer 1

您颠倒了参数：

在result_RHS中你调用：RHS(T_a,theta)，但是RHS的参数定义是def RHS(theta,T_a)

交换定义中的那些，错误不再发生。您的定义应如下所示：

def RHS(T_a, theta)

Answer 2

您的 RHS 函数的定义与所有其他函数略有不同，因为它的第一个参数是 theta，第二个参数是 T_a：

def RHS(theta,T_a):
    return sig*T_g(T_a,theta)**4*np.cos(theta)

我认为这就是为什么你在这里以错误的顺序传递参数：

lambda theta: RHS(T_a,theta)

按照正确的顺序获取它们，您应该没问题。

附带说明一下，您的某些导入看起来可能会导致奇怪的错误：

from numpy import *
from math import *

Numpy 和 math 模块至少有几个函数名相同，例如 sqrt。只做 import math 和 import numpy as np 并通过模块名称访问函数更安全。否则，当您调用 sqrt 时发生的情况可能会根据您执行导入的顺序而改变。

被零除？（以牛顿迭代法）

Division by zero ? (In newton iteration method)

python

optimization

zero

被零除 ？ （以牛顿迭代法）

Division by zero ? (In newton iteration method)

python

optimization

zero

被零除？（以牛顿迭代法）