Python中是否有多个集成商提供可变的集成限制例如scipy和高精度例如mpmath？

问题描述

我可以将scipy quad和nquad用于涉及可变积分限制的四重积分。问题是，当无法达到要求的公差时，使用的默认精度会引发错误。使用mpmath积分器，我可以通过设置mp.dps =任意来定义任意精度，但是我看不到限制是否以及如何像nquad一样变得可变。 Mpmath还可以通过Quadgl中的Gauss-Legendre方法提供非常快速的执行，这是非常可取的，因为我的函数很流畅，但是要花费大量的时间才能完成四个集成。请帮忙。下面只是一个无法实现我的目标的简单函数：

from datetime import datetime
import scipy
from scipy.special import jn,jn_zeros
import numpy as np
import matplotlib.pyplot as plt
from mpmath import *
from mpmath import mp
from numpy import *
from scipy.optimize import *

# Set the precision
mp.dps = 15#; mp.pretty = True

# Setup shortcuts,so we can just write exp() instead of mp.exp(),etc.
F = mp.mpf
exp = mp.exp
sin = mp.sin
cos = mp.cos
asin = mp.asin
acos = mp.acos
sqrt = mp.sqrt
pi = mp.pi
tan = mp.tan

start = datetime.Now()
print(start)

#optionsy={'limit':100,'epsabs':1.49e-1,'epsrel':1.49e-01}
#optionsx={'limit':100,'epsrel':1.49e-01}

def f(x,y,z):
    return 2*sqrt(1-x**2) + y**2.0 + z

def rangex(y,z):
    return [-1,1]

def rangey(z):
    return [1,2]

def rangez():
    return [2,3]


def result():
    return quadgl(f,rangex,rangey,rangez)

"""
#The below works:

def result():
    return quadgl(f,[-1,1],[1,2],[2,3])
"""

print(result())

end = datetime.Now()
print(end-start)

解决方法

下面是一个简单的示例，说明我如何仅使用mpmath进行三重集成。这不能解决四个集成的高精度问题。无论如何，执行时间甚至是一个更大的问题。欢迎任何帮助。

from datetime import datetime
import scipy
import numpy as np
from mpmath import *
from mpmath import mp
from numpy import *

# Set the precision
mp.dps = 20#; mp.pretty = True

# Setup shortcuts,so we can just write exp() instead of mp.exp(),etc.
F = mp.mpf
exp = mp.exp
sin = mp.sin
cos = mp.cos
asin = mp.asin
acos = mp.acos
sqrt = mp.sqrt
pi = mp.pi
tan = mp.tan

start = datetime.now()
print('start: ',start)

def f3():
    def f2(x):
        def f1(x,y):
            def f(x,y,z):
                return 1.0 + x*y + y**2.0 + 3.0*z
            return quadgl(f,[-1.0,1],[1.2*x,1.0],[y/4,x**2.0])
        return quadgl(f1,[-1,1.0])
    return quadgl(f2,1.0])

print('result =',f3())

end = datetime.now()
print('duration in mins:',end-start)

#start:  2020-08-19 17:05:06.984375
#result = 5.0122222222222221749
#duration: 0:01:35.275956

此外，尝试将一个（第一个）scipy积分与一个三重mpmath积分器相结合的尝试，即使使用最简单的功能，似乎也不会产生超过24小时的任何输出。以下代码有什么问题？

from datetime import datetime
import scipy
import numpy as np
from mpmath import *
from mpmath import mp
from numpy import *

from scipy import integrate

# Set the precision
mp.dps = 15#; mp.pretty = True

# Setup shortcuts,start)

#Function to be integrated
def f(x,z,w):
    return 1.0 + x + y + z + w 
    
#Scipy integration:FIRST INTEGRAL
def f0(x,z):
    return integrate.quad(f,-20,10,args=(x,z),epsabs=1.49e-12,epsrel=1.4e-8)[0]


#Mpmath integrator of function f0(x,z): THREE OUTER INTEGRALS
def f3():
    def f2(x):
        def f1(x,y):
            return quadgl(f0,[-2,x],[-10,y])
        return quadgl(f1,x])
    return quadgl(f2,f3())

end = datetime.now()
print('duration:',end-start)

下面是完整的代码，提出了原始问题。它包含使用scipy进行四个集成：


# Imports
from datetime import datetime
import scipy.integrate as si
import scipy
from scipy.special import jn,jn_zeros
from scipy.integrate import quad
from scipy.integrate import nquad
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import fixed_quad
from scipy.integrate import quadrature
from mpmath import mp

from numpy import *
from scipy.optimize import *

# Set the precision
mp.dps = 30

# Setup shortcuts,etc.
F = mp.mpf
exp = mp.exp
sin = mp.sin
cos = mp.cos
asin = mp.asin
acos = mp.acos
sqrt = mp.sqrt
pi = mp.pi
tan = mp.tan

start = datetime.now()
print(start)

R1 = F(6.37100000000000e6)
k1 = F(8.56677817058932e-8)
R2 = F(1.0)
k2 = F(5.45789437248245e-01)
r = F(12742000.0)

#Replace computed initial constants with values presuming is is faster,like below:
#a2 = R2/r
#print(a2) 
a2 = F(0.0000000784806152880238581070475592529)

def u1(phi2):
    return r*cos(phi2)-r*sqrt(a2**2.0-(sin(phi2))**2.0)
def u2(phi2):
    return r*cos(phi2)+r*sqrt(a2**2.0-(sin(phi2))**2.0)

def om(u,phi2):
    return u-r*cos(phi2)
def mp2(phi2):
    return r*sin(phi2)

def a1(u):
    return R1/u

optionsx={'limit':100,'epsabs':1.49e-14,'epsrel':1.49e-11}
optionsy={'limit':100,'epsrel':1.49e-10}

#---- in direction u
def a1b1_u(x,u):
    return 2.0*u*sqrt(a1(u)**2.0-(sin(y))**2.0)

def oa2_u(x,u,phi2):
    return (mp2(phi2)*sin(y)*cos(x)+om(u,phi2)*cos(y) 
                    - sqrt((mp2(phi2)*sin(y)*cos(x)+om(u,phi2)*(cos(y)))**2.0 
                           + R2**2.0-om(u,phi2)**2.0-mp2(phi2)**2.0))

def ob2_u(x,phi2)*cos(y) 
                    + sqrt((mp2(phi2)*sin(y)*cos(x)+om(u,phi2)**2.0-mp2(phi2)**2.0))

def func1_u(x,phi2):
    return (-exp(-k1*a1b1_u(x,u)-k2*ob2_u(x,phi2))+exp(+k2*oa2_u(x,phi2)))*sin(y)*cos(y)
 
#--------joint_coaxial integration: u1
def fg_u1(u,phi2):
    return nquad(func1_u,[[-pi,pi],[0,asin(a1(u))]],args=(u,phi2),opts=[optionsx,optionsy])[0]

#Constants to be used for normalization at the end or in the interim inegrals if this helps adjust values for speed of execution
piA1 = pi*(R1**2.0-1.0/(2.0*k1**2.0)+exp(-2.0*k1*R1)*(2.0*k1*R1+1.0)/(2.0*k1**2.0))
piA2 = pi*(R2**2.0-1.0/(2.0*k2**2.0)+exp(-2.0*k2*R2)*(2.0*k2*R2+1.0)/(2.0*k2**2.0))

#----THIRD integral of u1
def third_u1(u,phi2):
    return fg_u1(u,phi2)*u**2.0
def third_u1_I(phi2):
    return quad(third_u1,u1(phi2),u2(phi2),args = (phi2),epsabs=1.49e-20,epsrel=1.49e-09)[0]
    
#----FOURTH integral of u1
def fourth_u1(phi2):
    return third_u1_I(phi2)*sin(phi2)*cos(phi2)
def force_u1():
    return quad(fourth_u1,0.0,asin(a2),args = (),epsrel=1.49e-08)[0]


force_u1 = force_u1()*r**2.0*2.0*pi*k2/piA1/piA2

print('r = ',r,'force_u1 =',force_u1)

end = datetime.now()
print(end)

args = {
            'p':r,'q':force_u1,'r':start,'s':end
        }   

#to txt file
f=open('Sphere-test-force-u-joint.txt','a')

f.write('\n{p},{q},{r},{s}'.format(**args))
#f.flush()
f.close()

我想根据情况将epsrel设置得足够低。 epsabs通常是先验未知的，因此我知道我应该将其降低以免它占用输出，在这种情况下，它会引入计算假象。当我降低该值时，会出现错误警告，指出舍入误差很大，总误差可能被低估了，以实现所需的公差。

好的，我来回答一下，很难在注释中添加代码

MP数学的简单优化遵循以下简单规则：

y ^2.0非常昂贵（log，exp，...），请替换为y * y
y ²仍然很昂贵，请替换为y * y
乘法比求和要贵得多，用（x + y）* y代替x * y + y ** 2.0
除法比乘法更昂贵，将y / 4替换为0.25 * y

代码，Win 10 x64，Python 3.8

def f3():
    def f2(x):
        def f1(x,z):
                return 1.0 + (x+y)*y + 3.0*z
            return mpmath.quadgl(f,[0.25*y,x*x])
        return mpmath.quadgl(f1,1.0])
    return mpmath.quadgl(f2,1.0])

我的计算机上的时间从12.9秒减少到10.6秒，大约降低了20％

虽然问题不关乎速度，但后者与在询问精度和公差之前实际执行四重积分紧密相关。为了测试速度，我设置（增加了）所有四个epsrel = 1e-02，这将原始代码的时间减少到2:14（小时）。然后，我简化了每个Severin的功能，并实现了一些memoization。这些将时间累计减少到1:29（小时）。此处提供了代码的编辑行：

from memoization import cached

@cached(ttl=10)
def u1(phi2):
    return r*cos(phi2)-r*sqrt(a2*a2-sin(phi2)*sin(phi2))
@cached(ttl=10)
def u2(phi2):
    return r*cos(phi2)+r*sqrt(a2*a2-sin(phi2)*sin(phi2))
@cached(ttl=10)
def om(u,phi2):
    return u-r*cos(phi2)
@cached(ttl=10)
def mp2(phi2):
    return r*sin(phi2)
@cached(ttl=10)
def a1(u):
    return R1/u

optionsx={'limit':100,'epsrel':1.49e-02}
optionsy={'limit':100,'epsrel':1.49e-02}

def a1b1_u(x,u):
    return 2.0*u*sqrt(a1(u)*a1(u)-sin(y)*sin(y))

def oa2_u(x,phi2)*(cos(y)))**2.0 
                           + 1.0-om(u,phi2)*om(u,phi2)-mp2(phi2)*mp2(phi2)))

def ob2_u(x,phi2)-mp2(phi2)*mp2(phi2)))

def third_u1(u,phi2)*u*u

def third_u1_I(phi2):
    return quad(third_u1,epsrel=1.49e-02)[0]
    
def force_u1():
    return quad(fourth_u1,epsrel=1.49e-02)[0]

但是，输出是由于引入的公差不足引起的伪像。我可以逐步将epsrel设置为较低的值，并查看结果是否在实际时间内以可用的scipy精度收敛到实际值。希望这能更好地说明原始问题。

mpmath numpy python scipy scipy