[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]
..
1. Regla de Simpson 3/8
Referencia: Chapra 21.2.1 p.631, Burden 4.2 p147, Rodríguez 7.1.8 p288
I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)]
Es el resultado cuando para el integral se utiliza el resultado de una interpolación con polinomio de tercer grado.
I = \int_a^b f(x) \delta x I \cong \int_a^b f_3 (x) \delta xAl desarrollar la fórmula de la Regla de Simpson de 3/8 para un segmento con tres tramos con distancia h:
I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)]siendo el tramo de un segmento [a,b]
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]
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]
2. Ejercicio
Para integrar la función en el intervalo [1,3] con 6, 18, 24 y 72 tramos,
f(x)= \sqrt {(x)} \sin(x) 1 \leq x \leq 3Para el ejercicio planteado en la regla de trapecio, usando seis tramos, se aplica el método cada tres tramos.
tramos = 6
h = \frac{3-1}{6} = \frac{1}{3} I\cong \frac{3}{8}\Big(\frac{1}{3}\Big)[f(1)+3f\Big(1+\frac{1}{3}\Big) +3f\Big(1+\frac{2}{3}\Big) + f(2)] + + \frac{3}{8}\Big(\frac{1}{3}\Big)[f(2)+3f\Big(2+\frac{1}{3}\Big) +3f\Big(2+\frac{2}{3}\Big)+ f(3)] f(1)= \sqrt {(1)} \sin(1) = 0.8414 f(4/3)= \sqrt {(4/3)} \sin(4/3) = 1.1222 f(5/3)= \sqrt {(5/3)} \sin(5/3) = 1.2850 f(2)= \sqrt {(2)} \sin(2) = 1.2859 f(7/3)= \sqrt {(7/3)} \sin(7/3) = 1.1045 f(8/3)= \sqrt {(8/3)} \sin(8/3) = 0.7467 f(3)= \sqrt {(3)} \sin(3) = 0.2444 I\cong \frac{3}{24}[0.8414+3(1.1222)+3(1.2850)+1.2859] + \frac{3}{24}[1.2859+3(1.1045)+3(0.7467)+0.2444] I\cong 2.0542las muestras para el ejercicio son:
>>> import numpy as np
>>> fx = lambda x: np.sqrt(x)*np.sin(x)
>>> xi = [1, 1+1/3, 1+2/3, 2, 1+4/3, 1+5/3, 3]
>>> fx(xi)
array([0.84147098, 1.12229722, 1.28506615, 1.28594075,
1.10453193, 0.74672307, 0.24442702]
>>> I1=(3/8)*(1/3)*(fx(1)+3*fx(1+1/3)+3*fx(1+2/3)+fx(2))
>>> I2=(3/8)*(1/3)*(fx(2)+3*fx(1+4/3)+3*fx(1+5/3)+fx(3))
>>> I1+I2
2.0542043270535757
>>>
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]
3. Algoritmo en Python para f(x)
A partir del algoritmo para Simpson de 1/3, se realizan modificaciones para obtener el algoritmo de Simpson de 3/8:
Resultados de algoritmo
tramos: 3 Integral fx con Simpson 3/8: 2.0636730597016992
tramos: 6 Integral con Simpson 3/8: 2.0542043270535757
tramos: 9 Integral fx con Simpson 3/8: 2.053741914538051
tramos: 12 Integral fx con Simpson 3/8: 2.0536653010019474
Caso: f(x) es una expresión matemática
# Regla Simpson 3/8 para f(x) entre [a,b],tramos import numpy as np # INGRESO fx = lambda x: np.sqrt(x)*np.sin(x) a = 1 # intervalo de integración b = 3 tramos = 12 # multiplo de 3 # validar: tramos debe ser múltiplo de 3 while tramos%3 >0: print('tramos: ',tramos) txt = 'tramos debe ser múltiplo de 3:' tramos = int(input(txt)) # PROCEDIMIENTO muestras = tramos + 1 # Regla de Simpson 3/8 h = (b-a)/tramos xi = a suma = 0 # integral numérico for i in range(0,tramos-2,3): #muestras-3 S38 = (3/8)*h*(fx(xi)+3*fx(xi+h)+3*fx(xi+2*h)+fx(xi+3*h)) suma = suma + S38 xi = xi + 3*h # avanza 3 tramos # SALIDA print('tramos:', tramos) print('Integral fx con Simpson 3/8: ', suma)
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]
..
3.1 Gráfica
Se puede observar mejor lo expresado usando la gráfica para observar los tramos, muestras, por segmento.
Se puede cambiar el número de tramos en el algoritmo anterior, considerando que deben ser múltiplo de 3.
Instrucciones en Python adicionales al algoritmo anterior
# GRAFICA --------------------- import matplotlib.pyplot as plt # puntos de muestras muestras = tramos + 1 xi = np.linspace(a,b,muestras) fi = fx(xi) # fx suave aumentando muestras muestrasfxSuave = tramos*20 + 1 xk = np.linspace(a,b,muestrasfxSuave) fk = fx(xk) # Simpson 3/8 relleno y bordes, cada 3 tramos for i in range(0,muestras-2,3): x_tramo = xi[i:(i+3)+1] f_tramo = fi[i:(i+3)+1] # interpolación polinomica a*(x**3)+b*(x**2)+c*x+d coef = np.polyfit(x_tramo, f_tramo, 3) # [a,b,c,d] px = lambda x: coef[0]*(x**3)+coef[1]*(x**2)+coef[2]*x+coef[3] xp = np.linspace(x_tramo[0],x_tramo[-1],21) fp = px(xp) plt.plot(xp,fp,linestyle='dashed',color='orange') relleno = 'lightgreen' if (i/3)%2==0: # bloque 3 tramos, es par relleno ='lightblue' plt.fill_between(xp,fp,fp*0,color=relleno) # Divisiones entre Simpson 3/8 for i in range(0,muestras,1): tipolinea = 'dotted' if i%3==0: # i es multiplo de 3 tipolinea = 'dashed' plt.vlines(xi[i],0,fi[i], linestyle=tipolinea, color='orange') # Graficar f(x), puntos plt.plot(xk,fk,label='f(x)') plt.plot(xi,fi,'o',color='orange',label ='muestras') plt.xlabel('x') plt.ylabel('f(x)') titulo = 'Regla de Simpson 3/8' titulo = titulo + ', tramos:'+str(tramos) titulo = titulo + ', Area:'+str(suma) plt.title(titulo) plt.legend() plt.tight_layout() plt.show()
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]
..
4. Algoritmo para x[i], f[i] como muestras en Python
Cuando se integra sobre muestras dadas por los vectores xi, fi. Se revisa que los tramos entre muestras sean iguales para un Simpson 3/8 y que existan suficientes puntos para completar el integral.
xi = [1. , 1.33333333, 1.66666667, 2. ,
2.33333333, 2.66666667, 3. ]
fi = [0.84147098, 1.12229722, 1.28506615, 1.28594075,
1.10453193, 0.74672307, 0.24442702]
tolera = 1e-7
Resultados con el algoritmo
tramos: 6 Integral con Simpson 3/8: 2.054204305707957
Instrucciones en Python
# Integración Simpson 3/8 para muestras xi,fi import numpy as np # INGRESO xi = [1. , 1.33333333, 1.66666667, 2. , 2.33333333, 2.66666667, 3. ] fi = [0.84147098, 1.12229722, 1.28506615, 1.28594075, 1.10453193, 0.74672307, 0.24442702] # PROCEDIMIENTO casicero=1e-7 # decimales para truncar # vectores como arreglo, numeros reales xi = np.array(xi,dtype=float) fi = np.array(fi,dtype=float) n = len(xi) # revisa tramos iguales o equidistantes errado = -1 if (n-1)>=3: # len(tramos)>=3 d2 = np.diff(xi,2) # diferencias entre tramos errado = np.max(np.abs(d2)) if errado<casicero: # tramos equidistantes suma = 0 # integral numerico i = 0 h = xi[i+1]-xi[i] while (i+3)<n: S38 = (3/8)*h*(fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]) suma = suma + S38 i = i+3 msj = [] # mensajes de error if errado>=casicero: donde = np.argmax(d2)+2 msj.append(['tramos no equidistantes desde i:',donde]) sobra = (n-1)%3 # tramos multiplo 3 if sobra>0: # tramos no multiplo 3 msj.append(['tramos sobran:',sobra]) # SALIDA print('tramos: ',len(xi)-1) if len(msj)==0: print('Integral con Simpson 3/8: ',suma) else: print(' Revisar errores...') for unmsj in msj: print(unmsj[0],unmsj[1])
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]
..
5. Algoritmo: Integración Formulas Compuestas con Simpson 3/8, Simpson 1/3 y Trapecio. Tramos no iguales
Cuando los tramos no son múltiplos de 3, faltan tramos para completar la regla de Simpson 3/8. Se integra el último tramo con Simpson 1/3 o trapecios.
Los datos de ejemplo son:
xi = [1. , 1.33333333, 1.66666667, 2. ,
2.33333333, 2.66666667, 3.,
3.2, 3.4, 3.5 ]
fi = [0.84147098, 1.12229722, 1.28506615, 1.28594075,
1.10453193, 0.74672307, 0.24442702,
-0.10442284,-0.47119451,-0.65625532]
resultados:
Fórmulas compuestas, tramos: 9 métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado tramos iguales en xi: [3 0 0 3 0 0 2 0 1 0] ['S38', 1.168687718313123] ['S38', 0.8855165873948337] ['S13', -0.04296392333333337] ['trp', -0.05637249150000005] Integral: 1.9548678908746233
instrucciones en Python
# Integración Simpson 3/8, Simpson 1/3 y trapecio # tramos no iguales, formulas compuestas import numpy as np # INGRESO xi = [1. , 1.33333333, 1.66666667, 2. , 2.33333333, 2.66666667, 3., 3.2, 3.4, 3.5 ] fi = [0.84147098, 1.12229722, 1.28506615, 1.28594075, 1.10453193, 0.74672307, 0.24442702, -0.10442284,-0.47119451,-0.65625532] # PROCEDIMIENTO casicero = 1e-7 # digitos en muestras # 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 tramo_igual = -1*np.ones(n,dtype=int) integra_parcial =[] suma = 0 # integral total i = 0 while i<(n-1): #tramos tramo_integrado = False h = xi[i+1]-xi[i] # tres tramos iguales if (tramo_integrado==False) and (i+3)<n: d2 = np.diff(xi[i:i+4],2) # diferencias entre tramos errado = np.max(np.abs(d2)) 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 tramo_igual[i:i+3] = [3,0,0] # Simpson 3/8 integra_parcial.append(['S38',S38]) i = i+3 tramo_integrado = True # dos tramos iguales if (tramo_integrado==False) and (i+2)<n: d2 = np.diff(xi[i:i+3],2) # diferencias entre tramos errado = np.max(np.abs(d2)) if errado<casicero: # Simpson 1/3 S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2]) suma = suma + S13 tramo_igual[i:i+2] = [2,0] integra_parcial.append(['S13',S13]) i = i+2 tramo_integrado = True # un tramo igual if (tramo_integrado == False) and (i+1)<n: trapecio =(h/2)*(fi[i]+fi[i+1]) suma = suma + trapecio tramo_igual[i:i+1] = [1] # usar trapecio integra_parcial.append(['trp',trapecio]) i = i+1 tramo_integrado = True tramo_igual[n-1]=0 # ultima casilla # SALIDA 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:',tramo_igual) for unparcial in integra_parcial: print(unparcial) print('Integral:',suma)
5.1 Gráfica – Integración Compuesta, Simpson 1/3 y Trapecio
Se pueden añadir instrucciones al algoritmo de integral con tramos no equidistantes, para observar la secuencia de fórmulas apliacadas en la gráfica.
# GRAFICA --------------------- import matplotlib.pyplot as plt metodotipo = [['na','grey'], ['Trapecio','lightblue'], ['Simpson1/3','lightgreen'], ['Simpson3/8','cyan']] tramos = n-1 # etiquetas linea plt.plot(xi[0],fi[0],'o', label=metodotipo[1][0], color=metodotipo[1][1]) plt.plot(xi[0],fi[0],'o', label=metodotipo[2][0], color=metodotipo[2][1]) plt.plot(xi[0],fi[0],'o', label=metodotipo[3][0], color=metodotipo[3][1]) # relleno y bordes tipo = 0 for i in range(0,tramos,1): if tramo_igual[i]>0: tipo = tramo_igual[i] plt.fill_between([xi[i],xi[i+1]], [fi[i],fi[i+1]],[0,0], color=metodotipo[tipo][1]) # Divisiones entre tramos for i in range(0,n,1): tipolinea = 'dashed' # inicia un método if tramo_igual[i]==0: # dentro de un método tipolinea = 'dotted' plt.vlines(xi[i],0,fi[i],linestyle=tipolinea, color='orange') # puntos de xi,fi plt.plot(xi,fi,marker='o',linestyle='dashed', color='orange',label='f(xi)') plt.axhline(0,color='black') # eje x plt.xlabel('x') plt.ylabel('f(x)') plt.title('Integral numérico xi,fi con Fórmulas Compuestas') plt.legend() plt.tight_layout() plt.show()
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ Algoritmo xi,fi ]
[ Integral Compuesto xi,fi ]