Métodos de Integración Numérica
Métodos para Diferenciación Numérica
Tablas de Diferencias divididas centradas, hacia adelante y hacia atras
Curso con Python – MATG1052/MATG1013-FCNM-ESPOL
Métodos de Integración Numérica
Métodos para Diferenciación Numérica
Tablas de Diferencias divididas centradas, hacia adelante y hacia atras
Referencia: Chapra Fig.23.1 pag.669 pdf.693, Burden 4.1 p167, Rodriguez 8.2,3,4,6 p324
Primera derivada
f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h) f'(x_i) = \frac{-f(x_{i+2})+4f(x_{i+1})-3f(x_i)}{2h} + O(h^2)Segunda derivada
f''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+f(x_i)}{h^2} + O(h) f''(x_i) = \frac{-f(x_{i+3})+4f(x_{i+2})-5f(x_{i+1})+2f(x_i)}{h^2} + O(h^2)Tercera derivada
f'''(x_i) = \frac{f(x_{i+3})-3f(x_{i+2})+3f(x_{i+1})-f(x_i)}{h^3} + O(h) f'''(x_i) = \frac{-3f(x_{i+4})+14f(x_{i+3})-24f(x_{i+2})+18f(x_{i+1})-5f(x_i)}{2h^3} + O(h^2)Cuarta derivada
f''''(x_i) = \frac{f(x_{i+4})-4f(x_{i+3})+6f(x_{i+2})-4f(x_{i+1})+f(x_i)}{h^3} + O(h)Primera derivada
f'(x_i) = \frac{f(x_{i+1})-f(x_{i-1})}{2h} + O(h^2) f'(x_i) = \frac{-f(x_{i+2})+8f(x_{i+1})-8f(x_{i-1}) +f(x_{i-2})}{12h} + O(h^4)Segunda derivada
f''(x_i) = \frac{f(x_{i+1})-2f(x_{i})+f(x_{i-1})}{h^2} + O(h^2) f''(x_i) = \frac{-f(x_{i+2})+16f(x_{i+1})-30f(x_{i})+16f(x_{i-1})-f(x_{i-2})}{12h^2} + O(h^4)Tercera derivada
f'''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+2f(x_{i-1})-f(x_{i-2})}{2h^3} + O(h^2) f'''(x_i) = \frac{-f(x_{i+3})+8f(x_{i+2})-13f(x_{i+1})+13f(x_{i-1})-8f(x_{i-2})+f(x_{i-3})}{8h^3} + O(h^4)Primera derivada
f'(x_i) = \frac{f(x_{i})-f(x_{i-1})}{h} + O(h) f'(x_i) = \frac{3f(x_{i})-4f(x_{i-1})+f(x_{i-2})}{2h} + O(h^2)Segunda derivada
f''(x_i) = \frac{f(x_{i})-2f(x_{i-1})+f(x_{i-2})}{h^2} + O(h) f''(x_i) = \frac{2f(x_{i})-5f(x_{i-1})+4f(x_{i-2})-f(x_{i-3})}{h^2} + O(h^2)Tercera derivada
f'''(x_i) = \frac{f(x_{i})-3f(x_{i-1})+3f(x_{i-2})-f(x_{i-3})}{h^3} + O(h) f'''(x_i) = \frac{5f(x_{i})-18f(x_{i-1})+24f(x_{i-2})-14f(x_{i-3})+3f(x_{i-4})}{2h^3} + O(h^2)Referencia: Chapra 23.1 p668 pdf692, Rodriguez 8.2 p324
Como referencia, el polinomio de Taylor muestra una aproximación de una función f(x):
P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 + ...Una aproximación a primera derivada, usa los primeros dos términos del polinomio de Taylor alrededor de xi en para un punto a la derecha xi+1 a una distancia h = xi+1–xi
f(x_{i+1}) = f(x_i)+\frac{f'(x_i)}{1!} (h) + \frac{f''(x_i)}{2!}(h)^2 + ...se puede simplificar en un polinomio de grado uno y un término de error:
f_{i+1} = f_i + f'_i (h) + O(h^2) ...Despejando la expresión para f’i
La expresión también es la primera diferencia finita dividida con un error del orden O(h). (tema usado en interpolación).
Revise que el término de la expresión queda O(h2)/h con lo que se disminuye el exponente en uno.
Se realiza el mismo procedimiento que el anterior, usando un punto xi+1 y xi-1 alrededor de xi. En el término xi-1 el valor de h es negativo al invertir el orden de la resta.
f_{i+1} = f_i+\frac{f'_i}{1!}(h) + \frac{f''_i}{2!}(h)^2 f_{i-1} = f_i-\frac{f'_i}{1!}(h) + \frac{f''_i}{2!}(h)^2restando la ecuaciones se tiene que
f_{i+1} - f_{i-1} = f'_i (h)+f'_i (h) f_{i+1} - f_{i-1} = 2h f'_iLa expresión de primera derivada usando un punto antes y otro después del punto central queda como:
con un error del orden O(h2)
Al continuar con el procedimiento mostrado se pueden obtener las fórmulas para segundas derivadas, las que se resumen en las tablas de Fórmulas de diferenciación por diferencias divididas.
Para visualizar el concepto de Cuadratura de Gauss de dos puntos considere lo siguiente:
Se tiene un corte transversal del un recipiente rectangular lleno de líquido limitado en x entre [-1,1], al que se le ha depositado encima otro recipiente con perfil f(x) = x2 hasta que reposan sus extremos en x[-1,1].
La altura de ambos recipientes es la misma.
La superficie entre f(x) y el eje x es el integral de f(x) en el intervalo.
Si suponemos que la figura es el corte transversal de una vasija y la parte en amarillo es líquida, la vasija ha desplazado el líquido que ocupa ahora el «área» mostrada en la gráfica que corresponde al integral de f(x)=x2. entre [-1,1].
Ahora, suponga que se perfora el perfil de f(x) en dos puntos equidistantes cercanos x=0. Los orificios permitirían el desplazamiento del liquido al interior de f(x) que dejando pasar suficiente tiempo, perimitiría tener todo el líquido en el recipiente rectangular entre [-1,1] como una línea horizontal.
Podría medir la altura que tiene el líquido y que tiene un equivalente en un punto f(x1). Debería encontrar el valor de x1 que permite disponer del mismo valor entre el área debajo de f(x) y el rectángulo del corte transversal amarillo ahora formado.
Se usa el resultado analítico del integral restando el área del rectángulo obtenido al evaluar la funcion f(x) entre [0,1], teniendo un problema de busqueda de raíces. Obtenemos el valor de x1.
Se muestra que el área bajo f(x) es equivalente al área del rectángulo conformado.
Si utilizamos el desplazamiento horizontal desde el centro para un punto encontrado como un «factor», tendremos que el área del rectángulo se mantendría equivalente, y el desplazamiento proporcional a la mitad del intervalo si se varía el intervalo de observacion. Este factor coincide con el factor de Cuadratura de Gauss de dos puntos.
funcion fx: x**2 Integral Fx: x**3/3 I analitico: 0.6666666666666666 I aproximado: 0.6666666666654123 desplaza centro: 0.5773502691890826 factor desplaza: 0.5773502691890826 Factor CuadGauss: 0.5773502691896258 erradoFactor: 1.2545520178264269e-12 error integral: 1.2543299732215019e-12
El error del integral es del orden de 10-12
Cambiamos la figura geométrica a un trapecio generado por la recta que pasa por los puntos xi desplazados desde el centro.
Usamos la función f(x) = x2 + x + 1, observaremos si los resultado son equivalentes.
La figura al inicio del experimento será:
Luego de realizar realizar el mismo cálculo anterior usando un equivalente a trapecio se tiene:
con valores numéricos:
funcion fx: x**2 + x + 1 Integral Fx: x**3/3 + x**2/2 + x I analitico: 2.6666666666666665 I aproximado: 2.6666666666654124 desplaza centro: 0.5773502691890826 factor desplaza: 0.5773502691890826 Factor CuadGauss: 0.5773502691896258 erradoFactor: 1.2545520178264269e-12 error integral: 1.2541079286165768e-12
El error del integral es también del orden de 10-12, además observe que el factor de cuadratura de Gauss se mantiene.
Realice el experimento usando un polinomio de grado superior y observe los errores para el integral y las diferencia con el coeficiente de Cuatratura de 2 puntos.
Para resumir la cantidad de instrucciones, se usa el método de la bisección desde la librería scipy y el subgrupo de funciones de optimización.
Los cálculos para realizar las gráficas se tratan en un bloque luego de mostrar los resultado principales.
# Integración: Cuadratura de Gauss de dos puntos # modelo con varios tramos entre [a,b] # para un solo segmento. import numpy as np import matplotlib.pyplot as plt import sympy as sym import scipy.optimize as op # INGRESO x = sym.Symbol('x') # fx = (x)**2 fx = x**2 + x + 1 # fx = 0.2 + 25.0*x-200*(x**2)+675.0*(x**3)-900.0*(x**4)+400.0*(x**5) a = -1 b = 1 muestras = 51 tolera = 1e-12 iteramax = 100 # PROCEDIMIENTO # Desarrollo analítico con Sympy Fx = sym.integrate(fx,x) Fxn = sym.lambdify('x',Fx,'numpy') Fiab = Fxn(b)-Fxn(a) # Busca igualar trapecio con Integral analitico fxn = sym.lambdify('x',fx,'numpy') base = b-a mitad = base/2 xc = (a+b)/2 # centro diferencia = lambda x: Fiab-base*(fxn(xc-x)+fxn(xc+x))/2 desplazado = op.bisect(diferencia,0,mitad, xtol=tolera,maxiter=iteramax) factor = desplazado/mitad # Integral aproximando con trapecio x0 = xc - factor*mitad x1 = xc + factor*mitad Faprox = base*(fxn(x0)+fxn(x1))/2 # Integral cuadratura Gauss xa = xc + mitad/np.sqrt(3) xb = xc - mitad/np.sqrt(3) FcuadG = base*(fxn(xa)+fxn(xb))/2 erradofactor = np.abs(FcuadG - Faprox) erradoIntegral = np.abs(Fiab-Faprox) # SALIDA print('funcion fx: ', fx) print('Integral Fx: ', Fx) print('I analitico: ', Fiab) print('I aproximado: ', Faprox) print('desplaza centro: ', desplazado) print('factor desplaza: ', factor) print('Factor CuadGauss: ', 1/np.sqrt(3)) print('erradoFactor: ', erradofactor) print('error integral: ', erradoIntegral) # Grafica # Para GRAFICAR # Para gráfica f(x) xi = np.linspace(a,b,muestras) fi = fxn(xi) # Para gráfica Trapecio m = (fxn(x1)-fxn(x0))/(x1-x0) trapeciof = lambda x: fxn(x0)+m*(x-x0) trapecioi = trapeciof(xi) # Areas Trapecio para cada punto que busca k = int(muestras/2) xicg = xi[k:muestras-1] Fcg = [base*(fxn(xi[k+0])+fxn(xi[k-0]))/2] for i in range(1,k,1): untrapecio = base*(fxn(xi[k+i])+fxn(xi[k-i]))/2 Fcg.append(untrapecio) # Punto buscado Fiaprox = base*(fxn(x1)+fxn(x0))/2 Fi = Fxn(xi)-Fxn(a) # Areas de curvas y trapecio plt.subplot(211) # Grafica superior plt.xlim(a,b) plt.plot(xi, fi, label='f(x)') # Solo fi # plt.fill_between(xi,0, fi, # label='integral fi', # color='yellow') # usando cuadratura plt.fill_between(xi,0, trapecioi, label='Cuadratura 2 puntos', color='yellow') plt.axvline(x0,color='white') plt.axvline(x1,color='white') plt.plot([x0,x1],[fxn(x0),fxn(x1)], 'ro', label='x0,x1') plt.axvline(0,color='black') plt.xlabel('x') plt.ylabel('f(x) y Cuadratura de 2 puntos') plt.legend() # Valores de integrales plt.subplot(212) # Grafica inferior plt.xlim(a,b) plt.axhline(Fiab, label='F[a,b]') # plt.plot(xi,Fi,label='F(x)') plt.plot(xicg,Fcg,color='orange',label='Aprox Trapecio') plt.axvline(x1,color='yellow') plt.axvline((1/np.sqrt(3))*(b-a)/2 + xc ,color='magenta') plt.plot(x1,Fiaprox,'ro', label='x0,x1') plt.axvline(0,color='black') plt.xlabel('x') plt.legend() plt.ylabel('Integrando') plt.show()
Referencia: Chapra 22.3 p655 pdf679, Rodriguez 7.3 p294, Burden 4.7 p220 pdf230
La cuadratura de Gauss aproxima el integral de una función en un intervalo [a,b] centrado en cero mediante un cálculo numérico con menos operaciones y evaluaciones de la función. Se representa como una suma ponderada:
para la fórmula de dos puntos se tiene obtiene:
c_0 = c_1 = 1 x_0 = -\frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}para un intervalo de evaluación desplazado en el eje x se requiere convertir los puntos al nuevo rango. Se desplaza el punto cero al centro del intervalo [a,b] y se obtiene:
x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1con lo que el resultado aproximado del integral se convierte en:
I \cong \frac{b-a}{2}(f(x_a) + f(x_b))cuya fórmula es semejante a una mejor aproximación de un trapecio, cuyos promedios de alturas son puntos internos de [a,b], concepto mostrado en la gráfica.
El algoritmo se desarrolla en un tramo en el intervalo [a,b] junto a la gráfica para mostrar el concepto. Para el ejemplo el integral buscado es:
\int_1^3 \sqrt{x} \sin(x) dxque usando cuadratura de Gauss con un solo intervalo[a,b] tiene como resultado:
Integral: 2.03821975687 >>>
el resultado se puede mejorar aumentando el número de tramos en el intervalo [a,b]. Por ejemplo, el resultado usando 4 tramos el resultado es semejante al usar el método del trapecio con 128 tramos, lo que muestra el ahorro en calculos entre los métodos
Integral: 2.05357719003 >>>
Instrucciones usando la cuadratura de Gauss como una función
# Integración: Cuadratura de Gauss de dos puntos # modelo con varios tramos entre [a,b] import numpy as np import matplotlib.pyplot as plt # cuadratura de Gauss de dos puntos def integraCuadGauss2p(funcionx,a,b): x0 = -1/np.sqrt(3) x1 = -x0 xa = (b+a)/2 + (b-a)/2*(x0) xb = (b+a)/2 + (b-a)/2*(x1) area = ((b-a)/2)*(funcionx(xa) + funcionx(xb)) return(area) # INGRESO fx = lambda x: np.sqrt(x)*np.sin(x) # intervalo de integración a = 1 b = 3 tramos = 4 # PROCEDIMIENTO muestras = tramos+1 xi = np.linspace(a,b,muestras) area = 0 for i in range(0,muestras-1,1): deltaA = integraCuadGauss2p(fx,xi[i],xi[i+1]) area = area + deltaA # SALIDA print('Integral: ', area)
La gráfica que complementa el resultado anterior, se realiza añadiendo las intrucciones presentadas a continuación.
Considere que la gráfica es útil con pocos tramos en el intervalo[a,b]
# GRAFICAR por cada Segmento # para concepto con pocos segmentos x0 = -1/np.sqrt(3) x1 = 1/np.sqrt(3) # arregos para gráficas xi = np.array([]) fi = np.array([]) xat = np.array([]) xbt = np.array([]) recta = np.array([]) muestrastramo = 10 subtramo = np.linspace(a,b,muestras) for i in range(0,tramos,1): at = subtramo[i] bt = subtramo[i+1] xit = np.linspace(at,bt,muestrastramo) fit = fx(xit) xi = np.concatenate((xi,xit)) fi = np.concatenate((fi,fit)) # puntos xa y xb por tramo xa = (bt+at)/2 + (bt-at)/2*(x0) xb = (bt+at)/2 + (bt-at)/2*(x1) xat = np.concatenate((xat,[xa])) xbt = np.concatenate((xbt,[xb])) # Recta entre puntos x0 y x1 por tramo m = (fx(xb)-fx(xa))/(xb-xa) b0 = fx(xa)- m*xa linea = b0 + m*xit recta = np.concatenate((recta,linea)) # Marcadores 'o' de xa y xb por tramos puntox = np.concatenate((xat,xbt)) puntoy = fx(puntox) # Graficando # Trazado de lineas plt.plot(xi,recta, label = 'grado 1', color = 'tab:orange') plt.fill_between(xi,0,recta, color='tab:olive') plt.plot(xi,fi, label='f(x)', color = 'blue') # Verticales para dividir los tramos for j in range(0,len(subtramo),1): plt.axvline(subtramo[j], color='tab:gray') # Marcadores de puntos xa y xb por tramos for j in range(0,len(xat),1): plt.axvline(xat[j], color='w') plt.axvline(xbt[j], color='w') plt.plot(puntox,puntoy, 'o', color='g') plt.title('Integral: Cuadratura Gauss') plt.xlabel('x') plt.ylabel('f(x)') plt.legend() plt.show()
Referencia: Chapra 21.2.1 p.631 pdf.655, Burden 4.2 p192 pdf202, Rodriguez 7.1.8 p288
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 \cong \int_a^b f_3 (x) \delta xAl desarrollar la fórmula se 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}donde z se encuentra entre[a,b]
Tarea: A partir del algoritmo para Simpson de 1/3, realizar las modificaciones para obtener el algoritmo para Simpson de 3/8
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 xUsando 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 xque 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}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.
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 3tramos = 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.054Note 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 >>>
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)
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
Referencia: Chapra 21.1 p621 pdf645, Burden 4.2 p192 pdf202, Rodriguez 7.1.1 p273
La integración numérica con la regla del trapecio usa un polinomio de primer grado.
Es la primera de las formulas cerradas de Newton-Cotes.
I = \int_a^b f(x) dx \cong \int_a^b f_1 (x) dxusando el rango entre [a,b] el polinomio se aproxima con una línea recta:
f_1 (x) = f(a) + \frac{f(b)-f(a)}{b-a} (x-a)el área bajo esta línea recta es una aproximación de la integral de f(x) entre los límites a y b.
El resultado del integral es la regla del trapecio:
I = (b-a) \frac{f(a)+f(b)}{2}que se interpreta como la multiplicación entre la base y altura promedio de un trapecio. También se llega al resultando sumando las áreas de sus componentes: rectángulo y triángulo.
Error de truncamiento se encuentra como el integral del término que le sigue al polinomio de Taylor en la aproximación, es decir el de grado 2, que al integrarlo fiene un orden de h3. (Rodriguez p275)
error_{truncar} = -\frac{h^3}{12}f''(z)a < z < b
Para integrar la función en el intervalo [1,3] con 4, 16, 32 ,64 y 128 tramos,
f(x)= \sqrt {(x)} \sin(x)Tramos = 4
La base de los trapecios (xi-xi+1), si esta base tiene valores equidistantes se sustituye por la variable h .
base =\frac{b-a}{tramos}=\frac{3-1}{4}= 0.5 = hpara cada tramo, el área de un trapecio es:
A_i = base\frac{f(x_i)+f(x_{i+1})}{2}por lo que cada trapecio en cada subintervalo tendrá un área de:
A_1 = (0.5)\frac{f(1)+f(1.5)}{2} A_2 = (0.5)\frac{f(1.5)+f(2)}{2} A_3 = (0.5)\frac{f(2)+f(2.5)}{2} A_4 = (0.5)\frac{f(2.5)+f(3)}{2}y el integral de todo el intervalo es la suma de todos los trapecios.
Integral = A_1+ A_2+ A_3+A_4si se quiere tomar el factor común entre las sumas de áreas de cada tramo, también se tiene la forma de la ecuación con un h equidistante entre cada tramo del intervalo.
Integral = \frac{0.5}{2}\Big(f(1)+2f(1.5) +2f(2)+2f(2.5)+f(3)\Big)Interpretando la fórmula, el integral es la suma de:
En caso que los puntos no sean equidistantes la simplificación NO procede, y se sumaran las áreas de los trapecios de cada tramo.
Siguiendo el procedimiento, para cada cantidad de tramos en el intervalo, los resultados serán:
tramos: 4 Integral: 1.99841708623
tramos: 16 Integral: 2.05019783717
tramos: 32 Integral: 2.05277225085
tramos: 64 Integral: 2.05341563776
tramos: 128 Integral: 2.05357647096
Algoritmo con fórmula simplificada para varios tramos, con puntos de muestras equidistantes h.
# Integración: Regla de los trapecios # Usando una función fx() 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 = 4 # PROCEDIMIENTO # Regla del Trapecio # Usando tramos equidistantes en intervalo h = (b-a)/tramos xi = a suma = fx(xi) for i in range(0,tramos-1,1): xi = xi + h suma = suma + 2*fx(xi) suma = suma + fx(b) area = h*(suma/2) # SALIDA print('tramos: ', tramos) print('Integral: ', area)
La gráfica como ayuda didáctica, si la cantidad de tramos sea «poca», muestra los trapecios, si aumenta la cantidad de tramos, no es necesario poner líneas verticales en blanco para separar los trapecios:
# GRAFICA # Puntos de muestra muestras = tramos + 1 xi = np.linspace(a,b,muestras) fi = fx(xi) # Linea suave muestraslinea = tramos*10 + 1 xk = np.linspace(a,b,muestraslinea) fk = fx(xk) # Graficando plt.plot(xk,fk, label ='f(x)') plt.plot(xi,fi, marker='o', color='orange', label ='muestras') plt.xlabel('x') plt.ylabel('f(x)') plt.title('Integral: Regla de Trapecios') plt.legend() # Trapecios plt.fill_between(xi,0,fi, color='g') for i in range(0,muestras,1): plt.axvline(xi[i], color='w') plt.show()
Para rellenar el área bajo la curva se usa la instrucción plt.fill_between(xi,0,fi, color=’g’). La división de los trapecios con una línea blanca (white) se realiza con la instrucción plt.axvline(xi[i], color=’w’).
Algoritmo usando las muestras y para tramos que pueden tener distancias arbitrarias. En este caso se usa la formula del área del trapecio para cada tramo.
Usado para ejercicios donde los datos proporcionados son muestras. Se adapta el bloque de ingreso y el procedimiento inicia desde «Regla de Trapecio». La cantidad de tramos será el numero de muestras menos uno tramos = len(xi) - 1
# Integración: Regla de los trapecios # Usando incluso muestras arbitrariamente espaciadas 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 = 4 # PROCEDIMIENTO # Puntos de muestra muestras = tramos + 1 xi = np.linspace(a,b,muestras) fi = fx(xi) # Regla del Trapecio # Usando puntos muestreados # incluso arbitrariamente espaciados suma = 0 for i in range(0,tramos,1): dx = xi[i+1]-xi[i] Atrapecio = dx*(fi[i]+fi[i+1])/2 suma = suma + Atrapecio integral = suma # SALIDA print('tramos: ', tramos) print('integral: ', integral)