Python Solving an Integral equation with uncertainties

[Arman777]

TL;DR Summary

Solving an Integral equation with uncertainties, by using fsolve and uncertainties packages in Python

I have some variables that are uncertain, these are

w_m = u.ufloat(0.1430, 0.0011)
z_rec = u.ufloat(1089.92, 0.25)
theta_srec = u.ufloat(0.0104110,  0.0000031)
r_srec = u.ufloat(144.43, 0.26)

and some constant values

c = 299792.458  # speed of light in [km/s]
N_eff = 3.046
w_r = 2.469 * 10**(-5) * (1 + (7/8)*(4/11)**(4/3) * N_eff)

My problem is I need to solve an integral defined by this functions

def D_zrec(z):
            return (c/100) / sqrt(w_m * (1+z)**3 + w_r * (1+z)**4 + (h**2 - w_m - w_r))

This function is evaluated for $dz$ but it also contains an unknown $h$ that we need to find, with corresponding uncertainty. So I need to write a code that solves h.

Here is my full code

from numpy import sqrt, vectorize
    from scipy.integrate import quad
    import uncertainties as u
    from uncertainties.umath import *
    from scipy.optimize import fsolve

    ####  Important Parameters #####
    c = 299792.458  # speed of light in [km/s]
    N_eff = 3.046
    w_r = 2.469 * 10**(-5) * (1 + (7/8)*(4/11)**(4/3) * N_eff)
    w_m = u.ufloat(0.1430, 0.0011)
    z_rec = u.ufloat(1089.92, 0.25)
    theta_srec = u.ufloat(0.0104110,  0.0000031)
    r_srec = u.ufloat(144.43, 0.26)

    D_zrec_true = r_srec / theta_srec

    @u.wrap
    def D_zrec_finder(h, w_m, z_rec):
        def D_zrec(z):
            return (c/100) / sqrt(w_m * (1+z)**3 + w_r * (1+z)**4 + (h**2 - w_m - w_r))
        D_zrec_test, error = quad(D_zrec, 0, z_rec)
        return D_zrec_true - D_zrec_test    def h0_finder(w_m, z_rec):
        sol = fsolve(D_zrec_finder, 0.67, args=(w_m, z_rec))[0]
        return sol

    print(h0_finder(w_m, z_rec))

So to summarize I have an integral named as D_zrec` that is a function of ` $z$ but also contains an unknown number `$h$ that we need to find by using fsolve.

I have found 3 sites that might be useful for the coder. Please look at them if you want to help

https://kitchingroup.cheme.cmu.edu/blog/2013/03/07/Another-approach-to-error-propagation/

https://kitchingroup.cheme.cmu.edu/blog/2013/07/10/Uncertainty-in-an-integral-equation/

https://kitchingroup.cheme.cmu.edu/blog/2013/01/23/Solving-integral-equations-with-fsolve/

I have looked at them to write my code but no luck.

Thanks for the help

 

[lofgran]

!

Hi there,

Thank you for sharing your problem and code. It seems like you have already done a lot of research and have a good understanding of the problem you are trying to solve. However, I have a few suggestions that may help you find a solution.

First, it may be helpful to break down your problem into smaller parts. For example, you could first try to solve the integral without the uncertainty in h, and then add in the uncertainty later. This will help you understand the code and make it easier to troubleshoot any errors.

Second, I noticed that you are using the `quad` function from `scipy.integrate` to solve the integral. However, this function is not compatible with uncertainties. Instead, you could try using the `uncertainties` package's `ufloat` and `umath` functions to handle the uncertainties in your calculations. You can also use the `vectorize` function to make your function compatible with the uncertainties package.

Finally, I would suggest trying to solve the integral using a different method, such as the `scipy.optimize` package's `fsolve` function. This function allows you to solve for a root of a function, which may be more appropriate for your problem.

I hope these suggestions are helpful. Good luck with your code!