# Program to calculate Double Integration

Write a program to calculate double integral numerically.

Example:

Input: Given the following integral. where Output: 3.915905


## Recommended: Please try your approach on {IDE} first, before moving on to the solution.

Explanation and Approach:

1. We need to decide what method we are going to use to solve the integral.
In this example, we are going to use Simpson 1/3 method for both x and y integration.
To do so, first, we need to decide the step size. Let h be the step size for integration with respect to x and k be the step size for integration with respect to y.
We are taking h=0.1 and k=0.15 in this example.
Refer for Simpson 1/3 rule
2. We need to create a table which consists of the value of function f(x, y) for all possible combination of all x and y points.
 x\y y0 y1 y2 …. ym x0 f(x0, y0) f(x0, y1) f(x0, y2) …. f(x0, ym) x1 f(x1, y0) f(x1, y1) f(x1, y2) …. f(x1, ym) x2 f(x2, y0) f(x2, y1) f(x2, y2) …. f(x2, ym) x3 f(x3, y0) f(x3, y1) f(x3, y2) …. f(x3, ym) …. …. …. …. …. …. …. …. …. …. …. …. xn f(xn, y0) f(xn, y1) f(xn, y2) …. f(xn, ym)

In the given problem,

x0=2.3
x2=2.4
x3=3.5

y0=3.7
y1=3.85
y2=4
y3=4.15
y4=4.3

3. After generating the table, we apply Simpson 1/3 rule (or whatever rule is asked in the problem) on each row of the table to find integral wrt y at each x and store the values in an array ax[].
4. We again apply Simpson 1/3 rule(or whatever rule asked) on the values of array ax[] to calculate the integral wrt x.

Below is the implementation of the above code:

## C++

 // C++ program to calculate  // double integral value     #include  using namespace std;     // Change the function according to your need  float givenFunction(float x, float y)  {      return pow(pow(x, 4) + pow(y, 5), 0.5);  }     // Function to find the double integral value  float doubleIntegral(float h, float k,                       float lx, float ux,                       float ly, float uy)  {      int nx, ny;         // z stores the table      // ax[] stores the integral wrt y      // for all x points considered      float z, ax, answer;         // Calculating the numner of points      // in x and y integral      nx = (ux - lx) / h + 1;      ny = (uy - ly) / k + 1;         // Calculating the values of the table      for (int i = 0; i < nx; ++i) {          for (int j = 0; j < ny; ++j) {              z[i][j] = givenFunction(lx + i * h,                                      ly + j * k);          }      }         // Calculating the integral value      // wrt y at each point for x      for (int i = 0; i < nx; ++i) {          ax[i] = 0;          for (int j = 0; j < ny; ++j) {              if (j == 0 || j == ny - 1)                  ax[i] += z[i][j];              else if (j % 2 == 0)                  ax[i] += 2 * z[i][j];              else                 ax[i] += 4 * z[i][j];          }          ax[i] *= (k / 3);      }         answer = 0;         // Calculating the final integral value      // using the integral obtained in the above step      for (int i = 0; i < nx; ++i) {          if (i == 0 || i == nx - 1)              answer += ax[i];          else if (i % 2 == 0)              answer += 2 * ax[i];          else             answer += 4 * ax[i];      }      answer *= (h / 3);         return answer;  }     // Driver Code  int main()  {      // lx and ux are upper and lower limit of x integral      // ly and uy are upper and lower limit of y integral      // h is the step size for integration wrt x      // k is the step size for integration wrt y      float h, k, lx, ux, ly, uy;         lx = 2.3, ux = 2.5, ly = 3.7,      uy = 4.3, h = 0.1, k = 0.15;         printf("%f", doubleIntegral(h, k, lx, ux, ly, uy));      return 0;  }

## Java

 // Java program to calculate  // double integral value  class GFG{        // Change the function according to your need  static float givenFunction(float x, float y)  {      return (float) Math.pow(Math.pow(x, 4) +                           Math.pow(y, 5), 0.5);  }     // Function to find the double integral value  static float doubleIntegral(float h, float k,                      float lx, float ux,                      float ly, float uy)  {      int nx, ny;         // z stores the table      // ax[] stores the integral wrt y      // for all x points considered      float z[][] = new float, ax[] = new float, answer;         // Calculating the numner of points      // in x and y integral      nx = (int) ((ux - lx) / h + 1);      ny = (int) ((uy - ly) / k + 1);         // Calculating the values of the table      for (int i = 0; i < nx; ++i)      {          for (int j = 0; j < ny; ++j)           {              z[i][j] = givenFunction(lx + i * h,                                      ly + j * k);          }      }         // Calculating the integral value      // wrt y at each point for x      for (int i = 0; i < nx; ++i)       {          ax[i] = 0;          for (int j = 0; j < ny; ++j)          {              if (j == 0 || j == ny - 1)                  ax[i] += z[i][j];              else if (j % 2 == 0)                  ax[i] += 2 * z[i][j];              else                 ax[i] += 4 * z[i][j];          }          ax[i] *= (k / 3);      }         answer = 0;         // Calculating the final integral value      // using the integral obtained in the above step      for (int i = 0; i < nx; ++i)      {          if (i == 0 || i == nx - 1)              answer += ax[i];          else if (i % 2 == 0)              answer += 2 * ax[i];          else             answer += 4 * ax[i];      }      answer *= (h / 3);         return answer;  }     // Driver Code  public static void main(String[] args)  {      // lx and ux are upper and lower limit of x integral      // ly and uy are upper and lower limit of y integral      // h is the step size for integration wrt x      // k is the step size for integration wrt y      float h, k, lx, ux, ly, uy;         lx = (float) 2.3; ux = (float) 2.5; ly = (float) 3.7;      uy = (float) 4.3; h = (float) 0.1; k = (float) 0.15;         System.out.printf("%f", doubleIntegral(h, k, lx, ux, ly, uy));  }  }     /* This code contributed by PrinciRaj1992 */

## Python3

 # Python3 program to calculate   # double integral value      # Change the function according  # to your need   def givenFunction(x, y):          return pow(pow(x, 4) + pow(y, 5), 0.5)      # Function to find the double integral value   def doubleIntegral(h, k, lx, ux, ly, uy):          # z stores the table       # ax[] stores the integral wrt y       # for all x points considered       z = [[None for i in range(50)]                 for j in range(50)]      ax = [None] * 50        # Calculating the numner of points       # in x and y integral       nx = round((ux - lx) / h + 1)       ny = round((uy - ly) / k + 1)         # Calculating the values of the table       for i in range(0, nx):           for j in range(0, ny):               z[i][j] = givenFunction(lx + i * h,                                       ly + j * k)                  # Calculating the integral value       # wrt y at each point for x       for i in range(0, nx):           ax[i] = 0         for j in range(0, ny):                              if j == 0 or j == ny - 1:                   ax[i] += z[i][j]               elif j % 2 == 0:                  ax[i] += 2 * z[i][j]               else:                  ax[i] += 4 * z[i][j]                      ax[i] *= (k / 3)              answer = 0        # Calculating the final integral       # value using the integral       # obtained in the above step       for i in range(0, nx):           if i == 0 or i == nx - 1:               answer += ax[i]           elif i % 2 == 0:              answer += 2 * ax[i]           else:              answer += 4 * ax[i]              answer *= (h / 3)          return answer      # Driver Code   if __name__ == "__main__":         # lx and ux are upper and lower limit of x integral       # ly and uy are upper and lower limit of y integral       # h is the step size for integration wrt x       # k is the step size for integration wrt y       lx, ux, ly = 2.3, 2.5, 3.7     uy, h, k = 4.3, 0.1, 0.15        print(round(doubleIntegral(h, k, lx, ux, ly, uy), 6))          # This code is contributed   # by Rituraj Jain

## C#

 // C# program to calculate  // double integral value  using System;     class GFG  {     // Change the function according to your need  static float givenFunction(float x, float y)  {          return (float) Math.Pow(Math.Pow(x, 4) +                               Math.Pow(y, 5), 0.5);      }             // Function to find the double integral value      static float doubleIntegral(float h, float k,                          float lx, float ux,                          float ly, float uy)      {          int nx, ny;                 // z stores the table          // ax[] stores the integral wrt y          // for all x points considered          float [, ] z = new float[50, 50];          float [] ax = new float;          float answer;                 // Calculating the numner of points          // in x and y integral          nx = (int) ((ux - lx) / h + 1);          ny = (int) ((uy - ly) / k + 1);                 // Calculating the values of the table          for (int i = 0; i < nx; ++i)          {              for (int j = 0; j < ny; ++j)               {                  z[i, j] = givenFunction(lx + i * h,                                          ly + j * k);              }          }                 // Calculating the integral value          // wrt y at each point for x          for (int i = 0; i < nx; ++i)           {              ax[i] = 0;              for (int j = 0; j < ny; ++j)              {                  if (j == 0 || j == ny - 1)                      ax[i] += z[i, j];                  else if (j % 2 == 0)                      ax[i] += 2 * z[i, j];                  else                     ax[i] += 4 * z[i, j];              }              ax[i] *= (k / 3);          }                 answer = 0;                 // Calculating the final integral value          // using the integral obtained in the above step          for (int i = 0; i < nx; ++i)          {              if (i == 0 || i == nx - 1)                  answer += ax[i];              else if (i % 2 == 0)                  answer += 2 * ax[i];              else                 answer += 4 * ax[i];          }          answer *= (h / 3);                 return answer;      }             // Driver Code      public static void Main()      {          // lx and ux are upper and lower limit of x integral          // ly and uy are upper and lower limit of y integral          // h is the step size for integration wrt x          // k is the step size for integration wrt y          float h, k, lx, ux, ly, uy;                 lx = (float) 2.3; ux = (float) 2.5; ly = (float) 3.7;          uy = (float) 4.3; h = (float) 0.1; k = (float) 0.15;                 Console.WriteLine(doubleIntegral(h, k, lx, ux, ly, uy));      }  }     // This code contributed by ihritik

Output:

3.915905


Attention reader! Don’t stop learning now. Get hold of all the important DSA concepts with the DSA Self Paced Course at a student-friendly price and become industry ready.

My Personal Notes arrow_drop_up Check out this Author's contributed articles.

If you like GeeksforGeeks and would like to contribute, you can also write an article using contribute.geeksforgeeks.org or mail your article to contribute@geeksforgeeks.org. See your article appearing on the GeeksforGeeks main page and help other Geeks.

Please Improve this article if you find anything incorrect by clicking on the "Improve Article" button below.