Program to calculate Double Integration
Last Updated :
18 Nov, 2022
Improve
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}\:dxdy wheref(x, y)=\sqrt{x^4+y^5}\ Output: 3.915905
Explanation and Approach:
- 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
- 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
- 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[].
- 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++ 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 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 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# 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 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
