[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ gráfica ] [ Algoritmo 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) ] [ gráfica ] [ Algoritmo 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) ] [ gráfica ] [ Algoritmo 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 = 6 # 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 xi = np.linspace(a,b,muestras) fi = fx(xi) # Regla de Simpson 3/8 h = (b-a)/tramos suma = 0 # integral numérico for i in range(0,tramos-2,3): #muestras-3 S38 = (3/8)*h*(fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]) suma = suma + S38 # SALIDA print('tramos:', tramos) print('Integral fx con Simpson 3/8: ', suma)
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ gráfica ] [ Algoritmo xi,fi ]
..
4. Gráfica para integral f(x) con Simpson 3/8
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 titulo = 'Regla de Simpson 3/8' titulo = titulo + ', tramos:'+str(tramos) titulo = titulo + ', Area:'+str(suma) fx_existe = True try: # fx suave aumentando muestras muestrasfxSuave = tramos*10 + 1 xk = np.linspace(a,b,muestrasfxSuave) fk = fx(xk) except NameError: fx_existe = False try: # existen mensajes de error msj_existe = len(msj) except NameError: msj = [] # 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' if len(msj)==0: # sin errores plt.vlines(xi[i],0,fi[i],linestyle=tipolinea, color='orange') # Graficar f(x), puntos if fx_existe==True: plt.plot(xk,fk,label='f(x)') plt.plot(xi,fi,'o',color='orange',label ='muestras') plt.xlabel('x') plt.ylabel('f(x)') plt.title(titulo) plt.legend() plt.tight_layout() plt.show()
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ gráfica ] [ Algoritmo xi,fi ]
..
5. 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]
Resultados con el algoritmo
tramos: 6 Integral fx con Simpson3/8: 2.0542043057079566
La gráfica se puede obtener usando el bloque correspondiente en la sección anterior.
# 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) muestras = len(xi) tramos = muestras-1 msj = [] # mensajes de error if tramos<3: # puntos insuficientes msj.append(['tramos insuficientes:',tramos]) suma = 0 # integral numerico i = 0 while i<(muestras-3) and len(msj)==0: # i<tramos, al menos dos tramos h = xi[i+1]-xi[i] # tamaño de paso, supone constante if (i+2)<muestras: # dos tramos d2x = np.diff(xi[i:i+4],2) # diferencias entre tramos errado = np.max(np.abs(d2x)) if errado<casicero: # 3 tramos equidistantes S38 = (3/8)*h*(fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]) suma = suma + S38 i = i+3 else: donde = np.argmax(np.abs(d2x))+i+1 msj.append(['no equidistantes en i',donde]) # SALIDA print('tramos: ',len(xi)-1) if len(msj)>0: # mensajes de error print('Revisar errores en xi:',xi) print('tramos:',np.diff(xi,1)) for unmsj in msj: print('',unmsj[0],':',unmsj[1]) else: print('Integral fx con Simpson3/8:',suma)
[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo f(x) ] [ gráfica ] [ Algoritmo xi,fi ]
..