CODEWITHSUNDEEP
pythonfunctionrecursion
user1162809
user1162809

Reputation: 67

Recursive integration

This is a simplified version of the actual question I am working on.

Suppose that we have a sequence of (mathematical) functions, f_0, f_1, f_2, ..., such that enter image description here

We fix a function f_0, say, f_0(x)=x. In this simple case, f_n(x) = 1/((n+1)!) x^{n+1}. Is there an elegant and efficient way to program this recursively, so that Python returns the function f_n given an arbitrary function f_0?

I started out trying to return f_2, but this already failed:

from scipy.integrate import quad

f_0=lambda x: x
f_1=lambda x: quad(f_0,0,x)
f_2=lambda x: quad(f_1,0,x)

returns the error

error: Supplied function does not return a valid float.

after trying to evaluate e.g. f_2(3).

Upvotes: 3

Views: 1609

Answers (2)

MB-F
MB-F

Reputation: 23637

There is nothing wrong with chaining numeric integrals. The problem is what quad returns:

Returns:
y : float The integral of func from a to b.

abserr : float An estimate of the absolute error in the result.

So you get two return values but only the first one is interesting (unless you want to propagate the error over all integrals, but I don't know how to do that). The second integral complains because the function returns a tuple of two values instead of only one scalar.

The following small modification will fix the error by selecting the first of quad's return values:

from scipy.integrate import quad

f_0=lambda x: x
f_1=lambda x: quad(f_0,0,x)[0]
f_2=lambda x: quad(f_1,0,x)[0]

For the sake of completeness, here is a recursive definition for iterating n times:

def integrate_a_lot(f, n):
    for _ in range(n):
        f = lambda x, f=f: quad(f, 0, x)[0]
    return f

f_2 = integrate_a_lot(f_0, 2)
f_42 = integrate_a_lot(f_0, 42)

Upvotes: 3

Mohammad Athar
Mohammad Athar

Reputation: 1980

you can use sympy: http://docs.sympy.org/dev/modules/integrals/integrals.html something like

x = sympy.Symbol('x')

def f_0(m_x):
    return m_x

def f_n(f_0,n,m_x):
    if n==0:
        return f_0(m_x)
    return sympy.integrate(f_n(f_0(m_x),n-1,m_x))

then call it with f_n(f_0,1,x)

Upvotes: 0

Related Questions