Recursive integration
This is a simplified version of the actual question I am working on.
Suppose we have a sequence of (mathematical) functions, f_0, f_1, f_2, ..., such that
We fix the 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 a function f_n given an arbitrary function f_0?
I started trying to return f_2, but that failed already:
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 an error
error: Supplied function does not return a valid float.
after trying to evaluate, eg. f_2 (3).
source to share
There is nothing wrong with a chain of numerical integrals. The problem is what returns quad
:
Returns:
y : float Integral of func from a to b.abserr : float Estimate of the absolute error of the result.
So, you get two return values, but only the first one is interesting (unless you want to propagate the error over all the integrals, but I don't know how). The second integral is complaining because the function returns a tuple of two values instead of one scalar.
The following minor modification will fix the error by choosing the first of the quad
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 completeness, here's a recursive definition to iterate 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)
source to share
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)
source to share