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]
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
>>>
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.0636730597016992tramos: 6
Integral con Simpson 3/8: 2.0542043270535757tramos: 9
Integral fx con Simpson 3/8: 2.053741914538051tramos: 12
Integral fx con Simpson 3/8: 2.0536653010019474Caso: 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)
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()
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.0542043057079566La 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)