5.3 Regla de Simpson 3/8 con Python

[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo /función Python]

..


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)]

regla de Simpson 3/8

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 x

Al 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 /función Python]

..


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 3

Para el ejercicio planteado en la regla de trapecio, usando seis tramos, se aplica el método cada tres tramos.

f(x)= \sqrt {(x)} \sin(x) 1 \leq x \leq 3

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.0542

las 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 /función Python]

..


Algoritmo en Python

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:  6
Integral con Simpson 3/8:  2.0542043270535757

Caso: f(x) es una expresión matemática

# Integración Simpson 3/8
import numpy as np
import matplotlib.pyplot as plt

def integrasimpson38(fx,a,b,tramos):
    ''' integral con método de Simpson 3/8
        tramos debe ser múltiplo de 3
    '''
    # Validar cantidad de tramos múltiplo de 3
    esmult3 = tramos%3
    if esmult3 != 0:
        suma = 'tramos debe ser múltiplo de 3'
    # Regla de Simpson 3/8
    if esmult3 == 0:
        h = (b-a)/tramos
        xi = a
        suma = 0
        for i in range(0,tramos,3):
            unS38 = (3/8)*h*(fx(xi)+3*fx(xi+h)+3*fx(xi+2*h)+fx(xi+3*h))
            suma = suma + unS38
            xi = xi + 3*h
    return(suma)

# PROGRAMA -----------------
# INGRESO
fx = lambda x: np.sqrt(x)*np.sin(x)

# intervalo de integración
a = 1
b = 3
tramos = 6

# PROCEDIMIENTO
Area = integrasimpson38(fx,a,b,tramos)

# SALIDA
print('tramos: ',tramos)
print('Integral con Simpson 3/8: ',Area)
if type(Area)==str:
    print('  Revisar errores...')

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.

reglas Simpson 3/8Instrucciones 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-2,3):
    relleno = 'lightgreen'
    if (i/3)%2==0: # i/3 multiplo 2
        relleno ='lightblue'
    xktramo = xk[i*4:(i+3)*4+1]
    fktramo = fk[i*4:(i+3)*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 3/8
for i in range(0,muestras,1):
    tipolinea = 'dotted'
    if i%3==0: # multiplo de 3
        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 3/8')
plt.legend()
plt.show()

Caso: x[i], f[i] son muestras

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

Instrucciones en Python

# Integración Simpson 1/3
# Usando una muestras xi,fi
import numpy as np
import matplotlib.pyplot as plt

def integrasimpson38_fi(xi,fi,tolera = 1e-10):
    ''' sobre muestras de fi para cada xi
        integral con método de Simpson 3/8
        respuesta es np.nan para tramos desiguales,
        no hay suficientes puntos.
    '''
    n = len(xi)
    i = 0
    suma = 0
    while not(i>=(n-3)):
        h  = xi[i+1]-xi[i]
        h1 = (xi[i+2]-xi[i+1])
        h2 = (xi[i+3]-xi[i+2])
        dh = abs(h-h1)+abs(h-h2)
        if dh<tolera:# tramos iguales
            unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
            unS38 = (3/8)*h*unS38
            suma = suma + unS38
        else:  # tramos desiguales
            suma = 'tramos desiguales'
        i = i + 3
    if (i+1)<n: # incompleto, tramos por calcular
        suma = 'tramos incompletos, faltan '
        suma = suma +str(n-(i+1))+' tramos'
    return(suma)

# PROGRAMA -----------------
# 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]
tolera = 1e-7

# PROCEDIMIENTO
n = len(xi)-1
Area = integrasimpson38_fi(xi,fi,tolera)

# SALIDA
print('tramos: ',n)
print('Integral con Simpson 3/8: ',Area)
if type(Area)==str:
    print('  Revisar errores')

[ Simpson 3/8 ] [ Ejercicio ] [Algoritmo /función Python]