<p><strong>tl;dr:</strong> Tiny numbers lead to <a href="https://en.wikipedia.org/wiki/Round-off_error" rel="nofollow noreferrer">roundoff errors</a> during integration?</p>
<p>It's your first question, <em>Welcome to Stack Exchange!</em></p>
<p>You didn't mention exactly what you mean by errors. If you had asked this in <a href="https://scicomp.stackexchange.com/">Scientific Computing SE</a> the readers there would have seen the line</p>
<p><code>integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))</code></p>
<p>and known that the errors returned by nquad are what you are talking about.</p>
<p>But here on Astronomy SE where there are few numerical questions, readers needed that to be explained clearly.</p>
<p>You are using <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.nquad.html" rel="nofollow noreferrer"><code>scipy.integrate.nquad</code></a> which &quot;Wraps <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html" rel="nofollow noreferrer"><code>quad</code></a> to enable integration over multiple variables.&quot;</p>
<p>I tried changing all of the parameters available that can potentially reduce the error;</p>
<pre><code>def j_nu_2(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,),
                       opts={'epsrel': 1E-16, 'epsabs': 1E-16, 'limit': 500,
                             'maxp1': 500, 'limlst': 500})
    return integrand,err
</code></pre>
<p>but I got <em>exactly</em> the same result. That means that the <code>quad</code> is seeing the simple, smooth shape and deciding internally that it doesn't need so many steps to obtain a good result.</p>
<p>In this case, I wonder if you are encountering some numerical errors due to the values being very very very small. Once your result drops below 1E-27, the relative error starts increasing.</p>
<p>I recommend that you pull out all of the constants you can from the integrand, and multiply them back in only after integration.</p>
<p>Try to keep your integrand closer to one, no 1E+30s or 1E-30s. Computers are finite things and they have limited digits.</p>
<hr />
<p><a href="https://i.sstatic.net/nKGtV.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/nKGtV.png" alt="output from script" /></a></p>
<p>my script (scroll to the bottom for the new stuff):</p>
<pre><code>import numpy as np
from scipy.integrate import nquad
import matplotlib.pyplot as plt


N0 = 1
h = 6.626e-27 #  Planck const in erg. s
c = 3e10  # Speed of light in cm/s
k = 1.38e-16  # Boltzmann constant in erg/K
r0 = 2.82e-13  # Classical electron radius in cm
## integating gamma--&gt; nu_i ,nu is parameter

def integrand_G(gamma,nu_i,nu):
    N_e = N0 * gamma**(-3) 
    N_ph = 8 * np.pi * nu_i**2 / (c**3 * (np.exp(h*nu_i/(k*T_ph))-1) ) 

    x = nu / (4 * gamma**2 * nu_i)
    const = 8 * np.pi * r0**2 * c * h
    if 0 &lt; x &lt; 1:
        fx = 2 * x**2 * np.log(x) + x  + x**2 - 2*x**3
        result = const * fx * N_e * N_ph
        return result
    else:
        return 0

def range_gamma(nu_i,nu):
    return [0.5* (nu / nu_i)**0.5,np.inf]

range_nui = [1e13,1e16]

def j_nu(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
    return integrand,err
    
T_ph =5e4
nu01 = np.logspace(13,16,50)
gamma1 = 10
nu_values =4/3*gamma1**2 * nu01  

j_nu_values = [j_nu(nu) for nu in nu_values]

j_nu_valuesN0 = 1
h = 6.626e-27 #  Planck const in erg. s
c = 3e10  # Speed of light in cm/s
k = 1.38e-16  # Boltzmann constant in erg/K
r0 = 2.82e-13  # Classical electron radius in cm
## integating gamma--&gt; nu_i ,nu is parameter

def integrand_G(gamma,nu_i,nu):
    N_e = N0 * gamma**(-3) 
    N_ph = 8 * np.pi * nu_i**2 / (c**3 * (np.exp(h*nu_i/(k*T_ph))-1) ) 

    x = nu / (4 * gamma**2 * nu_i)
    const = 8 * np.pi * r0**2 * c * h
    if 0 &lt; x &lt; 1:
        fx = 2 * x**2 * np.log(x) + x  + x**2 - 2*x**3
        result = const * fx * N_e * N_ph
        return result
    else:
        return 0

def range_gamma(nu_i,nu):
    return [0.5* (nu / nu_i)**0.5,np.inf]

range_nui = [1e13,1e16]

def j_nu(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
    return integrand,err
    
T_ph =5e4
nu01 = np.logspace(13,16,50)
gamma1 = 10
nu_values =4/3*gamma1**2 * nu01  

j_nu_values = [j_nu(nu) for nu in nu_values]

# new stuff

def plotit(jnuval):

    results, errors = zip(*jnuval)
    relative_errors = [e/r for (r, e) in jnuval]

    fig, (ax1, ax2) = plt.subplots(2, 1)
    ax1.plot(results, label='results')
    ax1.plot(errors, label='errors')
    ax1.set_yscale('log')
    ax1.legend()
    ax1.set_title('result and error')
    ax2.plot(relative_errors)
    ax2.set_yscale('log')
    ax2.set_title('relative errors')

    plt.show()

plotit(j_nu_values)

def j_nu_2(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,),
                           opts={'epsrel': 1E-16, 'epsabs': 1E-16, 'limit': 500,
                                 'maxp1': 500, 'limlst': 500})
    return integrand,err

j_nu_values_2 = [j_nu_2(nu) for nu in nu_values]

plotit(j_nu_values_2)
</code></pre>

CC BY-SA 4.0
2024-01-05T07:54:39.183Z
Your reply was very helpful to me, thank you very much for that.
Shayu xiao
CC BY-SA 4.0
2024-01-05T08:33:28.117Z
@Shayuxiao that's great! You don't need to accept this answer yet - why don't you wait until you've got your calculation running well and add an answer describing your final solution. This is just a guess - I don't know for sure if it's roundoff or not. Accepting quickly may discourage other answers that may turn out to be even more helpful!
uhoh
CC BY-SA 4.0
2024-01-05T01:43:42.773Z
2024-01-05T01:43:42.773Z
uhoh
<p><strong>NOTE:</strong> At the end of my post the output is shown as a list of tuples. The first value in each tuple is the result and the second value is the estimated error on that result. <em>Near the end, <strong>the estimated error is of order one part per thousand of the result, which is much too high for a numerical integration</strong> over such smoothly varying functions.</em></p>
<hr />
<h4>1. Inverse Compton radiation of electron collectives</h4>
<p><span class="math-container">$$j(\nu)=\int P(\nu)N(\gamma)\mathrm{d}\gamma=8\pi r_0^2chG$$</span>
thereinto</p>
<p><span class="math-container">$$G=\int_{\nu_1}^{\nu_2}\int_{\gamma_1=\frac12(\frac{\nu}{\nu_i})^{1/2}}^{\gamma_2} N(\gamma)f(\frac{\nu}{4\gamma^2\nu_i})n_{\mathrm{ph}}(\nu_i)\mathrm{d}\nu_i\mathrm{ d}\gamma $$</span>
thereinto
<span class="math-container">$$f(x)=2x^2\mathrm{ln}x+x+x^2-2x^3，0&lt;x&lt;1$$</span>
where <span class="math-container">$x = \nu/4\gamma^{2}\nu_{i}$</span></p>
<h4>2. Background photon number density distribution in a radiated field</h4>
<p><span class="math-container">$$n_{\mathrm{ph}}(\nu_i)=\frac{8\pi\nu_i^2}{c^3}\frac1{e^{\frac{h\nu_i}{kT}}-1}$$</span></p>
<h4>3. Electrons have a power-law distribution</h4>
<p><span class="math-container">$$N(\gamma) = N\gamma^{-n}$$</span></p>
<pre><code>N0 = 1
h = 6.626e-27 #  Planck const in erg. s
c = 3e10  # Speed of light in cm/s
k = 1.38e-16  # Boltzmann constant in erg/K
r0 = 2.82e-13  # Classical electron radius in cm
## integating gamma--&gt; nu_i ,nu is parameter

def integrand_G(gamma,nu_i,nu):
    N_e = N0 * gamma**(-3) 
    N_ph = 8 * np.pi * nu_i**2 / (c**3 * (np.exp(h*nu_i/(k*T_ph))-1) ) 

    x = nu / (4 * gamma**2 * nu_i)
    const = 8 * np.pi * r0**2 * c * h
    if 0 &lt; x &lt; 1:
        fx = 2 * x**2 * np.log(x) + x  + x**2 - 2*x**3
        result = const * fx * N_e * N_ph
        return result
    else:
        return 0

def range_gamma(nu_i,nu):
    return [0.5* (nu / nu_i)**0.5,np.inf]

range_nui = [1e13,1e16]

def j_nu(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
    return integrand,err
    
T_ph =5e4
nu01 = np.logspace(13,16,50)
gamma1 = 10
nu_values =4/3*gamma1**2 * nu01  

j_nu_values = [j_nu(nu) for nu in nu_values]

j_nu_values

</code></pre>
<p>the result is</p>
<pre><code>[(4.646891358736883e-25, 9.732928577801494e-34),
 (4.0358780167911975e-25, 3.7331558505127706e-34),
 (3.50520595581472e-25, 8.745269160624114e-37),
 (3.044311236175784e-25, 4.596524332722152e-34),
 (2.6440189335175633e-25, 6.9773996637028405e-37),
 (2.2963605183798276e-25, 2.32847015875059e-34),
 (1.9944152312040706e-25, 5.106011672520151e-35),
 (1.7321723128286085e-25, 6.167458221137299e-35),
 (1.5044113565389222e-25, 6.549329044601913e-35),
 (1.3065983749075148e-25, 6.1367624503818e-35),
 (1.1347955500240598e-25, 5.196262110593544e-35),
 (9.855828423034596e-26, 4.672407189709615e-35),
 (8.559899083742284e-26, 4.015760139155061e-35),
 (7.434369712819832e-26, 3.4425082332554275e-35),
 (6.456834652685337e-26, 2.957012268335199e-35),
 (5.607834334707997e-26, 2.5772446384380916e-35),
 (4.8704679021211965e-26, 2.3382198801399686e-35),
 (4.230056768532779e-26, 2.2967597901824962e-35),
 (3.6738524146513086e-26, 2.540934279247937e-35),
 (3.1907826064254274e-26, 3.2026146513980704e-35),
 (2.7712309826214774e-26, 4.4521946638201733e-35),
 (2.4068456241343245e-26, 6.432951098088548e-35),
 (2.0903728116056675e-26, 1.0795160235615729e-35),
 (1.8155125637870292e-26, 3.831168995920699e-36),
 (1.5767933111091125e-26, 2.6129469178833845e-36),
...
 (8.166934445646721e-28, 1.7567146260202947e-33),
 (7.0930754638972155e-28, 2.1515340634505886e-33),
 (6.160414678071309e-28, 6.411250698120658e-34),
 (5.350517710760285e-28, 2.308277273356424e-31),
 (4.647002306856651e-28, 2.4616285151340818e-31)]
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2024-01-03T21:16:06.303Z
"I expected [...], but instead I got [...]"
Darth Pseudonym
CC BY-SA 4.0
2024-01-04T01:39:59.843Z
I'm sorry I didn't express my problem, I mean when I do numerical integration, the result error of the integration is a bit large, how can I solve it?

My source code is attached to the issue
Shayu xiao
CC BY-SA 4.0
2024-01-05T00:27:54.093Z
What is "a bit large". Please follow the hint of DarthPseudonym and edit your question with contrasting result and expectation clearly and why you think that the difference matters if e.g. the expected result was obtained numerically differently.
planetmaker
CC BY-SA 4.0
2024-01-05T01:46:44.237Z
*Welcome to Stack Exchange!* You've done a good job asking your first question, but here in Astronomy folks will want you to explain even more clearly exactly what the problem is and won't see that your integration function returns errors and that you have shown the errors. I've made an edit to show you how to do that more clearly. I think I see the problem so I've posted an answer. Hopefully the question ***stays open so that further answers are not prevented.***
uhoh
<radiation>
The result of my calculation gives too much error, how to solve it?
Shayu xiao
# The result of my calculation gives too much error, how to solve it?

(Asked by: username_0 on 2024-01-03)

NOTE: At the end of my post the output is shown as a list of tuples. The first value in each tuple is the result and the second value is the estimated error on that result. Near the end, the estimated error is of order one part per thousand of the result, which is much too high for a numerical integration over such smoothly varying functions.

1. Inverse Compton radiation of electron collectives
$$j(\nu)=\int P(\nu)N(\gamma)\mathrm{d}\gamma=8\pi r_0^2chG$$
thereinto
$$G=\int_{\nu_1}^{\nu_2}\int_{\gamma_1=\frac12(\frac{\nu}{\nu_i})^{1/2}}^{\gamma_2} N(\gamma)f(\frac{\nu}{4\gamma^2\nu_i})n_{\mathrm{ph}}(\nu_i)\mathrm{d}\nu_i\mathrm{ d}\gamma $$
thereinto
$$f(x)=2x^2\mathrm{ln}x+x+x^2-2x^3，0
where $x = \nu/4\gamma^{2}\nu_{i}$
2. Background photon number density distribution in a radiated field
$$n_{\mathrm{ph}}(\nu_i)=\frac{8\pi\nu_i^2}{c^3}\frac1{e^{\frac{h\nu_i}{kT}}-1}$$
3. Electrons have a power-law distribution
$$N(\gamma) = N\gamma^{-n}$$
N0 = 1
h = 6.626e-27 #  Planck const in erg. s
c = 3e10  # Speed of light in cm/s
k = 1.38e-16  # Boltzmann constant in erg/K
r0 = 2.82e-13  # Classical electron radius in cm
## integating gamma--> nu_i ,nu is parameter

def integrand_G(gamma,nu_i,nu):
    N_e = N0 * gamma**(-3) 
    N_ph = 8 * np.pi * nu_i**2 / (c**3 * (np.exp(h*nu_i/(k*T_ph))-1) ) 

    x = nu / (4 * gamma**2 * nu_i)
    const = 8 * np.pi * r0**2 * c * h
    if 0 < x < 1:
        fx = 2 * x**2 * np.log(x) + x  + x**2 - 2*x**3
        result = const * fx * N_e * N_ph
        return result
    else:
        return 0

def range_gamma(nu_i,nu):
    return [0.5* (nu / nu_i)**0.5,np.inf]

range_nui = [1e13,1e16]

def j_nu(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
    return integrand,err
    
T_ph =5e4
nu01 = np.logspace(13,16,50)
gamma1 = 10
nu_values =4/3*gamma1**2 * nu01  

j_nu_values = [j_nu(nu) for nu in nu_values]

j_nu_values


the result is
[(4.646891358736883e-25, 9.732928577801494e-34),
 (4.0358780167911975e-25, 3.7331558505127706e-34),
 (3.50520595581472e-25, 8.745269160624114e-37),
 (3.044311236175784e-25, 4.596524332722152e-34),
 (2.6440189335175633e-25, 6.9773996637028405e-37),
 (2.2963605183798276e-25, 2.32847015875059e-34),
 (1.9944152312040706e-25, 5.106011672520151e-35),
 (1.7321723128286085e-25, 6.167458221137299e-35),
 (1.5044113565389222e-25, 6.549329044601913e-35),
 (1.3065983749075148e-25, 6.1367624503818e-35),
 (1.1347955500240598e-25, 5.196262110593544e-35),
 (9.855828423034596e-26, 4.672407189709615e-35),
 (8.559899083742284e-26, 4.015760139155061e-35),
 (7.434369712819832e-26, 3.4425082332554275e-35),
 (6.456834652685337e-26, 2.957012268335199e-35),
 (5.607834334707997e-26, 2.5772446384380916e-35),
 (4.8704679021211965e-26, 2.3382198801399686e-35),
 (4.230056768532779e-26, 2.2967597901824962e-35),
 (3.6738524146513086e-26, 2.540934279247937e-35),
 (3.1907826064254274e-26, 3.2026146513980704e-35),
 (2.7712309826214774e-26, 4.4521946638201733e-35),
 (2.4068456241343245e-26, 6.432951098088548e-35),
 (2.0903728116056675e-26, 1.0795160235615729e-35),
 (1.8155125637870292e-26, 3.831168995920699e-36),
 (1.5767933111091125e-26, 2.6129469178833845e-36),
...
 (8.166934445646721e-28, 1.7567146260202947e-33),
 (7.0930754638972155e-28, 2.1515340634505886e-33),
 (6.160414678071309e-28, 6.411250698120658e-34),
 (5.350517710760285e-28, 2.308277273356424e-31),
 (4.647002306856651e-28, 2.4616285151340818e-31)]
```

**Tags:** radiation
> **Comments:**
> **username_1**: "I expected [...], but instead I got [...]"  
> **username_0**: I'm sorry I didn't express my problem, I mean when I do numerical integration, the result error of the integration is a bit large, how can I solve it?  My source code is attached to the issue  
> **username_2**: What is "a bit large". Please follow the hint of DarthPseudonym and edit your question with contrasting result and expectation clearly and why you think that the difference matters if e.g. the expected result was obtained numerically differently.  
> **username_3**: *Welcome to Stack Exchange!* You've done a good job asking your first question, but here in Astronomy folks will want you to explain even more clearly exactly what the problem is and won't see that your integration function returns errors and that you have shown the errors. I've made an edit to show you how to do that more clearly. I think I see the problem so I've posted an answer. Hopefully the question ***stays open so that further answers are not prevented.***  

---

## 1 answers

### Answer by username_3 on 2024-01-05. Score: 1 | Accepted answer

tl;dr: Tiny numbers lead to roundoff errors during integration?
It's your first question, Welcome to Stack Exchange!
You didn't mention exactly what you mean by errors. If you had asked this in Scientific Computing SE the readers there would have seen the line
integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
and known that the errors returned by nquad are what you are talking about.
But here on Astronomy SE where there are few numerical questions, readers needed that to be explained clearly.
You are using scipy.integrate.nquad which "Wraps quad to enable integration over multiple variables."
I tried changing all of the parameters available that can potentially reduce the error;
def j_nu_2(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,),
                       opts={'epsrel': 1E-16, 'epsabs': 1E-16, 'limit': 500,
                             'maxp1': 500, 'limlst': 500})
    return integrand,err

but I got exactly the same result. That means that the quad is seeing the simple, smooth shape and deciding internally that it doesn't need so many steps to obtain a good result.
In this case, I wonder if you are encountering some numerical errors due to the values being very very very small. Once your result drops below 1E-27, the relative error starts increasing.
I recommend that you pull out all of the constants you can from the integrand, and multiply them back in only after integration.
Try to keep your integrand closer to one, no 1E+30s or 1E-30s. Computers are finite things and they have limited digits.


my script (scroll to the bottom for the new stuff):
import numpy as np
from scipy.integrate import nquad
import matplotlib.pyplot as plt


N0 = 1
h = 6.626e-27 #  Planck const in erg. s
c = 3e10  # Speed of light in cm/s
k = 1.38e-16  # Boltzmann constant in erg/K
r0 = 2.82e-13  # Classical electron radius in cm
## integating gamma--> nu_i ,nu is parameter

def integrand_G(gamma,nu_i,nu):
    N_e = N0 * gamma**(-3) 
    N_ph = 8 * np.pi * nu_i**2 / (c**3 * (np.exp(h*nu_i/(k*T_ph))-1) ) 

    x = nu / (4 * gamma**2 * nu_i)
    const = 8 * np.pi * r0**2 * c * h
    if 0 < x < 1:
        fx = 2 * x**2 * np.log(x) + x  + x**2 - 2*x**3
        result = const * fx * N_e * N_ph
        return result
    else:
        return 0

def range_gamma(nu_i,nu):
    return [0.5* (nu / nu_i)**0.5,np.inf]

range_nui = [1e13,1e16]

def j_nu(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
    return integrand,err
    
T_ph =5e4
nu01 = np.logspace(13,16,50)
gamma1 = 10
nu_values =4/3*gamma1**2 * nu01  

j_nu_values = [j_nu(nu) for nu in nu_values]

j_nu_valuesN0 = 1
h = 6.626e-27 #  Planck const in erg. s
c = 3e10  # Speed of light in cm/s
k = 1.38e-16  # Boltzmann constant in erg/K
r0 = 2.82e-13  # Classical electron radius in cm
## integating gamma--> nu_i ,nu is parameter

def integrand_G(gamma,nu_i,nu):
    N_e = N0 * gamma**(-3) 
    N_ph = 8 * np.pi * nu_i**2 / (c**3 * (np.exp(h*nu_i/(k*T_ph))-1) ) 

    x = nu / (4 * gamma**2 * nu_i)
    const = 8 * np.pi * r0**2 * c * h
    if 0 < x < 1:
        fx = 2 * x**2 * np.log(x) + x  + x**2 - 2*x**3
        result = const * fx * N_e * N_ph
        return result
    else:
        return 0

def range_gamma(nu_i,nu):
    return [0.5* (nu / nu_i)**0.5,np.inf]

range_nui = [1e13,1e16]

def j_nu(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,))
    return integrand,err
    
T_ph =5e4
nu01 = np.logspace(13,16,50)
gamma1 = 10
nu_values =4/3*gamma1**2 * nu01  

j_nu_values = [j_nu(nu) for nu in nu_values]

# new stuff

def plotit(jnuval):

    results, errors = zip(*jnuval)
    relative_errors = [e/r for (r, e) in jnuval]

    fig, (ax1, ax2) = plt.subplots(2, 1)
    ax1.plot(results, label='results')
    ax1.plot(errors, label='errors')
    ax1.set_yscale('log')
    ax1.legend()
    ax1.set_title('result and error')
    ax2.plot(relative_errors)
    ax2.set_yscale('log')
    ax2.set_title('relative errors')

    plt.show()

plotit(j_nu_values)

def j_nu_2(nu):
    integrand, err = nquad(integrand_G,[range_gamma,range_nui],args=(nu,),
                           opts={'epsrel': 1E-16, 'epsabs': 1E-16, 'limit': 500,
                                 'maxp1': 500, 'limlst': 500})
    return integrand,err

j_nu_values_2 = [j_nu_2(nu) for nu in nu_values]

plotit(j_nu_values_2)

> **Comments:**
> **username_0**: Your reply was very helpful to me, thank you very much for that.  
> **username_3**: @Shayuxiao that's great! You don't need to accept this answer yet - why don't you wait until you've got your calculation running well and add an answer describing your final solution. This is just a guess - I don't know for sure if it's roundoff or not. Accepting quickly may discourage other answers that may turn out to be even more helpful!  

---

