Program to calculate Double Integration (original) (raw)

Last Updated : 11 Jul, 2025

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\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)
  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.

`