Integración numérica¶
En análisis numérico, la integración numérica constituye una amplia gama de algoritmos para calcular el valor numérico de una integral definida y, por extensión, el término se usa a veces para describir algoritmos numéricos para resolver ecuaciones diferenciales. El término cuadratura numérica (a menudo abreviado a cuadratura) es más o menos sinónimo de integración numérica, especialmente si se aplica a integrales de una dimensión a pesar de que para el caso de dos o más dimensiones (integral múltiple) también se utiliza.
El problema básico considerado por la integración numérica es calcular una solución aproximada a la integral definida:
$$F(x)=\int_a^b f(x)\, dx$$Este problema también puede ser enunciado como un problema de valor inicial para una ecuación diferencial ordinaria, como sigue:
$$y'(x) = f(x), \quad y(a) = 0$$Encontrar $y(b)$ es equivalente a calcular la integral. Los métodos desarrollados para ecuaciones diferenciales ordinarias, como el método de Runge-Kutta, pueden ser aplicados al problema reformulado.
Razones para la integración numérica¶
Hay varias razones para llevar a cabo la integración numérica. La principal puede ser la imposibilidad de realizar la integración de forma analítica. Es decir, integrales que requerirían de un gran conocimiento y manejo de matemática avanzada pueden ser resueltas de una manera más sencilla mediante métodos numéricos. Incluso existen funciones integrables pero cuya primitiva no puede ser calculada, siendo la integración numérica de vital importancia. La solución analítica de una integral nos arrojaría una solución exacta, mientras que la solución numérica nos daría una solución aproximada. El error de la aproximación, depende del método que se utilice y de qué tan fino sea, puede llegar a ser tan pequeño que es posible obtener un resultado idéntico a la solución analítica en las primeras cifras decimales.
Métodos para integrales unidimensionales¶
Los métodos de integración numérica pueden ser descritos generalmente como combinación de evaluaciones del integrando para obtener una aproximación a la integral. Una parte importante del análisis de cualquier método de integración numérica es estudiar el comportamiento del error de aproximación como una función del número de evaluaciones del integrando. Un método que produce un pequeño error para un pequeño número de evaluaciones es normalmente considerado superior. Reduciendo el número de evaluaciones del integrando se reduce el número de operaciones aritméticas involucradas, y por tanto, se reduce el error de redondeo total. También, cada evaluación cuesta tiempo, y el integrando puede ser arbitrariamente complicado.
De todos modos, un modo de integración por fuerza bruta puede hacerse siempre, de un modo muy simplista, evaluando el integrando con incrementos muy pequeños.
Métodos basados en funciones de interpolación¶
Hay una extensa familia de métodos que se basan en aproximar la función a integrar $f(x)$ por otra función $g(x)$ de la cual se conoce la integral exacta. La función que sustituye la original se encuentra de forma que en un cierto número de puntos tenga el mismo valor que la original. Como los puntos extremos forman parte siempre de este conjunto de puntos, la nueva función se llama una interpolación de la función original. Cuando los puntos extremos no se utilizan para encontrar la función que sustituye a la original entonces se dice extrapolación. Típicamente estas funciones son polinomios.
Fórmulas de Newton-Cotes¶
La interpolación con polinomios evaluada en puntos igualmente separados en $[a, b]$ da las fórmulas de Newton-Cotes, de las que la regla del rectángulo, la del trapecio y la de Simpson son ejemplos. Si se escogen los nodos hasta $k=n+1$ será la fórmula de Newton-Cotes cerrada y si se escogen $k=n-1$ será la fórmula de Newton-Cotes abierta.
Regla del rectángulo¶
El método más simple de este tipo es hacer a la función interpoladora ser una función constante (un polinomio de orden cero) que pasa a través del punto $(a, f(a))$. Este método se llama la regla del rectángulo:
$$\int_a^b f(x) dx \sim (b-a) f(a)$$
Regla del trapecio¶
La regla se basa en aproximar el valor de la integral de $f(x)$ por el de la función lineal, que pasa a través de los puntos $(a, f(a))$ y $(b, f(b))$. La integral de ésta es igual al área del trapecio bajo la gráfica de la función lineal.
Regla del trapecio Simple¶
Para realizar la aproximación por esta regla es necesario usar un polinomio de primer orden, y esta es representada por: $$P_1(x) = f(a) + \frac{f(b) - f(a)}{b - a}(x - a)$$ Entonces al sustituir en la integral tenemos: $$\begin{align} I & = \int_a^bf(x)dx \approx \int_a^b P_1(x)dx \\ & \approx \int_a^b\left[f(a) + \frac{f(b) - f(a)}{b - a}(x - a)\right] \end{align}$$ Por último al resolver esa integral nos queda: $$\int_{a}^{b} f(x)dx \approx (b-a)\frac{f(a) + f(b)}{2}$$
Regla del trapecio compuesta¶
La regla del trapecio compuesta o regla de los trapecios es una forma de aproximar una integral definida utilizando $n$ trapecios. En la formulación de este método se supone que $f$ es continua y positiva en el intervalo $[a,b]$. De tal modo la integral definida $\int_a^b f(x)\, dx$ representa el área de la región delimitada por la gráfica de $f$ y el eje $x$, desde $x=a$ hasta $x=b$. Primero se divide el intervalo $[a,b]$ en $n$ subintervalos, cada uno de ancho $\Delta x=(b-a)/n$.
Después de realizar todo el proceso matemático se llega a la siguiente fórmula:
$$\int_a^b f(x)\, dx \sim \frac{h}{2} [f(a)+2f(a+h)+2f(a+2h)+...+f(b)]$$Donde $\textstyle h= \frac{b-a}{n}$ y $n$ es el número de divisiones.
La expresión anterior también se puede escribir como:
$$\int_a^b f(x) dx \sim \frac{b-a}{n} \left[ \frac{f(a) + f(b)}{2} + \sum_{k=1}^{n-1} f\left(a + k \frac{b-a}{n}\right) \right]$$
Regla de Simpson¶
La función interpoladora puede ser un polinomio de grado 2 que pasa a través de los puntos $(a, f(a))$, $(\frac{a+b}{2}, f(\frac{a+b}{2}))$ y $(b, f(b))$. Este método se llama regla de Simpson: $$\int_{a}^{b} f(x) \, dx \approx \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right)+f(b)\right]$$.
Regla de Simpson 1/3 Compuesta¶
En el caso de que el intervalo $[a, b]$ no sea lo suficientemente pequeño, el error al calcular la integral puede ser muy grande. Para ello, se recurre a la fórmula compuesta de Simpson o de segmentos múltiples, que consiste en dividir el intervalo $[a, b]$ en $n$ subintervalos iguales (con n par), de manera que $h = \frac{(b-a)}{n}$. Entonces la regla compuesta viene dada por:
$$I \approx \frac{h}{3} [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \dots + f(x_n)]$$O en su forma general:
$$I \approx \frac{h}{3}\left[f(x_0) + 4\sum_{i=1,3,5,\dots}^{n-1}f(x_i) + 2\sum_{j=2,4,6,\dots}^{n-2}f(x_j) + f(x_n)\right]$$from numpy import sqrt,exp,pi,linspace
def f(x):
#funcion=3*x
funcion=x**2
#funcion=(1/sqrt(2*pi))*exp(-(x**2)/2.0)
return funcion
#limites
A=1
B=2
%matplotlib inline
import matplotlib.pyplot as plt
x=linspace(A,B)
fig = plt.figure()
axes1 = fig.add_axes([0.1, 0.1, 0.8, 0.9])
axes1.fill_between(x, 0, f(x), color="green", alpha=0.5)
axes1.set_xlim([A-1,B+1])
plt.show()
def trapecio_simple(f,A,B):
simple=(B-A)*((f(A)+f(B))/2)
return simple
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
print 'El resultado de la regla del trapecio simple es', trapecio_simple(f,A,B)
def trapecio_compuesto(f,A,B,n):
k=1
valor=0
while k <= n-1:
c=A+k*(B-A)/n
valor+=f(c)
k += 1
compuesta=((B-A)/n)*((f(A) + f(B))/2.0 + valor)
return compuesta
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
N=raw_input("En cuanto se se divide el intervalo [a,b]: ")
n=int(N)
print 'El resultado de la regla del trapecio compuesta es', trapecio_compuesto(f,A,B,n)
def simpson(f,A,B):
h=(B-A)/6
sim_simp=h*(f(A)+4*f((A+B)/2)+f(B))
return sim_simp
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
print 'El resultado de la regla simple de Simpson es', simpson(f,A,B)
from math import *
def simpson13(n,A,B,f):
#calculamos h
h = (B-A) / n
#Inicializamos nuestra varible donde se almacenara las sumas
suma = 0.0
#hacemos un ciclo para ir sumando las areas
for i in range(1, n):
#calculamos la x
#x = a - h + (2 * h * i)
x = A + i * h
# si es par se multiplica por 4
if(i % 2 == 0):
suma = suma + 2 * f(x)
#en caso contrario se multiplica por 2
else:
suma = suma + 4 * f(x)
#sumamos los el primer elemento y el ultimo
suma = suma + f(A) + f(B)
#Multiplicamos por h/3
rest = suma * (h / 3)
#Retornamos el resultado
return (rest)
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
N=raw_input("En cuanto se se divide el intervalo [a,b]")
n=int(N)
print 'El resultado de la regla de de Simpson', simpson13(n,A,B,f)
Ejercicio¶
Un cuerpo se mueve a lo largo de una línea recta de acuerdo a la ley $v=t^3-4t^2+10$ m/s. Calcular el desplazamiento total del móvil cuatro segundos despues que inició el movimiento.
Un objeto se mueve con movimiento rectilíneo de modo tal que su velocidad en el instante t es $v(t) = t^2 - 2t$ metros por segundo. Halle:
a) El desplazamiento del objeto durante los tres primeros segundos.
b) La distancia recorrida durante ese tiempo.
#ley v=t3−4t2+10 m/s
from numpy import sqrt,exp,pi,linspace
def v(t):
#funcion=3*x
velocidad = (t**3)-(4*(t**2))+10
#funcion=(1/sqrt(2*pi))*exp(-(x**2)/2.0)
return velocidad
#limites
A=0
B=4
v(2)
%matplotlib inline
import matplotlib.pyplot as plt
x=linspace(A,B)
fig = plt.figure()
axes1 = fig.add_axes([0.1, 0.1, 0.8, 0.9])
axes1.fill_between(x, 0, f(x), color="green", alpha=0.5)
axes1.set_xlim([A-1,B+1])
plt.show()
#Con regla del trapecio compuesta
a = float(input("Ingresa el valor de a que vamos a evaluar: "))
b = float(input("Ingresa el valor de b que vamos a evaluar: "))
n = input("Ingresa el numero de divisiones que relizaremos: ")
f = ((b-a)/n)
f1 = (v(a)+v(b))/2
e = []
for k in range (1,n-1):
c = v((a+(k*((b-a)/n))))
e.append(c)
z = sum(e)
print "El desplazamiento total del móvil en 4 segundos es de: ", f*(f1+z), "metros"
#Con regla de Simpson
def sim(v,A,B):
h=(B-A)/6
sim_simp = h*(v(A)+4*v((A+B)/2)+v(B))
return sim_simp
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
print 'El resultado de la regla simple de Simpson es', sim(v,A,B), "metros"
#ley v(t)=t2−2t
from numpy import sqrt,exp,pi,linspace
def v1(t):
#funcion=3*x
velocidad = (t**2)-(2*t)
#funcion=(1/sqrt(2*pi))*exp(-(x**2)/2.0)
return velocidad
#limites
A=0
B=3
v1(3)
%matplotlib inline
import matplotlib.pyplot as plt
x=linspace(A,B)
fig = plt.figure()
axes1 = fig.add_axes([0.1, 0.1, 0.8, 0.9])
axes1.fill_between(x, 0, v1(x), color="red", alpha=0.5)
axes1.set_xlim([A-1,B+1])
plt.show()
#Con regla del trapecio compuesta
a = float(input("Ingresa el valor de a que vamos a evaluar: "))
b = float(input("Ingresa el valor de b que vamos a evaluar: "))
n = input("Ingresa el numero de divisiones que relizaremos: ")
f = ((b-a)/n)
f1 = (v1(a)+v1(b))/2
e = []
for k in range (1,n-1):
c = v1((a+(k*((b-a)/n))))
e.append(c)
z = sum(e)
print "El desplazamiento total del móvil en 4 segundos es de: ", f*(f1+z), "metros"
#Con regla de Simpson
def sim(v1,A,B):
h=(B-A)/6
sim_simp = h*(v1(A)+4*v1((A+B)/2)+v1(B))
return sim_simp
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
print 'El resultado de la regla simple de Simpson es', sim(v1,A,B), "metros"
#Distancia recorrida en 3 segundos
def v2(t):
der = (2*t)-2
return der
v2(3)
#Con regla del trapecio compuesta
a = float(input("Ingresa el valor de a que vamos a evaluar: "))
b = float(input("Ingresa el valor de b que vamos a evaluar: "))
n = input("Ingresa el numero de divisiones que relizaremos: ")
f = ((b-a)/n)
f1 = (v2(a)+v2(b))/2
e = []
for k in range (1,n-1):
c = v2((a+(k*((b-a)/n))))
e.append(c)
z = sum(e)
print "El desplazamiento total del móvil en 4 segundos es de: ", f*(f1+z), "metros"
#Con regla de Simpson
def sim(v2,A,B):
h=(B-A)/6
sim_simp = h*(v2(A)+4*v2((A+B)/2)+v2(B))
return sim_simp
l=raw_input("Cuales son los limites de integracion [a,b]: ")
lim=l.split(',')
A=float(lim[0])
B=float(lim[1])
print 'El resultado de la regla simple de Simpson es', sim(v2,A,B), "metros"