0

I was trying to create a function that solves multiple integrals using Simpson's method and my idea was to integrate out variables one by one until I'm left with a single integral but I'm getting recursion depth exceeded and I'm unable to figure out why (I suspect that's because I'm using too many lambda functions).

def simpson(dim:int, f, dom:list, h=0.01): #dim:dimensions, f:function to be integrated, dom:domain, h:step size

    def integrate(f, interval, h):      #integrates single variable function f using simpson's method 
        integral = 0
        x=interval[0]
        while x<interval[1]:
            integral+=h/6*(f(x)+4*f(x+h/2)+f(x+h))
            x+=h
        return integral

    for d in range(dim-1,-1,-1):            #integrate out variables one by one
        if d==0:
            return integrate(lambda x:f([x]), dom[0], h)
        f = lambda x: integrate(lambda y:f([*x,y]), dom[d], h) #new function that takes one less variable (last variable integrated out)
#test function of two variables (takes input as a vector)   
def T(x):                        
    x,y = x
    return -x**2-2*y**2+2*x*y+2*x+72
print(simpson(2, T, [(0,8),(0,6)])/48)
Akshat Sharma
  • 121
  • 4
  • It's an infinite loop, right? By the time you call `integrate`, `f` is a function that calls `integrate` with a function that calls `f`, so when it calls `f` it goes right back to integrate. You don't have any way to escape from the recursion. – Tim Roberts Apr 03 '21 at 21:16
  • I've checked that it works if `T` is replaced by a function of a single variable, but in that case `f = ...` line is never called. – rpoleski Apr 03 '21 at 21:20
  • @TimRoberts perhaps there is some flaw in my reasoning then? I've been thinking about this for more than 2 hours now but I'm unable to figure out what's wrong. – Akshat Sharma Apr 03 '21 at 21:29
  • @rpoleski Ohh, so I think that suggests the integrate function is fine. – Akshat Sharma Apr 03 '21 at 21:30
  • I'll be honest, my head exploded when I tried to dig in to what you're doing. You mentioned a recursion issue, I tried to find where you break out of the recursion, and I don't see an escape. – Tim Roberts Apr 03 '21 at 21:32
  • 1
    There is one `f()` at a time, no closure is created for the "old" `f()`, what you may expect. It might not be that obvious since you have a single argument, though the endless recursion is a strong indicator. However if you try `f=lambda x,y:x+y` and then `f=lambda x:f(1,x)`, and try invoking the result as `f(2)`, you will immediately see that the "old" `f()` is simply gone, and the new one is just trying to call itself. You may want to make `simpson()` itself recursive, instead of that loop inside. Then there will be closures – tevemadar Apr 03 '21 at 22:06
  • @tevemadar Yes it works in 1D. Yes you're right the original f is never called, I was expecting the original f to be eventually called and all the lambda functions to 'wrap' around that (to create a composite function), but I think that's not what happens. – Akshat Sharma Apr 03 '21 at 22:13

0 Answers0