In [1]:
import matplotlib.pyplot as plt
import matplotlib as mp
import numpy as np
import math
# Esta linea hace que los graficos aparezcan en el notebook en vez de una ventana nueva.
%matplotlib inline
In [2]:
# Cambia el tamano de los ticks (los puntos de los ejes)
mp.rcParams['xtick.labelsize'] = 13
mp.rcParams['ytick.labelsize'] = 13
In [3]:
# Parametros iniciales.
e = 1.
k_factorial = 1.
N_max = 10
e_vs_n = [e] # Lista que va a contener los elementos de la serie.
for i in range(1,N_max): # Ciclo que calcula los elementos de la serie y los suma.
k_factorial *= i
e += 1. / k_factorial
e_vs_n.append(e)
In [4]:
# Instruccion para imprimir nuestra aproximacion a e en cada iteracion.
for i in range(N_max):
print(i, e_vs_n[i])
In [5]:
plt.plot(range(N_max), e_vs_n) # Genera un grafico de e_vs_n v
plt.axhline(math.e, color='0.5') # Genera una linea en 'e'
plt.ylim(0,3) # Cambia los limites del eje y del grafico
# Cambia los labels del grafico.
plt.xlabel('N', fontsize=20)
plt.ylabel('$e_{N}$', fontsize=20)
Out[5]:
In [6]:
# Se hace un ciclo for dentro de una lista. i.e. Comprehension list
diferencia = [math.fabs(e_i - math.e) for e_i in e_vs_n] # Nota. fabs convierte a float y luego calcular abs
# Se hace un grafico con escala logaritmica en el eje y.
plt.plot(range(N_max), diferencia)
plt.yscale('log')
plt.xlabel('N', fontsize=20)
plt.ylabel('$e_N - e_{real}$', fontsize=20)
"""
NOTA:
Lo anterior es equivalente a hacer
plt.semilogy(range(N_max), diferencia))
plt.xlabel('N', fontsize=20)
plt.ylabel('$e_N - e_{real}$', fontsize=20)
"""
Out[6]:
In [7]:
import numpy as np
In [8]:
epsilon = np.logspace(-1, -15, 14, base=10.)
print(epsilon)
In [9]:
print(type(e_vs_n))
print(type(epsilon))
In [10]:
dsindx = (np.sin(1.+epsilon) - np.sin(1.)) / epsilon
print(dsindx)
In [11]:
plt.semilogx(epsilon, dsindx - np.cos(1.), 'o')
plt.axhline(0, color='0.8')
plt.xlabel('epsilon', fontsize=20)
plt.ylabel('$\\frac{d}{dx}\\sin(x) - \\cos(x)$', fontsize=20)
_ = plt.xticks(epsilon[::2])
Haciendo una expansion de taylor sobre la integral se obtiene:
$$ \int_{x_{0}}^{x_{0}+\Delta x}\left[ f(x_{0}) + f'(x_{o})(x'-x_{0}) + \frac{1}{2}f''(x_{0})(x'-x_{0})^{2} + \ldots \right]dx' $$Note que el error de truncacion es del mismo orden en $\Delta x$ que el ultimo termino previo a truncar.
Donde $O(\Delta x^3) \sim -f''(x_{0})\frac{\Delta x^3}{12}$
(graficos van aca)
dividir [a,b] en $N = \frac{b-a}{\Delta x}$ tramos.
$$ \int_{a}^{b}f(x)dx \sim \sum_{i=0}^{N-1}\left[ \frac{f(a+\Delta xi)+f(a+\Delta x(i+1))}{2}\Delta x + O^{*}(\Delta x^3) \right] $$Nos interesa encontrar una cota superior para nuestro error. Como estamos sumando N tramos, estamos incluyendo N errores, y dado que el error de cada tramo va como $\Delta x^3$ y N va como $\Delta x$, el error finalmente queda como $O(\Delta x^2)$
$$ \sim \left[ \sum_{i=0}^{N-1}\frac{f(a+\Delta xi)+f(a+\Delta x(i+1))}{2}\Delta x \right] + NO^{*}(\Delta x^3) \rightarrow O(\Delta x^2) $$La regla compuesta asumiendo $\Delta x$ constante
In [12]:
def integral_trap(f, intervalo):
'''
Integracion numerica utilizando el metodo o regla trapezoidal.
Note que el intervalo debe estar equiespaciado para este metodo.
Parameters
----------
f : numpy.ndarray
Integrando.
intervalo : numpy.ndarray
El intervalo de integracion.
Returns
-------
res : double
Resultado de la integracion.
'''
res = f[0]+f[len(intervalo-1)] # Se inicializa la variable que contendra el resultado final.
dx = intervalo[1] - intervalo[0] # Se escribe el delta x.
for i in range(len(intervalo)-1): # Se escribe el ciclo que calculara la integral numericamente.
res += 2*(f[i]+f[i+1])
res *= dx/2.
return res
Se hace una suerte de "doble Taylor"
$$ \int_{x_{0}}^{x_{0}+2\Delta x}f(x)dx = \int_{x_{0}}^{x_{0}+2\Delta x}\left[ f(x_{0} + f'(x_{0})(x-x_{0}) + f''(x_{0})\frac{(x-x_{0})^{2}}{2} + f'''(x_{0})\frac{(x-x_{0})^{3}}{6} + f^{\text{iv}}(x_{0})\frac{(x-x_{0})^{4}}{24} + \ldots \right] $$Note que para la regla compuesta el orden del error va como $\Delta x^4$
In [13]:
def simpson(f, a, b, n):
'''
Integracion numerica utilizando el metodo o regla trapezoidal.
Note que el intervalo debe estar equiespaciado para este metodo.
Parameters
----------
f : function
Integrando.
a : int or double
Valor inical del intervalo sobre el cual se va a integrar.
b : int or double
Valor final del intervalo sobre el cual se va a integrar.
n : int
Numero de particiones del intervalo. Debe ser par.
Returns
-------
res : double
Resultado de la integracion.
Raises
------
ValueError
Cuando n es impar.
'''
if n % 2: # El algoritmo funciona solo para n par. Aca nos aseguramos de recibir un n par.
raise ValueError("n debe ser par (se recibio n=%d)." % n)
dx = float(b - a) / n # Se inicializa el dx.
i = a # Se inicializa en a un iterador
res = 0. # Se inicializa res, en donde se guardara el resultado de la integral.
sumaPar = sumaImp = 0. # Se inicializan las sumas parciales.
k = 0
while k < n/2-1: # El ciclo que calcula la integral.
i += dx
sumaImp += f(i)
i += dx
sumaPar += f(i)
k += 1
sumaImp += f(i+dx) # Se agrega el ultimo termino de la suma impar.
res += dx/3.*(f(a) + 4 * sumaImp + 2 * sumaPar + f(b))
return res
Cuando $\Delta x$ no es constante, la expresion anterior pierde validez. En este caso se debe implementar un $\Delta x$ variable y la expresion queda de la siguiente manera:
$$ \int_{a}^{b}f(x)dx \sim \frac{1}{2}\sum_{i=0}^{N-1}(x_{i+1}-x_{i})(f(x_{i+1})+f(x_{i})) $$