CODEWITHSUNDEEP
pythonscipyequationnumerical-integration
Illa Rivero Losada
Illa Rivero Losada

Reputation: 151

How to solve an integral equation in python?

I'm trying to solve this integral equation using Python:

enter image description here

where z ranges from 0 to 1.

The scipy.quad function only provides the numerical solution for a certain interval, but it doesn't provide the solution over the interval.

def f(z,Om,Ol): return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)
quad(lambda r:f(r,Om,Ol),0,1)
(0.77142706642781111, 8.5645609096719596e-15)

But I don't know how to get a full vector in this interval, as you get when using scipy.odeint to solve a differential equation.

In the other hand, sympy.integrate can't do it. I get a stack overflow. Plus, I can't figure out how to substitute the symbols by a list,i.e.:

sy.integrate(x**2,x).subs(x,1)
1/3
sy.integrate(x**2,x).subs(x,[1,2])
TypeError: unhashable type: 'list'

So the question is: does anyone know how to solve an integral equation using python?

Upvotes: 3

Views: 22257

Answers (4)

Illa Rivero Losada
Illa Rivero Losada

Reputation: 151

I finally used the "quad" function with a for statement to solve my problem: 

import pylab as p
import scipy as s
from scipy.integrate import odeint,quad

def Dl_lflrw(z,Om,Ol):
    c = 1.
    H0 = 1.
    y = []
    for i in z:
        y1 = (c/H0)*(1+i)*quad(lambda r:f(r,Om,Ol),0,i)[0]
        y.append(y1)
    return p.asarray(y)

def f(z,Om,Ol):
    return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)

Thank you all for your ideas, they was really helpful.

Upvotes: 0

rocksportrocker
rocksportrocker

Reputation: 7419

I suppose that the z before the integral is a typo which should be z1, and you are looking for z1 given DL.

First you have to implement the right hand side (rhs):

def f(z,Om,Ol): 
        return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)
def rhs(z1, Om, Ol, c, H0):
        return c/H0*(1+z1)*quad(lambda r:f(r, Om, Ol), 0, z1)[0]

Now you have to find a z0 such that rhs(z1, ...) = DL, which is the same as

rhs(z1, ...) - DL = 0

Which means that your problem is reduced to finding the zero (there is only one, because rhs is monotone in), of

f(z1) = rhs(z1, ...) - DL

Here you can apply many methods for root finding (see http://en.wikipedia.org/wiki/Root-finding_algorithm) from http://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding

Upvotes: 1

Eric O. Lebigot
Eric O. Lebigot

Reputation: 94575

I understand that you want to solve a differential equation dF/dz1 = f(z1, Om, Ol) and want F(z1) at different locations. If this is the case, then the Ordinary Differential Equation (ODE) routines of SciPy are the way to go. You might want to check odeint(), in particular, as it can give you the values of your integral at locations that you specify.

Upvotes: 4

Profane
Profane

Reputation: 1128

In the sympy examples section at, http://docs.sympy.org/0.7.1/modules/integrals.html, they show solutions to nearly identical problems. Can you post your sympy code?

For scipy, did you try using a tuple, which is hashable, instead of a list? e.g.:

sy.integrate(x**2,x).subs(x,(1,2,))

Upvotes: 0

Related Questions