¿Existe un integrador múltiple en Python que proporcione límites de integración variable (como scipy) y alta precisión (como mpmath)?
Puedo usar scipy quad y nquad para una integración cuádruple que implique límites de integración variable. El problema es que la precisión predeterminada utilizada genera un error cuando no se puede alcanzar la tolerancia solicitada. Con el integrador mpmath, puedo definir cualquier precisión arbitraria con la configuración de mp.dps = arbitrary, pero no puedo ver si los límites pueden volverse variables como con nquad y cómo. Mpmath también proporciona una ejecución muy rápida con el método Gauss-Legendre en quadgl, lo cual es muy deseable, porque mi función es fluida, pero toma una cantidad exorbitante de tiempo con scipy para completar cuatro integraciones. Por favor ayuda. La siguiente es solo una función simple que falla en mi objetivo:
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)
Respuestas
Ok, déjame responder algo, es difícil poner código en los comentarios
La optimización simple con matemáticas MP consiste en seguir reglas simples:
- y 2.0 es MUY caro (log, exp, ...), reemplácelo con y * y
- y 2 sigue siendo caro, reemplácelo con y * y
- la multiplicación es mucho más cara que la suma, reemplace x * y + y ** 2.0 con (x + y) * y
- La división es más cara que la multiplicación, reemplace y / 4 con 0.25 * y
Código, 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])
en mi computadora pasó de 12,9 segundos a 10,6 segundos, aproximadamente un 20% de descuento
A continuación se muestra un ejemplo simple de cómo puedo hacer solo una integración triple con mpmath. Esto no aborda la alta precisión con cuatro integraciones. En cualquier caso, el tiempo de ejecución es un problema aún mayor. Cualquier ayuda bienvenida.
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
Además, un intento de combinar una (primera) integración scipy seguida de un integrador triple mpmath no parece producir ningún resultado durante más de 24 horas, incluso con una función más simple. ¿Qué hay de malo en el siguiente código?
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)
A continuación se muestra el código completo, para el cual se planteó la pregunta original. Contiene el uso de scipy para realizar cuatro integraciones:
# 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()
Estoy interesado en establecer el epsrel lo suficientemente bajo, según el caso. El epsabs es generalmente desconocido a priori, por lo que entiendo que debería hacerlo muy bajo para evitar que se apodere de la salida, en cuyo caso introduce un articacto computacional. Cuando lo hago más bajo, se genera una advertencia de error de que los errores de redondeo son significativos y que el error total puede subestimarse para lograr la tolerancia deseada.
Si bien la cuestión no es sobre velocidad, esta última está íntimamente relacionada con hacer práctica la ejecución de una integración cuádruple antes de la indagación sobre precisión y tolerancia. Para probar la velocidad, configuré (aumenté) los cuatro epsrel = 1e-02, lo que redujo el tiempo del código original a 2:14 (horas). Luego simplifiqué los poderes según Severin e implementé algunas memorias . Estos redujeron el tiempo acumulativamente a 1:29 (horas). Las líneas editadas del código se proporcionan aquí:
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]
Sin embargo, la salida es un artefacto causado por la tolerancia inadecuada introducida. Puedo establecer progresivamente el epsrel en valores más bajos y ver si el resultado converge a un valor realista en tiempo realista con la precisión scipy disponible. Espero que esto ilustre mucho mejor la pregunta original.