7.4 Cuadratura de Gauss

Referencia: Chapra 22.3 p655 pdf679, Rodriguez 7.3 p294, Burden 4.7 p220 pdf230

La cuadratura de Gauss aproxima el integral de una función en un intervalo [a,b] centrado en cero mediante un cálculo numérico con menos operaciones y evaluaciones de la función. Se representa como una suma ponderada:

I \cong c_0f(x_0) + c_1f(x_1)

para la fórmula de dos puntos se tiene obtiene:

c_0 = c_1 = 1 x_0 = -\frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}

para un intervalo de evaluación desplazado en el eje x se requiere convertir los puntos al nuevo rango. Se desplaza el punto cero al centro del intervalo [a,b] y se obtiene:

x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{b-a}{2}(f(x_a) + f(x_b))

cuya fórmula es semejante a una mejor aproximación de un trapecio, cuyos promedios de alturas son puntos internos de [a,b], concepto mostrado en la gráfica.


Algoritmo

El algoritmo se desarrolla en un tramo en el intervalo [a,b] junto a la gráfica para mostar el concepto. Para el ejemplo el integral buscado es:

\int_1^3 \sqrt{x} \sin(x) dx

que usando cuadratura de Gauss con un solo intervalo[a,b] tiene como resultado:

Integral:  2.03821975687
>>>

el resultado se puede mejorar aumentando el número de tramos en el intervalo [a,b]. Por ejemplo, el resultado usando 4 tramos el resultado es semejante al usar el método del trapecio con 128 tramos, lo que muestra el ahorro en calculos entre los métodos

Integral:  2.05357719003
>>> 

Las intrucciones en Python son:

# Integración: Cuadratura de Gauss de dos puntos
# modelo con varios tramos entre [a,b]
import numpy as np
import matplotlib.pyplot as plt

# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    return(area)

#Función a evaluar
def funcionx(x):
    fxi = np.sqrt(x)*np.sin(x)
    return(fxi)

# PROGRAMA
a = 1 # float(input('a: '))
b = 3 # float(input('b: '))
tramos = 4 # int(input('tramos: '))

# PROCEDIMIENTO
muestras = tramos+1
xi = np.linspace(a,b,muestras)
area = 0
for i in range(0,muestras-1,1):
    deltaA = integraCuadGauss2p(funcionx,xi[i],xi[i+1])
    area = area + deltaA
# SALIDA
print('Integral: ', area)

Gráfica por tramos

La gráfica que complementa el resultado anterior se realiza añadiendo las intrucciones presentadas a continuación.

Considere que la gráfica es útil con pocos tramos en el intervalo[a,b]

# Gráfica por cada subtramo
x0 = -1/np.sqrt(3) 
x1 = 1/np.sqrt(3)
muestrastramo = 20
subtramo = np.linspace(a,b,tramos+1)
xi = np.array([])
fi = np.array([])
xat = np.array([])
xbt = np.array([])
recta = np.array([])
for i in range(0,tramos,1):
    at = subtramo[i]
    bt = subtramo[i+1]
    xit = np.linspace(at,bt,muestrastramo)
    fit = funcionx(xit)
    xi = np.concatenate((xi,xit))
    fi = np.concatenate((fi,fit))
    # puntos xa y xb por tramo
    xa = (bt+at)/2 + (bt-at)/2*(x0)
    xb = (bt+at)/2 + (bt-at)/2*(x1)
    xat = np.concatenate((xat,[xa]))
    xbt = np.concatenate((xbt,[xb]))
    # Recta entre puntos x0 y x1 por tramo
    pendiente = (funcionx(xb)-funcionx(xa))/(xb-xa)
    b0 = funcionx(xa)-pendiente*xa
    linea = b0 + pendiente*xit
    recta = np.concatenate((recta,linea))
# Marcadores 'o' de xa y xb por tramos
puntox = np.concatenate((xat,xbt))
puntoy = funcionx(puntox)
    
# Trazado de lineas
plt.plot(xi,recta, label = 'grado 1', color = 'tab:orange')
plt.fill_between(xi,0,recta, color='tab:olive')
plt.plot(xi,fi, label='f(x)', color = 'b')
# Verticales para dividir los tramos
for j in range(0,len(subtramo),1):
    plt.axvline(subtramo[j], color='tab:gray')
# Marcadores de puntos xa y xb por tramos
for j in range(0,len(xat),1):
    plt.axvline(xat[j], color='w')
    plt.axvline(xbt[j], color='w')
plt.plot(puntox,puntoy, 'o', color='g')
plt.title('Integral: Cuadratura Gauss')
plt.xlabel('x')
plt.legend()
plt.show()

7.3 Regla de Simpson 3/8

Referencia: Chapra 21.2.1 p.631 pdf.655, Burden 4.2 p192 pdf202, Rodriguez 7.1.8 p288

Es el resultado cuando se realiza una interpolación con polinomio de tercer grado.

I = \int_a^b f(x) \delta x \cong \int_a^b f_3 (x) \delta x

Usando un polinomio de Lagrange de tercer grado:

I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)]

siendo
h=\frac{b-a}{3}

Error de Truncamiento

error_{truncamiento} = -\frac{3}{80} h^5 f^{(4)} (z)

donde z se encuentra entre[a,b]

Tarea: Desarrollar el algoritmo en python

7.2 Regla de Simpson 1/3

Referencia: Chapra 21.2.1 p631 pdf655, Burden 4.2 p192 pdf202, Rodriguez 7.1.4 p281

Es el resultado cuando se realiza una interpolación con polinomio de segundo grado.

I = \int_a^b f(x) \delta x \cong \int_a^b f_2 (x) \delta x

Usando un polinomio de Lagrange de segundo grado:

I = \int_{x_0}^{x_2} \Big[ \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} f(x_0) + + \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} f(x_1) + + \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} f(x_2) \Big] \delta x

que tiene como resultado para un solo tramo:

I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]

siendo

h=\frac{b-a}{2}

Error de truncamiento

la cota del error de truncamiento se estima como O(h5)

error_{trunca} = -\frac{h^5}{90} f^{(4)}(z)

para un valor de z entre [a,b]

para cuantificar el valor, se puede usar la diferencia finita Δ4f, pues con la derivada sería muy laborioso.

Algoritmo en Python

Del ejercicio con trapecios, se repite el ejercicio con n tramos; usando dos tramos o tres puntos por cada iteración. Cada iteración se procesa avanzando dos puntos xi. Ejemplo:

tramos:  8
Integral:  2.05363501328
>>> 

A continuación se presentan dos formas de algoritmos que entregan los mismos resultados. Como tarea puede revisar la diferencia de tiempos de ejecución de cada uno.

Algoritmo conceptual

Se realiza mediante la aplicación directa de la forma conceptual de la fórmula para cada segmento conformado de dos tramos.

# Integración: Regla Simpson 1/3
import numpy as np
import matplotlib.pyplot as plt

#Función a evaluar
def funcionx(x):
    f = np.sqrt(x)*np.sin(x)
    return(f)

# PROGRAMA
a = 1 # float(input('a: '))
b = 3 # float(input('b: '))
tramos = 8 # int(input('tramos: '))
# validar tramos par

# PROCEDIMIENTO
h = (b-a)/tramos
x = a
area = 0
for i in range(0,tramos,2):
    deltaA = (h/3)*(funcionx(x)+4*funcionx(x+h)+funcionx(x+2*h))
    area = area + deltaA
    x = x + 2*h

# SALIDA
print('Integral: ', area)

Algoritmo mejorado

Presenta una forma meojrada del algoritmo para minimizar las operaciones entre cada tramo. El algoritmo toma en cuenta que en el primer segmento, el extremo izquierdo se suma una vez, pero el extremo derecho se vuelve a repetir la suma con el siguiente segmento, es decir se suman dos veces los segmentos pares en el interior del intervalo [a,b].

Se puede observar mejor lo expresado usando la gráfica, en la que cada segmento se destaca con diferente color.

# Integración: Regla Simpson 1/3
# Usar cantidad de tramos pares
import numpy as np
import matplotlib.pyplot as plt

#Función a evaluar
def funcionx(x):
    f = np.sqrt(x)*np.sin(x)
    return(f)

# PROGRAMA
a = float(input('a: '))
b = float(input('b: '))
tramos = int(input('tramos: '))
esimpar = tramos%2
while (esimpar == 1):
    tramos = int(input('tramos es par: '))

# PROCEDIMIENTO
h = (b-a)/tramos
x = a
# segmento: por cada dos tramos
suma = funcionx(x)
for i in range(0,tramos-2,2):
    x = x + h
    suma = suma + 4*funcionx(x)
    x = x + h
    suma = suma + 2*funcionx(x)
# último segmento 
x= x + h
suma = suma + 4*funcionx(x)
suma = suma + funcionx(b)
area = h*(suma/3)

# SALIDA
print('Integral: ', area)

Tarea: realizar la gráfica de la función y las áreas bajo las curvas

7.1 Regla del trapecio

Referencia: Chapra 21.1 p621 pdf645, Burden 4.2 p192 pdf202, Rodriguez 7.1.1 p273

Es la primera de las formulas cerradas de Newton-Cotes, usando un polinomio de primer grado.

I = \int_a^b f(x) dx \cong \int_a^b f_1 (x) dx

usando el rango entre [a,b] el polinomio se aproxima con una línea recta:

f_1 (x) = f(a) + \frac{f(b)-f(a)}{b-a} (x-a)

el área bajo esta línea recta es una aproximación de la integral de f(x) entre los límites a y b.

El resultado del integral es la regla del trapecio:

I = (b-a) \frac{f(a)+f(b)}{2}

que se interpreta como la multiplicación entre la base y altura promedio de un trapecio.

Algoritmo en Python

Para integrar la función \sqrt {(x)} \sin(x) en el intervalo [1,3] con 4, 16, 32 ,64 y 128 tramos, los resultados serán:

tramos: 4
Integral: 1.99841708623
tramos:  16
Integral:  2.05019783717
tramos:  32
Integral:  2.05277225085
tramos:  64
Integral:  2.05341563776
tramos: 128
Integral:  2.05357647096

instrucciones en python:

# Integración: Regla de los trapecios
import numpy as np
import matplotlib.pyplot as plt

#Función a evaluar
def funcionx(x):
    fxi = np.sqrt(x)*np.sin(x)
    return(fxi)

# PROGRAMA
a = 1 # float(input('a: '))
b = 3 # float(input('b: '))
tramos = 4 # int(input('tramos: '))

# PROCEDIMIENTO
h = (b-a)/tramos
x = a
suma = funcionx(x)
for i in range(0,tramos-1,1):
    x = x+h
    suma = suma + 2*funcionx(x)
suma = suma + funcionx(b)
area = h*(suma/2)

# SALIDA
print('tramos: ', tramos)
print('Integral: ', area)

La gráfica como ayuda didáctica, si la cantidad de tramos sea «poca», muestra los trapecios, si aumenta la cantidad de tramos, no es necesario poner líneas verticales en blanco para separar los trapecios:

# Puntos para la gráfica
muestras = tramos+1
xi = np.linspace(a,b,muestras)
fi = funcionx(xi)

# Gráfica
# Referencia función contínua
xig = np.linspace(a,b,muestras*10)
fig = funcionx(xig)
plt.plot(xig,fig)
# Trapecios
plt.fill_between(xi,0,fi, color='g')
plt.title('Integral: Regla de Trapecios')
for i in range(0,muestras,1):
    plt.axvline(xi[i], color='w')
plt.plot(xi,fi, 'o') # puntos muestra
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()

Para rellenar el área bajo la curva se usa la instrucción plt.fill_between(xi,0,fi, color=’g’). La división de los trapecios con una línea blanca (white) se realiza con la instrucción plt.axvline(xi[i], color=’w’).