Evaluate double integrals using Simpsons composite rule

Question

I have previously created this code to evaluate an integral using Simpsons composite rule, and now need to modify it to evaluate double integrals. Does anyone have any idea how to go about this? I would be grateful for any help..

#include <iostream>
#include <cmath>

using namespace std;

float f(float x); //returns the function we are integrating
void simpsons_rule(float a, float b, int n); // Carries out Simpsons composite rule


int main()
{
    cout << "This program finds an approximation of the integral of a function under two limits using Simpsons composite rule." << endl;
    // User is introduced to program and given explanation of what it does
    float a,b;
    int n;
    cout << "Enter the lower limit " << endl;
    cin >> a;
    cout << "Enter the upper limit " << endl;
    cin >> b; // User has freedom to specify the upper and lower limits under which the function will be integrated
    cout << "Enter the number of intervals (must be an even number) " << endl; 

    do
    {
     cin >> n; // User can also choose the number of intervals

    if(n%2 == 0 && n!=0) // If an even number of intervals was entered, the program continues and simpsons rule is carried out
    {
        simpsons_rule(a,b,n);
    }
    else
    {
        cout << "Invalid number of intervals, please re-enter a value" << endl; //If an odd number of intervals was entered, the user is prompted to enter a valid number
    }
    }while(cin); // do-while loop used to ensure even number of intervals was entered

    return 0;
}

float f(float x)
{
    return(pow(x,4) - 3*pow(x,3) + pow(x,2) + x + 1); // function to be integrated
}

void simpsons_rule(float a, float b, int n)
{
    float s2 = 0;
    float s3 = 0;
    int i;
    float h = (b-a)/n; // step length is calculated

      for(i=1; i <= (n/2)-1; i++)
      {

          if((i%2) == 0) // checks if i is even
          {
            s2 = s2 + f(a+(i*h)); // Finds the sum of all the values of f(x) where x is even

          }
      }

      for(i=1; i<=(n/2); i++)
      {
         if((i%2) == 1) // checks if i is odd
         {
            s3 = s3 + f(a+2*(i*h)); // Finds the sum of all the values of f(x) where x is odd
         }
      }

    float s = (h/3)*(f(a)+ f(b) + (2*s2) + (4*s3)); // This variable contains the final approximaton of the integral

    cout << "The value of the integral under the specified limits is: " << s << endl;

Answer 1

Unfortunately Simpson's rule can't be applied directly to multiple integrals. What you need to do is derive interpolant surfaces or hypersurfaces for double or triple integrals, respectively. For a double integral you end up evaluating the function on a grid of nine points instead of the three points you use in the single integral case. For a triple integral you use a 3-D lattice of 27 points. Needless to say, it gets pretty complicated.

An easier approach is some sort of Monte Carlo integration, in which you randomly sample the function many times, take the average of all the function samples, and multiply by the area of integration. The drawback here is that the error is inversely proportional to the square root of the number of samples, so 4 times as many samples only halves your expected error. If you have plenty of time and your accuracy requirements aren't great, this is probably the approach you should try.