Referencia: Chapra 25.3.3 p746 pdf 770, Rodriguez 9.1.8 p358
Para una ecuación diferencial de primer orden con una condición de inicio, la fórmula de Runge-Kutta de 4to orden se obtiene de la expresión con cinco términos:
y_{i+1} = y_i + aK_1 + bK_2 + cK_3 + dK_4
siendo:
y'(x) = f(x_i,y_i) y(x_0) = y_0
debe ser equivalente a la serie de Taylor de 5 términos:
que usando aproximaciones de derivadas, se obtienen:
# Runge Kutta de 4do ordendefrungekutta4(d1y,x0,y0,h,muestras):
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,2),dtype=float)
# incluye el punto [x0,y0]
estimado[0] = [x0,y0]
xi = x0
yi = y0
for i inrange(1,tamano,1):
K1 = h * d1y(xi,yi)
K2 = h * d1y(xi+h/2, yi + K1/2)
K3 = h * d1y(xi+h/2, yi + K2/2)
K4 = h * d1y(xi+h, yi + K3)
yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
xi = xi + h
estimado[i] = [xi,yi]
return(estimado)
Note que el método de Runge-Kutta de 4to orden es similar a la regla de Simpson 1/3. La ecuación representa un promedio ponderado para establecer la mejor pendiente.
La segunda parte corresponde a Runge-Kutta de 4to Orden
Como tema de introducción observar dos minutos del video sugerido a partir de donde se encuentra marcado el enlace. En este caso, en combate aereo, las armas se encuentran fijas en las alas.
Video Revisar:
Luego de observar el video de introducción conteste las siguientes preguntas:
¿ Que trayectoria siguen los proyectiles al salir del cañon?
¿ Que trayectoria siguen los aviones, el perseguido y el que caza?
¿ Cuál es la relación entre las trayectorias de los dos aviones?
Los métodos de Runge-Kutta mejoran la aproximación a la respuesta sin requerir determinar las expresiones de las derivadas de orden superior. Los métodos usan una corrección a la derivada tomando valores de puntos alrededor referenciado al tamaño de paso h.
Por ejemplo, Runge-Kutta de 2do Orden usa el promedio entre los incrementos xi y xi+h, calculados como términos K1 y K2.
# EDO. Método de RungeKutta 2do Orden # estima la solucion para muestras espaciadas h en eje x# valores iniciales x0,y0# entrega arreglo [[x,y]]import numpy as np
defrungekutta2(d1y,x0,y0,h,muestras):
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,2),dtype=float)
# incluye el punto [x0,y0]
estimado[0] = [x0,y0]
xi = x0
yi = y0
for i inrange(1,tamano,1):
K1 = h * d1y(xi,yi)
K2 = h * d1y(xi+h, yi + K1)
yi = yi + (K1+K2)/2
xi = xi + h
estimado[i] = [xi,yi]
return(estimado)
Ejercicio
Para probar el algoritmo se usa la ecuación del problema presentado en ‘EDO con Taylor‘ :
En los métodos con Taylor para Ecuaciones Diferenciales Ordinarias (EDO), se aproxima el resultado a n términos de la serie, para lo cual se ajusta la expresión del problema a la derivada correspondiente.
Ejemplo:
Se pide encontrar puntos de la solución en la ecuación diferencial usando los tres primeros términos de la serie de Taylor con h=0.1 y punto inicial x0=0, y0=1
y'-y -x +x^2 -1 = 0
La solución empieza usando la Serie de Taylor para tres términos:
y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h E = \frac{h^3}{3!} y'''(z) = O(h^3)
Luego, al despejar de la fórmula planteada el valor de y’, se puede obtener y" a combinar las ecuaciones.
ecuaciones que permiten estimar nuevos valores yi+1 para nuevos puntos muestra distanciados en i*h desde el punto inicial.
Se empieza evaluando el nuevo punto a una distancia x1= x0+h del punto de origen con lo que se obtiene y1 , repitiendo el proceso para el siguiente punto en forma sucesiva
Algoritmo con Python
Para simplificar los calculos se crea una función edo_taylor3t() para encontrar los valores para una cantidad de muestras distanciadas entre si h veces del punto inicial [x0,y0]
# EDO. Método de Taylor 3 términos # estima solucion para muestras separadas h en eje x# valores iniciales x0,y0# entrega arreglo [[x,y]]import numpy as np
defedo_taylor3t(d1y,d2y,x0,y0,h,muestras):
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,4),dtype=float)
# incluye el punto [x0,y0]
estimado[0] = [x0,y0,0,0]
x = x0
y = y0
for i inrange(1,tamano,1):
estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
x = x+h
estimado[i,0:2] = [x,y]
return(estimado)
La función se puede usar en un programa que plantee resolver el problema propuesto:
# PROGRAMA PRUEBA# Ref Rodriguez 9.1.1 p335 ejemplo.# prueba y'-y-x+(x**2)-1 =0, y(0)=1# INGRESO.# d1y = y', d2y = y''
d1y = lambda x,y: y -x**2 + x + 1
d2y = lambda x,y: y -x**2 - x + 2
x0 = 0
y0 = 1
h = 0.1
muestras = 5
# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]
# SALIDAprint('estimado[xi, yi, d1yi, d2yi]')
print(puntos)
La ecuación diferencial ordinaria del ejercicio tiene una solución conocida, lo que permite encontrar el error real en cada punto respecto a la aproximación estimada.
y = e^x + x + x^2
Note que el error crece al distanciarse del punto inicial
# ERROR vs solución conocida
y_sol = lambda x: ((np.e)**x) + x + x**2
yi_psol = y_sol(xi)
errores = yi_psol - yi
errormax = np.max(np.abs(errores))
# SALIDAprint('Error máximo estimado: ',errormax)
print('entre puntos: ')
print(errores)
Para visualizar los resultados se usan los puntos estimados y muchos más puntos para la solución conocida en el intervalo de estudio. Con la gráfica se podrán observar las diferencias entre las soluciones encontradas.
# GRAFICA [a,b+2*h]
a = x0
b = h*muestras+2*h
muestreo = 10*muestras+2
xis = np.linspace(a,b,muestreo)
yis = y_sol(xis)
# Gráficaimport matplotlib.pyplot as plt
plt.plot(xis,yis, label='y conocida')
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()
Tarea: Realizar el ejercicio con más puntos muestra, donde se visualice que el error aumenta al aumentar la distancia del punto inicial [x0,y0]
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.
Tarea
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.
Algoritmo con Python
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)
# SALIDAprint('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 inrange(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()
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 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:
con 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.
Ejercicio
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) dx
que 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
>>>
Algoritmo con Python
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 puntosdefintegraCuadGauss2p(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 inrange(0,muestras-1,1):
deltaA = integraCuadGauss2p(fx,xi[i],xi[i+1])
area = area + deltaA
# SALIDAprint('Integral: ', area)
Gráfica por tramos
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 inrange(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 tramosfor j inrange(0,len(subtramo),1):
plt.axvline(subtramo[j], color='tab:gray')
# Marcadores de puntos xa y xb por tramosfor j inrange(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 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.
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/3import 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 inrange(0,tramos,2):
deltaA = (h/3)*(fx(xi)+4*fx(xi+h)+fx(xi+2*h))
area = area + deltaA
xi = xi + 2*h
# SALIDAprint('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.
# Integración: Regla Simpson 1/3# Validar cantidad de tramos paresimport 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 inrange(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
# SALIDAprint('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
Algoritmo con muestras
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,fiimport numpy as np
import matplotlib.pyplot as plt
defintegrasimpson13_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
whilenot(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'
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)
# SALIDAprint('tramos: ',len(xi)-1)
print('Integral con Simpson 1/3: ',Area)
iftype(Area)==str:
print(' Revisar errores...')
resultados
tramos: 4
Integral con Simpson 1/3: 2.0549261966666665
>>>
La integración numérica con la regla del trapecio usa un polinomio de primer grado como aproximación de f(x) entre los extremos a y b del intervalo.
Es la primera de las formulas cerradas de Newton-Cotes.
I = \int_a^b f(x) dx \cong \int_a^b f_1 (x) dx
usando 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.
El 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 tiene un orden de h3. (Rodríguez p275)
error_{truncar} = -\frac{h^3}{12}f''(z)
a < z < b
Algoritmo en Python
Ejemplo de Integración con el método del trapecio
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 = h
para 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:
y el integral de todo el intervalo es la suma de todos los trapecios.
Integral = A_1+ A_2+ A_3+A_4
si 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:
una vez la función evaluada en cada extremo,
mas dos veces cada valor de la función evaluada en los puntos intermedios.
El resultado de la suma se multiplica por h/2.
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 en Python – Regla o Método del trapecio
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 inrange(0,tramos-1,1):
xi = xi + h
suma = suma + 2*fx(xi)
suma = suma + fx(b)
area = h*(suma/2)
# SALIDAprint('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 inrange(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 con muestras
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 una muestras xi,fiimport numpy as np
import matplotlib.pyplot as plt
defintegratrapecio_fi(xi,fi):
''' sobre muestras de fi para cada xi
integral con método de trapecio
'''
n = len(xi)
suma = 0
for i inrange(0,n-1,1):
dx = xi[i+1]-xi[i]
untrapecio = dx*(fi[i+1]+fi[i])/2
suma = suma + untrapecio
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 = integratrapecio_fi(xi,fi)
# SALIDAprint('tramos: ',len(xi)-1)
print('Integral con trapecio: ',Area)
En algunos casos, los datos (x,y) no tienen una relación de tipo funcional y(x), entonces no se pueden aplicar directamente los métodos de interpolación revisados.
Por ejemplo, en la trayectoria del balón en el «gol imposible», la gráfica de la trayectoria en el espacio o sus proyecciones en los planos dependen del parámetro tiempo «t»en lugar de una relación de x,y,z
Sin embargo si las coordenadas (x,y) se expresan como funciones de otra variable t denomiada parámetro, entonces los puntos x(t), y(t) tienen relación funcional, y se pueden construir polinomios de interpolación.
Ejemplo
Las coordinadas x(t) y y(t) del recorrido de un cohete registradas en los instantes t fueron:
ti = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]
Usaremos un algoritmo en Python para mostrar la trayectoria x,y para el problema planteado.
Al realizar la interpolación de los puntos para obtener polinomios que dependen de «t» se obtiene:
polinomios con los que se puede realizar la gráfica px(t), py(t) en forma separada. Pero para comprender mejor la trayectoria del cohete, se utiliza la gráfica px(t) vs py(t) en el intervalo t entre[0,3]
Tiene como objetivo incorporar condiciones adicionales para la función en los extremos del intervalo donde se encuentran los puntos o «nodos».
Por ejemplo, si los datos son los de una trayectoria en un experimento de física, se podría disponer de la aceleración el punto inicial de las muestras (izquierda) y de salida (derecha) del intervalo de las muestras.
El método busca obtener un polinómio de tercer grado para cada sub-intervalo o tramo entre «nodos» consecutivos (j, j+1), de la forma:
Para n datos existen n-1 tramos, cuatro incógnitas (coeficientes) que evaluar por cada tramo, como resultado 4(n-1) incógnitas para todo el intervalo .
Para los términos (xj+1– xj) de los tramos que se usan varias veces en el desarrollo, se simplifican como hj:
Para generar el sistema de ecuaciones, se siguen los siguientes criteros:
1. Los valores de la función deben ser iguales en los nodos interiores
S_j(x_{j+1}) = f(x_{j+1}) = S_{j+1}(x_{j+1})
2. Las primeras derivadas en los nodos interiores deben ser iguales
S'_j(x_{j+1}) = S'_{j+1}(x_{j+1})
3. Las segundas derivadas en los nodos interiores deben ser iguales
S''_j(x_{j+1}) = S''_{j+1}(x_{j+1})
4. El primer y último polinomio deben pasar a través de los puntos extremos
f(x_0) = S_0(x_0) = a_0 f(x_n) = S_n(x_n) = a_n
5. Una de las condiciones de frontera se satisface:
5a. frontera libre o natural: Las segundas derivadas en los nodos extremos son cero
S''(x_0) = S''(x_n) = 0
5b. frontera sujeta: las primeras derivadas en los nodos extremos son conocidas
S'(x_0) = f'(x_0)
S'(x_n) = f'(x_n)
Tarea: Revisar en los textos el planteamiento de las ecuaciones para resolver el sistema que se genera al plantear el polinomio.
Ubicar en el texto las ecuaciones resultantes, que son las que se aplicarán en el algoritmo.
Algoritmo en Python
El algoritmo parte de lo desarrollado para «trazadores lineales», donde se presentaron los bloques para:
construir el trazador en una tabla de polinomios por tramos,
graficar los trazadores en cada tramo, requientdo,
evaluar cada uno de ellos en cada tramo con muestras suficientes para presentar una buena precisión en la gráfica.
Del algoritmo básico se modifica entonces el bloque del cálculo de los polinomios de acuerdo al planteamiento de formulas y procedimientos para trazadores cúbicos naturales (enviado a revisar como tarea).
# Trazador cúbico natural# Condición: S''(x_0) = S''(x_n) = 0import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
deftraza3natural(xi,yi):
n = len(xi)
# Valores h
h = np.zeros(n-1, dtype = float)
for j inrange(0,n-1,1):
h[j] = xi[j+1] - xi[j]
# Sistema de ecuaciones
A = np.zeros(shape=(n-2,n-2), dtype = float)
B = np.zeros(n-2, dtype = float)
S = np.zeros(n, dtype = float)
A[0,0] = 2*(h[0]+h[1])
A[0,1] = h[1]
B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])
for i inrange(1,n-3,1):
A[i,i-1] = h[i]
A[i,i] = 2*(h[i]+h[i+1])
A[i,i+1] = h[i+1]
factor21 = (yi[i+2]-yi[i+1])/h[i+1]
factor10 = (yi[i+1]-yi[i])/h[i]
B[i] = 6*(factor21 - factor10)
A[n-3,n-4] = h[n-3]
A[n-3,n-3] = 2*(h[n-3]+h[n-2])
factor12 = (yi[n-1]-yi[n-2])/h[n-2]
factor23 = (yi[n-2]-yi[n-3])/h[n-3]
B[n-3] = 6*(factor12 - factor23)
# Resolver sistema de ecuaciones S
r = np.linalg.solve(A,B)
for j inrange(1,n-1,1):
S[j] = r[j-1]
S[0] = 0
S[n-1] = 0
# Coeficientes
a = np.zeros(n-1, dtype = float)
b = np.zeros(n-1, dtype = float)
c = np.zeros(n-1, dtype = float)
d = np.zeros(n-1, dtype = float)
for j inrange(0,n-1,1):
a[j] = (S[j+1]-S[j])/(6*h[j])
b[j] = S[j]/2
factor10 = (yi[j+1]-yi[j])/h[j]
c[j] = factor10 - (2*h[j]*S[j]+h[j]*S[j+1])/6
d[j] = yi[j]
# Polinomio trazador
x = sym.Symbol('x')
px_tabla = []
for j inrange(0,n-1,1):
pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2
pxtramo = pxtramo + c[j]*(x-xi[j])+ d[j]
pxtramo = pxtramo.expand()
px_tabla.append(pxtramo)
return(px_tabla)
# PROGRAMA -----------------------# INGRESO , Datos de prueba
xi = np.array([0.1 , 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.8, 1.7, 2.0])
muestras = 10 # entre cada par de puntos# PROCEDIMIENTO# Tabla de polinomios por tramos
n = len(xi)
px_tabla = traza3natural(xi,fi)
# SALIDAprint('Polinomios por tramos: ')
for tramo inrange(1,n,1):
print(' x = ['+str(xi[tramo-1])
+','+str(xi[tramo])+']')
print(str(px_tabla[tramo-1]))
con lo que se obtiene el resultado por cada tramo
Polinomios por tramos:
x = [0.1,0.2]
-146.666666666667*x**3 + 44.0*x**2 + 0.566666666666666*x + 1.1
x = [0.2,0.3]
283.333333333333*x**3 - 214.0*x**2 + 52.1666666666667*x - 2.34
x = [0.3,0.4]
-136.666666666667*x**3 + 164.0*x**2 - 61.2333333333333*x + 9.0
>>>
Para el trazado de la gráfica se añade al final del algoritmo:
# GRAFICA# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
whilenot(tramo>=n):
a = xi[tramo-1]
b = xi[tramo]
xtramo = np.linspace(a,b,muestras)
# evalua polinomio del tramo
pxtramo = px_tabla[tramo-1]
pxt = sym.lambdify('x',pxtramo)
ytramo = pxt(xtramo)
# vectores de trazador en x,y
xtraza = np.concatenate((xtraza,xtramo))
ytraza = np.concatenate((ytraza,ytramo))
tramo = tramo + 1
# Gráfica
plt.plot(xi,fi,'ro', label='puntos')
plt.plot(xtraza,ytraza, label='trazador'
, color='blue')
plt.title('Trazadores Cúbicos Naturales')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()
Si los polinomios no se igualan entre los nodos, tampoco sus velocidades y aceleraciones; podría considerar la experiencia semejante a variaciones de velocidad y aceleración en una trayectoria:
Car sales woman scares customers. maxman.tv. 4 enero 2016
Pilot shows off the impact of flying at 9.5G.@loteofficial. 2 Marzo2023
El concepto de trazador se originó en la técnica de dibujo que usa una cinta delgada y flexible (spline) para dibujar curvas suaves a través de un conjunto de puntos.
La unión más simple entre dos puntos es una línea recta. El método crea un polinomio para cada par de puntos consecutivos en el intervalo, por lo que el resultado será una tabla de polinomios.
Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto de funciones lineales.
Observe que la expresión de f(x) para un tramo entre dos puntos es el polinomio de grado 1 realizado con diferencia finita avanzadas o las diferencias divididas.
Las ecuaciones se pueden usar para evaluar la función en cualquier punto entre x0 y xn. Al localizar primero el intervalo dentro del cual está el punto, puede seleccionar el polinomio que corresponde a ese tramo.
Polinomios por tramos:
x = [0.1,0.2] p(x) = 3.5*x + 1.1
x = [0.2,0.3] p(x) = -1.0*x + 2.0
x = [0.3,0.4] p(x) = 3.0*x + 0.8
Algoritmo en Python
Datos de los puntos como ejemplo para el algoritmo
xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]
El método con trazadores lineales, permite plantear los bloques necesarios para manejar varios polinomios, uno por cada tramo entre dos puntos dentro del intervalo del problema
El método permitirá disponer de un punto de partida para trazadores de mayor grado, por ejemplo los cúbicos.
Los polinomios de cada tramo se almacenan en una tabla, cada uno puede ser utilizado individualmente en su respectivo tramo, por ejemplo para realizar la gráfica de línea entre tramos.
# Trazador (spline) lineal, grado 1import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
deftrazalineal(xi,fi):
n = len(xi)
x = sym.Symbol('x')
px_tabla = []
tramo = 1
whilenot(tramo>=n):
# con 1ra diferencia finita avanzada
numerador = fi[tramo]-fi[tramo-1]
denominador = xi[tramo]-xi[tramo-1]
m = numerador/denominador
pxtramo = fi[tramo-1] + m*(x-xi[tramo-1])
px_tabla.append(pxtramo)
tramo = tramo + 1
return(px_tabla)
# PROGRAMA# INGRESO , Datos de prueba
xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]
muestras = 10 # entre cada par de puntos# PROCEDIMIENTO# Tabla de polinomios por tramos
n = len(xi)
px_tabla = trazalineal(xi,fi)
# SALIDAprint('Polinomios por tramos: ')
for tramo inrange(1,n,1):
print(' x = ['+str(xi[tramo-1])
+','+str(xi[tramo])+']')
print(str(px_tabla[tramo-1]))
Se obtiene como resultado:
Polinomios por tramos:
x = [0.1,0.2]
3.5*x + 1.1
x = [0.2,0.3]
-1.0*x + 2.0
x = [0.3,0.4]
3.0*x + 0.8
>>>
Para añadir la gráfica se añaden las instrucciones para:
evaluar el polinomio en cada tramo
concatenar los resultados de todos los tramos en los vectores xtraza, ytraza.
poner en la gráfica los puntos del problema y las líneas que genera cada polinomio
# GRAFICA# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
whilenot(tramo>=n):
a = xi[tramo-1]
b = xi[tramo]
xtramo = np.linspace(a,b,muestras)
# evalua polinomio del tramo
pxtramo = px_tabla[tramo-1]
pxt = sym.lambdify('x',pxtramo)
ytramo = pxt(xtramo)
# vectores de trazador en x,y
xtraza = np.concatenate((xtraza,xtramo))
ytraza = np.concatenate((ytraza,ytramo))
tramo = tramo + 1
# Gráfica
plt.plot(xi,fi,'o', label='puntos')
plt.plot(xtraza,ytraza, label='trazador')
plt.title('Trazadores lineales (splines)')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()