Related Articles

# Program to calculate Double Integration

• Last Updated : 18 Nov, 2019

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 needfloat givenFunction(float x, float y){    return pow(pow(x, 4) + pow(y, 5), 0.5);}  // Function to find the double integral valuefloat 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 Codeint 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 valueclass GFG{    // Change the function according to your needstatic 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 valuestatic 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 Codepublic 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 valueusing System;  class GFG{  // Change the function according to your needstatic 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. Participate in the Scholorship Test for First-Step-to-DSA Course for Class 9 to 12 students.

My Personal Notes arrow_drop_up