क्या पाइथन में एक मल्टीपल इंटीग्रेटर है जो दोनों वैरिएबल इंटीग्रेशन लिमिट्स (जैसे स्कैपी) और हाई प्रिसिजन (mpmath की तरह) प्रदान करता है?
मैं चर एकीकरण के लिए एक चौगुनी एकीकरण के लिए घसीट क्वाड और न्यूक्ड का उपयोग कर सकता हूं। समस्या यह है कि उपयोग की गई डिफ़ॉल्ट परिशुद्धता एक त्रुटि उठाती है जब सहनशीलता का अनुरोध प्राप्त नहीं किया जा सकता है। Mpmath इंटीग्रेटर के साथ, मैं mp.dps = मनमानी सेटिंग के साथ किसी भी मनमानी परिशुद्धता को परिभाषित कर सकता हूं, लेकिन मैं यह नहीं देख सकता कि क्या और कैसे सीमाएं nquad की तरह परिवर्तनशील हो सकती हैं। Mpmath भी क्वाडल में गॉस-लीजेंड विधि के साथ एक बहुत तेज़ निष्पादन प्रदान करता है, जो अत्यधिक वांछनीय है, क्योंकि मेरा कार्य चिकनी है, लेकिन चार एकीकरण को पूरा करने के लिए डरपोक के साथ अत्यधिक मात्रा में समय लगता है। कृपया मदद करे। नीचे केवल एक सरल कार्य है जो मेरे लक्ष्य को विफल करता है:
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)
जवाब
ठीक है, मुझे जवाब में कुछ कहना है, टिप्पणी में कोड डालना मुश्किल है
एमपी गणित के साथ सरल अनुकूलन सरल नियमों का पालन करना है:
- y 2.0 बहुत महंगा है (लॉग, ऍक्स्प, ...), y * y के साथ बदलें
- y 2 अभी भी महंगा है, y * y के साथ बदलें
- गुणन योग की तुलना में बहुत अधिक महंगा है, x * y + y ** 2.0 को (x + y) * y से प्रतिस्थापित करें
- विभाजन गुणन की तुलना में अधिक महंगा है, y / 4 को 0.25 * y से बदलें
कोड, विन 10 x64, पायथन 3.8
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
इसके अलावा, ट्रिपल 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)
नीचे पूर्ण कोड है, जिसके लिए मूल प्रश्न उठाया गया था। इसमें चार एकीकरण करने के लिए स्किपी का उपयोग होता है:
# 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()
मैं मामले के आधार पर एप्सेल को पर्याप्त रूप से कम करने में रुचि रखता हूं। एप्सबस आम तौर पर अज्ञात एप्रीओरी होता है, इसलिए मैं समझता हूं कि आउटपुट को पकड़ने से बचने के लिए मुझे इसे बहुत कम करना चाहिए, इस मामले में यह एक कम्प्यूटेशनल आर्टिकैक्ट का परिचय देता है। जब मैं इसे कम करता हूं, तो एक त्रुटि चेतावनी उठाई जाती है कि राउंड-ऑफ त्रुटियां महत्वपूर्ण हैं और प्राप्त की जाने वाली वांछित सहिष्णुता के लिए कुल त्रुटि को कम करके आंका जा सकता है।
जबकि प्रश्न गति के बारे में नहीं है, उत्तरार्द्ध सटीकता और सहिष्णुता के बारे में पूछताछ से पहले एक चौगुनी एकीकरण के निष्पादन को व्यावहारिक बनाने के साथ जुड़ा हुआ है। गति का परीक्षण करने के लिए, मैंने सभी चार एप्सेल = 1e-02 को बढ़ाया (बढ़ाया), जिससे मूल कोड का समय घटकर 2:14 (घंटे) हो गया। फिर मैंने सेवेरिन के प्रति शक्तियों को सरल किया और कुछ संस्मरण लागू किए । इनका समय कम से कम 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]
हालांकि, आउटपुट एक विरूपण साक्ष्य है जो अपर्याप्त सहिष्णुता के कारण होता है। मैं उत्तरोत्तर निचले मानों के लिए एप्सेल सेट कर सकता हूं और देख सकता हूं कि क्या परिणाम उपलब्ध स्काइप सटीकता के साथ यथार्थवादी समय में वास्तविक मूल्य में परिवर्तित हो जाता है। आशा है कि यह मूल प्रश्न को बेहतर ढंग से प्रदर्शित करता है।