Programa para calcular Doble Integración

Posted on by Rudeus Greyrat

Escriba un programa para calcular integrales dobles numéricamente.

Ejemplo:

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

Explicación y enfoque:

  1. Necesitamos decidir qué método vamos a usar para resolver la integral.
    En este ejemplo, vamos a utilizar el método de Simpson 1/3 para la integración de x e y.
    Para hacerlo, primero, debemos decidir el tamaño del paso. Sea h el tamaño de paso para la integración con respecto a x y k sea el tamaño de paso para la integración con respecto a y.
    Tomamos h=0.1 y k=0.15 en este ejemplo.
    Referirse a la regla de Simpson 1/3
  2. Necesitamos crear una tabla que consista en el valor de la función f(x, y) para todas las combinaciones posibles de todos los puntos x e y.
    x\y y0 y1 y2 …. mmm
    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)

    En el problema dado,

    x0=2.3
x2=2.4
x3=3.5

y0=3.7
y1=3.85
y2=4
y3=4.15
y4=4.3
  3. Después de generar la tabla, aplicamos la regla de Simpson 1/3 (o cualquier regla que se solicite en el problema) en cada fila de la tabla para encontrar la integral wrt y en cada x y almacenar los valores en una array ax[].
  4. Nuevamente aplicamos la regla de Simpson 1/3 (o cualquier regla que se solicite) en los valores de la array ax[] para calcular la integral wrt x.

A continuación se muestra la implementación del código anterior:

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 
Producción:

3.915905

Publicación traducida automáticamente

Artículo escrito por Archaic y traducido por Barcelona Geeks. The original can be accessed here. Licence: CCBY-SA

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *