[ 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 xSe 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 xque 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 3Para 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.054Note 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]