Python에 변수 통합 한계 (scipy와 같은)와 높은 정밀도 (mpmath와 같은)를 모두 제공하는 다중 통합자가 있습니까?

Aug 17 2020

가변 적분 한계를 포함하는 사중 적분에 scipy quad 및 nquad를 사용할 수 있습니다. 문제는 요청 된 허용 오차를 달성 할 수 없을 때 사용 된 기본 정밀도가 오류를 발생 시킨다는 것입니다. mpmath integrator를 사용하면 mp.dps = 임의 설정으로 임의의 정밀도를 정의 할 수 있지만 nquad와 같이 제한이 어떻게 변할 수 있는지 여부와 방법을 알 수 없습니다. Mpmath는 또한 quadgl에서 Gauss-Legendre 방법으로 매우 빠른 실행을 제공하는데, 이는 내 기능이 부드럽기 때문에 매우 바람직하지만 4 개의 통합을 완료하는 데 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)

답변

1 SeverinPappadeux Aug 20 2020 at 00:54

좋아, 대답에 뭔가를 넣어 보자. 주석에 코드를 넣기가 어렵습니다.

MP 수학을 사용한 간단한 최적화는 간단한 규칙을 따르는 것입니다.

y ^2.0 은 매우 비쌉니다 (log, exp, ...), y * y로 대체
y ² 는 여전히 비싸므로 y * y로 대체
곱셈은 합산보다 훨씬 비쌉니다. x * y + y ** 2.0을 (x + y) * y로 바꿉니다.
나눗셈은 곱셈보다 비싸므로 y / 4를 0.25 * y로 바꿉니다.

코드, 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])

내 컴퓨터에서 12.9 초에서 10.6 초로 약 20 % 할인되었습니다.

1 gerryD Aug 19 2020 at 14:08

다음은 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이므로 출력을 확보하지 않으려면 매우 낮게 설정해야하며,이 경우 계산 기술이 도입됩니다. 더 낮게 설정하면 반올림 오류가 중요하고 원하는 허용 오차를 달성하기 위해 총 오류가 과소 평가 될 수 있다는 오류 경고가 발생합니다.

gerryD Aug 24 2020 at 06:40

문제는 속도에 관한 것이 아니지만 후자는 정밀도와 허용 오차에 대한 질문에 앞서 4 중 통합을 실행하는 것과 밀접하게 관련되어 있습니다. 속도를 테스트하기 위해 epsrel = 1e-02 4 개를 모두 설정 (증가)하여 원래 코드의 시간을 2:14 (시간)로 줄였습니다. 그런 다음 Severin 당 권한을 단순화하고 일부 메모를 구현했습니다 . 이로 인해 누적 시간이 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을 더 낮은 값으로 설정하고 결과가 사용 가능한 scipy 정밀도로 현실적인 시간에 현실적인 값으로 수렴하는지 확인할 수 있습니다. 이것이 원래 질문을 훨씬 더 잘 설명하기를 바랍니다.