5.2 Regla de Simpson 1/3 con Python

[ 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 x$

Se 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 x$

que 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(h⁵)

$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 Δ⁴f, 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 3$

Para 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.054$

Note 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 x_i, x_i+h, x_i+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]