1. Integrales Dobles sobre área rectangular
Referencia: Chapra 21.5 p643, Burden 4.9 p174
Un ejemplo sencillo de integral doble es una función sobre un área rectangular. Las técnicas revisadas en las secciones anteriores de pueden reutilizar para la aproximación de integrales múltiples.
Considere la integral doble
\int_R \int f(x,y) \delta Asiendo
R = {(x,y) |
x0 ≤ x ≤ xn,
y0 ≤ y ≤ yn }
que es una región rectangular en el plano.
\int_{x_o}^{x_n}\int_{y_o}^{y_n} f(x,y) \delta y \delta x \int_{y_o}^{y_n} \int_{x_o}^{x_n} f(x,y) \delta x \delta y2. Ejercicio
Referencia: 3Eva2018TII_T1 Integral doble con Cuadratura de Gauss
Aproxime el resultado de la integral doble (en área rectangular)
\displaystyle \int_{0}^{\pi/4} \displaystyle\int_{0}^{\pi/4} \Big( 2y \sin(x) + \cos ^2 (x) \Big) \delta y \delta xConsidere los tramosx = 5 como tamaños de paso en eje x , tramosy = 4 como tamaños de paso en eje y.
3. Desarrollo Analítico
Use las fórmulas compuestas de Simpson, haciendo primero constante un valor en y, luego propagando el proceso. Al obtener los integrales en un sentido, repetir el proceso con los resultados en el otro eje.
j: 0 , yj[0]: 0.0
Fórmulas compuestas, tramos: 5
métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
tramos iguales en xi: [3 0 0 2 0 0]
Simpson 3/8 : 0.4378989163640264
Simpson 1/3 : 0.20482799857900125
j: 1 , yj[1]: 0.19634954084936207
Fórmulas compuestas, tramos: 5
métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
tramos iguales en xi: [3 0 0 2 0 0]
Simpson 3/8 : 0.4807008818748735
Simpson 1/3 : 0.2770455037573468
j: 2 , yj[2]: 0.39269908169872414
Fórmulas compuestas, tramos: 5
métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
tramos iguales en xi: [3 0 0 2 0 0]
Simpson 3/8 : 0.5235028473857204
Simpson 1/3 : 0.34926300893569256
j: 3 , yj[3]: 0.5890486225480862
Fórmulas compuestas, tramos: 5
métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
tramos iguales en xi: [3 0 0 2 0 0]
Simpson 3/8 : 0.5663048128965673
Simpson 1/3 : 0.42148051411403814
j: 4 , yj[4]: 0.7853981633974483
Fórmulas compuestas, tramos: 5
métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
tramos iguales en xi: [3 0 0 2 0 0]
Simpson 3/8 : 0.6091067784074142
Simpson 1/3 : 0.4936980192923838
Fórmulas compuestas, tramos: 4
métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
tramos iguales en xi: [3 0 0 1 0]
Simpson 3/8 : 0.4802254950852897
trapecio : 0.20524320554554915
xi [0. 0.1571 0.3142 0.4712 0.6283 0.7854]
yj [0. 0.1963 0.3927 0.589 0.7854]
f(x,y)
[[1. 0.9755 0.9045 0.7939 0.6545 0.5 ]
[1. 1.037 1.0259 0.9722 0.8853 0.7777]
[1. 1.0984 1.1472 1.1505 1.1162 1.0554]
[1. 1.1598 1.2686 1.3287 1.347 1.333 ]
[1. 1.2213 1.3899 1.507 1.5778 1.6107]]
<strong>Integral Volumen: 0.6854</strong>687006308389
4. Algoritmo en Python
El punto de partida es usar las Fórmulas compuestas con Δx segmentos desiguales, donde el algoritmo selecciona aplicar Simpson 3/8, Simpson 1/3 y trapecios según se cumplan los requisitos.
Instrucciones en Python
# 3Eva_IIT2018_T1 Integral doble con Cuadratura de Gauss
# Integrales de volumen o dobles
import numpy as np
# INGRESO
fxy = lambda x,y: 2*y*np.sin(x)+(np.cos(x))**2
x0 = 0 # intervalo x
xn = np.pi/4
y0 = 0 # intervalo y
yn = np.pi/4
tramosx = 5 # tramos por eje
tramosy = 4
titulo = 'f(x,y)'
precision = 4 # decimales a mostrar
# PROCEDIMIENTO
xi = np.linspace(x0,xn,tramosx+1)
yj = np.linspace(y0,yn,tramosy+1)
n = len(xi)
m = len(yj)
hx = (xn-x0)/tramosx # tamaños de paso
hy = (yn-y0)/tramosy
Xi,Yj = np.meshgrid(xi,yj)
# Matriz u[xi,yj] para f(x,y)
F = np.zeros(shape=(n,m),dtype=float)
F = fxy(Xi,Yj)
def simpson_compuesto(xi,fi,vertabla=False,casicero=1e-7):
'''Método compuesto Simpson 3/8, Simpson 1/3 y trapecio
salida: integral,cuenta_h
'''
# vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# Método compuesto Simpson 3/8, Simpson 1/3 y trapecio
cuenta_h = (-1)*np.ones(n,dtype=int) # sin usar
suma = 0 # integral total
suma_parcial =[] # integral por cada método
i = 0
while i<(n-1): # i<tramos, al menos un tramo
tramo_usado = False
h = xi[i+1]-xi[i] # tamaño de paso, supone constante
# tres tramos iguales
if (tramo_usado==False) and (i+3)<n:
d2x = np.diff(xi[i:i+4],2) # diferencias entre tramos
errado = np.max(np.abs(d2x))
if errado<casicero: # Simpson 3/8
S38 = (3/8)*h*(fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3])
suma = suma + S38
cuenta_h[i:i+3] = [3,0,0] # Simpson 3/8
suma_parcial.append(['Simpson 3/8',S38])
i = i+3
tramo_usado = True
# dos tramos iguales
if (tramo_usado==False) and (i+2)<n:
d2x = np.diff(xi[i:i+3],2) # diferencias entre tramos
errado = np.max(np.abs(d2x))
if errado<casicero: # Simpson 1/3
S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
suma = suma + S13
cuenta_h[i:i+2] = [2,0]
suma_parcial.append(['Simpson 1/3',S13])
i = i+2
tramo_usado = True
# un tramo igual
if (tramo_usado == False) and (i+1)<n:
trapecio = (h/2)*(fi[i]+fi[i+1])
suma = suma + trapecio
cuenta_h[i:i+1] = [1] # usar trapecio
suma_parcial.append(['trapecio',trapecio])
i = i+1
tramo_usado = True
cuenta_h[n-1] = 0 # ultima casilla
if vertabla==True: #mostrar datos parciales
print('Fórmulas compuestas, tramos:',n-1)
print('métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado')
print('tramos iguales en xi:',cuenta_h)
for unparcial in suma_parcial:
print('',unparcial[0],':',unparcial[1])
return(suma,cuenta_h)
# PROGRAMA --------------
areaj = []
cuenta_hj = []
j = 0 # y[0] constante
for j in range(0,tramosy+1,1): # cada yj
print('j:',j,', yj['+str(j)+']:',yj[j])
xj = Xi[j]
fj = F[j]
area,cuenta_h=simpson_compuesto(xj,fj,
vertabla=True,
casicero=1e-7)
areaj.append(area)
cuenta_hj.append(cuenta_h)
area,cuenta_h = simpson_compuesto(yj,areaj,
vertabla=True,
casicero=1e-7)
# SALIDA
np.set_printoptions(precision)
print('xi',xi)
print('yj',yj)
print(titulo)
print(F)
print('Integral Volumen:', area)
4. Gráfica de malla
Permite observar la superficie de f(x,y) sobre la que se calculará el volumen respecto al plano x,y en cero.
Instrucciones adicionales al algoritmo
# GRAFICA 3D de malla ------
import matplotlib.pyplot as plt
# Malla para cada eje X,Y
fig3D = plt.figure()
graf3D = fig3D.add_subplot(111, projection='3d')
graf3D.plot_wireframe(Xi,Yj,F,label='f(x,y)')
graf3D.plot(0,0,0,'ro',label='[0,0,0]')
graf3D.set_title(titulo)
graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_title(titulo)
graf3D.legend()
graf3D.view_init(35, -45)
plt.tight_layout()
#plt.show()
5. Gráfica de barras
Realizado para observar las barras con lados de segmentos dx=hx, dy=hy.
Instrucciones complementarias a la gráfica anterior.
# Grafica barras 3D ------
factor = 0.95 # barra proporcion dx, dy
fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111,projection='3d')
F_min = np.min([np.min(F),0])
F_max = np.max([np.max(F)])
color_F = F[:m-1,:n-1].ravel()/F_max
color_ij = plt.cm.rainbow(color_F)
graf_3D.bar3d(Xi[:m-1,:n-1].ravel(),
Yj[:m-1,:n-1].ravel(),
np.zeros((n-1)*(m-1)),
hx*factor,hy*factor,
F[:m-1,:n-1].ravel(),
color=color_ij)
graf_3D.plot(0,0,0,'ro',label='[0,0,0]')
# escala eje z
graf_3D.set_zlim(F_min,F_max)
# etiqueta
graf_3D.set_xlabel('xi')
graf_3D.set_ylabel('yj')
graf_3D.set_zlabel('z')
graf_3D.set_title(titulo)
graf_3D.legend()
# Barra de color
color_escala = plt.Normalize(F_min,F_max)
color_barra = plt.cm.ScalarMappable(norm=color_escala,
cmap=plt.colormaps["rainbow"])
fig_3D.colorbar(color_barra,ax=graf_3D,label="F[xi,yj]")
graf_3D.view_init(35, -45)
plt.tight_layout()
plt.show()