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()

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()

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]

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/m³ y ρ_agua    = 1000 kg/m³ 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 = h₀ = 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()

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. gráfica 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(h_x²) y el de primera derivada O(h_t), al reemplazar los valores será la suma}

O(h_x²) + O(h_t) = (1/3)² + 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()

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(h³)

por lo que el error está en el orden de (0.01)³ = 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

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( f(0.3169)+f(1.1830)\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( f(1.8169)+f(2.6830)\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()

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 math
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 = 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()

Ejercicio: 1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos

a) Realice el planteamiento del sistema de ecuaciones que permita determinar la cantidad máxima de pacientes de cada grupo que podrían ser atendidos usando todos los productos disponibles. Una vez planteadas las ecuaciones, se le indica que la capacidad de atención para emergencia sea fija en K = 10 pacientes (variable libre).

Producto\ Paciente	Maternidad	Pos - operatorio	Covid_19	emergencia	Suministro diario
Producto A	0.2	0.1	1.7	0.25	135
Producto B	0.5	2	0.05	0.4	320
Producto C	1.5	0.2	0.75	1.4	410

Sistema de ecuaciones usando la columna de emergencias como variable libre y valor constante para la variable igual a K=10

0.2 x_1 + 0.1 x_2 +1.7 x_3 = 135-0.25 K 0.5 x_1 + 2 x_2 +0.05 x_3 = 320-0.4 K 1.5 x_1 + 0.2 x_2 +0.75 x_3 = 410-1.4 K

b) Muestre los pasos detallados para la matriz aumentada y pivoteo parcial por filas.

matriz aumentada

\begin{pmatrix} 0.2 & 0.1 & 1.7 & 135-0.25*10 \\ 0.5 & 2 &0.05 & 320-0.4*10 \\ 1.5 & 0.2 & 0.75 &410-1.4*10 \end{pmatrix}

pivoteo parcial por filas

\begin{pmatrix} 1.5 & 0.2 & 0.75 & 396 \\ 0.5 & 2 &0.05 & 316 \\ 0.2 & 0.1 & 1.7 & 132.5 \end{pmatrix}

c) Desarrolle al menos 3 iteraciones para el método requerido, con expresiones completas.
1.5 x_1 + 0.2 x_2 +0.75 x_3 = 396

0.5 x_1 + 2 x_2 +0.05 x_3 = 316 0.2 x_1 + 0.1 x_2 +1.7 x_3 = 132.5

ecuaciones a usar

x_1 = \frac{1}{1.5}(396 - 0.2 x_2 - 0.75 x_3) x_2 = \frac{1}{2}(316 - 0.5 x_1 - 0.05 x_3) x_3 = \frac{1}{1.7}(132.5 - 0.2 x_1 - 0.1 x_2)

Dado que el valor de la variable libre se establecía con K=10, el vector inicial podría ser el doble de este valor

X_0 = [20,20,20]

itera=1

X_0 = [20,20,20] x_1 = \frac{1}{1.5}(396 - 0.2(20) - 0.75(20)) =251.33 x_2 = \frac{1}{2}(316 - 0.5(20) - 0.05 (20)) = 152.5 x_3 = \frac{1}{1.7}(132.5 - 0.2 (20) - 0.1 (20)) =74.41 diferencia = [231.33, 132.5, 54.41] errado = max |[231.33, 132.5, 54.41]| = 231.33

itera=2

X_1 = [251.33, 152.5, 74.41] x_1 = \frac{1}{1.5}(396 - 0.2(152.5) - 0.75(74.41)) =206.46 x_2 = \frac{1}{2}(316 - 0.5(251.33) - 0.05 (74.41)) = 96.30 x_3 = \frac{1}{1.7}(132.5 - 0.2 (251.33) - 0.1 (152.5)) =39.40 diferencia = [-44.86, -59.27, -211.92] errado = max |[-44.86, -59.27, -211.92]| = 211.92

itera=3

X_2 = [231.85, 105.39, 48.16] x_1 = \frac{1}{1.5}(396 - 0.2(105.39) - 0.75(48.16)) =231.85 x_2 = \frac{1}{2}(316 - 0.5(231.85) - 0.05 (48.16)) = 105.39 x_3 = \frac{1}{1.7}(132.5 - 0.2 ( 231.85) - 0.1 (105.39)) =48.16 diferencia = [25.39, 121.09, -158.29] errado = max |[25.39, 121.09, -158.29]| = 158.29 X_3 = [231.85, 105.39, 48.16]

d) Realice las observaciones necesarias sobre los errores entre iteraciones y la convergencia.

El error entre iteraciones disminuye, el método converge a:

respuesta X: 
[[228.22]
 [ 99.81]
 [ 44.92]]

e) Si se decide no atender a los pacientes del grupo emergencias, ¿Qué aumento individual de cada una de otros grupos de pacientes podría soportarse con la cantidad diaria de alimento disponible? (use el algoritmo.py).

Se interpreta como K=0 y usando el algoritmo se obtiene:

respuesta X: 
[[237.23]
 [ 99.54]
 [ 45.64]]

el aumento de pacientes entre grupos es

diferencia = [9.01, -0.27, 0.72]

en términos de pacientes, como número entero, solo se gana 9 pacientes para el primer grupo. Se descarta la parte decimal si los pacientes no se cuentan en miles, por lo que las otras variaciones se interpretan como cero.

Algoritmo en Python

Datos tomados luego de pivoteo parcial por filas

Resultados:

0 [20. 20. 20.]
0 [251.33333333 152.5         74.11764706]
1 [206.60784314  93.31372549  39.10784314]
2 [232.00424837 105.37034314  47.85121107]
3 [226.02501538  98.80265763  44.15418589]
4 [228.74921937 100.38989151  45.24395951]
5 [227.99270138  99.68159617  44.83009822]
6 [228.2940714   99.8810722   44.96076477]
7 [228.20214132  99.80246303  44.91357559]
8 [228.23621714  99.82662528  44.92901496]
9 [228.22527582  99.81772034  44.92358473]
10 [228.22917825  99.82059143  44.92539577]
11 [228.22788993  99.81957054  44.92476777]
12 [228.22834004  99.81990832  44.92497939]
13 [228.2281892   99.8197905   44.92490656]
14 [228.22824132  99.81983004  44.92493124]
numero de condicion: 2.166985328561448
respuesta con Jacobi
[[228.22824132]
 [ 99.81983004]
 [ 44.92493124]]
verificando:
[[396.00002641]
 [316.00002729]
 [132.00001438]]

Algoritmo en Python

# 1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos
import numpy as np

def jacobi(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    xnuevo = np.copy(X)

    itera = 0
    print(itera, X)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        print(itera, xnuevo)
        X = np.copy(xnuevo)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
    return(X)


# INGRESO
A = np.array([[1.5, 0.2, 0.75],
              [0.5, 2.0, 0.05],
              [0.2, 0.1, 1.7 ]],
             dtype=float)

B = np.array([396.,316.,132.],
             dtype=float)
tolera = 1e-4

X = np.array([20.,20.,20.],
             dtype=float)

# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

respuesta = jacobi(A,B,tolera,X)

verifica = np.dot(A,respuesta)

# SALIDA
print('numero de condicion:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('verificando:')
print(verifica)

Ejercicio: 1Eva_2022PAOI_T1 Impacto en trayectoria del drone

Desarrollo analítico

a) Realice el planteamiento del problema usando inicialmente las trayectorias en el eje x, donde para el intervalo de operación del misil antidrone, se observa más de un impacto.

x₁(t) = x₂(t)

f(t) = cos(t) -  sin(0.75 t) =0

y₁(t) = y₂(t)

sin(2 t) =kt

k = \frac{sin(2 t)}{t}

b) Usando el método de Newton-Raphson encuentre el valor de t en el cual se pretende realizar el impacto al drone. Realice al menos 3 iteraciones de forma analítica, use tolerancia de 10^-4,

f(t) = cos(t) - sin(0.75 t) f'(t) = - sin(t) - 0.75 cos(0.75 t)

Como punto inicial para encontrar la raíz de f(t) podría ser t₀=4 para el punto marcado en rojo. Para el método de Newton-Raphson se tiene que

t_{i+1} = t_i - \frac{f(t_i)}{f'(t_i)} error = |t_{i+1} - t_i|

iteración 1 t₀=4

t_1 = 4 - \frac{cos(4) - sin(0.75*4)}{- sin(4) - 0.75 cos(0.75*4)} t_1 = 4 - \frac{-0.7947}{1.4992} = 4.5300 tramo = |4.5300-4| = 0.53

iteración 2 t₁=4.53

t_{2} = 4.53 - \frac{cos(4.53) - sin(0.75*4.53)}{- sin(4.53) - 0.75 cos(0.75*4.53)} t_2 = 4.53 - \frac{-0.0717}{1.7089} = 4.4880 tramo= |4.4880-4.53| = 0.042

iteración 3 t₂=4.4880

t_{3} = 4.4880 - \frac{cos(4.4880) - sin(0.75*4.4880)}{ - sin(4.4880) - 0.75 cos(0.75*4.4880)} t_3 = 4.4880 - \frac{-0.0000179}{1.7061} = 4.4879 tramo= |4.4879-4.4880| = 0.0001

c) Realice el análisis de la convergencia del método.

El error disminuye, el método converge. La raiz se encuentra en t=4.487989505154422

d) Con el resultado de t anterior, determine el valor de la constante k para la expresión de y₂(t) que asegura el impacto contra el drone.

y_1(t) = y_2(t) sin(2 t) =kt k = \frac{sin(2 t)}{t} = \frac{sin(2*4.487989505154422)}{4.487989505154422} = 0.096714

Desarrollo con Algoritmo

Resultados

['xi', 'xnuevo', 'tramo']
[[4.0000e+00 4.5301e+00 5.3009e-01]
 [4.5301e+00 4.4880e+00 4.2071e-02]
 [4.4880e+00 4.4880e+00 3.0277e-05]]
raiz en:  4.487989505154422
con error de:  3.0276981949128867e-05

Algoritmo en Python

# Método de Newton-Raphson
import numpy as np

# INGRESO
fx  = lambda t: np.cos(1*t) - np.sin(0.75*t)
dfx = lambda t: -np.sin(t) - np.cos(0.75*t)*0.75

x0 = 4
tolera = 0.0001

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)