Tarea 3 José Luis Aguilar ipynb - Colaboratory

Vista previa del material en texto

3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 1/8
Tarea 3 Metodos Numéricos
José Luis Aguilar Estrada
import numpy as np
# Sympy libreria para manejar álgebra simbólica
import sympy as sym
import matplotlib.pyplot as plt
# Método de Lagrange
def lagrange(xn: np.ndarray,fn: np.ndarray):
 n = len(xn) # longitud de los datos
 x = sym.Symbol('x') # x variable con Sympy
 sumf = 0 # valor inicial de la suma 
 # índice k para suma y j para poducto
 for k in range(0,n):
 dividendo = 1 # valor inicial del producto numerador
 divisor = 1 # valor inicial del producto denominador
 for j in range(0,n):
 if (j!=k): # condición de no igual
 dividendo = dividendo*(x-xn[j])
 divisor = divisor*(xn[k]-xn[j])
 Lk = dividendo/divisor # L_k en la fórmula
 sumf = sumf + Lk*fn[k] # sumf =f(x) en la fórmula
 # simplifica el polinomio f(x)=sumf obtenido usando Sympy
 polinomio = sumf.expand()
 # Evaluación numérica 
 px = sym.lambdify(x,polinomio)
 return px,polinomio
#Regla de Simpson
def simpson(f,a,b,n):
 h = (b-a)/n # tamaño del intervalo
 sumpar = 0 # inicio suma par
 sumimp = 0 # inicio suma impar
 n0 = int(n/2+1) # para el ciclo for
 # Método de integración
 for i in range(1,n0): # rango de 1 hasta n/2
 ximp = a + (2*i-1)*h # 
 sumimp = sumimp + f(ximp)
 if i < n/2:
 xpar = a + 2*i*h
 sumpar = sumpar + f(xpar)
 # valor de la integral 
 I = (f(a) + f(b) + 4*sumimp + 2*sumpar )*h/3. 
 return I
#Regla del trapecio
def trapecio(f,a,b,n):
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 2/8
def trapecio(f,a,b,n):
 h = (b-a)/n # tamaño del intervalo
 suma = 0 # inicio de suma
 # Método de integración
 for i in range(1,n): # rango de 1 hasta n-1
 xi = a + i*h
 suma = suma + f(xi)
 # valor de la integral 
 I = (f(a) + f(b) + 2*suma )*h/2.
 return I
1
# Datos usando arreglos numpy
xn = np.array([0,np.pi/8,np.pi/4,3*np.pi/8,np.pi/2])
fn = xn*np.sin(xn)
 
 
px,polinomio = lagrange(xn,fn)
 
# Puntos para la gráfica
puntos = 100
a = np.min(xn) # valor minímo de los datos 
b = np.max(xn) # valor minímo de los datos
# se hace la división de puntos desde a hasta b 
pxn = np.linspace(a,b,puntos)
# se evalua en el polinomio simplificado
pfn = px(pxn)
# Salida
 
print('Polinomio de Lagrange: ')
print(polinomio)
Polinomio de Lagrange: 
-0.0908249619904673*x**4 - 0.135031875647655*x**3 + 1.08326240948747*x**2 - 0.0163888
# Gráfica
plt.figure(figsize=(9,6))
#plt.legend()
plt.xlabel('x')
plt.ylabel('f(x)')
plt.xlim(a - abs(0.05*(a-b)),b + abs(0.05*(a-b)))
plt.plot(pxn,pfn,'b',label="Lagrange")
plt.legend(loc='best')
plt.show()
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 3/8
2
Primeramente tenemos que donde 
Entonces
Con lo cual
Y ahora
a + b = mav2 a = = = vdv
dt
dv
dx
dx
dt
dv
dx
a + b = mv ⇒ 1 =v2 dv
dx
mv
a +bv2
dv
dx
Δx = dx = dx = dv∫ xf
x0
∫ xf
x0
mv
a +bv2
dv
dx
∫ 30
15
mv
a +bv2
a = 8.27
b = 2000
m = 5400
 
F = lambda v: (m*v)/(a*v**2 + b)
 
# intervalo de la integral
vi = 15
vf = 30
n = 100 # divisiones
 
 
I = simpson(F,vi,vf,n)
print('El auto recorrera ',I,'m, hasta pararce')
El auto recorrera 292.0088397507007 m, hasta pararce 
3
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 4/8
a)
g = lambda x: ((x**4)*(np.exp(x)))/(np.exp(x) - 1)**(2)
 
 
a = 0.01
b_1 = 1/0.02
b_2 = 1/0.04
b_3 = 1/0.06
b_4 = 1/0.08
n = 100 # divisiones
 
 
#las integrales son
I1 = trapecio(g,a,b_1,n)
I2 = trapecio(g,a,b_2,n)
I3 = trapecio(g,a,b_3,n)
I4 = trapecio(g,a,b_4,n)
 
#Los Cv/3R son
C1 = 3*(1/b1)**3*I1
C2 = 3*(1/b2)**3*I2
C3 = 3*(1/b3)**3*I3
C4 = 3*(1/b4)**3*I4
#Vector de T/T_D
T = [1/b1,1/b2,1/b3,1/b4]
#Vector de c_V/3R para cada T/T_D
C = [C1,C2,C3,C4]
b)
print('Tabla:')
for i in range(0,len(T)):
 print('T/T_D =',T[i],'\tc_V/3R =',C[i])
Tabla: 
T/T_D = 0.02 c_V/3R = 0.0006234081367747244 
T/T_D = 0.04 c_V/3R = 0.004987324157961326 
T/T_D = 0.06 c_V/3R = 0.01682853278424321 
T/T_D = 0.08 c_V/3R = 0.0397015423718304 
px,polinomio = lagrange(T,C)
 
# Salida
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 5/8
# Salida
 
print('Polinomio de Lagrange: ')
print(polinomio)
Polinomio de Lagrange: 
74.0522574210418*x**3 + 0.460344865844083*x**2 - 0.016771211670232*x + 0.000182276364
c)
a = 0.
b = 0.1
n = 100 # divisiones 
 
X = np.linspace(a,b,n)
Y = px(X)
# Gráfica
plt.figure(figsize=(9,6))
#plt.legend()
plt.xlabel('$T/T_D$')
plt.ylabel('$C_v/3R$')
plt.xlim(a - abs(0.05*(a-b)),b + abs(0.05*(a-b)))
plt.plot(X,Y,'b',label="Lagrange")
plt.legend(loc='best')
plt.show()
d)
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 6/8
a = 0.
b = 0.1
n = 100 # divisiones 
 
X_teorico = np.linspace(a,b,n)
# Gráfica
plt.figure(figsize=(12,8))
#plt.legend()
plt.xlabel('$T/T_D$')
plt.ylabel('$C_v/3R$')
plt.xlim(a - abs(0.05*(a-b)),b + abs(0.05*(a-b)))
plt.plot(X,Y,'b',label="Lagrange")
plt.plot(X_teorico,(((4*(np.pi)**(4))/(5))*(X_teorico)**(3)),'r',label="Ecuación teorica")
plt.legend(loc='best')
plt.show()
e)
g = lambda x: ((x**4)*(np.exp(x)))/(np.exp(x) - 1)**(2)
 
 
a = 0.01
n = 100 # divisiones
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 7/8
n = 100 # divisiones 
T_a = np.linspace(0.02, 0.26, 13)
 
 
B_a = 1/T_a
I_a = simpson(g,a,B_a,n)
C_a = 3*(T_a)**3*I_a
#Definimos el integrando
g = lambda x: ((x**4)*(np.exp(x)))/(np.exp(x) - 1)**(2)
 
 
a_I = 0.01
n_I = 100 # divisiones 
n_T = 13# division de los puntos
 
#Para 1/b = 0.2,0.4,...,2.6
T_b = np.linspace(0.2,2.6 , n_T)
 
# Vector de constantes b
B_b = 1./T_b
I_b = simpson(g,a,B_b,n)# Evalúa las integrales
C_b = 3*(T_b)**3*I_b# Evalúa c_V/3R
f)
 
a = np.min(T) 
b = np.max(T_b) 
n = 100 # divisiones 
 
Xt = np.linspace(a,b,puntos)
Yt = (4*(np.pi**4)/5)*Xt**3
# Gráfica
plt.figure(figsize=(12,8))
#plt.legend()
plt.xlabel('$T/T_D$')
plt.ylabel('$C_v/3R$')
plt.xlim(a - abs(0.05*(a-b)),b + abs(0.05*(a-b)))
plt.plot(T_a,C_a,'b',label="1")
plt.plot(T_b,C_b,'r',label="2")
plt.plot(T,C,'k',label="3")
plt.legend(loc='best')
plt.show()
3/12/2020 Tarea 3 José Luis Aguilar.ipynb - Colaboratory
https://colab.research.google.com/drive/17DZC59zS1eMCpCRQUVmYblKgbFp9ayi4#scrollTo=AIYH8AVRCjA2&printMode=true 8/8