s2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

Ejercicio: 2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

Resultados con el algoritmo

Tabla de resultados
[[ 1.    1.    1.    1.    1.    1.    1.    1.    1.    1.  ]
 [ 0.59  0.58  0.56  0.55  0.54  0.53  0.53  0.52  0.51  0.5 ]
 [-0.31 -0.3  -0.3  -0.29 -0.28 -0.28 -0.27 -0.27 -0.26 -0.26]
 [-0.95 -0.93 -0.91 -0.89 -0.88 -0.86 -0.84 -0.82 -0.81 -0.79]
 [-0.81 -0.79 -0.78 -0.76 -0.74 -0.73 -0.71 -0.7  -0.68 -0.67]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]]

2Eva2022PAOII_T3 EDP Parabolica 02

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 1
Tb = 0
#T0 = 25
fx = lambda x: np.cos(3*np.pi/2*x)
# longitud en x
a = 0
b = 1
# Constante K
K = 2
# Tamaño de paso
dx = 0.2
dt = dx/100
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot): # igual con lazo for
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
    
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

s2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

Ejercicio: 2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

literal a

literal b

Resultado con el algoritmo

 [ ti, xi, yi]
[[0.000000e+00 1.000000e+00 1.000000e-02]
 [5.000000e-01 1.002494e+00 1.109532e-02]
 [1.000000e+00 1.004975e+00 1.220285e-02]
 ...
 [1.995000e+02 1.325161e+00 1.156081e+00]
 [2.000000e+02 1.325242e+00 1.157780e+00]
 [2.005000e+02 1.325322e+00 1.159465e+00]]

2Eva2022PAOII_T2 poblacion protestantes
Instrucciones en Python

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,3),dtype=float)
    tabla[0] = [t0,x0,y0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi]
    tabla = np.array(tabla)
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
b = 0.02
d = 0.015
r = 0.1

# Ecuaciones
f = lambda t,x,y : (b-d*x)*x
g = lambda t,x,y : (b-d*y)*y + r*b*(x-y)

# Condiciones iniciales
t0 = 0
x0 = 1
y0 = 0.01

# parámetros del algoritmo
h = 0.5
muestras = 401

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt

plt.plot(ti,xi, label='xi poblacion')
plt.plot(ti,yi, label='yi protestante')

plt.title('población y protestantes')
plt.xlabel('t años')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

# gráfica xi vs yi
plt.plot(xi,yi)

plt.title('población y protestantes [xi,yi]')
plt.xlabel('x población')
plt.ylabel('y protestantes')
plt.grid()
plt.show()

s2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

Ejercicio: 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

literal a

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt v = 1800 \ln\Big(\frac{160000}{160000-2500t}\Big) - 9.8t

Seleccionando el método de Simpson de 3/8, se requieren al menos 3 tramos o segmentos para usarlo, que generan 4 muestras. El vector de tiempo se obtiene como:

v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3
h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)

siendo los vectores:

ti = [ 0. 10. 20. 30.]
vi = [ 0. 207.81826623 478.44820899 844.54060574]

la aplicación del método de Simpson de 3/8 es:

I = \frac{3}{8}(10) \Bigg(1800 \ln\Big(\frac{160000}{160000-2500(0)}\Big) - 9.8(0) +3(1800 \ln\Big(\frac{160000}{160000-2500(10)}\Big) - 9.8(10)) +3(1800 \ln\Big(\frac{160000}{160000-2500(20)}\Big) - 9.8(20)) +1800 \ln\Big(\frac{160000}{160000-2500(30)}\Big) - 9.8(30) \Bigg) = I = \frac{3}{8}(10) \Big(v(0)+3(v(10))+3(v(20))+v(30) \Big) I = \frac{3}{8}(10) \Big(0+3(207.81)+3(478.44)+844.54 \Big) I = 10887.52

literal b

para el primer segmento se usa t entre [0,10]

x_a = \frac{0+10}{2} + \frac{1}{\sqrt{3}}\frac{10-0}{2} = 7.88 x_b = \frac{0+10}{2} - \frac{1}{\sqrt{3}}\frac{10-0}{2} = 2.11 I = \frac{10-0}{2}\Big(v(7.88)+v(2.11)\Big)=995.79

para el 2do segmento se usa t entre [10,20]

x_a = \frac{10+20}{2} + \frac{1}{\sqrt{3}}\frac{20-10}{2} = 17.88 x_b = \frac{10+20}{2} - \frac{1}{\sqrt{3}}\frac{20-10}{2} = 12.11 I = \frac{20-10}{2}\Big(v(17.88)+v(12.11)\Big) =3368.42

para el 3er segmento se usa t entre [20,30]

x_a = \frac{20+30}{2} + \frac{1}{\sqrt{3}}\frac{30-20}{2} = 27.88 x_b = \frac{20+30}{2} - \frac{1}{\sqrt{3}}\frac{30-20}{2} = 22.11 I = \frac{30-20}{2}\Big(v(27.88)+v(22.11)\Big) = 6515.23 Altura = 995.79+ 3368.42 + 6515.23 = 10879.44

literal c

el error es la diferencia entre los métodos
error_entre = |10887.52-10879.44| = 8.079

Resultados con algoritmo

Método de Simpon 3/8
ti
[ 0. 10. 20. 30.]
vi
[ 0. 207.81826623 478.44820899 844.54060574]
Altura con Simpson 3/8 : 10887.52511781406
segmento Cuad_Gauss :    [995.792, 3368.421, 6515.231]
Altura Cuadratura Gauss: 10879.445437288954
diferencia s3/8 y Cuad_Gauss: 8.079680525106596
>>>

Instrucciones en Python

# 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos
import numpy as np

# INGRESO
v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3

# PROCEDIMIENTO literal a
def integrasimpson38_fi(xi,fi,tolera = 1e-10):
    ''' sobre muestras de fi para cada xi
        integral con método de Simpson 3/8
        respuesta es np.nan para tramos desiguales,
        no hay suficientes puntos.
    '''
    n = len(xi)
    i = 0
    suma = 0
    while not(i>=(n-3)):
        h  = xi[i+1]-xi[i]
        h1 = (xi[i+2]-xi[i+1])
        h2 = (xi[i+3]-xi[i+2])
        dh = abs(h-h1)+abs(h-h2)
        if dh<tolera:# tramos iguales
            unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
            unS38 = (3/8)*h*unS38
            suma = suma + unS38
        else:  # tramos desiguales
            suma = 'tramos desiguales'
        i = i + 3
    if (i+1)<n: # incompleto, tramos por calcular
        suma = 'tramos incompletos, faltan '
        suma = suma +str(n-(i+1))+' tramos'
    return(suma)

h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)
altura = integrasimpson38_fi(ti,vi)

# SALIDA
print('Método de Simpon 3/8')
print('ti')
print(ti)
print('vi')
print(vi)
print('Altura con Simpson 3/8 :',altura)

# PROCEDIMIENTO literal b
# 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)

area = 0
area_i =[]
for i in range(0,tramos,1):
    deltaA = integraCuadGauss2p(v,ti[i],ti[i+1])
    area = area + deltaA
    area_i.append(deltaA)
# SALIDA
print('segmento Cuad_Gauss :   ', area_i)
print('Altura Cuadratura Gauss:', area)

print('diferencia s3/8 y Cuad_Gauss:',altura-area)

import matplotlib.pyplot as plt
plt.plot(ti,vi)
plt.plot(ti,vi,'o')
plt.title('v(t)')
plt.xlabel('t (s)')
plt.ylabel('v (m/s)')
plt.grid()
plt.show()

s1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

Ejercicio: 1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

La variable independiente para la trayectoria es tiempo, con datos en el vector de ti.

ti  = [0, 1, 2, 3, 4]
xti = [2, 1, 3, 4, 2]
yti = [0, 1, 5, 1, 0]

Considerando que los puntos marcan posiciones por donde debe pasar el dron y se define la trayectoria, se usarán todos los puntos. Cada polinomio será de grado 4 al incluir los 5 puntos disponibles para cada eje.

Nota: podría usar polinomios de menor grado, siempre que considere que se debe completar la trayectoria y regresar al punto de salida.

px(t) = 2\frac{(t-1)(t-2)(t-3)(t-4)}{(0-1)(0-2)(0-3)(0-4)} + 1 \frac{(t-0)(t-2)(t-3)(t-4)}{(1-0)(1-2)(1-3)(1-4)} + 3 \frac{(t-0)(t-1)(t-3)(t-4)}{(2-0)(2-1)(2-3)(2-4)} + 4 \frac{(t-0)(t-1)(t-2)(t-4)}{(3-0)(3-1)(3-2)(3-4)} + 2 \frac{(t-0)(t-1)(t-2)(t-3)}{(4-0)(4-1)(4-2)(4-3)}

simplificando con el algoritmo:

px(t) = \frac{1}{12}t^4 - \frac{7}{6}t^3 + \frac{53}{12}t^2 - \frac{13}{3}t + 2

Realizando lo mismo con el algoritmo para polinomio de Lagrange se obtiene:

py(t) = \frac{11}{12}t^4 - \frac{22}{3}t^3 + \frac{205}{12}t^2 - \frac{29}{3}t

se muestra la gráfica de trayectorias por cada eje vs tiempo

Observaciones: La trayectoria usada tiene el mismo punto de salida como de retorno. La trayectoria presenta lóbulos que podrían ser reducidos y minimizar uso de recursos como bateria. Considere usar trazadores cúbicos y observe la misma gráfica de trayectorias x(t) vs y(t).

Resultado con el algoritmo:

Polinomio de Lagrange x: 
x**4/12 - 7*x**3/6 + 53*x**2/12 - 13*x/3 + 2
Polinomio de Lagrange y: 
11*x**4/12 - 22*x**3/3 + 205*x**2/12 - 29*x/3

Algoritmo en Python

# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
ti  = [0,1,2,3,4]
xti = [2,1,3,4,2]
yti = [0,1,5,1,0]

# PROCEDIMIENTO
x = sym.Symbol('x')

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

# para ejex
polinomiox = interpola_lagrange(ti,xti)
polisimplex = polinomiox.expand()
px = sym.lambdify(x,polisimplex)

# para ejey
polinomioy = interpola_lagrange(ti,yti)
polisimpley = polinomioy.expand()
py = sym.lambdify(x,polisimpley)

# Puntos para la gráfica
muestras = 101
a = np.min(ti)
b = np.max(ti)
ti = np.linspace(a,b,muestras)
pxi = px(ti)
pyi = py(ti)

# SALIDA
print('Polinomio de Lagrange x: ')
print(polisimplex)
print('Polinomio de Lagrange y: ')
print(polisimpley)

# Gráfica
figura, enplano = plt.subplots()
plt.scatter(xti,yti, color='red')
plt.plot(pxi,pyi)
plt.ylabel('y(t)')
plt.xlabel('x(t)')
plt.title('trayectoria 2D')
plt.grid()

figura, entiempo = plt.subplots()
plt.plot(ti,pxi, label = 'px')
plt.plot(ti,pyi, label = 'py')
plt.legend()
plt.title('posicion en tiempo')
plt.xlabel('t')
plt.ylabel('p(t)')
plt.grid()

plt.show()

s1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Ejercicios: 1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Se requiere determinar la distribución de cupos en base a los costos relativos al promedio por estudiante para docencia, infraestructura y servicios mostrados en la tabla.

Costo referencial /carrera Mecatrónica Computación Civil Matemáticas
Docencia 1.5 0.9 0.6 0.7
Infraestructura 0.8 1.4 0.4 0.5
Servicios 0.45 0.55 1.1 0.5

Con los datos del total de recursos relativos al promedio por estudiante disponibles son docencia 271, infraestructura 250 y servicios 230.

1.5 a + 0.9 b + 0.6 c + 0.7 d = 271 0.8 a + 1.4 b + 0.4 c + 0.5 d = 250 0.45 a + 0.55 b + 1.1 c + 0.5 d = 230

se indica que en carreras como matemáticas de baja demanda, se establece el cupo de 10,

1.5 a + 0.9 b + 0.6 c + 0.7 (10) = 271 0.8 a + 1.4 b + 0.4 c + 0.5(10) = 250 0.45 a + 0.55 b + 1.1 c + 0.5(10) = 230

el sistema se convierte en:

1.5 a + 0.9 b + 0.6 c = 271 - 0.7 (10) 0.8 a + 1.4 b + 0.4 c = 250 - 0.5(10) 0.45 a + 0.55 b + 1.1 c = 230 - 0.5(10)

Para usar un método iterativo se convierte a matriz aumentada:

\begin{pmatrix} 1.5 & 0.9 & 0.6 & \Big| & 264 \\ 0.8 & 1.4 & 0.4 & \Big| & 245 \\ 0.45 & 0.55 & 1.1 &\Big| & 225 \end{pmatrix}

con pivoteo parcial por filas, la matriz aumentada se mantiene igual, pues los valores de la diagonal ya son los mayores posibles según el algoritmo.

Para un método iterativo se despeja una ecuación por cada incógnita.

a = \frac{1}{1.5}(264 - 0.9 b - 0.6 c) b = \frac{1}{1.4}(245 - 0.8 a - 0.4 c) c = \frac{}{1.1}(225 -0.45 a - 0.55 b)

Para los valores iniciales se consideran números mayores que cero, pues existen recursos para los cupos. No se admiten cupos negativos.

X_0 = [50,50,50]

Las iteraciones para el método iterativo de Gauss-Seidel

itera = 0

a = \frac{1}{1.5}(264 - 0.9 (50) - 0.6 (50)) = 126 b = \frac{1}{1.4}(245 - 0.8 (126) - 0.4 (50)) = 88.714 c = \frac{}{1.1}(225 -0.45 (126) - 0.55 (88.714)) =108.642 diferencia = [126-50, 88.714-50, 108.42-50] diferencia = [76, 38.714, 58.642] errado = max|[76, 38.714, 58.642]| =76 X = [126, 88.714, 108.42]

itera = 1

a = \frac{1}{1.5}(264 - 0.9 (88.714) - 0.6 (108.42)) = 79.314 b = \frac{1}{1.4}(245 - 0.8 (79.314) - 0.4 (108.42)) = 98.637 c = \frac{}{1.1}(225 -0.45 (79.314) - 0.55 (98.637)) =122.780 diferencia = [79.314-126, 88, 98.637-88.714, 122.780-108.42] diferencia = [46.685,9.922, 14.137] errado = max| [46.685,9.922, 14.137] | = 46.685

el error disminuye en la iteración

X = [79.314 , 98.637, 122.780]

itera = 2

a = \frac{1}{1.5}(264 - 0.9 (79.314) - 0.6 (122.780)) = 67.705 b = \frac{1}{1.4}(245 - 0.8 (67.705) - 0.4 (122.780)) = 101.230 c = \frac{}{1.1}(225 -0.45 (67.705) - 0.55 (101.230)) =126.232 diferencia = [67.705-79.314, 101.230-98.637, 126.232-122.780] diferencia = [-11.608, 2.594, 3.451] errado = max| [-11.608, 2.594, 3.451] | = 11.608

el error disminuye en la iteración, se considera que el método converge

X = [67.705 , 101.230, 126.232]

con el algoritmo se tiene como resultado:

[126.          88.71428571 108.64285714]
[76.         38.71428571 58.64285714]

[ 79.31428571  98.63673469 122.78033395]
[46.68571429  9.92244898 14.13747681]

[ 67.7058256  101.23086138 126.23218611]
[11.60846011  2.59412669  3.45185216]

[ 64.76860873 101.92302755 127.08769174]
[2.93721688 0.69216617 0.85550564]

[ 64.01110677 102.11145563 127.30336487]
[0.75750196 0.18842808 0.21567312]

[ 63.81178067 102.16373537 127.3587675 ]
[0.1993261  0.05227973 0.05540263]

[ 63.75825178 102.17849398 127.37328637]
[0.05352889 0.01475862 0.01451887]

[ 63.74358906 102.18272443 127.37716953]
[0.01466272 0.00423045 0.00388316]

[ 63.73949753 102.18395297 127.37822907]
[0.00409153 0.00122854 0.00105954]

[ 63.73833659 102.18431364 127.37852366]
[0.00116094 0.00036067 0.0002946 ]

[ 63.73800235 102.18442047 127.37860699]
[3.34240232e-04 1.06824300e-04 8.33224905e-05]

respuesta X: 
[[ 63.73800235]
 [102.18442047]
 [127.37860699]]
verificar A.X=B: 
[[264.00014614]
 [245.00003333]
 [225.        ]]
>>>

se interpreta la respuesta como la parte entera de la solución:

cupos = [ 63, 102 , 127]

s1Eva_2022PAOII_T1 Esfera flotando en agua

Ejercicio: 1Eva_2022PAOII_T1 Esfera flotando en agua

Según el principio de Arquímedes, la fuerza de flotación o empuje es igual al peso de el fluido desplazado por la porción sumergida de un objeto.

F_{empuje} = F_{peso} \rho_{agua} V_{sumergido} \text{ } g = \rho_{esfera}V_{esfera} \text{ } g V_{sumergido} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{esfera} - V_{sobreagua} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{sobreagua} = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) V_{esfera} V_{esfera} = \frac{4}{3}\pi r^3 \frac{\pi h^2}{3}(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) \frac{4}{3}\pi r^3 h^2(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3

El planteamiento para la búsqueda de raíces es f(x) = 0, que para este caso será:

f(h) = h^2(3r-h) - \Big( 1 - \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3 = 0

usando los valores dados para el ejercicio, r=1 y ρesfera = 200 Kg/m3 y ρagua    = 1000 kg/m3 se tiene que:

f(h) = h^2(3-h) - \Big( 1 - \frac{200}{1000}\Big) 4 f(h) = h^2(3-h) - \frac{16}{5}

Se observa la gráfica de f(h) en el intervalo de h entre[0,2] interpretado como totalmente sumergida y totalmente flotando sobre el agua, confirmando que existe una raíz

Para el caso de aplicar el método del punto fijo se plantea que x=g(x),

h = g(h) h^2(3-h) = \frac{16}{5}

con lo que se puede plantear dos ecuaciones al despejar h

h = \sqrt{ \frac{16}{5(3-h)}} h = 3-\frac{16}{5 h^2}

Iteraciones de la primera ecuación

itera = 0 ; h = h0 = 0.5 ;

g(h) = \sqrt{ \frac{16}{5(3-0.5)}} = 1.1313 tramo = |1.1313-0.5|=0.6313

itera = 1 ; h = 1.1313 ;

g(h) = \sqrt{ \frac{16}{5(3-1.1313)}} = 1.3086 tramo = |1.3086-1.1313| = 0.1772

itera = 2 ; h = 1.3086 ;

g(h) = \sqrt{ \frac{16}{5(3-1.3086)}} = 1.3754 tramo = |1.3754-1.3086| = 0.0668

Observando los errores o tramos en cada iteración se tiene que se reduce, el método converge.


resultados.txt

x,g(x),tramo
0.5 1.131370849898476 0.631370849898476
1.131370849898476 1.308619626317284 0.17724877641880799
1.308619626317284 1.3754802083033437 0.06686058198605971
1.3754802083033437 1.4035002223557855 0.02802001405244181
1.4035002223557855 1.4157629993958152 0.012262777040029649
1.4157629993958152 1.4212317895316 0.005468790135784829
1.4212317895316 1.4236912066694054 0.0024594171378053975
1.4236912066694054 1.424801422465215 0.0011102157958096104
1.424801422465215 1.4253034412081806 0.0005020187429656264
raiz: 1.4253034412081806

Algoritmo en Python

# Algoritmo de punto fijo
# [a,b] intervalo de búsqueda
# error = tolera

import numpy as np
import matplotlib.pyplot as plt

def puntofijo(gx,a,tolera, iteramax = 15):
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    print('x,g(x),tramo')
    print(a,b,tramo)
    while(tramo>=tolera and i<=iteramax ):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        print(a,b,tramo)
        i = i + 1
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax ):
        respuesta = np.nan
    return(respuesta)

# PROGRAMA ---------
# INGRESO
fx = lambda h: h**2*(3-h)-16/5
gx = lambda h: np.sqrt(16/(5*(3-h)))

#fx = lambda h: h**2*(3-h)-16/5
#gx = lambda h: 3-16/(5*(h**2))

x0 = 0.5
tolera = 0.001
iteramax = 50  # itera máximo
a = 0     # intervalo
b = 2
muestras = 51  # gráfico

# PROCEDIMIENTO
respuesta = puntofijo(gx,x0,tolera)

# SALIDA
print('raiz:',respuesta)

hi = np.linspace(a,b,muestras)
fi = fx(hi)
gi = gx(hi)
plt.plot(hi,fi,label='f(h)')
plt.plot(hi,gi,label='g(h)')
plt.plot(hi,hi,label='Identidad')
plt.axhline(0,color='grey')
plt.grid()
plt.xlabel('h')
plt.ylabel('f(h)')
plt.title('esfera sumergida')
plt.legend()
plt.show()

s2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Ejercicio: 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Para la ecuación dada con Δx = 1/3, Δt = 0.02, en una revisíón rápida para cumplir la convergencia dt<dx/10, condición que debe verificarse con la expresión obtenida para λ al desarrollar el ejercicio.

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 0 \leq x \leq 2, t>0

literal a. grafica de malla

literal b. Ecuaciones de diferencias divididas a usar

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 \frac{\partial ^2 U}{\partial x^2} = 9 \frac{\partial U}{\partial t} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = 9 \frac{u_{i,j+1}-u_{i,j}}{\Delta t}

se agrupan las constantes,

\frac{\Delta t}{9(\Delta x)^2} \Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) = u[i,j+1]-u[i,j]

literal d Determine el valor de λ

\lambda = \frac{\Delta t}{9(\Delta x)^2} =\frac{0.02}{9(1/3)^2} = 0.02

valor de λ que es menor que 1/2, por lo que el método converge.

continuando luego con la ecuación general,

\lambda \Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) = u[i,j+1]-u[i,j] \lambda u[i-1,j]-2 \lambda u[i,j] + \lambda u[i+1,j] \Big) = u[i,j+1]-u[i,j]

literal c. Encuentre las ecuaciones considerando las condiciones dadas en el problema.

\lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j] = u[i,j+1]

el punto que no se conoce su valor es u[i,j+1] que es la ecuación buscada.

u[i,j+1] = \lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j]

literal e iteraciones

iteración  i=1, j=0

u[1,1] = \lambda u[0,0]+(1-2 \lambda ) u[1,0] + \lambda u[2,0] u[1,1] =0.02 \cos \Big( \frac{\pi}{2}(0-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}\big(\frac{1}{3}-3\big)\Big) + 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big) \Big) u[1,1] =0.02(0)+(0.96)(-0.5)+0.02(-0.8660)=-0.4973

iteración  i=2, j=0

u[2,1] = \lambda u[1,0]+(1-2 \lambda ) u[2,0] + \lambda u[3,0] u[2,1] = 0.02 \cos \Big( \frac{\pi}{2}(\frac{1}{3}-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}(\frac{2}{3}-3)\Big)+ + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{3}{3}-3\big)\Big) u[2,1] = 0.02 (-0.5) + (0.96 ) (-0.866025) + 0.02 (-1) =-0.8614

iteración  i=3, j=0

u[3,1] = \lambda u[2,0]+(1-2 \lambda ) u[3,0] + \lambda u[4,0] u[3,1] = 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big)\Big)+(1-2 (0.02) ) \cos \Big( \frac{\pi}{2}(1-3)\Big) + + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{4}{3}-3\big)\Big) u[3,1] = 0.02 (-0.866025)+(0.96 ) (-1) + 0.02 (-0,866025) = -0,9946

literal f

la cotas de errores de truncamiento en la ecuación corresponden a segunda derivada O(hx2) y el de primera derivada O(ht), al reemplazar los valores será la suma}

O(hx2) + O(ht) = (1/3)2 + 0.02 = 0,1311

literal g

Resultados usando el algoritmo en Python

Tabla de resultados
[[ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-1.     -0.9946 -0.9893 -0.984  -0.9787 -0.9735 -0.9683 -0.9631 -0.9579  -0.9528]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]]

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito, usando diferencias finitas
# 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 0
Tb = 0
T0 = lambda x: np.cos((np.pi/2)*(x-3))
# longitud en x
a = 0.0
b = 2.0
# Constante K
K = 9
# Tamaño de paso
dx = 1/3
dt = 0.02
tramos = int(np.round((b-a)/dx,0))
muestras = tramos + 1
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.linspace(a,b,muestras)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0(xi[1:ultimox])
u[ultimox,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
    u[0,j+1] = Ta
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    u[m-1,j+1] = Tb
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

s2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

Ejercicio: 2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

La corriente del inductor y(t) para t≥0 se deriva para tener la expresión solo derivadas:

\frac{\delta}{\delta t}y(t) + 2 y(t) + 5 \int_{-\infty}^t y(\tau) \delta \tau = 10 \mu(t)

Para t>0 que es donde transcurre el experimento, el escalón es una constante, se tiene que:

\frac{\delta ^2}{\delta t^2}y(t) + 2 \frac{\delta}{\delta t}y(t) + 5 y(t) = 0

tomando las condiciones iniciales dadas para t=0, y(0)=2, y'(0)=–4

literal a

EL resultadoes perado es el planteamiento del problema. Se reescribe la ecuación con la nomenclatura simplificada y se resordena segun el modelo del método:

y'' = - 2y' - 5 y

luego se sustituye la variable y se convierte a las ecuaciones:

z =y' = f_x(t,y,z) z' = - 2z - 5 y = g_z(t,y,z)

se usa una tabla para llevar el registro de operaciones:

Se plantea las operaciones:

K1y = h * f(ti,yi,zi)
K1z = h * g(ti,yi,zi)

K2y = h * f(ti+h, yi + K1y, zi + K1z)
K2z = h * g(ti+h, yi + K1y, zi + K1z)

yi = yi + (K1y+K2y)/2
zi = zi + (K1z+K2z)/2
ti = ti + h

literal b

El resultado esperado es la aplicación correcta de los valores en las expresiones para al menos tres iteraciones usando h=0.01

itera = 0

K1y = 0.01 y'(0) = 0.01(-4) = -0.04 K1z = 0.01 (- 2z(0) - 5 y(0)) = 0.01(- 2(-4) - 5 (2)) = -0.02 K2y = 0.01 (-4-0.02) = -0.0402 K2z = 0.01 (-2(-4-0.02)-5(2-0.04)) = -0.0176 yi = yi + \frac{K1y+K2y}{2} = 2+\frac{-0.04-0.0402} {2} = 1.9599 zi = zi + \frac{K1z+K2z}{2} = -4 +\frac{-0.02-0.0176}{2} = -4.0188 ti = ti + h = 0+0.01 = 0.01

itera = 1

K1y = 0.01(-4.0188) = -0.040188 K1z = 0.01(- 2(-4.0188) - 5 (1.9599)) = -0.0176 K2y = 0.01 (-4.0188-0.0176) = -0.0403 K2z = 0.01 (-2(-4.0188-0.0176)-5(1.9599-0.040188)) = -0.0152 yi = 1.9599 +\frac{-0.040188-0.0403} {2} = 1.9196 zi = -4.0188 +\frac{-0.0176-0.0152}{2} = -4.0352 ti = ti + h = 0.01+0.01 = 0.02

itera = 2

K1y = 0.01(-4.0352) = -0.040352 K1z = 0.01(- 2(-4.0352) - 5 (1.9196)) = -0.0152 K2y = 0.01 (-4.0352-0.0152) = -0.0405 K2z = 0.01 (-2(-4.0352-0.0152)-5(1.9196-0.040352)) = -0.0129 yi = 1.9196 +\frac{-0.040352-0.0405} {2} =1.8792 zi = -4.0352 +\frac{-0.0152-0.0129}{2} = -4.0494 ti = ti + h = 0.02+0.01 = 0.03

Resultados con el algoritmo en Python

   ti,   yi,    zi,      K1y,    K1z,    K2y,     K2z
[[ 0.00  2.0000 -4.0000  0.0000  0.0000  0.0000   0.0000]
 [ 0.01  1.9599 -4.0188 -0.0400 -0.0200 -0.0402  -0.0176]
 [ 0.02  1.9196 -4.0352 -0.0401 -0.0176 -0.0403  -0.0152]
 [ 0.03  1.8792 -4.0494 -0.0403 -0.0152 -0.0405  -0.0129]
...

Literal c

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

por lo que el error está en el orden de (0.01)3 = 0.000001


Literal d

Se requiere presentar el resultado para el intervalo t entre [0,5]. Siendo el tamaño de paso h=0.01 que es pequeño, se requieren realizar (5-0)/0.01=500 iteraciones, que es más práctico realizarlas usando el algoritmo.

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,7),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -2*z -5*y + 0

t0 = 0
y0 = 2
z0 = -4

h = 0.01
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=4)
print('ti, yi, zi, K1y, K1z, K2y, K2z')
print(tabla)

# GRAFICA
plt.plot(ti,yi, color = 'orange', label='y_RK(t)')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('y(t) con Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Nota: En el curso TELG1001 Señales y Sistemas, la solución se realiza con Transformadas de Laplace

s2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Ejercicio: 2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Literal a. Integral con Simpson 1/3

Para la ecuación en el intervalo x entre [0,3] aplicando dos veces la fórmula en el intervalo se requieren al menos dos tramos cada Simpson de 1/3. Por lo que la cantidad de tramos es 4 (dos por cada formula, y dos fórmulas) que corresponden a 5 muestras.

El tamaño de paso se calcula como:

h = \frac{b-a}{tramos}=\frac{3-0}{4} = \frac{3}{4} = 0.75

representados en una gráfica como

A = \int_0^3 \frac{e^x \sin(x)}{1+x^2} \delta x

con lo que se define la función del integral f(x)

f(x) = \frac{e^x \sin(x)}{1+x^2}

Con lo que aplicando la fórmula se puede obtener los valores de las muestras:

xi= [0. 0.75   1.5    2.25   3.    ]
fi= [0. 0.9235 1.3755 1.2176 0.2834]

Nota: realizar las expresiones completas para las fórmulas si no adjunta el algoritmo en Python

Aplicando Simpson de 1/3 en cada tramo se tiene:

A_s= \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 0 + 4(0.9235) + 1.3755 \Big) + + \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 1.3755 + 4(1.2176) +0.2834 \Big) A_s = 2.8998

Literal b. Integral con Cuadratura de Gauss

Para usar dos veces el método de Cuadratura de Gauss se usan dos intervalos, con lo que las muestras en x serán:

xj= [0. 1.5 3. ]

se calculan los valores para el tramo [0, 1.5]:

x_{a1} = \frac{0+1.5}{2} - \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 0.3169 x_{b1} = \frac{0+1.5}{2} + \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 1.1830 A_{g1} =\frac{1.5-0}{2} \Big( \frac{f(0.3169)+f(1.1830)}{2} \Big) = 1.2361

se calculan los valores para el tramo [1.5, 3]:

x_{a2} = \frac{1.5+3}{2} - \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 1.8169 x_{b2} = \frac{1.5+3}{2} + \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 2.6830 A_{g2} =\frac{3-1.5}{2} \Big( \frac{f(1.8169)+f(2.6830)}{2} \Big) = 1.6329

El total del integral para el intervalo  [0,3]

A_g = A_{g1} + A_{g2} = 2.8691

Al comparar los resultados entre los métodos del literal a y b

errado = |A_s - A_g| = 2.8998 - 2.8691 = 0.0307

Instrucciones integradas en Python

# 2Eva_2022PAOI_T1
# Comparar integrales numéricos Simpson
#   y Cuadratura de Gauss
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: (np.exp(x)*np.sin(x))/(1+x**2)
a = 0
b = 3

# PROCEDIMIENTO
# Aplicando Simpson 1/3
tramos = 4
muestras = tramos+1
xi = np.linspace(a,b,muestras)
fi = fx(xi)
hs = xi[1]-xi[0]
As = (1/3)*(3/4)*(fi[0]+4*fi[1]+fi[2])
As = As + (1/3)*(3/4)*(fi[2]+4*fi[3]+fi[4])
erradoS = 2*(hs**5)

# Aplicando Cuadratura de Gauss
tramosG = 2
muestrasG = tramosG+1
xj = np.linspace(a,b,muestrasG)
hg = xj[1]-xj[0]

xa1 = (xj[0]+xj[1])/2 - (1/np.sqrt(3))*(xj[1]-xj[0])/2
xb1 = (xj[0]+xj[1])/2 + (1/np.sqrt(3))*(xj[1]-xj[0])/2
Ag1 = (hg/2)*(fx(xa1)+fx(xb1))

xa2 = (xj[1]+xj[2])/2 - (1/np.sqrt(3))*(xj[2]-xj[1])/2
xb2 = (xj[1]+xj[2])/2 + (1/np.sqrt(3))*(xj[2]-xj[1])/2
Ag2 = (hg/2)*(fx(xa2)+fx(xb2))

Ag = Ag1 + Ag2

# error entre métodos
errado = np.abs(As-Ag)

# SALIDA
print('xi=',xi)
print('fi=',fi)
print('A Simpson =', As)
print('error Truncamiento Simpson 2*(h^5):',
      erradoS)
print('Cuadratura de Gauss xa,xb por tramos:',
      [xa1,xb1,xa2,xb2])
print('  fx(xa),fx(xb) por tramos:',
      [fx(xa1),fx(xb1),fx(xa2),fx(xb2)])
print('  integral de cada tramo:', [Ag1,Ag2])
print('A Gauss =', Ag)
print('errado entre Simpson y Gauss',errado)

# Grafica con mejor resolucion
xk = np.linspace(a,b,5*tramos+1)
fk = fx(xk)

plt.plot(xk,fk)
plt.plot(xi,fi,'o')
for unx in xi:
    plt.axvline(unx,color='red',
                linestyle='dotted')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('(np.exp(x)*np.sin(x))/(1+x**2)')
plt.grid()
plt.show()

s1Eva_2022PAOI_T3 Interpolar crecimiento de contagios

Ejercicio: 1Eva_2022PAOI_T3 Interpolar crecimiento de contagios

Día del mes 1 8 15 22
Contagios 1 5.6 27 43.5

a) Realice el planteamiento del sistema de ecuaciones que se usaría usando el método de interpolación polinómica.

El modelo de polinomio de grado máximo que se puede obtener es grado 3:

p_3(t) = a_0 + a_1 t + a_2 t^2 + a_3 t^3

por lo que usando los valores de los puntos dados en la tabla:

p_3(1) = a_0 + a_1 (1) + a_2 (1)^2 + a_3 (1)^3 = 1 p_3(8) = a_0 + a_1 (8) + a_2 (8)^2 + a_3 (8)^3 = 5.6 p_3(15) = a_0 + a_1 (15) + a_2 (15)^2 + a_3 (15)^3 = 27 p_3(22) = a_0 + a_1 (22) + a_2 (22)^2 + a_3 (22)^3 = 43.5

b) Realice el planteamiento del sistema de ecuaciones en su forma matricial y muestre la matriz aumentada.

\begin{pmatrix} 1 & 1 & 1^2 & 1^2\\ 1 & 8 & 8^2 & 8^3 \\ 1 & 15 & 15^2 & 15^3 \\ 1 & 22 & 22^2 & 22^3 \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} \begin{pmatrix} 1 \\ 56 \\ 27 \\43.5 \end{pmatrix}

matriz aumentada,

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1 & 8 & 8^2 & 8^3 & 56 \\ 1 & 15 & 15^2 & 15^3 & 27\\ 1 & 22 & 22^2 & 22^3 & 43.5\end{pmatrix}

c) Desarrolle el pivoteo parcial por filas, indicando las operaciones realizadas en éste proceso

pivoteo parcial por filas

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1 & 22 & 22^2 & 22^3 & 43.5 \\ 1 & 15 & 15^2 & 15^3 & 27 \\ 1 & 8 & 8^2 & 8^3 & 56\end{pmatrix}

d) Usando el método directo de Gauss-Jordan, muestre las expresiones necesarias para el algoritmo.

eliminación hacia adelante

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1-1 & 22-1 & 22^2-1^2 & 22^3 -1^3& 43.5 - 1\\ 1-1 & 15-1 & 15^2 -1^2& 15^3 -1^3& 27 -1\\ 1-1 & 8-1 & 8^2 -1^2& 8^3 -1^3& 56-1\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 14 & 224 & 3376& 26\\ 0 & 7 & 63& 511 & 55\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 14-\frac{14}{21} 21 & 224-\frac{14}{21}483& 3376 -\frac{14}{21}10647& 26-\frac{14}{21}42.5\\ 0 & 7-\frac{7}{21}21 & 63-\frac{7}{21}483& 511-\frac{7}{21}10647 & 55-\frac{7}{21}42.5\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & -98 & -3038 & 40.83 \end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & -98-\frac{98}{-98}98 & -3038 -\frac{98}{-98}3722& 40.83 -\frac{98}{-98}(-2.33)\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & 0 & 686 & -7.23\end{pmatrix}

realizando el proceso de eliminación hacia atrás, semejante al método anterior se obtiene

\begin{pmatrix} 1 & 0 & 0 & 0 & 2.98\\ 0 & 1 & 0 & 0& -2.39 \\ 0 & 0 & 1 & 0 & 0.424 \\ 0 & 0 & 0 &1 & -0.0105\end{pmatrix}

con lo que el vector resultado es:

X= [2.98, -2.39, 0.42, -0.0105]

El polinomio de interpolación resultante es:

p(t)= 2.98 -2.39 t + 0.42 t^2 -0.0105 t^3

e) Para el día 19 se encuentra que el valor correspondiente a contagios es de 37%. Estime el error presentado del modelo para ese día.

p(19)= 2.98 -2.39 (19) + 0.42 (19)^2 -0.0105 (19)^3 = 38.42 error = |38.42-37| = 1.42

f) Desarrolle el ejercicio usando otro método para encontrar el polinomio de interpolación.

usando diferencias finitas

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[  0.    1.    1.    4.6  16.8 -21.7   0. ]
 [  1.    8.    5.6  21.4  -4.9   0.    0. ]
 [  2.   15.   27.   16.5   0.    0.    0. ]
 [  3.   22.   43.5   0.    0.    0.    0. ]]

polinomio:

p(t) = 1+\frac{4.6}{1! (7)}(t-1) + + \frac{16.8}{2!(7^2)}(t-1)(t-8) + +\frac{-21.7}{3!(7^3}(t-1)(t-8)(t-15)

Algoritmo en Python

Para literal f

# Polinomio interpolación
# Diferencias finitas avanzadas
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([1,8,15,22],dtype=float)
fi = np.array([1,5.6,27,43.5],dtype=float)

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()