What a workout, early in the morning!
Overview of the integral
I will use a few screenshots because stackoverflow don't support equations.

This can be rearranged as

Looking to the integrand there is some hope that the integration on z
can be solved analytically. The most difficult part of the integral is when both s
and y
are close to zero. both the logarithm terms on the right most side of the expression and the denominator in the square root in the left will be close to zero.

Numerical complexity experiments
Since you mentioned that it is taking forever, and we identified a some difficulties when integrating in s
and y
the first test is to do a double integral on F(z, y)
integrate.nquad(F, [bounds_s, bounds_y])
Wall time: 128 ms
(15.340467397261975, 6.86478152545078e-07)
It was easily integrated, so let's the term G(w, y)
and solve the triple integral

integrate.nquad(
lambda w, s, y: w**2 * (1-fd1(w)) * (1-fd123(w,y)) * F(s,y) ,
[bounds_w, bounds_s, bounds_y])
Wall time: 1min 5s
(5107.5853057484455, 0.00022839419762021862)
It was slower, it is accurate to about 7 decimal places, and runs in 1 minute.
So let's try to solve the integration in z
analytically.
The integration in z

I will take advantage of sympy for this piece
import sympy
from sympy.abc import a,b,z
Iz = sympy.integrate(1 / (1 + sympy.exp(2*a) + sympy.cosh(b*z)), z)
Iz
That gives

In our case b
is always non-zero.
Final calculation
The above expression gives us the indefinite integral, now calculate the difference of this expression replacing z=y
and z=-y
to compute the definite integral. This gives us one more function.
def Iwyz(w, y):
a = (w*y + w - mu)/TTF
b = 1/TTF;
tz = np.tanh(b*y/2);
ea = np.exp(-2*a)
c = np.sqrt(1 + 2*ea);
# Logarithm may be problematic non positive numbers are given
# write the sum of logarithms as logarithm of product
# since we know that the integral will be positive, the value
# argument must be postive as well
#
# also takes advantage of tanh(-b*y/2) = -tanh(b*y/2) and reuse tz
return (ea / c) * (np.log((-c-tz)*(c+tz)/((-c+tz)*(c-tz))));
And finally we can evaluate the initial integral as
integrate.nquad(lambda w, s, y: w**2 * (1-fd1(w)) * (1-fd123(w,y)) * F(s,y) * Iwyz(w, y) , [bounds_w, bounds_s, bounds_y])
Numerical considerations
The integration above will be problematic where w
is big, because if a
and b*y
will be big, thus c
and tz
will be very close to one. Then you can apply the approximations
Consider the case for large positive a
and b*y
you have

Thus the term log(c - tz)
can be well approximated by log(exp(-2*a) + 2*exp(-b*y))
. You could write similar expressions for the four combinations of sign for a
and b*y
. But since the integration takes place for non-negative w
and y
The updated function could be
def Iwyz(w, y):
a = (w*y + w - mu)/TTF
if w == 0:
return z/(np.exp(2*a) + 2)
else:
b = w/TTF;
tz = np.tanh(b*y/2);
ea = np.exp(-2*a)
c = np.sqrt(1 + 2*ea);
# Logarithm may be problematic non positive numbers are given
# write the sum of logarithms as logarithm of product
# since we know that the integral will be positive, the value
# argument must be postive as well
#
# also takes advantage of tanh(-b*y/2) = -tanh(b*y/2) and reuse tz
eby = np.exp(-b*y)
# both w and y are non-negative, we only have to handle this case
if eby + ea < 1e-6:
return (ea / c) * (2*log(2) - log(ea + 2 * eby))
else:
return (ea / c) * (np.log(abs(c+tz)/(abs(c-tz))));
You may want to choose a different condition to switch to the asymptotic approximation, you could also chose a different approximation.
The integration results using asymptotic approximation if eby + ea < 1e-6
3min 41s
(0.038781260929779036, 1.4898372223919173e-08)
To evaluate the impact of the approximation I will recalculate applying the approximation in based on a different condition.
Results using asymptotic approximation if eby + ea < 1e-9
Wall time: 3min 44s
(0.038781260929776414, 1.4898372223919173e-08)
The results match for 13 decimal places, meaning that the error of the integration is higher than the error introduced by the asymptotic approximation.