Python 中是否有多个积分器提供可变积分限制(如 scipy)和高精度(如 mpmath)? [英] Is there a multiple integrator in Python providing both variable integration limits (like scipy) and high precision (like mpmath)?

问题描述

我可以使用 scipy quad 和 nquad 进行涉及可变积分限制的四重积分.问题是当无法达到要求的容差时，使用的默认精度会引发错误.使用 mpmath 积分器，我可以通过设置 mp.dps = 任意来定义任意精度，但是我看不到限制是否以及如何像 nquad 一样变得可变.Mpmath 还在 quadgl 中使用 Gauss-Legendre 方法提供了非常快速的执行，这是非常可取的，因为我的函数是平滑的，但是使用 scipy 完成四个积分需要大量时间.请帮忙.下面只是一个没有达到我的目标的简单函数:

from datetime import datetime进口scipy从 scipy.special 导入 jn, jn_zeros将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt从 mpmath 导入 *从 mpmath 导入 mp从 numpy 导入 *从 scipy.optimize 导入 *# 设置精度mp.dps = 15#;mp.pretty = 真# 设置快捷方式，所以我们可以只写 exp() 而不是 mp.exp() 等.F = mp.mpfexp = mp.exp罪 = mp.sincos = mp.cosasin = mp.asinacos = mp.acossqrt = mp.sqrtpi = mp.pitan = mp.tan开始 = datetime.now()打印(开始)#optionsy={'limit':100,'epsabs':1.49e-1,'epsrel':1.49e-01}#optionsx={'limit':100, 'epsabs':1.49e-1, 'epsrel':1.49e-01}定义 f(x,y,z):返回 2*sqrt(1-x**2) + y**2.0 + z定义范围x(y，z):返回 [-1,1]定义范围(z):返回 [1,2]定义范围():返回 [2,3]定义结果():返回 quadgl(f, rangex, rangey, rangez)"#下面的作品:定义结果():返回 quadgl(f, [-1,1], [1,2], [2,3])"打印(结果())结束 = datetime.now()打印(结束开始)

解决方案

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

from datetime import datetime进口scipy将 numpy 导入为 np从 mpmath 导入 *从 mpmath 导入 mp从 numpy 导入 *# 设置精度mp.dps = 20#;mp.pretty = 真# 设置快捷方式，所以我们可以只写 exp() 而不是 mp.exp() 等.F = mp.mpfexp = mp.exp罪 = mp.sincos = mp.cosasin = mp.asinacos = mp.acossqrt = mp.sqrtpi = mp.pitan = mp.tan开始 = datetime.now()打印('开始:'，开始)定义 f3():定义 f2(x):def f1(x,y):定义 f(x,y,z):返回 1.0 + x*y + y**2.0 + 3.0*z返回 quadgl(f, [-1.0, 1], [1.2*x, 1.0], [y/4, x**2.0])返回 quadgl(f1, [-1, 1.0], [1.2*x, 1.0])返回 quadgl(f2, [-1.0, 1.0])打印('结果='，f3())结束 = datetime.now()print('持续时间以分钟为单位:',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进口scipy将 numpy 导入为 np从 mpmath 导入 *从 mpmath 导入 mp从 numpy 导入 *从 scipy 导入集成# 设置精度mp.dps = 15#;mp.pretty = 真# 设置快捷方式，所以我们可以只写 exp() 而不是 mp.exp() 等.F = mp.mpfexp = mp.exp罪 = mp.sincos = mp.cosasin = mp.asinacos = mp.acossqrt = mp.sqrtpi = mp.pitan = mp.tan开始 = datetime.now()打印('开始:'，开始)#需要集成的功能def f(x,y,z,w):返回 1.0 + x + y + z + w#Scipy 集成:第一个集成def f0(x,y,z):返回integrate.quad(f, -20, 10, args=(x,y,z), epsabs=1.49e-12, epsrel=1.4e-8)[0]#Mpmath 函数 f0(x,y,z) 积分器:三个外积分定义 f3():定义 f2(x):def f1(x,y):返回 quadgl(f0, [-1.0, 1], [-2, x], [-10, y])返回 quadgl(f1, [-1, 1.0], [-2, x])返回 quadgl(f2, [-1.0, 1.0])打印('结果='，f3())结束 = datetime.now()打印('持续时间:'，结束开始)

下面是完整的代码，针对它提出了最初的问题.包含使用scipy进行四次集成:

<预><代码># 进口从日期时间导入日期时间导入 scipy.integrate 作为 si进口scipy从 scipy.special 导入 jn, jn_zeros从 scipy.integrate 导入四边形从 scipy.integrate 导入 nquad将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt从 scipy.integrate 导入 fixed_quad从 scipy.integrate 进口正交从 mpmath 导入 mp从 numpy 导入 *从 scipy.optimize 导入 *# 设置精度mp.dps = 30# 设置快捷方式，所以我们可以只写 exp() 而不是 mp.exp() 等.F = mp.mpfexp = mp.exp罪 = mp.sincos = mp.cosasin = mp.asinacos = mp.acossqrt = mp.sqrtpi = mp.pitan = mp.tan开始 = datetime.now()打印(开始)R1 = F(6.37100000000000e6)k1 = F(8.56677817058932e-8)R2 = F(1.0)k2 = F(5.45789437248245e-01)r = F(12742000.0)#Replace 计算的初始常数与值假设更快，如下所示:#a2 = R2/r#打印(a2)a2 = F(0.0000000784806152880238581070475592529)定义 u1(phi2):返回 r*cos(phi2)-r*sqrt(a2**2.0-(sin(phi2))**2.0)定义 u2(phi2):返回 r*cos(phi2)+r*sqrt(a2**2.0-(sin(phi2))**2.0)def om(u,phi2):返回 u-r*cos(phi2)def mp2(phi2):返回 r*sin(phi2)定义 a1(u):返回 R1/uoptionsx={'limit':100,'epsabs':1.49e-14,'epsrel':1.49e-11}选项={'限制':100，'epsabs':1.49e-14，'epsrel':1.49e-10}#---- 在你的方向定义 a1b1_u(x,y,u):返回 2.0*u*sqrt(a1(u)**2.0-(sin(y))**2.0)def oa2_u(x,y,u,phi2):返回 (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,y,u,phi2):返回 (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 func1_u(x,y,u,phi2):返回 (-exp(-k1*a1b1_u(x,y,u)-k2*ob2_u(x,y,u,phi2))+exp(+k2*oa2_u(x,y,u,phi2)))*sin(y)*cos(y)#--------joint_coaxial 集成:u1def fg_u1(u,phi2):返回 nquad(func1_u, [[-pi, pi], [0, asin(a1(u))]], args=(u,phi2), opts=[optionsx,optionsy])[0]#Constants 用于在最后或中间 inegrals 中进行规范化，如果这有助于调整执行速度的值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))#----u1的三次积分def third_u1(u,phi2):返回 fg_u1(u,phi2)*u**2.0def third_u1_I(phi2):返回 quad(third_u1, u1(phi2), u2(phi2), args = (phi2), epsabs=1.49e-20, epsrel=1.49e-09)[0]#----u1的第四次积分定义第四_u1(phi2):返回third_u1_I(phi2)*sin(phi2)*cos(phi2)def force_u1():返回 quad(fourth_u1, 0.0, asin(a2), args = (), epsabs=1.49e-20, epsrel=1.49e-08)[0]force_u1 = force_u1()*r**2.0*2.0*pi*k2/piA1/piA2打印('r = ', r, 'force_u1 =', force_u1)结束 = datetime.now()打印(结束)参数 = {'p':r,'q':force_u1,'r':开始，'发送}#到txt文件f=open('Sphere-test-force-u-joint.txt', 'a')f.write('\n{p},{q},{r},{s}'.format(**args))#f.flush()f.close()

我有兴趣根据情况将 epsrel 设置得足够低.epsabs 通常是先验未知的，因此我知道我应该将其设置得非常低以避免它占用输出，在这种情况下它会引入计算问题.当我降低它时，会出现错误警告，指出舍入误差很大，并且可能低估了要实现的所需容差的总误差.

I can use scipy quad and nquad for a quadruple integration involving variable integration limits. The problem is that the default precision used raises an Error when the tolerance requested cannot be achieved. With mpmath integrator, I can define any arbitrary precision with setting mp.dps = arbitrary, but I can't see if and how the limits can become variable like with nquad. Mpmath also provides a very fast execution with Gauss-Legendre method in quadgl, which is highly desirable, because my function is smooth, but takes an exorbitant amount of time with scipy to complete four integrations. Please help. The below is only a simple function that fails my goal:

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, 'epsabs':1.49e-1, '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)

解决方案

Below is a simple example of how I can do only triple integration with mpmath. This does not address high precision with four integrations. In any case, execution time is even a bigger problem. Any help welcome.

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], [1.2*x, 1.0])
    return quadgl(f2, [-1.0, 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

Furthermore, an attempt to combine one (first) scipy integration followed by a triple mpmath integrator does not seem to produce any output for more than 24 hours even with a simplest function. What is wrong with the following code?

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, 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)

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


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

print('result =', f3())

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

Below is the full code, for which the original question was raised. It contains the use of scipy to carry out four integrations:

# 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, 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)

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, 'epsabs':1.49e-14, 'epsrel':1.49e-10}

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

def oa2_u(x,y,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,y,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 func1_u(x,y,u,phi2):
    return (-exp(-k1*a1b1_u(x,y,u)-k2*ob2_u(x,y,u,phi2))+exp(+k2*oa2_u(x,y,u,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 = (), epsabs=1.49e-20, 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()

I am interested in setting the epsrel sufficiently low, depending on the case. The epsabs is generally unknown apriori, so I understand that I should make it very low to avoid it taking hold of the output, in which case it introduces an computational articact. When I make it lower, an Error warning is raised that the round-off errors are significant and the total error may be underestimated for the desired tolerance to be achieved.

这篇关于Python 中是否有多个积分器提供可变积分限制(如 scipy)和高精度(如 mpmath)?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！