Adaptive Simpsons Quadrature Algorithm for Double Integrals?

Question

I'm currently using Numerical Analysis 10th edition by Richard L Burden as a reference for approximate Integration techniques. In there it describes the Adaptive Simpsons Quadrature rule that inputs only the bounds and an error tolerance, and spits out the approximate integral within precision of the error tolerance. This method is much more effective than the standard Simpsons rule where you have to input number of iterations and not know how close it is to the actual solution. However, the book goes on to describe a method for Double Integrals using Simpson's rule, but not an algorithm Adaptive Simpsons Quadrature rule for double integrals. Does anyone know a pseudo algorithm for an Adaptive Simpsons rule for double integrals??

For reference, this is the pseudo algorithm for Composite Simpsons rule for single integrals: Inputs bounds (a, b) and n # of iterations

 `NAME: compositeSimpsons(a, b, n):
h=(b-a)/n
first = f(a)
last = f(b)
sum=0
x = a+h
for(i=2:n-1)
   if(i%2==0) // even
      sum+=4*(x)
   else // odd
      sum+=2*f(x) 
x+=h
end for
return (h/3) * (first+sum+last)`

And here is the pseudo-algorithm for Adaptive Simpsons Quadrature for single integrals: (Input bounds a, b) and tolerance (tol)

`NAME: adaptiveQuadratureSimspons(a, b, tol):
myStack.push(a)
myStack.push(b)
I=0
while(myStack is not empty)
    bb = myStack.pop()
    aa = myStack.pop()
    I1 = compositeSimpsons(aa, bb, 2)
    m = (aa+bb)/2
    I2 = compositeSimpsons(aa, mm, 2) + compositeSimspons(mm, bb, 2)
    if(|I2-I1|/15 < (bb-aa)*tol)
        I += I2
    else
        myStack.push(m)
        myStack.push(bb)
        myStack.push(aa)
    myStackl.push(m)
end while
return I`

The algorithm for Simpsons rule for two integrals gets very complex fast as you're replacing the x variable with each iteration with a different subdivision, so I won't detail it here unless necessary. However, I know that the problem isn't that algorithm as I've tried it many times and works fine for many different double integral problems. I tried to use the same logic found in the adaptive Simpsons rule my double integral adaptive Simpsons rule by replacing compositeSimpsons() with my compositeSimpsonsDouble(), but it entered an infinite loop as the difference between I2 and I1 was always less than the tolerance. Any help? Coding this in Java

Nico Schlömer Nico Schlömer · Accepted Answer · 2019-12-24T10:38:21

In the lingo of numerical quadrature, "double integrals" don't play as big as a role as the domain you want to integrate your function over. In 1D it's always an interval, in 2D it can be a disk, a rectangle, a triangle, the plane with weight function exp(-r**2) etc. Perhaps your double integral is one of these. For all these different domains, you have different integration techniques. See https://github.com/nschloe/quadpy for some examples.

For adaptive quadrature in 2D, my first impulse would be to check if the domain can be approximated well by a number of triangles. Like intervals in 1D, those can be easily split into smaller triangles if the error estimator recommends so.

Check https://github.com/nschloe/quadpy/wiki/Adaptive-quadrature for how to do this with quadpy.

Adaptive Simpsons Quadrature Algorithm for Double Integrals?

1 Answers