Program to calculate Double Integration
Write a program to calculate double integral numerically.
Example:
Input: Given the following integral. whereOutput: 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++
// 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
Contact Us