OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Different answer for numerical integral from nquad and Mathematica

  • Thread starter Thread starter QFTheorist
  • Start date Start date
Q

QFTheorist

Guest
I am trying to write a function that returns the numerical answer of a very complicated integral. I have the following code

Code:
import numpy as np
import scipy as sp
import sympy as smp
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.integrate import nquad
import cmath
from scipy.special import gamma #need the gamma function from scipy
#loading Gegenbauer from scipy.special
from scipy.special import eval_gegenbauer as G #gegenbauer

options={'limit':1000, 'epsabs':1e-12, 'epsrel':1e-12}

def IntZ(z,sp,s,l,lp):
    return (1-z**2)**((d-4)/2)/2*G(l,(d-3)/2,z)*G(lp,(d-3)/2,cmath.sqrt(-(3*s*(-1+z**2)+sp*(3+z**2))/((s*(-1+z**2)-sp*(3+z**2)))))
def ker(z,sp,s):
#     return -2+s/(sp-s)+(4*sp*(s+2*sp))/((s+2*sp)**2-s**2*z**2)
    return s*(1/(complex(0,-10**(-20))-s+sp)+(-1+z)/(s+2*sp-s*z)-(1+z)/(s+2*sp+s*z))

#These are all the prefactors of the integrand.

def prefactor(sp,s,sigma,l,lp):
    return ((gamma(l+1)*gamma((d-3)/2)*2**(2*d-3)*(np.pi)**((d-3)/2)*gamma((d-3)/2)*(2*lp+d-3))/(2**(d+1)*(np.pi)**((d-1)/2)*gamma(l+d-3)))*(((s-mu)**((d-3)/2)*2*cmath.sqrt((sp-1)*sigma)*cmath.sqrt(sp)*(sp-1)**(2*lp+1/2))/(cmath.sqrt(s)*(1-sp-sigma)*(sp-mu)**((d-3)/2)*(sp)**(2*lp+1/2)))

#this is the final integrand (sp and z) with all the factors coded in.
def integrand(z,sp,s,sigma,l,lp):
    return IntZ(z,sp,s,l,lp)*ker(z,sp,s)*prefactor(sp,s,sigma,l,lp)

def prefactor2(sp,s,sigma,lp):
    return ((s-mu)**((d-3)/2)*2*cmath.sqrt((sp-1)*sigma)*cmath.sqrt(sp)*(sp-1)**(2*lp+1/2))/(cmath.sqrt(s)*(1-sp-sigma)*(sp-mu)**((d-3)/2)*(sp)**(2*lp+1/2))

def fell(s,sigma,l,lp):
    return nquad(integrand,[[-1,1],[1,np.inf]],args=(s,sigma,lp,l),opts=[options,options])

Note that this only gives me the real part of the integral, which is fine. I also have the integraion defined in Mathematica with the following function:

Code:
d = 5; 
\[Mu] = 0; 
(*Here the z-integral is done numerically*)
Integral2[s_, \[Sigma]_, l_, lp_, MaxRes_, WorkPres_] := 
 NIntegrate[(l! Gamma[(d - 3)/2])/(
   2^(d + 1) \[Pi]^((d - 1)/2) Gamma[l + d - 3]) (s - \[Mu])^((d - 3)/
    2)/s^(1/2) 2^(2 d - 3) \[Pi]^((d - 3)/2) Gamma[(d - 3)/2] sp^(
   1/2)/(sp - \[Mu])^((d - 3)/2) (2 lp + d - 3) ((sp - 1)/sp)^(
   2 lp + 1/2) (2 Sqrt[(sp - 1) \[Sigma]])/(
   1 - sp - \[Sigma]) (s/(sp - (s + I 10^-20)) + t/(sp - t) + u/(
       sp - u) /. {u -> -s - t} /. {t -> (s (z - 1))/2}) (1 - z^2)^((
   d - 4)/2)/
   2 GegenbauerC[l, (d - 3)/2, 
    z] (GegenbauerC[lp, (d - 3)/2, Sqrt[
        1 + (4 a)/(-a + sp)]] /. {a -> (s t u)/(
         s t + t u + u s)} /. {u -> -s - t} /. {t -> (s (z - 1))/
       2}), {z, -1, 1}, {sp, 1, Infinity}, MaxRecursion -> MaxRes, 
  WorkingPrecision -> WorkPres, Method -> "GlobalAdaptive"]

When I run the code fell(200,45,4,4), I get the answer (915.9500431066836, 0.0008916854858398438). This looks reasonable enough since the error is very small. However, when I run the same integral in Mathematica, with the same values and MaxRecursion and WorkingPrecision set to 100 each, I get (Mathematica returns both the real and the imaginary part)

Code:
946.243126048692722003932251899650713550717875143021622444266774813501\
3068351836006101088945258421851 - 
 288.75854540972949504131273023372598463063894553182017391668357286673\
70457344996286394039275246094406 I

The reason for doing this integral using Python instead of Mathematica was because in Python it's much faster. Since both the answers are different, I am not sure which one I should believe in.

For the integral in Python, I have tried to change values of limit, epsabs and epsrel but it does not help a lot in changing the answer. The discrepancy becomes worse for larger values of the first two arguments of the function, where I get for fell(1800,1800,4,4) the answer (5582.867565858222, 90.11617121750788) but Mathematica returns

Code:
6600.16959988860178799109412672538944693500880870136958557262039403043\
5553207340978200772084077156640 - 
 5628.2185973484253452277428879848996442046986842344701373348873276687\
65333223655602231358008282651070 I

What is going wrong?
<p>I am trying to write a function that returns the numerical answer of a very complicated integral. I have the following code</p>
<pre><code>import numpy as np
import scipy as sp
import sympy as smp
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.integrate import nquad
import cmath
from scipy.special import gamma #need the gamma function from scipy
#loading Gegenbauer from scipy.special
from scipy.special import eval_gegenbauer as G #gegenbauer

options={'limit':1000, 'epsabs':1e-12, 'epsrel':1e-12}

def IntZ(z,sp,s,l,lp):
return (1-z**2)**((d-4)/2)/2*G(l,(d-3)/2,z)*G(lp,(d-3)/2,cmath.sqrt(-(3*s*(-1+z**2)+sp*(3+z**2))/((s*(-1+z**2)-sp*(3+z**2)))))
def ker(z,sp,s):
# return -2+s/(sp-s)+(4*sp*(s+2*sp))/((s+2*sp)**2-s**2*z**2)
return s*(1/(complex(0,-10**(-20))-s+sp)+(-1+z)/(s+2*sp-s*z)-(1+z)/(s+2*sp+s*z))

#These are all the prefactors of the integrand.

def prefactor(sp,s,sigma,l,lp):
return ((gamma(l+1)*gamma((d-3)/2)*2**(2*d-3)*(np.pi)**((d-3)/2)*gamma((d-3)/2)*(2*lp+d-3))/(2**(d+1)*(np.pi)**((d-1)/2)*gamma(l+d-3)))*(((s-mu)**((d-3)/2)*2*cmath.sqrt((sp-1)*sigma)*cmath.sqrt(sp)*(sp-1)**(2*lp+1/2))/(cmath.sqrt(s)*(1-sp-sigma)*(sp-mu)**((d-3)/2)*(sp)**(2*lp+1/2)))

#this is the final integrand (sp and z) with all the factors coded in.
def integrand(z,sp,s,sigma,l,lp):
return IntZ(z,sp,s,l,lp)*ker(z,sp,s)*prefactor(sp,s,sigma,l,lp)

def prefactor2(sp,s,sigma,lp):
return ((s-mu)**((d-3)/2)*2*cmath.sqrt((sp-1)*sigma)*cmath.sqrt(sp)*(sp-1)**(2*lp+1/2))/(cmath.sqrt(s)*(1-sp-sigma)*(sp-mu)**((d-3)/2)*(sp)**(2*lp+1/2))

def fell(s,sigma,l,lp):
return nquad(integrand,[[-1,1],[1,np.inf]],args=(s,sigma,lp,l),opts=[options,options])
</code></pre>
<p>Note that this only gives me the real part of the integral, which is fine. I also have the integraion defined in Mathematica with the following function:</p>
<pre><code>d = 5;
\[Mu] = 0;
(*Here the z-integral is done numerically*)
Integral2[s_, \[Sigma]_, l_, lp_, MaxRes_, WorkPres_] :=
NIntegrate[(l! Gamma[(d - 3)/2])/(
2^(d + 1) \[Pi]^((d - 1)/2) Gamma[l + d - 3]) (s - \[Mu])^((d - 3)/
2)/s^(1/2) 2^(2 d - 3) \[Pi]^((d - 3)/2) Gamma[(d - 3)/2] sp^(
1/2)/(sp - \[Mu])^((d - 3)/2) (2 lp + d - 3) ((sp - 1)/sp)^(
2 lp + 1/2) (2 Sqrt[(sp - 1) \[Sigma]])/(
1 - sp - \[Sigma]) (s/(sp - (s + I 10^-20)) + t/(sp - t) + u/(
sp - u) /. {u -> -s - t} /. {t -> (s (z - 1))/2}) (1 - z^2)^((
d - 4)/2)/
2 GegenbauerC[l, (d - 3)/2,
z] (GegenbauerC[lp, (d - 3)/2, Sqrt[
1 + (4 a)/(-a + sp)]] /. {a -> (s t u)/(
s t + t u + u s)} /. {u -> -s - t} /. {t -> (s (z - 1))/
2}), {z, -1, 1}, {sp, 1, Infinity}, MaxRecursion -> MaxRes,
WorkingPrecision -> WorkPres, Method -> "GlobalAdaptive"]
</code></pre>
<p>When I run the code <code>fell(200,45,4,4)</code>, I get the answer <code>(915.9500431066836, 0.0008916854858398438)</code>. This looks reasonable enough since the error is very small. However, when I run the same integral in Mathematica, with the same values and <code>MaxRecursion</code> and <code>WorkingPrecision</code> set to 100 each, I get (Mathematica returns both the real and the imaginary part)</p>
<pre class="lang-none prettyprint-override"><code>946.243126048692722003932251899650713550717875143021622444266774813501\
3068351836006101088945258421851 -
288.75854540972949504131273023372598463063894553182017391668357286673\
70457344996286394039275246094406 I
</code></pre>
<p>The reason for doing this integral using Python instead of Mathematica was because in Python it's much faster. Since both the answers are different, I am not sure which one I should believe in.</p>
<p>For the integral in Python, I have tried to change values of <code>limit</code>, <code>epsabs</code> and <code>epsrel</code> but it does not help a lot in changing the answer. The discrepancy becomes worse for larger values of the first two arguments of the function, where I get for <code>fell(1800,1800,4,4)</code> the answer <code>(5582.867565858222, 90.11617121750788)</code> but Mathematica returns</p>
<pre class="lang-none prettyprint-override"><code>6600.16959988860178799109412672538944693500880870136958557262039403043\
5553207340978200772084077156640 -
5628.2185973484253452277428879848996442046986842344701373348873276687\
65333223655602231358008282651070 I
</code></pre>
<p>What is going wrong?</p>
 

Online statistics

Members online
0
Guests online
5
Total visitors
5
Top