5.2 Regla de Simpson 1/3

Referencia: Chapra 21.2.1 p631 pdf655, Burden 4.2 p192 pdf202, Rodriguez 7.1.4 p281

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

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 tiene como resultado para un solo tramo:

I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]

siendo

h=\frac{b-a}{2}

Error de truncamiento

la cota del error de truncamiento se estima como O(h5)

error_{trunca} = -\frac{h^5}{90} f^{(4)}(z)

para un valor de z entre [a,b]

para cuantificar el valor, se puede usar la diferencia finita Δ4f, pues con la derivada sería muy laborioso.


Ejercicio

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

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

tramos = 4

h = \frac{3-1}{4} = 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)] 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.

tramos: 4
Integral:  2.0549261957703937
>>>

Algoritmo en Python

Del ejercicio con trapecios, se repite el ejercicio con n tramos; usando dos tramos o tres puntos por cada iteración. Cada iteración se procesa avanzando dos puntos xi. Ejemplo:

tramos:  8
Integral:  2.05363501328
>>> 

A continuación se presentan dos formas de algoritmos que entregan los mismos resultados. Como tarea puede revisar la diferencia de tiempos de ejecución de cada uno.

Algoritmo conceptual

Se realiza mediante la aplicación directa de la forma conceptual de la fórmula para cada segmento conformado de dos tramos.

# 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

# PROCEDIMIENTO
# Tarea: validar tramos par

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

Algoritmo 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)]
# Integración: Regla Simpson 1/3
# Validar cantidad de tramos pares
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):
    tramos = int(input('tramos es par: '))
    esimpar = tramos%2

# PROCEDIMIENTO
# Regla de Simpson 1/3, varios tramos
h = (b-a)/tramos
xi = a
# segmento por cada dos tramos
suma = fx(xi)
for i in range(0,tramos-2,2):
    xi = xi + h
    suma = suma + 4*fx(xi)
    xi = xi + h
    suma = suma + 2*fx(xi)
# último segmento
xi = xi + h
suma = suma + 4*fx(xi)
suma = suma + fx(b)
area = (h/3)*suma

# SALIDA
print('tramos: ', tramos)
print('Integral: ', area)

Gráfica

Se puede observar mejor lo expresado usando la gráfica, en la que cada segmento se destaca con diferente color.

Tarea: realizar la gráfica de la función y las áreas bajo las curvas