Integración Numérica

Utilizando la interpretación geométrica de la Integración de Riemann, definida como el área bajo la curva de la función, es posible aproximar una integral definida utilizando una sumatoria de áreas de rectángulos de la siguiente forma:

\[ \int_{a}^{b} f(x)\,dx \approx \sum_{i=0}^{N-1}f(x_i)\Delta x \approx \sum_{i=1}^{N}f(x_i)\Delta x, \]

donde \(x_i\) corresponden a los puntos dentro del intervalo \([a, b]\) donde se define la integral y \(\Delta x = (b - a)/N\) el ancho de cada rectángulo. Es decir, se calcula el área de cada rectángulo como \(f(x_i)\,\Delta x\) y luego se suman. Un esquema de la idea se presenta a continuación

Se puede notar que entre mayor sea el valor de \(N\) el cálculo de la integral se hace de forma más precisa.

Ejercicios

Parte 1

Utilizando la información entregada anteriormente desarrolle un programa en Python que calcule el valor de la integral definida para las siguientes funciones:

• \(f(x)=\sin(x) + x^2\)
• \(f(x)=x^5 + x^2 + e^{x}\)
• \(f(x)=\sin(x) / x\)

donde \(a\), \(b\) y \(N\) deben ser ingresados por el usuario.

Solución

from math import sin

N = int(input("Ingrese N: "))
a = float(input("Ingrese a: "))
b = float(input("Ingrese b: "))

dx = (b - a) / N
x = a
integral = 0
while x < b:
    integral += (sin(x) + x ** 2) * dx
    x += dx

print(integral)

x = a + dx
integral = 0
while x <= b:
    integral += (sin(x) + x ** 2) * dx
    x += dx

print(integral)

Parte 2

Extienda la misma idea pero para aproximar la siguiente integral doble:

\[ \int_{a}^{b}\int_{c}^{d} f(x,y)\,dy\,dx \approx \sum_{i=0}^{N_x-1}\sum_{j=0}^{N_y-1}f(x_i, y_j)\Delta y \Delta x,\quad \Delta x = \frac{b - a}{N_x}, \quad \Delta y = \frac{d - c}{N_y}. \]

Pruebe su programa calculando la siguiente integral doble

\[\int_{0}^{1}\int_{0}^{1} -(x^2 + y^2)\,dy\,dx\]

sabiendo que el resultado es \(-2/3\).

Pruebe otras funciones y verifique con la siguiente aplicación: Wolfram Alpha

Solución

Nx = int(input("Ingrese Nx: "))
Ny = int(input("Ingrese Ny: "))
a = float(input("Ingrese a: "))
b = float(input("Ingrese b: "))
c = float(input("Ingrese c: "))
d = float(input("Ingrese d: "))

dx = (b - a) / Nx
dy = (d - c) / Ny
x = a
integral = 0
while x < b:
    y = c
    while y < d:
        integral += -(x ** 2 + y ** 2) * dy * dx
        y += dy
    x += dx

print(integral)

integral = 0
x = a + dx
while x <= b:
    y = c + dy
    while y <= d:
        integral += -(x ** 2 + y ** 2) * dy * dx
        y += dy
    x += dx

print(integral)