มีตัวรวมหลายตัวใน Python ที่ให้ทั้งขีด จำกัด การรวมตัวแปร (เช่น scipy) และความแม่นยำสูง (เช่น mpmath) หรือไม่
ฉันสามารถใช้ scipy quad และ nquad สำหรับการรวมสี่เท่าที่เกี่ยวข้องกับขีด จำกัด การรวมตัวแปร ปัญหาคือความแม่นยำเริ่มต้นที่ใช้ทำให้เกิดข้อผิดพลาดเมื่อไม่สามารถบรรลุความคลาดเคลื่อนที่ร้องขอได้ ด้วยตัวรวม mpmath ฉันสามารถกำหนดความแม่นยำใด ๆ โดยพลการด้วยการตั้งค่า mp.dps = โดยพลการ แต่ฉันไม่สามารถดูได้ว่าขีด จำกัด สามารถเปลี่ยนแปลงได้อย่างไรและอย่างไรเช่นกับ nquad Mpmath ยังให้การดำเนินการที่รวดเร็วมากด้วยเมธอด Gauss-Legendre ในรูปสี่เหลี่ยมจัตุรัสซึ่งเป็นที่ต้องการอย่างมากเนื่องจากฟังก์ชันของฉันราบรื่น แต่ต้องใช้เวลามากเกินไปในการทำ scipy เพื่อดำเนินการรวมสี่อย่าง กรุณาช่วย. ด้านล่างนี้เป็นเพียงฟังก์ชันง่ายๆที่ทำให้เป้าหมายของฉันล้มเหลว:
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)
คำตอบ
โอเคให้ฉันใส่คำตอบยากที่จะใส่รหัสในความคิดเห็น
การเพิ่มประสิทธิภาพอย่างง่ายด้วย MP math คือการปฏิบัติตามกฎง่ายๆ:
- y 2.0มีราคาแพงมาก (log, exp, ... ) แทนที่ด้วย y * y
- y 2ยังคงแพงแทนที่ด้วย y * y
- การคูณมีราคาแพงกว่าการสรุปมากแทนที่ x * y + y ** 2.0 ด้วย (x + y) * y
- การหารมีราคาแพงกว่าการคูณแทนที่ y / 4 ด้วย 0.25 * y
รหัส, Win 10 x64, Python 3.8.0
def f3():
def f2(x):
def f1(x,y):
def f(x,y,z):
return 1.0 + (x+y)*y + 3.0*z
return mpmath.quadgl(f, [-1.0, 1], [1.2*x, 1.0], [0.25*y, x*x])
return mpmath.quadgl(f1, [-1, 1.0], [1.2*x, 1.0])
return mpmath.quadgl(f2, [-1.0, 1.0])
บนคอมพิวเตอร์ของฉันไปจาก 12.9 วินาทีเป็น 10.6 วินาทีลดลงประมาณ 20%
ด้านล่างนี้เป็นตัวอย่างง่ายๆของวิธีที่ฉันสามารถผสานรวมกับ 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], [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
นอกจากนี้ความพยายามที่จะรวมการรวม 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, 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)
ด้านล่างนี้คือรหัสแบบเต็มซึ่งเป็นคำถามเดิมที่เกิดขึ้น ประกอบด้วยการใช้ 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, 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()
ฉันสนใจที่จะตั้งค่า epsrel ให้ต่ำพอสมควรขึ้นอยู่กับกรณี โดยทั่วไปแล้ว epsabs นั้นไม่รู้จัก apriori ดังนั้นฉันจึงเข้าใจว่าฉันควรทำให้ต่ำมากเพื่อหลีกเลี่ยงไม่ให้มันค้างเอาท์พุทซึ่งในกรณีนี้จะแนะนำข้อต่อเชิงคำนวณ เมื่อฉันทำให้ต่ำลงจะมีการแจ้งเตือนข้อผิดพลาดขึ้นว่าข้อผิดพลาดในการปัดเศษมีความสำคัญและอาจมีการประเมินข้อผิดพลาดทั้งหมดต่ำเกินไปสำหรับค่าความคลาดเคลื่อนที่ต้องการ
ในขณะที่คำถามไม่ได้เกี่ยวกับความเร็ว แต่คำถามหลังมีความเชื่อมโยงอย่างใกล้ชิดกับการทำให้การดำเนินการรวมสี่เท่าในทางปฏิบัติก่อนที่จะมีการสอบถามเกี่ยวกับความแม่นยำและความอดทน ในการทดสอบความเร็วฉันตั้งค่า (เพิ่มขึ้น) 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, 'epsabs':1.49e-14, 'epsrel':1.49e-02}
optionsy={'limit':100, 'epsabs':1.49e-14, 'epsrel':1.49e-02}
def a1b1_u(x,y,u):
return 2.0*u*sqrt(a1(u)*a1(u)-sin(y)*sin(y))
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
+ 1.0-om(u,phi2)*om(u,phi2)-mp2(phi2)*mp2(phi2)))
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
+ 1.0-om(u,phi2)*om(u,phi2)-mp2(phi2)*mp2(phi2)))
def third_u1(u,phi2):
return fg_u1(u,phi2)*u*u
def third_u1_I(phi2):
return quad(third_u1, u1(phi2), u2(phi2), args = (phi2), epsabs=1.49e-20, epsrel=1.49e-02)[0]
def force_u1():
return quad(fourth_u1, 0.0, asin(a2), args = (), epsabs=1.49e-20, epsrel=1.49e-02)[0]
อย่างไรก็ตามผลลัพธ์เป็นสิ่งประดิษฐ์ที่เกิดจากความอดทนไม่เพียงพอที่นำมาใช้ ฉันสามารถตั้งค่า epsrel เป็นค่าที่ต่ำลงเรื่อย ๆ และดูว่าผลลัพธ์จะแปลงเป็นค่าจริงในเวลาจริงหรือไม่ด้วยความแม่นยำในการสแกนที่มีอยู่ หวังว่านี่จะช่วยอธิบายคำถามเดิมได้ดีขึ้นมาก