Search code examples
pythonvariablesscipyintegralquad

Integral with variable limits in python


I'm trying to compute the following integral

enter image description here

using scipy, with the following program:

def E(z):
    result = 1/np.sqrt(Om*(1 + z)**3 + Ode + Ox*(1 + z)**2)
    return result 

def r(z, E): 
    result, error  = quad(E, 0, z) # integrate E(z) from 0 to z
    return result

z being the independent variable, while Om Ode and Ox are simply constants (previously assigned). When I then try to call the function:

z = np.linspace(1e-3, 4, 300)
plt.plot(z, r(z))

I get the error

flip, a, b = b < a, min(a, b), max(a, b)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() 
or a.all()

What is the problem? Is scipy.quad unable to integrate up to a variable? Thank you so much for your help


Solution

  • You could use a combination of Python's map(function, iterable, ...) function,

    Return an iterator that applies function to every item of iterable, yielding the results.

    and functools partial(func[,*args][, **keywords]) method:

    Return a new partial object which when called will behave like func called with the positional arguments args and keyword arguments keywords.

    import numpy as np
    from scipy.integrate import quad
    import matplotlib.pyplot as plt
    from functools import partial
    
    
    def E(z):
        Om = 0.32
        Ode = 0.68
        Ox = 0
        result = 1/np.sqrt(Om*(1 + z)**3 + Ode + Ox*(1 + z)**2)
    
        return result
    
    
    def r(z):
        result = np.array(
            list(map(partial(quad, E, 0), z))
        )[:, 0]  # integrate E(z) from 0 to z
    
        return result
    
    
    z = np.linspace(1e-3, 4, 300)
    fig, ax = plt.subplots()
    ax.plot(z, r(z))
    fig.show()
    

    enter image description here