[ Simpson 1/3 ] [ Ejercicio ] [Algoritmo Python] [ función Python]
..
La regla de Simpson 1/3
Referencia: Chapra 21.2.1 p631, Burden 4.3 p144, Rodríguez 7.1.4 p281
I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]
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
Se puede obtener usando un polinomio de Lagrange de segundo grado:
I = \int_{x_0}^{x_2} \bigg[ \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) \bigg] \delta x
que simplificando tiene como resultado para un solo tramo:
I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]
siendo h el tamaño de paso, donde para la expresión el divisor debe ser par, múltiplo de 2, al cubrir todo el intervalo [a,b]. En caso de usar valores de muestras xi, fi, el valor de h debe ser constante.
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 dentro del intervalo [a,b] de integración.
para cuantificar el valor, se puede usar la diferencia finita Δ4f, pues con la derivada sería muy laborioso.
[ Simpson 1/3 ] [ Ejercicio ] [Algoritmo Python] [ función Python]
..
Ejercicio
Para integrar la función en el intervalo [1,3] con 4, 16, 32 ,64 y 128 tramos,
f(x)= \sqrt {(x)} \sin(x)
1 \leq x \leq 3
Para el ejercicio planteado en la regla de trapecio, usando cuatro tramos, se aplica el método cada dos tramos.
tramos = 4
h = \frac{3-1}{4} = \frac{1}{2} = 0.5
I\cong \frac{0.5}{3}[f(1)+4f(1.5) + f(2)] +
+ \frac{0.5}{3}[f(2)+4f(2.5) + f(3)]
f(1)= \sqrt {(1)} \sin(1) = 0.8414
f(1.5)= \sqrt {(1.5)} \sin(1.5) = 1.2216
f(2)= \sqrt {(2)} \sin(2) = 1.2859
f(2.5)= \sqrt {(2.5)} \sin(2.5) = 0.9462
f(3)= \sqrt {(3)} \sin(3) = 0.2444
I\cong \frac{0.5}{3}[0.8414+4(1.2216) + 1.2859] +
+ \frac{0.5}{3}[1.2859+4(0.9462) + 0.2444]
I\cong 2.054
Note que al usar Simpson 1/3 con 4 tramos el resultado tiene los 2 primeros decimales iguales a usar Trapecio con 16 tramos.
>>> import numpy as np
>>> fx = lambda x: np.sqrt(x)*np.sin(x)
>>> xi = [1, 1+1/2, 1+2/2, 1+3/2,3]
>>> xi
[1, 1.5, 2.0, 2.5, 3]
>>> fx(xi)
array([0.84147098, 1.22167687, 1.28594075,
0.94626755, 0.24442702])
>>> (0.5/3)*(fx(1)+4*fx(1.5)+fx(2)) + (0.5/3)*(fx(2)+4*fx(2.5)+fx(3))
2.0549261957703937
>>>
Para más de 4 tramos es preferible realizar las operaciones con un algoritmo.
[ Simpson 1/3 ] [ Ejercicio ] [Algoritmo Python] [ función Python]
..
Algoritmo en Python
Del ejercicio con trapecios, se repite el ejercicio con n tramos; usando dos tramos (tres puntos) en cada iteración.
Cada iteración se procesa avanzando dos puntos xi, xi+h, xi+2h . Ejemplo:
tramos: 8
Integral: 2.053709383061734
>>>
Instrucciones en Python
Se realiza mediante la aplicación directa de la fórmula para cada segmento conformado de dos tramos. Se verifica que el valor de tramos sea par.
# Integración: Regla Simpson 1/3
import numpy as np
import matplotlib.pyplot as plt
# INGRESO:
fx = lambda x: np.sqrt(x)*np.sin(x)
# intervalo de integración
a = 1
b = 3
tramos = 8
# Validar cantidad de tramos pares
esimpar = tramos%2
while (esimpar == 1):
print('tramos: ',tramos)
tramos = int(input('tramos debe ser par: '))
esimpar = tramos%2
# PROCEDIMIENTO
# Regla de Simpson 1/3
h = (b-a)/tramos
xi = a
area = 0
for i in range(0,tramos,2):
deltaA = (h/3)*(fx(xi)+4*fx(xi+h)+fx(xi+2*h))
area = area + deltaA
xi = xi + 2*h
# SALIDA
print('tramos:', tramos)
print('Integral: ', area)
Gráfica
Se puede observar mejor lo expresado usando la gráfica para observar los tramos, muestras, por segmento. Para este caso se ejecuta el algoritmo anterior usando tramos =8
Instrucciones en Python adicionales al algoritmo anterior
# GRAFICA
# fx muestras por tramo
muestras = tramos + 1
xi = np.linspace(a,b,muestras)
fi = fx(xi)
fi0 = np.zeros(muestras) # linea base
# fx suave aumentando muestras
muestrasfxSuave = 4*tramos + 1
xk = np.linspace(a,b,muestrasfxSuave)
fk = fx(xk)
# Relleno
for i in range(0,muestras-1,2):
relleno = 'lightgreen'
if (i/2)%2==0: # i/2 par
relleno ='lightblue'
xktramo = xk[i*4:(i+2)*4+1]
fktramo = fk[i*4:(i+2)*4+1]
plt.fill_between(xktramo,fktramo,fktramo*0,color=relleno)
# Funcion f(x)
plt.plot(xk,fk, label='f(x)')
plt.plot(xi,fi,'o', label='f(xi)')
# Divisiones entre Simpson 1/3
for i in range(0,muestras,1):
tipolinea = 'dotted'
if i%2==0: # i par
tipolinea = 'dashed'
plt.vlines(xi[i],fi0[i],fi[i],
linestyle=tipolinea)
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f')
plt.title('Integral: Regla de Simpson 1/3')
plt.legend()
plt.show()
[ Simpson 1/3 ] [ Ejercicio ] [Algoritmo Python] [ función Python]
..
Algoritmo como función de Python
Caso: f(x) es una expresión matemática
A partir del ejercicio anterior, se resume en una función de Python si la entrada f(x) es una expresión matemática:
def integrasimpson13(fx,a,b,tramos):
''' integral con método de Simpson 1/3
tramos debe debe ser múltiplo de 2
'''
# Validar cantidad de tramos pares
esmult2 = tramos%2
if esmult2 != 0:
suma = 'tramos debe ser múltiplo de 2'
# Regla de Simpson 1/3
if esmult2 == 0:
h = (b-a)/tramos
xi = a
suma = 0
for i in range(0,tramos,2):
unS13 = (h/3)*(fx(xi)+4*fx(xi+h)+fx(xi+2*h))
suma = suma + unS13
xi = xi + 2*h
return(suma)
Caso: x[i], f[i] son muestras
xi = [1. , 1.5, 2. , 2.5, 3.]
fi = [0.84147098, 1.22167687, 1.28594075,
0.94626755, 0.24442702]
Cuando se integra sobre muestras dadas por los vectores xi, fi. Se revisa que los tramos entre muestras sean iguales para un Simpson 1/3 y que existan suficientes puntos para completar el integral.
# Integración Simpson 1/3
# Usando una muestras xi,fi
import numpy as np
import matplotlib.pyplot as plt
def integrasimpson13_fi(xi,fi,tolera = 1e-10):
''' sobre muestras de fi para cada xi
integral con método de Simpson 1/3
respuesta es np.nan para tramos desiguales,
no hay suficientes puntos.
'''
n = len(xi)
i = 0
suma = 0
while not(i>=(n-2)):
h = xi[i+1]-xi[i]
dh = abs(h - (xi[i+2]-xi[i+1]))
if dh<tolera:# tramos iguales
unS13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
suma = suma + unS13
else: # tramos desiguales
suma = 'tramos desiguales en i: '+str(i)
suma = suma ' , '+str(i+2)
i = i + 2
if i<(n-1): # incompleto, faltan tramos por calcular
suma = 'tramos incompletos, faltan '
suma = suma ++str((n-1)-i)+' tramos'
return(suma)
# PROGRAMA -----------------
# INGRESO
xi = [1. , 1.5, 2. , 2.5, 3.]
fi = [0.84147098, 1.22167687, 1.28594075,
0.94626755, 0.24442702]
# PROCEDIMIENTO
Area = integrasimpson13_fi(xi,fi)
# SALIDA
print('tramos: ',len(xi)-1)
print('Integral con Simpson 1/3: ',Area)
if type(Area)==str:
print(' Revisar errores...')
resultados
tramos: 4
Integral con Simpson 1/3: 2.0549261966666665
>>>
[ Simpson 1/3 ] [ Ejercicio ] [Algoritmo Python] [ función Python]
Fórmula con varios segmentos y h constante
Usado cuando el intervalo a integrar tiene varios segmentos, cada segmento tiene dos tramos. Ejemplo para dos segmentos, cuatro tramos, semejante al usado en la gráfica. La simplificación es válida si h es constante.
I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)] +
+ \frac{h}{3}[f(x_2)+4f(x_3) + f(x_4)]
tomando factor común h/3
I\cong \frac{h}{3}[f(x_0)+4f(x_1) + 2f(x_2) +
+4f(x_3) + f(x_4)]
[ Simpson 1/3 ] [ Ejercicio ] [Algoritmo Python] [ función Python]
