1
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy import exp
import scipy as sp
import mpmath as mp
from sympy import polylog
from math import *

def F(s,y):
  return (16/pi) * (sqrt((1-s)/(s+y))) * (log(s+y/2+sqrt(s*(s+y)))-log(y/2))
def fd1(w):
  gw=np.exp((w-mu)/TTF)
  return 1/(gw+1)
def fd123(w,y):
  gwy=np.exp((w*(1+2*y)-mu)/TTF)
  return 1/(gwy+1)
def fd12(w,y,z):
  gwyz=np.exp((w*(1+y-z)-mu)/TTF)
  return 1/(gwyz+1)
def fd13(w,y,z):
  gwyz=np.exp((w*(1+y+z)-mu)/TTF)
  return 1/(gwyz+1)
def bounds_y():
  return [0,10]
def bounds_s(y):
  return [0,1]
def bounds_w(y,z):
  return [0,10]
def bounds_z(w,y,z):
  return [-y,y]
def integrand1(z,w,s,y):
  return 144 * w**(3) *  fd13(w,y,z) * fd12(w,y,z) * (1-fd1(w)) * (1-fd123(w,y)) * F(s,y)
mu=0.22
TTF=0.5;
integrate.nquad(integrand1,[bounds_z,bounds_w,bounds_s,bounds_y])

Basically i am trying to calculate tetra-integral which depends on four variables but it kept taking too much time and eventually gave wrong answer. It also gave IntegrationWarning: The integral is probably divergent, or slowly convergent. I dont know where and what i am doing wrong.

4
  • 1
    A few comments: (i) please always post code that is well formed and enough to run without too much tinkering. (ii) your bounds include infinity, is that correct? They don't seem to work with F (iii) the parameter names on the bounds function are misleading (e.g. bounds_z(w,y,z) does not receive z as argument but w, s and y), (iv) your function F seems to be undefined for the point F(0, 0), which is part of your integration domain, if I understood correctly (v) Is your formula correct at all? integrand1(0, 0, 0, 0) gives a division by 0 (and it is within your bounds).
    – nonDucor
    Apr 18, 2021 at 23:46
  • (i) Right I ll take care next time. (ii) Yes bounds are correct but you can replace 10 or 20 instead of infinity.(iii) Parameters of bounds (e.g. bounds_z(w,y,z) can be changed with any other alphabets. The basic purpose of parameters is our positional arguments reside(Rule for integration in python). (iv) You may be right that it's undefined but quesion is how to remove this so i can get desired results (v) Formula is 100% correct.
    – Awais
    Apr 19, 2021 at 8:53
  • Posting a bounty will draw attention to your question, but that might just result in downvotes if your question isn't detailed enough; probably review How to ask and check whether you could still improve this. (I don't understand the topic well enough to decide whether this is complete or even a programming question at all.)
    – tripleee
    Apr 28, 2021 at 8:14
  • i am trying to solve my integral with nquad method. yes it is complete and a programming question.
    – Awais
    Apr 28, 2021 at 8:17

1 Answer 1

7
+25

What a workout, early in the morning!

Overview of the integral

I will use a few screenshots because stackoverflow don't support equations.

$$
\int_{y=0}^{10}\int_{s=0}^{1} \int_{w=0}^{10} \int_{z=-y}^{y} 144 \cdot w^3 \cdot  fd13(w,y,z) \cdot fd12(w,y,z) \cdot (1-fd1(w)) \cdot (1-fd123(w,y)) \cdot F(s,y) dz dw ds dy
$$

This can be rearranged as

144  \cdot \int_{y=0}^{10}\int_{s=0}^{1} F(s,y) \int_{w=0}^{10} \underbrace{(1-fd123(w,y)) \cdot (1-fd1(w)) \cdot w^3}{G(w, y)} \underbrace{\int{z=-y}^{y} fd13(w,y,z) \cdot fd12(w,y,z) \cdot  dz}_{I_z(w, y)} dw ds dy

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.

integrand

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 three level

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

z params substitution

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 enter image description here

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

assymptotic

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.

10
  • Thank you for such a great effort. But The integration in Z portion is wrong b/c w is also being multiplied with z substitution won't be that simple.
    – Awais
    Apr 30, 2021 at 17:36
  • The worst problem when substituting b = w/TTF is that in some cases b*y and a are large, and you will have do some special approximaton in this case.
    – Bob
    May 2, 2021 at 20:25
  • Your numerical consideration approximation worked quite well. But actually, I want to pass a list to mu. But when i change value of mu from 0.22 to something else then it gives math range error.
    – Awais
    May 4, 2021 at 13:45
  • Well, there is nothing special about 0.22, the curve changes continuously for real mu. For some values the convergence may become more problematic. I tested here for mu=0.5, and it gives (0.0815841273785013, 1.4897562550925374e-08) after 6 min 31s.
    – Bob
    May 4, 2021 at 16:04
  • Exactly but when i put mu 0.8 and TTF 0.25 then it gave math range error, Do you know why it's like that? Or what should i do?
    – Awais
    May 4, 2021 at 16:24

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.