Cálculo de volúmenes, mediante sumas de Riemann, en un script de Python

Los ejercicios 15 y 16, pp. 995 del libro de Stewart (Multivariable Calculus), sugieren el uso de un calculador programable, un computador o el comando suma sobre un CAS para estimar dos integrales dobles mediante sumas de Riemann. En ambos casos se pide el uso de la regla del punto medio para muestrear los valores de la función en el caso de 1, 4, 16, 64, 256 y 1024 cuadrados (series 1, 2, 3, 4, 5 y 6; respectivamente) de igual tamaño para la región R = [0,1]x[0,1]. Para realizar esta tarea decidí usar Python.

En términos generales, el código lo que hace es determinar los valores de (x, y) en el centro de cada sub cuadrado, estimar el valor de la función para cada par (x,y), multiplicarlos por su respectivo elemento de área (que sería equivalente a tener el volumen de un prisma de base cuadrada) y sumar luego todos estos elementos para obtener el volumen total en la región R. Las integrales consideradas se encuentran a continuación:

doble2

El código Python para evaluar el primer caso es el siguiente:

#!/usr/bin/env python
# -*- coding: utf-8

from os import system
from math import sqrt,exp

system("clear")

n = int(raw_input("número de serie = ? "))

#Número de filas o columnas
rows = 2**(n - 1)

#Para calcular el valor inicial de x o y
nxy = n

#Número de cuadrados a generar
num_cuad = 2**(2*n-2)

print "número de cuadrados a generar = ", num_cuad
print "número de filas o columnas = ", rows

#Determina el elemento de área

delta_x = 1./rows
delta_y = 1./rows

delta_area = delta_x * delta_y

print "Elemento de área = ", delta_area

#Determina valor inicial de y
y = 1./2**nxy

#Determinación del incremento para x e y
inc= 2*y

volumen = 0

for i in range(rows):

	x = 1./2**nxy #valor inicial de x

	for j in range(rows):

		funcion = sqrt(1+x*exp(-y)) * delta_area

		volumen += funcion

		#print x,y,funcion

		x +=inc 

	y +=inc		

print "volumen = %.5f" % volumen

Su ejecución produce los siguientes valores:

    1.14161, 1.14319, 1.14353, 1.14362, 1.14364 y 1.14364

para 1, 4, 16, 64, 256 y 1024 cuadrados, respectivamente.

Para verificar si el código estaba funcionando adecuadamente se empleó el software Wolfram Mathematica 9 (para Linux) que tiene un comando, NIntegrate, el cual permite la evaluación de integrales múltiples. El resultado es idéntico al del script de Python; tal como se observa en la siguiente imagen:

doble1

A continuación, se da una imagen de la superficie que se está integrando:

integral_doble

Para evaluar la segunda integral, en la misma región, sólo hay que modificar el código anterior para incluir la función seno de la librería math y modificar la línea:

    funcion = sqrt(1+x*exp(-y)) * delta_area

del bucle de j por ésta:

    funcion = sin(x + sqrt(y)) * delta_area

Los valores respectivos de 0.93459, 0.88199, 0.86575, 0.86049, 0.85875 y 0.85816 se producen para 1, 4, 16, 64, 256 y 1024 cuadrados. Sin embargo, la concordancia entre dos valores consecutivos (diferencia < 1×10-5) sólo se alcanza cuando se generan 1.048.576 cuadrados (1024×1024) para un resultado de 0.85785. No obstante, este cálculo se produce en mi sistema en poco más de 2 segundos. El valor de 0.85785 es idéntico al producido por Mathematica y coincidencialmente es también idéntico al obtenido con 256 cuadrados (16×16).

La superficie de integración es la siguiente:

integral_doble2

Se intuye porqué las sumas Riemann requieren una mayor sub división del área para una mejor aproximación. El caso anterior, a pesar de su complicada expresión funcional, da una superficie que se asemeja en mucho a un plano y, además, poco inclinado con relación al plano z = 0.

Esta entrada fue publicada en Código Python, Software Libre. Guarda el enlace permanente.

4 respuestas a Cálculo de volúmenes, mediante sumas de Riemann, en un script de Python

  1. Juanlu001 dijo:

    Esto en la próxima versión de SciPy (al caer) será mucho más fácil:

    >>> from scipy.integrate import nquad
    >>> from numpy import sqrt, exp
    >>> nquad(lambda x, y: sqrt(1 + x * exp(-y)), [(0, 1), (0, 1)])
    (1.1436438408450194, 1.3528376126647705e-14)
    

    el segundo valor es el error cometido. Más sencillo imposible😉

  2. No lo creo. Ve esto:

    http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-4-Matplotlib.ipynb

    Por otra parte, ni siquiera puedo correr tu código de 3 líneas anteriores. Usas Python 3.x porque yo todavía estoy con el 2.6.6?.

    En mi versión de python (puedo correr la 3 pero no reconoce la librería scipy de python 2 que está por defecto) para hacer lo que tu haces el código es el siguiente:

    >>> from scipy.integrate import dblquad 
    >>> from math import sqrt,exp
    >>> dblquad(lambda x,y: sqrt(1 + x * exp(-y)), 0, 1, lambda x:0, lambda x:1)
    (1.1436438408450194, 1.2696997240768885e-14)
    

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s