Có trình tích hợp nhiều trong Python cung cấp cả giới hạn tích hợp biến (như scipy) và độ chính xác cao (như mpmath) không?
Tôi có thể sử dụng scipy quad và nquad để tích hợp bốn liên quan đến các giới hạn tích hợp biến. Vấn đề là độ chính xác mặc định được sử dụng gây ra Lỗi khi không thể đạt được dung sai được yêu cầu. Với công cụ tích phân mpmath, tôi có thể xác định bất kỳ độ chính xác tùy ý nào bằng cách đặt mp.dps = tùy ý, nhưng tôi không thể biết liệu các giới hạn có thể biến đổi như thế nào với nquad hay không. Mpmath cũng cung cấp một quá trình thực thi rất nhanh với phương pháp Gauss-Legendre trong quadgl, điều này rất đáng mong đợi, bởi vì hàm của tôi rất mượt, nhưng mất một khoảng thời gian cắt cổ với scipy để hoàn thành bốn tích hợp. Hãy giúp tôi. Dưới đây chỉ là một chức năng đơn giản không đạt được mục tiêu của tôi:
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)
Trả lời
Được rồi, hãy để tôi giải đáp thắc mắc, khó điền mã vào phần nhận xét
Tối ưu hóa đơn giản với toán học MP là tuân theo các quy tắc đơn giản:
- y 2.0 là RẤT đắt (log, exp, ...), hãy thay thế bằng y * y
- y 2 vẫn đắt, thay bằng y * y
- phép nhân đắt hơn nhiều so với phép tính tổng, hãy thay x * y + y ** 2.0 bằng (x + y) * y
- Phép chia đắt hơn phép nhân, hãy thay y / 4 bằng 0,25 * y
Code, Win 10 x64, Python 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])
trên máy tính của tôi đã tăng từ 12,9 giây xuống 10,6 giây, giảm khoảng 20%
Dưới đây là một ví dụ đơn giản về cách tôi có thể thực hiện tích hợp ba lần với mpmath. Điều này không giải quyết được độ chính xác cao với bốn tích hợp. Trong mọi trường hợp, thời gian thực hiện thậm chí còn là một vấn đề lớn hơn. Mọi sự giúp đỡ hoan nghênh.
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
Hơn nữa, nỗ lực kết hợp một tích hợp scipy (đầu tiên) sau đó là tích hợp ba mpmath dường như không tạo ra bất kỳ đầu ra nào trong hơn 24 giờ ngay cả với một chức năng đơn giản nhất. Điều gì sai với đoạn mã sau?
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)
Dưới đây là mã đầy đủ, mà câu hỏi ban đầu đã được đưa ra. Nó bao gồm việc sử dụng scipy để thực hiện bốn tích hợp:
# 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()
Tôi quan tâm đến việc đặt epsrel đủ thấp, tùy thuộc vào từng trường hợp. Các epsabs nói chung là apriori không xác định, vì vậy tôi hiểu rằng tôi nên đặt nó rất thấp để tránh nó chiếm đầu ra, trong trường hợp đó, nó giới thiệu một khớp tính toán. Khi tôi làm cho nó thấp hơn, cảnh báo Lỗi được đưa ra rằng các lỗi làm tròn là đáng kể và tổng sai số có thể bị đánh giá thấp hơn so với dung sai mong muốn đạt được.
Trong khi câu hỏi không phải là về tốc độ, thì câu hỏi sau có liên quan mật thiết đến việc thực hiện tích hợp bộ bốn trong thực tế trước khi đặt câu hỏi về độ chính xác và dung sai. Để kiểm tra tốc độ, tôi đặt (tăng) cả bốn epsrel = 1e-02, điều này làm giảm thời gian của mã gốc xuống còn 2:14 (giờ). Sau đó, tôi đơn giản hóa quyền hạn cho mỗi Severin và thực hiện một số bản ghi nhớ . Những điều này đã giảm thời gian tích lũy xuống còn 1:29 (giờ). Các dòng đã chỉnh sửa của mã được cung cấp ở đây:
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]
Tuy nhiên, kết quả đầu ra là một tác phẩm gây ra bởi dung sai không đầy đủ được đưa vào. Tôi có thể dần dần đặt epsrel thành các giá trị thấp hơn và xem liệu kết quả có hội tụ thành giá trị thực trong thời gian thực hay không với độ chính xác scipy có sẵn. Hy vọng điều này minh họa câu hỏi ban đầu tốt hơn nhiều.