Error propagation of a variable for an integral

Question

I have an integral that depends on two parameters $a\pm\delta a$ and $b\pm \delta b$. I am doing this integral numerically and no python function can calculate the integral with uncertainties.

So I have calculated the integral for each min, max values of a and b. As a result I have obtained 4 values, such that;

$$(a + \delta a, b + \delta b) = 13827.450210 \pm 0.000015~~(1)$$ $$(a + \delta a, b - \delta b) = 13827.354688 \pm 0.000015~~(2)$$ $$(a - \delta a, b + \delta b) = 13912.521548 \pm 0.000010~~(3)$$ $$(a - \delta a, b - \delta b) = 13912.425467 \pm 0.000010~~(4)$$

So it is clear that $(2)$ gives the min and $(3)$ gives the max. Let us show the result of the integral as $c \pm \delta c$. So my problem is what is $c$ and $\delta c$ here?

The integral is something like this

$$I(a,b,x) =C\int_0^b \frac{dx}{\sqrt{a(1+x)^3 + \eta(1+x)^4 + (\gamma^2 - a - \eta)}}$$

where $\eta$ and $\gamma$ are constant.

Note: You guys can also generalize it by taking $\eta \pm \delta \eta$ but it is not necessary for now.

I have to take derivatives or integrals numerically. There's no known analytical solution for the integral.

$\eta = 4.177 \times 10^{-5}$, $a = 0.1430 \pm 0.0011$, $b = 1089.92 \pm 0.25$, $\gamma = 0.6736 \pm 0.0054$, $C = 2997.92458$

Answer 1

What is inside the square root is $$\gamma ^2+ (3 a+4 \eta )x+ 3( a+2 \eta )x^2+ (a+4 \eta )+\eta x^4\tag 1$$ Write it as $$\eta\, (x-r_1 ) (x-r_2 ) (x-r_3 ) (x-r_4)$$ where the $r_i$ are the roots of the quartic polynomial given in $(1)$ .

So, we need to compute $$I(a,b)=\frac C {\sqrt \eta}\,\int_0^b\frac{dx}{\sqrt{(x-r_1 ) (x-r_2) (x-r_3 ) (x-r_4 )}}$$ and we have an elliptic integral of the first kind (have a look here).

So, now, we can compute all the partial derivatives with repect to $(\eta,r_1,r_2,r_3,r_4)$ and use the chain rule.

So , assuming no cross terms, thz final result write $$I = I_0 +\frac {\partial I}{\partial a} ( a-a_0)+\frac {\partial I}{\partial b} (b-b_0)+\frac {\partial I}{\partial \eta} (\eta-\eta_0)+\frac {\partial I}{\partial \gamma} (\gamma-\gamma_0)$$ with $$I_0=13869.7187382790600280056975524$$ $$\frac {\partial I}{\partial a}=-38667.5002882782982646434723$$ $$\frac {\partial I}{\partial b}=0.1916010843310452774261082$$ $$\frac {\partial I}{\partial \eta}=-1517907.851327789447487779$$ $$\frac {\partial I}{\partial \gamma}=-3984.5811163972118547061439$$

Answer 2

from numpy import sqrt
from scipy import integrate
import uncertainties as u
from uncertainties.umath import *

#Important Parameters
C = 2997.92458  # speed of light in [km/s]
eta = 4.177 * 10**(-5)
a = u.ufloat(0.1430, 0.0011)
b = u.ufloat(1089.92, 0.25)
gama = u.ufloat(0.6736, 0.0054)

@u.wrap
def D_zrec_finder(gama, a, b):
    def D_zrec(z):
        return C / sqrt(a * (1+z)**3 + eta * (1+z)**4 + (gama**2 - a - eta))
    result, error = integrate.quad(D_zrec, 0, b)
    return result


print((D_zrec_finder(gama, a, b)).n)
print((D_zrec_finder(gama, a, b)).s)

This works