如何解决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 2 仍然很昂贵，请替换为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精度收敛到实际值。希望这能更好地说明原始问题。
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
