scipy 执行多重积分时,'采用0个位置参数,但给定了1个'

ajsxfq5m  于 2022-11-10  发布在  其他
关注(0)|答案(2)|浏览(164)

我正在做一个4个变量的多重积分,其中2个有极限函数。但是错误出现在我的一个常极限变量上。真的不知道为什么。非常感谢你的建议!

from numpy import sqrt, sin, cos, pi, arcsin, maximum
from sympy.functions.special.delta_functions import Heaviside
from scipy.integrate import nquad

def bmax(x):
    return 1.14*10**9/sin(x/2)**(9/7)

def thetal(x,y,z):
    return arcsin(3.7*10**15*sqrt(cos(x/2)**2/10**6-1.23*10**10/z+0.003*sin(x/2)**2*(2.51*10**63/sin(x/2)**9/y**7-1))/(z*sin(x/2)**2*cos(x/2)*(2.51*10**63/sin(x/2)**9/y**7-1)))

def rt(x,y):
    return 3.69*10**12/(2.5*10**63/sin(x/2)**7*y**7-sin(x/2)**2)

def rd(x,y):
    return maximum(1.23*10**10,rt(x,y))

def rl(x,y):
    return rd(x,y)*(sqrt(1+5.04*10**16/(rd(x,y)*cos(x/2)**2))-1)/2

def wbound():
    return [1.23*10**10,3.1*10**16]

def zbound():
    return [10**(-10),pi-10**(-10)]

def ybound(z):
    return [0,bmax(z)-10**(-10)]

def xbound(z,y,w):
    return [thetal(z,y,w),pi-thetal(z,y,w)]

def f(x,y,z,w):
    return [5.77/10**30*sin(z)*sin(z/2)*y*sin(x)*Heaviside(w-rl(z,y))*Heaviside(w-rd(z,y))/w**2]

result = nquad(f, [xbound, ybound,zbound,wbound])
xfyts7mz

xfyts7mz1#

出现这个错误的原因是,尽管你不希望这些绑定依赖于变量,但nquad仍然会将变量传递给你提供给它的函数,所以绑定函数必须接受正确数量的变量:

def wbound():
    return [1.23*10**10,3.1*10**16]

def zbound(w_foo):
    return [10**(-10),pi-10**(-10)]

def ybound(z, w_foo):
    return [0,bmax(z)-10**(-10)]

def xbound(z,y,w): 
    return [thetal(z,y,w),pi-thetal(z,y,w)]

现在,函数zboundybound接受额外的变量,但只是忽略它们。
我不确定最后一次绑定xbound(...):是否要翻转变量yz?根据scipy.integrate.nquad的定义,正确的顺序应该是

def xbound(y,z,w):
    ...

编辑:正如kazemakase所指出的,函数f应该返回一个float而不是一个列表,所以return语句中的括号[...]应该被删除。

wwwo4jvm

wwwo4jvm2#

nquad的第二个参数需要一个bounds序列,语法相当严格。
如果被积函数f依赖于x, y, z, w,并且这是定义的顺序,则bounds中的项必须按顺序为xbybzbwb,其中每个边界可以是2元组,例如xb = (xmin, xmax),或者返回2元组的函数。
关键点是,这些函数的参数......当我们执行内部积分时,在dx中,我们有yzw来计算x中的边界,所以它必须是
def xb(y,z,w): return(..., ...)-同样,def yb(z,w): return (..., ...)
def zb(w): return (..., ...) .
关于积分的最后一个变量的界限必须是常数。
总结一下


# DEFINITIONS

def f(x, y, z, w): return ..  . # x inner integration, ..., w outer integration
def xb(y,z,w): return (...,...) # or simply xb=(...,...) if it's a constant
def yb(z,w): return (...,...)   # or yb=(...,...)
def zb(w): return (...,...)     # or zb=(...,...)
wb = (...,...)

# INTEGRATION

result, _ = nquad(f, [xb, yb, zb, wb])

相关问题