Open In App

Program to calculate Double Integration

Improve
Improve
Improve
Like Article
Like
Save Article
Save
Share
Report issue
Report

Write a program to calculate double integral numerically.

Example:

Input: Given the following integral.
\int _{3.7}^{4.3}\int _{2.3}^{2.5}\sqrt{x^4+y^5}\:dxdywheref(x, y)=\sqrt{x^4+y^5}\Output: 3.915905

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\yy0y1y2….ym
x0f(x0, y0)f(x0, y1)f(x0, y2)….f(x0, ym)
x1f(x1, y0)f(x1, y1)f(x1, y2)….f(x1, ym)
x2f(x2, y0)f(x2, y1)f(x2, y2)….f(x2, ym)
x3f(x3, y0)f(x3, y1)f(x3, y2)….f(x3, ym)
….….….….….….
….….….….….….
xnf(xn, y0)f(xn, y1)f(xn, y2)….f(xn, ym)
  1. 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
  1. 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[].
  2. 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 <bits/stdc++.h>
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[50][50], ax[50], answer;
 
    // Calculating the number 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[50][50], ax[] = new float[50], answer;
 
    // Calculating the number 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 number 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[50];
        float answer;
     
        // Calculating the number 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

                    

Javascript

// JavaScript program to calculate
// double integral value
 
// Change the function according to your need
function givenFunction(x, y){
    return Math.pow(Math.pow(x, 4) + Math.pow(y, 5), 0.5);
}
 
// Function to find the double integral value
function doubleIntegral(h, k, lx, ux, ly, uy){
     
    // z stores the table
    // ax[] stores the integral wrt y
    // for all x points considered
    var z = new Array(50);
    for(var i = 0; i < z.length; i++){
        z[i] = new Array(50);
    }
    var ax = new Array(50);
    let answer;
     
    // Calculating the number of points
    // in x and y integral
    let nx = Math.round((ux - lx) / h + 1);
    let ny = Math.round((uy - ly) / k + 1);
     
    // Calculating the values of the table
    for(let i = 0; i < nx; i++){
        for(let 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 (let i = 0; i < nx; ++i)
    {
        ax[i] = 0;
        for (let 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 (let 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;
}
 
// 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
let lx = 2.3, ux = 2.5, ly = 3.7;
let uy = 4.3, h = 0.1, k = 0.15;
 
console.log(doubleIntegral(h, k, lx, ux, ly, uy).toFixed(6));
 
// This code is contributed by lokeshmvs21.

                    
Output:
3.915905


Last Updated : 18 Nov, 2022
Like Article
Save Article
Previous
Next
Share your thoughts in the comments
Similar Reads